Hello,
I use the LU decomposition in fortran 77 to inverse matrix (typically n<100, declared as real(8)), it works really fine for almost all matrix but it fails in some case (when significant number is too high ). How can I improved the routine to work fine in all cases?Is Cholesky decomposition more stable? How to use it to inverse matrix?
Thank you.
kelment

Have you considered using singular value decomposition instead of LU decomposition?

Yes i've tried but the result is the same.... The svd decomposition work well (I've tested it by verifying that M=U*W*Vt, where U W and V are output of svdcmp) but the inversion doesn't work...
So if someone is able to inverse this 5*5 matrix (nevermind the method) :

54987.82 27668.27 -61321.53 -137689.0 148731.3
-121870.2 -61321.53 135907.7 305161.9 -329635.0
-273642.7 -137689.0 305161.9 685198.8 -740149.8
295588.0 148731.2 -329635.0 -740149.8 799507.6

Please tell me how you've done it !!
Thanks,
kelment

A bit disconcerted by a grumble where the example supplied is not even square? Is there another row somewhere?

Hi,

1) Cholesky decomposition is more stable but it can be used only for symmetric positive-defined matrices.

2) In any case LU and Cholesky is used for inverse of nxn square not-singular matrices.

3) If you have a nxm non-square rectangular matrix (n>m) you should calculate a pseudoinverse matrix for example by SVD or solve a linear system in terms of LMS via QR-decomposition or SVD.

Evgeny

Indeed, there's a row missing! (Is there numerical routine for copy/paste, i've got some difficulty....). Of course the 'real' matrix is square... I don't know what to do... maybe it's a numerical problem linked to my compiler or something else... no idea about this...
I've got this problem in 1-D simulation, when I applied this on 3-D, everything works fine with much larger matrix... So, let's go for 3-D !!!!

Hi friends,
I am new to NRC and I am having trouble with the ludcmp() function.... I suppose there is somethign wrong in my usage rather than the function itself.... The code snippet is given below..

*******************************************
main(){
float **a,**y,d,*col;
int i,j,*indx;
a=matrix(1,ROWS,1,COLS);
y=matrix(1,ROWS,1,COLS);
col=vector(1,ROWS);
indx=ivector(1,ROWS);


/*here i fill values into the input matrix*/
for(i=1;i<=ROWS;i++){
for(j=1;j<=COLS;j++){
*(*(a+i)+j)=3*i+j;
printf("%2.3f\t",a[i][j]);
}
printf("\n");
}


/*here starts the routine to find the inverse*/
ludcmp(a,ROWS,indx,&d);
for(j=1;j<=ROWS;j++){
for(i=1;i<=ROWS;i++) col[i]=0.0;
col[j]=1.0;
lubksb(a,ROWS,indx,col);
for(i=1;i<=ROWS;i++) y[i][j]=col[i];
}


/*here i print the result*/
for(i=1;i<=ROWS;i++){
for(j=1;j<=COLS;j++){
printf("%2.3f\t",y[i][j]);
}
printf("\n");
}
*******************************************

Is there anything really wrong?
I dont understand why I am getting garbage values in the output.By the way, I run this in Microsoft VisualStudio.
Thanks in advance,
CW

ROWS and COLS are #define'd to the same value 2.
Sorry for missing this detail in the previous post..
Once again,
Thanking you in advance,
CW.