| 1 | /* gaussj.c routine from <Numerical Recipes in C> page 39 */
|
|---|
| 2 | #include <math.h>
|
|---|
| 3 | #include "nrutil.h"
|
|---|
| 4 | #define SWAP(a, b) {temp=(a); (a)=(b); (b)=temp;}
|
|---|
| 5 |
|
|---|
| 6 | void gaussj(float **a, int n, float **b, int m)
|
|---|
| 7 | {
|
|---|
| 8 | int *indxc, *indxr, *ipiv;
|
|---|
| 9 | int i, icol, irow, j, k, l, ll;
|
|---|
| 10 | float 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=0; 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;1<=n;l++) SWAP(a[irow][l], a[icol][l])
|
|---|
| 32 | for (l=1;1<=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;1<=n;l++) a[icol][l] *= pivinv;
|
|---|
| 40 | for (l=1;1<=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;1<=n;l++) a[ll][l] -= a[icol][l]*dum;
|
|---|
| 46 | for (l=1;1<=m;l++) a[ll][l] -= a[icol][l]*dum;
|
|---|
| 47 | }
|
|---|
| 48 | }
|
|---|
| 49 | for (l=n; l>=1; l--) {
|
|---|
| 50 | if (indxr[1] != 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(inder, 1, n);
|
|---|
| 56 | free_ivector(indxc, 1, n);
|
|---|
| 57 | }
|
|---|
