ROOT  6.06/09
Reference Guide
MatrixInversion.icc
Go to the documentation of this file.
1 // @(#)root/smatrix:$Id$
2 // Authors: CLHEP authors, L. Moneta 2006
3 
4 #ifndef ROOT_Math_MatrixInversion_icc
5 #define ROOT_Math_MatrixInversion_icc
6 
7 
8 #include "Math/SVector.h"
9 #include "Math/Math.h"
10 #include <limits>
11 
12 
13 // inversion algorithms for matrices
14 // taken from CLHEP (L. Moneta May 2006)
15 
16 namespace ROOT {
17 
18  namespace Math {
19 
20 
21  /** General Inversion for a symmetric matrix
22  Bunch-Kaufman diagonal pivoting method
23  It is decribed in J.R. Bunch, L. Kaufman (1977).
24  "Some Stable Methods for Calculating Inertia and Solving Symmetric
25  Linear Systems", Math. Comp. 31, p. 162-179. or in Gene H. Golub,
26  /Charles F. van Loan, "Matrix Computations" (the second edition
27  has a bug.) and implemented in "lapack"
28  Mario Stanke, 09/97
29 
30  */
31 
32 
33 
34 template <unsigned int idim, unsigned int N>
35 template<class T>
37 
38  typedef T value_type;
39  //typedef double value_type; // for making inversions in double's
40 
41 
42  int i, j, k, s;
43  int pivrow;
44 
45  const int nrow = MatRepSym<T,idim>::kRows;
46 
47  // Establish the two working-space arrays needed: x and piv are
48  // used as pointers to arrays of doubles and ints respectively, each
49  // of length nrow. We do not want to reallocate each time through
50  // unless the size needs to grow. We do not want to leak memory, even
51  // by having a new without a delete that is only done once.
52 
53 
56 
57  typedef int* pivIter;
58  typedef T* mIter;
59 
60 
61  // Note - resize shuld do nothing if the size is already larger than nrow,
62  // but on VC++ there are indications that it does so we check.
63  // Note - the data elements in a vector are guaranteed to be contiguous,
64  // so x[i] and piv[i] are optimally fast.
65  mIter x = xvec.begin();
66  // x[i] is used as helper storage, needs to have at least size nrow.
67  pivIter piv = pivv.begin();
68  // piv[i] is used to store details of exchanges
69 
70  value_type temp1, temp2;
71  mIter ip, mjj, iq;
72  value_type lambda, sigma;
73  const value_type alpha = .6404; // = (1+sqrt(17))/8
74  // LM (04/2009) remove this useless check (it is not in LAPACK) which fails inversion of
75  // a matrix with values < epsilon in the diagonal
76  //
77  //const double epsilon = 32*std::numeric_limits<T>::epsilon();
78  // whenever a sum of two doubles is below or equal to epsilon
79  // it is set to zero.
80  // this constant could be set to zero but then the algorithm
81  // doesn't neccessarily detect that a matrix is singular
82 
83  for (i = 0; i < nrow; i++)
84  piv[i] = i+1;
85 
86  ifail = 0;
87 
88  // compute the factorization P*A*P^T = L * D * L^T
89  // L is unit lower triangular, D is direct sum of 1x1 and 2x2 matrices
90  // L and D^-1 are stored in A = *this, P is stored in piv[]
91 
92  for (j=1; j < nrow; j+=s) // main loop over columns
93  {
94  mjj = rhs.Array() + j*(j-1)/2 + j-1;
95  lambda = 0; // compute lambda = max of A(j+1:n,j)
96  pivrow = j+1;
97  ip = rhs.Array() + (j+1)*j/2 + j-1;
98  for (i=j+1; i <= nrow ; ip += i++)
99  if (std::abs(*ip) > lambda)
100  {
101  lambda = std::abs(*ip);
102  pivrow = i;
103  }
104 
105  if (lambda == 0 )
106  {
107  if (*mjj == 0)
108  {
109  ifail = 1;
110  return;
111  }
112  s=1;
113  *mjj = 1.0f / *mjj;
114  }
115  else
116  {
117  if (std::abs(*mjj) >= lambda*alpha)
118  {
119  s=1;
120  pivrow=j;
121  }
122  else
123  {
124  sigma = 0; // compute sigma = max A(pivrow, j:pivrow-1)
125  ip = rhs.Array() + pivrow*(pivrow-1)/2+j-1;
126  for (k=j; k < pivrow; k++)
127  {
128  if (std::abs(*ip) > sigma)
129  sigma = std::abs(*ip);
130  ip++;
131  }
132  // sigma cannot be zero because it is at least lambda which is not zero
133  if ( std::abs(*mjj) >= alpha * lambda * (lambda/ sigma) )
134  {
135  s=1;
136  pivrow = j;
137  }
138  else if (std::abs(*(rhs.Array()+pivrow*(pivrow-1)/2+pivrow-1))
139  >= alpha * sigma)
140  s=1;
141  else
142  s=2;
143  }
144  if (pivrow == j) // no permutation neccessary
145  {
146  piv[j-1] = pivrow;
147  if (*mjj == 0)
148  {
149  ifail=1;
150  return;
151  }
152  temp2 = *mjj = 1.0f/ *mjj; // invert D(j,j)
153 
154  // update A(j+1:n, j+1,n)
155  for (i=j+1; i <= nrow; i++)
156  {
157  temp1 = *(rhs.Array() + i*(i-1)/2 + j-1) * temp2;
158  ip = rhs.Array()+i*(i-1)/2+j;
159  for (k=j+1; k<=i; k++)
160  {
161  *ip -= static_cast<T> ( temp1 * *(rhs.Array() + k*(k-1)/2 + j-1) );
162 // if (std::abs(*ip) <= epsilon)
163 // *ip=0;
164  ip++;
165  }
166  }
167  // update L
168  ip = rhs.Array() + (j+1)*j/2 + j-1;
169  for (i=j+1; i <= nrow; ip += i++)
170  *ip *= static_cast<T> ( temp2 );
171  }
172  else if (s==1) // 1x1 pivot
173  {
174  piv[j-1] = pivrow;
175 
176  // interchange rows and columns j and pivrow in
177  // submatrix (j:n,j:n)
178  ip = rhs.Array() + pivrow*(pivrow-1)/2 + j;
179  for (i=j+1; i < pivrow; i++, ip++)
180  {
181  temp1 = *(rhs.Array() + i*(i-1)/2 + j-1);
182  *(rhs.Array() + i*(i-1)/2 + j-1)= *ip;
183  *ip = static_cast<T> ( temp1 );
184  }
185  temp1 = *mjj;
186  *mjj = *(rhs.Array()+pivrow*(pivrow-1)/2+pivrow-1);
187  *(rhs.Array()+pivrow*(pivrow-1)/2+pivrow-1) = static_cast<T> (temp1 );
188  ip = rhs.Array() + (pivrow+1)*pivrow/2 + j-1;
189  iq = ip + pivrow-j;
190  for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
191  {
192  temp1 = *iq;
193  *iq = *ip;
194  *ip = static_cast<T>( temp1 );
195  }
196 
197  if (*mjj == 0)
198  {
199  ifail = 1;
200  return;
201  }
202  temp2 = *mjj = 1.0f / *mjj; // invert D(j,j)
203 
204  // update A(j+1:n, j+1:n)
205  for (i = j+1; i <= nrow; i++)
206  {
207  temp1 = *(rhs.Array() + i*(i-1)/2 + j-1) * temp2;
208  ip = rhs.Array()+i*(i-1)/2+j;
209  for (k=j+1; k<=i; k++)
210  {
211  *ip -= static_cast<T> (temp1 * *(rhs.Array() + k*(k-1)/2 + j-1) );
212 // if (std::abs(*ip) <= epsilon)
213 // *ip=0;
214  ip++;
215  }
216  }
217  // update L
218  ip = rhs.Array() + (j+1)*j/2 + j-1;
219  for (i=j+1; i<=nrow; ip += i++)
220  *ip *= static_cast<T>( temp2 );
221  }
222  else // s=2, ie use a 2x2 pivot
223  {
224  piv[j-1] = -pivrow;
225  piv[j] = 0; // that means this is the second row of a 2x2 pivot
226 
227  if (j+1 != pivrow)
228  {
229  // interchange rows and columns j+1 and pivrow in
230  // submatrix (j:n,j:n)
231  ip = rhs.Array() + pivrow*(pivrow-1)/2 + j+1;
232  for (i=j+2; i < pivrow; i++, ip++)
233  {
234  temp1 = *(rhs.Array() + i*(i-1)/2 + j);
235  *(rhs.Array() + i*(i-1)/2 + j) = *ip;
236  *ip = static_cast<T>( temp1 );
237  }
238  temp1 = *(mjj + j + 1);
239  *(mjj + j + 1) =
240  *(rhs.Array() + pivrow*(pivrow-1)/2 + pivrow-1);
241  *(rhs.Array() + pivrow*(pivrow-1)/2 + pivrow-1) = static_cast<T>( temp1 );
242  temp1 = *(mjj + j);
243  *(mjj + j) = *(rhs.Array() + pivrow*(pivrow-1)/2 + j-1);
244  *(rhs.Array() + pivrow*(pivrow-1)/2 + j-1) = static_cast<T>( temp1 );
245  ip = rhs.Array() + (pivrow+1)*pivrow/2 + j;
246  iq = ip + pivrow-(j+1);
247  for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
248  {
249  temp1 = *iq;
250  *iq = *ip;
251  *ip = static_cast<T>( temp1 );
252  }
253  }
254  // invert D(j:j+1,j:j+1)
255  temp2 = *mjj * *(mjj + j + 1) - *(mjj + j) * *(mjj + j);
256  if (temp2 == 0)
257  std::cerr
258  << "SymMatrix::bunch_invert: error in pivot choice"
259  << std::endl;
260  temp2 = 1. / temp2;
261  // this quotient is guaranteed to exist by the choice
262  // of the pivot
263  temp1 = *mjj;
264  *mjj = static_cast<T>( *(mjj + j + 1) * temp2 );
265  *(mjj + j + 1) = static_cast<T>( temp1 * temp2 );
266  *(mjj + j) = static_cast<T>( - *(mjj + j) * temp2 );
267 
268  if (j < nrow-1) // otherwise do nothing
269  {
270  // update A(j+2:n, j+2:n)
271  for (i=j+2; i <= nrow ; i++)
272  {
273  ip = rhs.Array() + i*(i-1)/2 + j-1;
274  temp1 = *ip * *mjj + *(ip + 1) * *(mjj + j);
275 // if (std::abs(temp1 ) <= epsilon)
276 // temp1 = 0;
277  temp2 = *ip * *(mjj + j) + *(ip + 1) * *(mjj + j + 1);
278 // if (std::abs(temp2 ) <= epsilon)
279 // temp2 = 0;
280  for (k = j+2; k <= i ; k++)
281  {
282  ip = rhs.Array() + i*(i-1)/2 + k-1;
283  iq = rhs.Array() + k*(k-1)/2 + j-1;
284  *ip -= static_cast<T>( temp1 * *iq + temp2 * *(iq+1) );
285 // if (std::abs(*ip) <= epsilon)
286 // *ip = 0;
287  }
288  }
289  // update L
290  for (i=j+2; i <= nrow ; i++)
291  {
292  ip = rhs.Array() + i*(i-1)/2 + j-1;
293  temp1 = *ip * *mjj + *(ip+1) * *(mjj + j);
294 // if (std::abs(temp1) <= epsilon)
295 // temp1 = 0;
296  *(ip+1) = *ip * *(mjj + j)
297  + *(ip+1) * *(mjj + j + 1);
298 // if (std::abs(*(ip+1)) <= epsilon)
299 // *(ip+1) = 0;
300  *ip = static_cast<T>( temp1 );
301  }
302  }
303  }
304  }
305  } // end of main loop over columns
306 
307  if (j == nrow) // the the last pivot is 1x1
308  {
309  mjj = rhs.Array() + j*(j-1)/2 + j-1;
310  if (*mjj == 0)
311  {
312  ifail = 1;
313  return;
314  }
315  else
316  *mjj = 1.0f / *mjj;
317  } // end of last pivot code
318 
319  // computing the inverse from the factorization
320 
321  for (j = nrow ; j >= 1 ; j -= s) // loop over columns
322  {
323  mjj = rhs.Array() + j*(j-1)/2 + j-1;
324  if (piv[j-1] > 0) // 1x1 pivot, compute column j of inverse
325  {
326  s = 1;
327  if (j < nrow)
328  {
329  ip = rhs.Array() + (j+1)*j/2 + j-1;
330  for (i=0; i < nrow-j; ip += 1+j+i++)
331  x[i] = *ip;
332  for (i=j+1; i<=nrow ; i++)
333  {
334  temp2=0;
335  ip = rhs.Array() + i*(i-1)/2 + j;
336  for (k=0; k <= i-j-1; k++)
337  temp2 += *ip++ * x[k];
338  for (ip += i-1; k < nrow-j; ip += 1+j+k++)
339  temp2 += *ip * x[k];
340  *(rhs.Array()+ i*(i-1)/2 + j-1) = static_cast<T>( -temp2 );
341  }
342  temp2 = 0;
343  ip = rhs.Array() + (j+1)*j/2 + j-1;
344  for (k=0; k < nrow-j; ip += 1+j+k++)
345  temp2 += x[k] * *ip;
346  *mjj -= static_cast<T>( temp2 );
347  }
348  }
349  else //2x2 pivot, compute columns j and j-1 of the inverse
350  {
351  if (piv[j-1] != 0)
352  std::cerr << "error in piv" << piv[j-1] << std::endl;
353  s=2;
354  if (j < nrow)
355  {
356  ip = rhs.Array() + (j+1)*j/2 + j-1;
357  for (i=0; i < nrow-j; ip += 1+j+i++)
358  x[i] = *ip;
359  for (i=j+1; i<=nrow ; i++)
360  {
361  temp2 = 0;
362  ip = rhs.Array() + i*(i-1)/2 + j;
363  for (k=0; k <= i-j-1; k++)
364  temp2 += *ip++ * x[k];
365  for (ip += i-1; k < nrow-j; ip += 1+j+k++)
366  temp2 += *ip * x[k];
367  *(rhs.Array()+ i*(i-1)/2 + j-1) = static_cast<T>( -temp2 );
368  }
369  temp2 = 0;
370  ip = rhs.Array() + (j+1)*j/2 + j-1;
371  for (k=0; k < nrow-j; ip += 1+j+k++)
372  temp2 += x[k] * *ip;
373  *mjj -= static_cast<T>( temp2 );
374  temp2 = 0;
375  ip = rhs.Array() + (j+1)*j/2 + j-2;
376  for (i=j+1; i <= nrow; ip += i++)
377  temp2 += *ip * *(ip+1);
378  *(mjj-1) -= static_cast<T>( temp2 );
379  ip = rhs.Array() + (j+1)*j/2 + j-2;
380  for (i=0; i < nrow-j; ip += 1+j+i++)
381  x[i] = *ip;
382  for (i=j+1; i <= nrow ; i++)
383  {
384  temp2 = 0;
385  ip = rhs.Array() + i*(i-1)/2 + j;
386  for (k=0; k <= i-j-1; k++)
387  temp2 += *ip++ * x[k];
388  for (ip += i-1; k < nrow-j; ip += 1+j+k++)
389  temp2 += *ip * x[k];
390  *(rhs.Array()+ i*(i-1)/2 + j-2)= static_cast<T>( -temp2 );
391  }
392  temp2 = 0;
393  ip = rhs.Array() + (j+1)*j/2 + j-2;
394  for (k=0; k < nrow-j; ip += 1+j+k++)
395  temp2 += x[k] * *ip;
396  *(mjj-j) -= static_cast<T>( temp2 );
397  }
398  }
399 
400  // interchange rows and columns j and piv[j-1]
401  // or rows and columns j and -piv[j-2]
402 
403  pivrow = (piv[j-1]==0)? -piv[j-2] : piv[j-1];
404  ip = rhs.Array() + pivrow*(pivrow-1)/2 + j;
405  for (i=j+1;i < pivrow; i++, ip++)
406  {
407  temp1 = *(rhs.Array() + i*(i-1)/2 + j-1);
408  *(rhs.Array() + i*(i-1)/2 + j-1) = *ip;
409  *ip = static_cast<T>( temp1 );
410  }
411  temp1 = *mjj;
412  *mjj = *(rhs.Array() + pivrow*(pivrow-1)/2 + pivrow-1);
413  *(rhs.Array() + pivrow*(pivrow-1)/2 + pivrow-1) = static_cast<T>( temp1 );
414  if (s==2)
415  {
416  temp1 = *(mjj-1);
417  *(mjj-1) = *( rhs.Array() + pivrow*(pivrow-1)/2 + j-2);
418  *( rhs.Array() + pivrow*(pivrow-1)/2 + j-2) = static_cast<T>( temp1 );
419  }
420 
421  ip = rhs.Array() + (pivrow+1)*pivrow/2 + j-1; // &A(i,j)
422  iq = ip + pivrow-j;
423  for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
424  {
425  temp1 = *iq;
426  *iq = *ip;
427  *ip = static_cast<T>(temp1);
428  }
429  } // end of loop over columns (in computing inverse from factorization)
430 
431  return; // inversion successful
432 
433 }
434 
435 
436 
437 /**
438  LU factorization : code originally from CERNLIB dfact routine and ported in C++ for CLHEP
439 */
440 
441 template <unsigned int idim, unsigned int n>
442 template<class T>
443 int Inverter<idim,n>::DfactMatrix(MatRepStd<T,idim,n> & rhs, T &det, unsigned int *ir) {
444 
445  if (idim != n) return -1;
446 
447  int ifail, jfail;
448 
449  typedef T* mIter;
450 
451  typedef T value_type;
452  //typedef double value_type; // for making inversions in double's
453 
454 
455  value_type tf;
456  value_type g1 = 1.0e-19, g2 = 1.0e19;
457 
458  value_type p, q, t;
459  value_type s11, s12;
460 
461  // LM (04.09) : remove useless check on epsilon and set it to zero
462  const value_type epsilon = 0.0;
463  //double epsilon = 8*std::numeric_limits<T>::epsilon();
464  // could be set to zero (like it was before)
465  // but then the algorithm often doesn't detect
466  // that a matrix is singular
467 
468  int normal = 0, imposs = -1;
469  int jrange = 0, jover = 1, junder = -1;
470  ifail = normal;
471  jfail = jrange;
472  int nxch = 0;
473  det = 1.0;
474  mIter mj = rhs.Array();
475  mIter mjj = mj;
476  for (unsigned int j=1;j<=n;j++) {
477  unsigned int k = j;
478  p = (std::abs(*mjj));
479  if (j!=n) {
480  mIter mij = mj + n + j - 1;
481  for (unsigned int i=j+1;i<=n;i++) {
482  q = (std::abs(*(mij)));
483  if (q > p) {
484  k = i;
485  p = q;
486  }
487  mij += n;
488  }
489  if (k==j) {
490  if (p <= epsilon) {
491  det = 0;
492  ifail = imposs;
493  jfail = jrange;
494  return ifail;
495  }
496  det = -det; // in this case the sign of the determinant
497  // must not change. So I change it twice.
498  }
499  mIter mjl = mj;
500  mIter mkl = rhs.Array() + (k-1)*n;
501  for (unsigned int l=1;l<=n;l++) {
502  tf = *mjl;
503  *(mjl++) = *mkl;
504  *(mkl++) = static_cast<T>(tf);
505  }
506  nxch = nxch + 1; // this makes the determinant change its sign
507  ir[nxch] = (((j)<<12)+(k));
508  } else {
509  if (p <= epsilon) {
510  det = 0.0;
511  ifail = imposs;
512  jfail = jrange;
513  return ifail;
514  }
515  }
516  det *= *mjj;
517  *mjj = 1.0f / *mjj;
518  t = (std::abs(det));
519  if (t < g1) {
520  det = 0.0;
521  if (jfail == jrange) jfail = junder;
522  } else if (t > g2) {
523  det = 1.0;
524  if (jfail==jrange) jfail = jover;
525  }
526  if (j!=n) {
527  mIter mk = mj + n;
528  mIter mkjp = mk + j;
529  mIter mjk = mj + j;
530  for (k=j+1;k<=n;k++) {
531  s11 = - (*mjk);
532  s12 = - (*mkjp);
533  if (j!=1) {
534  mIter mik = rhs.Array() + k - 1;
535  mIter mijp = rhs.Array() + j;
536  mIter mki = mk;
537  mIter mji = mj;
538  for (unsigned int i=1;i<j;i++) {
539  s11 += (*mik) * (*(mji++));
540  s12 += (*mijp) * (*(mki++));
541  mik += n;
542  mijp += n;
543  }
544  }
545  // cast to avoid warnings from double to float conversions
546  *(mjk++) = static_cast<T>( - s11 * (*mjj) );
547  *(mkjp) = static_cast<T> ( -(((*(mjj+1)))*((*(mkjp-1)))+(s12)) );
548  mk += n;
549  mkjp += n;
550  }
551  }
552  mj += n;
553  mjj += (n+1);
554  }
555  if (nxch%2==1) det = -det;
556  if (jfail !=jrange) det = 0.0;
557  ir[n] = nxch;
558  return 0;
559 }
560 
561 
562 
563  /**
564  Inversion for General square matrices.
565  Code from dfinv routine from CERNLIB
566  Assumed first the LU decomposition via DfactMatrix function
567 
568  taken from CLHEP : L. Moneta May 2006
569  */
570 
571 template <unsigned int idim, unsigned int n>
572 template<class T>
574 
575 
576  typedef T* mIter;
577 
578  typedef T value_type;
579  //typedef double value_type; // for making inversions in double's
580 
581 
582  if (idim != n) return -1;
583 
584 
585  value_type s31, s32;
586  value_type s33, s34;
587 
588  mIter m11 = rhs.Array();
589  mIter m12 = m11 + 1;
590  mIter m21 = m11 + n;
591  mIter m22 = m12 + n;
592  *m21 = -(*m22) * (*m11) * (*m21);
593  *m12 = -(*m12);
594  if (n>2) {
595  mIter mi = rhs.Array() + 2 * n;
596  mIter mii= rhs.Array() + 2 * n + 2;
597  mIter mimim = rhs.Array() + n + 1;
598  for (unsigned int i=3;i<=n;i++) {
599  unsigned int im2 = i - 2;
600  mIter mj = rhs.Array();
601  mIter mji = mj + i - 1;
602  mIter mij = mi;
603  for (unsigned int j=1;j<=im2;j++) {
604  s31 = 0.0;
605  s32 = *mji;
606  mIter mkj = mj + j - 1;
607  mIter mik = mi + j - 1;
608  mIter mjkp = mj + j;
609  mIter mkpi = mj + n + i - 1;
610  for (unsigned int k=j;k<=im2;k++) {
611  s31 += (*mkj) * (*(mik++));
612  s32 += (*(mjkp++)) * (*mkpi);
613  mkj += n;
614  mkpi += n;
615  }
616  *mij = static_cast<T>( -(*mii) * (((*(mij-n)))*( (*(mii-1)))+(s31)) );
617  *mji = static_cast<T> ( -s32 );
618  mj += n;
619  mji += n;
620  mij++;
621  }
622  *(mii-1) = -(*mii) * (*mimim) * (*(mii-1));
623  *(mimim+1) = -(*(mimim+1));
624  mi += n;
625  mimim += (n+1);
626  mii += (n+1);
627  }
628  }
629  mIter mi = rhs.Array();
630  mIter mii = rhs.Array();
631  for (unsigned int i=1;i<n;i++) {
632  unsigned int ni = n - i;
633  mIter mij = mi;
634  //int j;
635  for (unsigned j=1; j<=i;j++) {
636  s33 = *mij;
637  mIter mikj = mi + n + j - 1;
638  mIter miik = mii + 1;
639  mIter min_end = mi + n;
640  for (;miik<min_end;) {
641  s33 += (*mikj) * (*(miik++));
642  mikj += n;
643  }
644  *(mij++) = static_cast<T> ( s33 );
645  }
646  for (unsigned j=1;j<=ni;j++) {
647  s34 = 0.0;
648  mIter miik = mii + j;
649  mIter mikij = mii + j * n + j;
650  for (unsigned int k=j;k<=ni;k++) {
651  s34 += *mikij * (*(miik++));
652  mikij += n;
653  }
654  *(mii+j) = s34;
655  }
656  mi += n;
657  mii += (n+1);
658  }
659  unsigned int nxch = ir[n];
660  if (nxch==0) return 0;
661  for (unsigned int mm=1;mm<=nxch;mm++) {
662  unsigned int k = nxch - mm + 1;
663  int ij = ir[k];
664  int i = ij >> 12;
665  int j = ij%4096;
666  mIter mki = rhs.Array() + i - 1;
667  mIter mkj = rhs.Array() + j - 1;
668  for (k=1; k<=n;k++) {
669  // 2/24/05 David Sachs fix of improper swap bug that was present
670  // for many years:
671  T ti = *mki; // 2/24/05
672  *mki = *mkj;
673  *mkj = ti; // 2/24/05
674  mki += n;
675  mkj += n;
676  }
677  }
678  return 0;
679 }
680 
681 
682  } // end namespace Math
683 } // end namespace ROOT
684 
685 #endif
static int DfinvMatrix(MatRepStd< T, idim, n > &rhs, unsigned int *work)
LU inversion of general square matrices.
Namespace for new ROOT classes and functions.
Definition: ROOT.py:1
double T(double x)
Definition: ChebyshevPol.h:34
MatRepSym Matrix storage representation for a symmetric matrix of dimension NxN This class is a templ...
Definition: HelperOps.h:35
static int DfactMatrix(MatRepStd< T, idim, n > &rhs, T &det, unsigned int *work)
LU Factorization method for inversion of general square matrices (see implementation in Math/MatrixIn...
Double_t x[n]
Definition: legend1.C:17
static void InvertBunchKaufman(MatRepSym< T, idim > &rhs, int &ifail)
Bunch-Kaufman method for inversion of symmetric matrices.
static Vc_ALWAYS_INLINE Vector< T > abs(const Vector< T > &x)
Definition: vector.h:450
iterator begin()
STL iterator interface.
Definition: SVector.icc:543
Expression wrapper class for Matrix objects.
Definition: Expression.h:134
TLine * l
Definition: textangle.C:4
REAL epsilon
Definition: triangle.c:617
Namespace for new Math classes and functions.
int * iq
Definition: THbookFile.cxx:86
float * q
Definition: THbookFile.cxx:87
const Int_t n
Definition: legend1.C:16
SVector: a generic fixed size Vector class.