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