LU Decomposition


Rariusz
08-23-2009, 05:03 AM
Hallo

I have the code:

#include <stdio.h>
#include <stdlib.h>
#define NN 4
#include "nrutil.h"
#include "nr.h"

int main(int argc, char *argv[])
{
float **M1=matrix(1,NN,1,NN);
int *M2=ivector(1,NN);
float A;
int j,i;

for(i=1;i<NN;i++){
for(j=1;j< NN;j++){
printf("Wprowadz element [%d %d]macierzy ",i,j);
scanf("%f",&(M1[i][j]));
}
}
for(i=1;i<NN;i++){
for(j=1;j<NN;j++){
printf("%f\t",M1[i][j]);
}
printf("\n");
}


ludcmp(M1,NN,M2,&A);
printf("\n");


for(i=1;i<NN;i++){
for(j=1;j<NN;j++){
printf("%f\t",M1[i][j]);
}
printf("\n");
}






free_matrix(M1,1,NN,1,NN);
free_ivector(M2,1,NN);


system("PAUSE");
return 0;
}

matrix:

2.000000 34.000000 3.000000
23.000000 44.000000 55.000000
23.000000 98.000000 34.000000

then:

ludcmp(M1,NN,M2,&A);

I have:

462002186468332547000000000000.000000 0.00000 697762776322846200000000000.00000
0.000000 44.00000 54.999653
0.000000 0.772727 -39.4999763


Where is the mistake?

=======================
Rariusz PL

davekw7x
08-23-2009, 09:12 AM
Where is the mistake?
Here's what I see: (I added a couple of comments)

#define NN 4
.
.
.
/* NN = 4, so you are reading a 3x3 matrix */
for (i = 1; i < NN; i++) {
for (j = 1; j < NN; j++) {
printf("Wprowadz element [%d %d]macierzy ", i, j);
scanf("%f", &(M1[i][j]));
}
}
.
.
.
ludcmp(M1, NN, M2, &A);/* NN = 4 so you are telling ludcmp it is a 4x4 matrix. */


You are reading a 3x3 matrix but you tell ludcmp that it's a 4x4.

Quickest fix with minimum changes: You can change the call to ludcmp to

ludcmp(M1, NN-1, M2, &A);


But I would probably do something like the following to make the intent more clear. (At least it's more clear to me)

.
.
.
int j, i;

/* Number of equations. Must be less than or equal to NN */
int num_equations = 3; /* for a 3x3 Matrix */

for (i = 1; i <= num_equations; i++) {
for (j = 1; j <= num_equations; j++) {
printf("Wprowadz element [%d %d]macierzy ", i, j);
scanf("%f", &(M1[i][j]));
printf("\n");
}
}

printf("\nInitial M Matrix:\n");
for (i = 1; i <= num_equations; i++) {
for (j = 1; j <= num_equations; j++) {
printf("%f\t", M1[i][j]);
}
printf("\n");
}

ludcmp(M1, num_equations, M2, &A);

printf("\nM Matrix after ludcmp:\n");
for (i = 1; i <= num_equations; i++) {
for (j = 1; j <= num_equations; j++) {
printf("%f\t", M1[i][j]);
}
printf("\n");
}
.
.
.


Here is my output with the matrix you used:

Initial M Matrix:
2.000000 34.000000 3.000000
23.000000 44.000000 55.000000
23.000000 98.000000 34.000000

M Matrix after ludcmp:
23.000000 44.000000 55.000000
0.086957 30.173912 -1.782609
1.000000 1.789625 -17.809797


Regards,

Dave

Rariusz
08-28-2009, 04:36 AM
Hallo

Thanks. All works:-D