source: Sophya/trunk/Poubelle/archTOI.old/gaussj.c@ 4010

Last change on this file since 4010 was 577, checked in by ansari, 26 years ago

SST
File size: 1.4 KB
Line 
1#include <math.h>
2#define NRANSI
3#include "nrutil.h"
4#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
5
6void gaussj(double **a, int n, double **b, int m)
7{
8 int *indxc,*indxr,*ipiv;
9 int i,icol,irow,j,k,l,ll;
10 double big,dum,pivinv,temp;
11
12 indxc=ivector(1,n);
13 indxr=ivector(1,n);
14 ipiv=ivector(1,n);
15 for (j=1;j<=n;j++) ipiv[j]=0;
16 for (i=1;i<=n;i++) {
17 big=0.0;
18 for (j=1;j<=n;j++)
19 if (ipiv[j] != 1)
20 for (k=1;k<=n;k++) {
21 if (ipiv[k] == 0) {
22 if (fabs(a[j][k]) >= big) {
23 big=fabs(a[j][k]);
24 irow=j;
25 icol=k;
26 }
27 } else if (ipiv[k] > 1) nrerror("gaussj: Singular Matrix-1");
28 }
29 ++(ipiv[icol]);
30 if (irow != icol) {
31 for (l=1;l<=n;l++) SWAP(a[irow][l],a[icol][l])
32 for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l])
33 }
34 indxr[i]=irow;
35 indxc[i]=icol;
36 if (a[icol][icol] == 0.0) nrerror("gaussj: Singular Matrix-2");
37 pivinv=1.0/a[icol][icol];
38 a[icol][icol]=1.0;
39 for (l=1;l<=n;l++) a[icol][l] *= pivinv;
40 for (l=1;l<=m;l++) b[icol][l] *= pivinv;
41 for (ll=1;ll<=n;ll++)
42 if (ll != icol) {
43 dum=a[ll][icol];
44 a[ll][icol]=0.0;
45 for (l=1;l<=n;l++) a[ll][l] -= a[icol][l]*dum;
46 for (l=1;l<=m;l++) b[ll][l] -= b[icol][l]*dum;
47 }
48 }
49 for (l=n;l>=1;l--) {
50 if (indxr[l] != indxc[l])
51 for (k=1;k<=n;k++)
52 SWAP(a[k][indxr[l]],a[k][indxc[l]]);
53 }
54 free_ivector(ipiv,1,n);
55 free_ivector(indxr,1,n);
56 free_ivector(indxc,1,n);
57}
58#undef SWAP
59#undef NRANSI
Note: See TracBrowser for help on using the repository browser.