home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / lsqrft15.zip / solve_da.c < prev    next >
Text File  |  1993-07-21  |  2KB  |  92 lines
  1. #include <math.h>
  2. #include "fit.h"
  3.  
  4. void gaussj();
  5.  
  6. /* This program solves the equation alpha*da = beta for da. */
  7. /* alpha is a matrix and da and beta are vectors */
  8. /* covar is used for temporary space */
  9.  
  10. void solve_for_da(alpha, covar, beta, da, mfit)
  11. double **alpha, **covar, *beta, *da;
  12. int mfit;
  13. {
  14. int i, j;
  15. double **mbeta;
  16.  
  17. mbeta = dmatrix(mfit,1);
  18.  
  19. for(i=0;i<mfit;i++){
  20.     mbeta[i][0] = beta[i];
  21.     for(j=0;j<mfit;j++)    covar[i][j] = alpha[i][j];
  22. }
  23. gaussj(covar,mfit,mbeta,1);
  24. for(i=0; i<mfit; i++) da[i] = mbeta[i][0];
  25. free_dmatrix(mbeta,mfit,1);
  26. }
  27.  
  28. #define SWAP(a,b) {double temp=(a);(a)=(b);(b)=temp;}
  29.  
  30. /* gaussj inverts the matrix a and performs */
  31. /* the same operations on the matrix b */
  32. /* This function is very similar to the one in */
  33. /* Numerical Recipes of the same name. */
  34.  
  35. void gaussj(a,n,b,m)
  36. double **a,**b;
  37. int n,m;
  38. {
  39.     int *indxc,*indxr,*ipiv;
  40.     int i,icol,irow,j,k,l,ll;
  41.     double big,dum,pivinv;
  42.  
  43.     indxc=ivector(n);
  44.     indxr=ivector(n);
  45.     ipiv=ivector(n);
  46.     for (j=0;j<n;j++) ipiv[j]=0;
  47.     for (i=0;i<n;i++) {
  48.         big=0.0;
  49.         for (j=0;j<n;j++)
  50.             if (ipiv[j] != 1)
  51.                 for (k=0;k<n;k++) {
  52.                     if (ipiv[k] == 0) {
  53.                         if (fabs(a[j][k]) >= big) {
  54.                             big=fabs(a[j][k]);
  55.                             irow=j;
  56.                             icol=k;
  57.                         }
  58.                     } else if (ipiv[k] > 1) myerror("GAUSSJ: Singular Matrix-1");
  59.                 }
  60.         ++(ipiv[icol]);
  61.         if (irow != icol) {
  62.             for (l=0;l<n;l++) SWAP(a[irow][l],a[icol][l])
  63.             for (l=0;l<m;l++) SWAP(b[irow][l],b[icol][l])
  64.         }
  65.         indxr[i]=irow;
  66.         indxc[i]=icol;
  67.         if (a[icol][icol] == 0.0) myerror("GAUSSJ: Singular Matrix-2");
  68.         pivinv=1.0/a[icol][icol];
  69.         a[icol][icol]=1.0;
  70.         for (l=0;l<n;l++) a[icol][l] *= pivinv;
  71.         for (l=0;l<m;l++) b[icol][l] *= pivinv;
  72.         for (ll=0;ll<n;ll++)
  73.             if (ll != icol) {
  74.                 dum=a[ll][icol];
  75.                 a[ll][icol]=0.0;
  76.                 for (l=0;l<n;l++) a[ll][l] -= a[icol][l]*dum;
  77.                 for (l=0;l<m;l++) b[ll][l] -= b[icol][l]*dum;
  78.             }
  79.     }
  80.     for (l=n-1;l>=0;l--) {
  81.         if (indxr[l] != indxc[l])
  82.             for (k=0;k<n;k++)
  83.                 SWAP(a[k][indxr[l]],a[k][indxc[l]]);
  84.     }
  85.     free(ipiv);
  86.     free(indxr);
  87.     free(indxc);
  88. }
  89.  
  90. #undef SWAP
  91.  
  92.