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