Hi, Your 4x4 matrix is near singular, its condition number is > 1e+16 ! Actually it is so badly conditioned that the matrix package of root version <=3.10.02 can not properly detect it and claims that the determinant = -157736 I have attached your macro modified/extended for the new matrix package in root 4.00.00: First you see me check the determinant. The "standard" Determinant, Inversion is delivered through a LU decomposition . The algorithm will get a zero on its diagonal and therefore think it is singular . Therefore the deteminant returned is zero . If I now would have performed an inversion, an invalid matrix would have been returned (check through A.Isvalid()) . Seeing that it is a difficult maybe singular case we proceed with a Singular Value Decomposition (SVD) , print the singular values and calculate the condition number . Even the inverse using svd is not able to get a matrix to so that A * A^-1 = E I wonder why you need the inverse ?? Are you trying to solve for x, A x = b . You do not want to do x = A^-1 b . It is not precise enough . Better inparticular your case is to do : TDecompSVD svd(A); svd.SetTol(1.0e-20); // explained in the attachment TVectorD x = b; svd.Solve(x); And now you can check A x, playing with the tolerance . Eddy --- Rene Brun <Rene.Brun@cern.ch> wrote: > Eddy, > > This code seems to run correctly with 3.10/02. > It crashes with 4.00/01 (CVS) > > Rene > ATTACHMENT part 2 message/rfc822 > Date: Wed, 28 Jan 2004 13:15:20 +0000 > From: Cozmin Timis <C.Timis@surrey.ac.uk> > To: roottalk@pcroot.cern.ch > Subject: [ROOT] TMatrix inversion results > > Hello, > > It seems like I have similar TMatrix inversion > problems as in: > http://root.cern.ch/root/roottalk/roottalk03/2876.html > > That was fixed for 3x3 matrix in the version I use > (3.05/05) but I have > problems in 4x4 case . > > > Here's the result of the macro matrix.C that runs > at the Root prompt. > First part is the 3x3 matrix used by Chris Jillings > but I'm puzzled by > the 4x4 case. > > Do I need a newer version than 3.05/05 ? > > Regards, > Cozmin > > > $ root > the current keyboard layout is 437 > ******************************************* > * * > * W E L C O M E to R O O T * > * * > * Version 3.05/05 8 May 2003 * > * * > * You are welcome to visit our Web site * > * http://root.cern.ch * > * * > ******************************************* > > Compiled for win32. > > CINT/ROOT C/C++ Interpreter version 5.15.86, Apr 29 > 2003 > Type ? for help. Commands must be C++ statements. > Enclose multiple statements between { }. > root [0] .x matrix.C > > Matrix 3x3 is as follows > > | 0 | 1 | 2 | > ------------------------------------------------------------------ > 0 |-8.761e-005 7.804e-005 -1.808e-005 > 1 | -0.0003604 0.0001103 2.77e-005 > 2 | 0.009239 -0.01581 -0.001075 > > > Matrix 3x3 is as follows > > | 0 | 1 | 2 | > ------------------------------------------------------------------ > 0 | -2601 -3009 -33.82 > 1 | 1069 -2126 -72.79 > 2 |-3.809e+004 5405 -150.3 > > > Matrix 3x3 is as follows > > | 0 | 1 | 2 | > ------------------------------------------------------------------ > 0 | 1 -1.388e-017 -8.674e-019 > 1 | 0 1 1.735e-018 > 2 | 7.105e-015 4.441e-015 1 > > > Matrix 3x3 is as follows > > | 0 | 1 | 2 | > ------------------------------------------------------------------ > 0 | 1 0 1.388e-017 > 1 | 0 1 0 > 2 | 0 0 1 > > > Matrix 4x4 is as follows > > | 0 | 1 | 2 | 3 > | > ------------------------------------------------------------------ > 0 |-4.288e+008 4.288e+008 4.288e+008 > -4.288e+008 > 1 | 2.115e+012 -2.115e+012 -2.115e+012 > 2.115e+012 > 2 |-5.269e+008 5.269e+008 5.269e+008 > -5.269e+008 > 3 | 2.513e+012 -2.513e+012 -2.513e+012 > 2.513e+012 > > > Matrix 4x4 is as follows > > | 0 | 1 | 2 | 3 > | > ------------------------------------------------------------------ > 0 | 2.863e+007 5351 1.756e+005 > 419 > 1 | 2.56e+005 506 2.627e+007 > 5125 > 2 | 2.863e+007 5351 -2.627e+007 > -5125 > 3 | 2.56e+005 506 -1.756e+005 > -419 > > > Matrix 4x4 is as follows > > | 0 | 1 | 2 | 3 > | > ------------------------------------------------------------------ > 0 | -0.2031 -0.0001831 -0.0625 > 0.0001831 > 1 | -2944 1.25 -2624 > 0 > 2 | -0.7969 -0.0001221 1.344 > 0.0004883 > 3 | 1152 1 -2304 > 0.625 > > > Matrix 4x4 is as follows > > | 0 | 1 | 2 | 3 > | > ------------------------------------------------------------------ > 0 | 0 0 0 > 0 > 1 | 2 0 -2 > 0 > 2 | -2 0 2 > 0 > 3 | 0 0 0 > 1 > > > > here's the macro "matrix.C". > > > { > TMatrixD m(3,3); > TMatrixD n(3,3); > > m[0][0] = -2600.61; > m[0][1] = -3009.4; > m[0][2] = -33.82; > > m[1][0] = 1068.99; > m[1][1] = -2125.77; > m[1][2] = -72.79; > > m[2][0] = -38087.3; > m[2][1] = 5404.86; > m[2][2] = -150.3; > > n = m; > Double_t det; > m.Invert(&det); > > m.Print(); > n.Print(); > (m*n).Print(); > (n*m).Print(); > > > // 4x4 matrix > Double_t x1=5351; > Double_t x2=506; > Double_t y1=419; > Double_t y2=5125; > > TMatrixD A (4,4); > TMatrixD B (4,4); > TMatrixD C (4,4); > > A[0][0]=x1*x1; > A[0][1]=x1; > A[0][2]=y1*y1; > A[0][3]=y1; > > A[1][0]=x2*x2; > A[1][1]=x2; > A[1][2]=y2*y2; > A[1][3]=y2; > > A[2][0]=x1*x1; > A[2][1]=x1; > A[2][2]=-y2*y2; > A[2][3]=-y2; > > A[3][0]=x2*x2; > A[3][1]=x2; > A[3][2]=-y1*y1; > A[3][3]=-y1; > > > === message truncated ===
This archive was generated by hypermail 2b29 : Sun Jan 02 2005 - 05:50:05 MET