Hello,

i've tried to find a inverse of a 4x4 matrix with the ludcmp and lubksb. I'm new in nr.
In the ludcmp routine i get an index vector with
0,1,3,3. Where is the index for the third row? What could be wrong in the routine? I`ve try to find the error since two days.

I know, that ludcmp changes the third and fourth row in its algorithm, but the vector of the pivoting indices is wrong (it should be 0,1,3,2).
In the inverse matrix the first colum is total wrong (0.8,0,0,0) and the third and fourth columns are interchanged.


Is there anybody who can help me or is there somebody else, who has this problem?

Could it be, that my indices are from 0 to n-1? But I don`t think so, because I had the indices in the code changed in this way.

PS: The arrays of the matrices and vectors are from 0 to n-1

void LUDCMP(double **a,int n,int *index, double *d)
{
int i,imax,j,k;
double big,dum,sum,temp;
double *vv;
vv=vector(n);
*d=1.0;
for(i=0;i<n;i++)
{
big=0.0;
for(j=0;j<n;j++)
if((temp=fabs(a[i][j]))>big) big=temp; //maximales Element in i-ter Spalte
if(big==0.0) printf("Singulaere Matrix");
vv[i]=1.0/big;
} //Skalierungsfaktor
for(j=0;j<n;j++)
{
for(i=0;i<j;i++)
{
sum=a[i][j];
for(k=0;k<i;k++) sum -= a[i][k]*a[k][j]; //Summe=Element a(i,j)-i-te Zeile*j-te Spalte
a[i][j]=sum; //a(i,j) neu (oberes Dreieck)
}
big=0.0;
for(i=j;i<n;i++) //Schleife f&#252;r unteres Dreieck
{
sum=a[i][j];
for(k=0;k<j;k++) sum -=a[i][k]*a[k][j];
a[i][j]=sum;
if((dum=vv[i]*fabs(sum))>=(big)) //Wenn Skalierung*Summe (Zeilenweise) gr&#246;&#223;er als big ist
{
big=dum; //big neu setzen
imax=i; //imax setzen
}
}
if(j!=imax) //wenn j ungleich imax ist
{
for(k=0;k<n;k++)
{
dum=a[imax][k]; //Zeile imax gegen Zeile j tauschen
a[imax][k]=a[j][k];
a[j][k]=dum;
}
*d=-(*d); // Vorzeichen von d &#228;ndern
vv[imax]=vv[j];
}
index[j]=imax; //an j-ter Stelle von index imax setzen
if(a[j][j]==0.0) a[j][j]=zero;
if(j!=(n-1))
{
dum=1.0/(a[j][j]);
for(i=j+1;i<n;i++) a[i][j] *= dum;
}
}
cout<<"LU-Zerlegung"<<endl;
for(i=0;i<n;i++)
cout<<index[i]<<" "<<vv[i]<<" "<<*d<<endl;
for(i=0;i<n;i++)
{
cout<<a[i][0]<<" "<<a[i][1]<<" "<<a[i][2]<<" "<<a[i][3]<<endl;
}
free_vector(vv,n);
}

I've checked your code against mine(I've also adjusted the indices to be zero based) and it looks exactly the same... could you please post the matrix that is giving you the wrong results? It could help me detect a bug in my code as well.
Thanks!
Fernando.