Hi,

(1)

I am using gaussj.c routine (chapter 2) to solve N simultaneous linear equations:

N=5;

double **coef_matrix, **sol_col;

coef_matrix=dmatrix(1, N, 1, N);
sol_col=dmatrix(1, N, 1, 1);

/* ************ creating coef_matrix[ ] [ ] matrix */
/* coef_matrix[ ] [ ] is two dimensional array of numbers */
/* it contains N rows and N columns */


for( j=-10 ; j<11 ; j++)
{
for(r=1 ; r<=N ; r++)
{
for(c=1 ; c<=N ; c++)
{
coef_matrix[r][c]= coef_matrix[r][c] + pow( j,(r+c-2)*1.0 );
}
}
}

/* ********* creating sol_col[ ][ ] matrix */
/* sol_col[ ] [ ] is a column matrix having N rows and 1 column */

for(n= 10; n< 300; n++)
{
for(r=1 ; r<=N ; r++)
{
c=1;
for(j=-10; j<11 ; j++)
{
sol_col[r][c]= sol_col[r][c] + (my_array [n+j] *pow(j,(r-1)*1.0));
}
}

gaussj(coef_matrix, N, sol_col, 1); /* call to NR recipe */

/* solution is sol_matrix[ ] [ ] */
}

:
:
:
:
free_dmatrix(coef_matrix, 1, N, 1, N);
free_dmatrix(sol_col, 1, N, 1, 1);



(a)
I do not get expected solution. I use sol_col[ ] [ ] further to calculate some other quantity; the values of this quantity are extremely large as opposed to expected small values.

(b) Is there any mistake in the above segment of the code?

Regards
smp

I am using gaussj.c routine
.
.
.
double **coef_matrix, **sol_col;

coef_matrix=dmatrix(1, N, 1, N);
sol_col=dmatrix(1, N, 1, 1);
.
.
.
.
.
gaussj(coef_matrix, N, sol_col, 1);
.
.
.



Assuming that you are using unmodified code from the book or the NR distribution, here is the signature for gaussj() in gaussj.c:

void gaussj(float **a, int n, float **b, int m)


To use this, your matrix variables would be floats, not doubles, and you would use matrix(), not dmatrix(), to allocate the storage.

Regards,

Dave