# dmatrix() function from NR

• 12-24-2009
swapnilp
dmatrix() function from NR
Hi,

Following is the segment of my main C code:
poly_fit () is my own function.
I am using Numerical Recipe subroutine gaussj.c and function dmatrix() from nrutil.c

Basically I want to solve linear equations.
I am fitting a 4th degree polynomial fucntion to some data points.

I get matrix equation: A*B=C
where A= coef_matrix in my code
B= set of 5 coefficients to be found e.g. {a0,a1,a2,a3,a4} = sol_col in my code
C= data points obtained from array old_est[] as argument to my function = sol_col in my code
(sol_col is replaced by the solution from gaussj.c )

N = 5 (N is global variable, there are other glbal bariables also in this segment).

My problems:
1. Is use of dmatrix() okay?
2. I get the output which contains very large unexpected numbers alternately i.e. 1 st no ok, 2nd no bad, 3 rd no ok,....so on.

I AM SURE THAT PROBLEM IS WITHIN MY poly_fit() FUNCTION REGARDING dmatrix() and gaussj.c

Any help will be appreciated.

Code:

```void poly_fit(double ip[], double old_est[], double err[])  {         double **coef_matrix, **sol_col;                      coef_matrix=dmatrix(1, N, 1, N);         sol_col=dmatrix(1, N, 1, 1);           int i,j,n,r,c;         nst_neg=-1*nst; //-------------resetting to zero ------------------------     for(r=1;r<=N;r++)         {           for(c=1;c<=N;c++)                 {                   coef_matrix[r][c]=0.0;                        }         } //-------------------------------------------------------------     for(j=nst_neg;j<nst+1;j++)  // j varies from -10 to +10       {         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 );                 }           }        }     for(n=nst+1;n<nend+2;n++)  // n varies from 10 to 290       {             for(r=1;r<=N;r++)                 {                     c=1;                     sol_col[r][c]=0;                 }                        for(j=nst_neg;j<nst+1;j++)      // j varies from -10 to +10             {                       for(r=1;r<=N;r++)                     {                     c=1;                     sol_col[r][c]= sol_col[r][c] + (old_est[n+j]*pow(j,(r-1)*1.0));                                           }                }                     //---------------------------------------------------------                       gaussj(coef_matrix, N, sol_col, 1);       /* numerical recipe for solving simultaneous equations */                 //-----------------------------------------------------------             err[n] =sol_col[1][1]- ip[n] ;                                  }             free_dmatrix(coef_matrix, 1, N, 1, N);                      free_dmatrix(sol_col, 1, N, 1, 1);                              return ;  }```
thanks