I am new to using NR. I am trying to solve a 100x100 system of equations using ludcmp in C. I use a=matrix(1,rows, 1, cols) and b=vector(1,rows) and then use ludcmp and lubksb. I have checked it for small systems (3 variables) and it works OK. But for the 100x100 system, it gives me strange answer with negative values and large numbers. What is wrong? Can anybody help? (Please send me email).

Thanks,
bladeru99er



#include <stdio.h>
#include <malloc.h>
#include <stdlib.h>
#include "nr.h"
#include "nrutil.h"

#define dataType float

int input(FILE *fA, FILE *fB, dataType** A, dataType* B, int r, int c) {
int rr, cc;
for (rr=1; rr<=r; rr++) {
for (cc=1; cc<=c; cc++) {
//fprintf(stderr, "reading (%d,%d)\t", rr,cc);
fscanf(fA, "%f", &((A[rr])[cc]));
}
fscanf(fB, "%f", &B[rr]);
//fprintf(stderr, "%f \t", B[rr]);
}

return 0;
}

int printVec(FILE *fVec, dataType* vec, int r, int c) {
int rr, cc;
for (rr=1; rr<=r; rr++) {
fprintf(fVec, "%f\n", vec[rr]);
}
return 0;
}



int main(int argc, char* argv[]) {
dataType **A, *B, d;
FILE *fA, *fX, *fB;
int *indx, i, j;
int rows, cols;

if (argc < 6) {
printf("usage: lusolve <vector A filename> <vector X output filename> <vector B filename> <# array rows> <# array cols>\n");
exit(0);
}

rows = atoi(argv[4]);
cols = atoi(argv[5]);
fprintf(stderr, "%dx%d\n", rows, cols);
if ((rows <=0) || (cols <=0)) {
printf("Either the number of rows or columns provided (%s, %s) is incorrect.\n", argv[4], argv[5]);
exit(0);
}

if (rows != cols) {fprintf(stderr, "No. equations not same as no. variables, can't solve, quitting now.\n"); exit(0);}

A = matrix(1, rows, 1, cols);
B = vector(1, rows);
indx = ivector(1, rows);

if ((A == NULL) || (B == NULL) || (indx == NULL)) {
printf("Failed to allocate required memory, something is terribly wrong.\n");
exit(0);
}

fA = fopen(argv[1], "rt");
fX = fopen(argv[2], "wt");
fB = fopen(argv[3], "rt");
if ((fA == NULL) || (fX == NULL) || (fB == NULL)) {
printf("Failed to open required files, something is terribly wrong.\n");
exit(0);
}

input(fA, fB, A, B, rows, cols);
ludcmp(A, rows, indx, &d);
lubksb(A, rows, indx, B);

printVec(stdout, B, rows, 1);

fclose(fA);fclose(fX);fclose(fB);
return 0;
}

I am using ludcmp and lubksb to solve a 87X87 system and I get the same problem. I've checked the solution vector and it seems to be memory adresses. I am using the way to solve the linear system as described on page 51 of the book. If anybody knows what the problem is please post.

Dear bladeru99er,

You are using the float datatype to solve a 100x100 dense
linear algebraic equation system. This is not a good idea.
The roundoff errors can (except lucky cases) destroy your solution even if the algorithm is not failing. Try double, but for
"bad" matrices couldn't be enough.
Consider that there are no "recipes" for ill conditioned matrices.

Regards

Max