Logo ROOT   6.18/05
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
20namespace 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
38template <unsigned int idim, unsigned int N>
39template<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
445template <unsigned int idim, unsigned int n>
446template<class T>
447int 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
575template <unsigned int idim, unsigned int n>
576template<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
float * q
Definition: THbookFile.cxx:87
int * iq
Definition: THbookFile.cxx:86
static void InvertBunchKaufman(MatRepSym< T, idim > &rhs, int &ifail)
Bunch-Kaufman method for inversion of symmetric matrices.
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 int DfinvMatrix(MatRepStd< T, idim, n > &rhs, unsigned int *work)
LU inversion of general square matrices.
Expression wrapper class for Matrix objects.
MatRepSym Matrix storage representation for a symmetric matrix of dimension NxN This class is a templ...
SVector: a generic fixed size Vector class.
Definition: SVector.h:75
iterator begin()
STL iterator interface.
Definition: SVector.icc:553
const Double_t sigma
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
double T(double x)
Definition: ChebyshevPol.h:34
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
static constexpr double s
static constexpr double mm
auto * l
Definition: textangle.C:4
REAL epsilon
Definition: triangle.c:617