//$Id: gaussj.C,v 1.1 2012/04/18 02:21:52 zjcao Exp $ #ifdef newc #include <iostream> #include <iomanip> #include <fstream> #include <cstdlib> #include <cstring> #include <cmath> using namespace std; #else #include <iostream.h> #include <iomanip.h> #include <fstream.h> #include <stdlib.h> #include <stdio.h> #include <string.h> #include <math.h> #endif /* Linear equation solution by Gauss-Jordan elimination. a[0..n-1][0..n-1] is the input matrix. b[0..n-1] is input containing the right-hand side vectors. On output a is replaced by its matrix inverse, and b is replaced by the corresponding set of solution vectors */ int gaussj(double *a, double *b, int n) { double swap; int *indxc, *indxr, *ipiv; indxc = new int[n]; indxr = new int[n]; ipiv = new int[n]; int i, icol, irow, j, k, l, ll; double big, dum, pivinv, temp; for (j = 0; j < n; j++) ipiv[j] = 0; for (i = 0; i < n; i++) { big = 0.0; for (j = 0; j < n; j++) if (ipiv[j] != 1) for (k = 0; k < n; k++) { if (ipiv[k] == 0) { if (fabs(a[j * n + k]) >= big) { big = fabs(a[j * n + k]); irow = j; icol = k; } } else if (ipiv[k] > 1) { cout << "gaussj: Singular Matrix-1" << endl; for (int ii = 0; ii < n; ii++) { for (int jj = 0; jj < n; jj++) cout << a[ii * n + jj] << " "; cout << endl; } return 1; // error return } } ipiv[icol] = ipiv[icol] + 1; if (irow != icol) { for (l = 0; l < n; l++) { swap = a[irow * n + l]; a[irow * n + l] = a[icol * n + l]; a[icol * n + l] = swap; } swap = b[irow]; b[irow] = b[icol]; b[icol] = swap; } indxr[i] = irow; indxc[i] = icol; if (a[icol * n + icol] == 0.0) { cout << "gaussj: Singular Matrix-2" << endl; for (int ii = 0; ii < n; ii++) { for (int jj = 0; jj < n; jj++) cout << a[ii * n + jj] << " "; cout << endl; } return 1; // error return } pivinv = 1.0 / a[icol * n + icol]; a[icol * n + icol] = 1.0; for (l = 0; l < n; l++) a[icol * n + l] *= pivinv; b[icol] *= pivinv; for (ll = 0; ll < n; ll++) if (ll != icol) { dum = a[ll * n + icol]; a[ll * n + icol] = 0.0; for (l = 0; l < n; l++) a[ll * n + l] -= a[icol * n + l] * dum; b[ll] -= b[icol] * dum; } } for (l = n - 1; l >= 0; l--) { if (indxr[l] != indxc[l]) for (k = 0; k < n; k++) { swap = a[k * n + indxr[l]]; a[k * n + indxr[l]] = a[k * n + indxc[l]]; a[k * n + indxc[l]] = swap; } } delete[] indxc; delete[] indxr; delete[] ipiv; return 0; } // for check usage /* int main() { double *A,*b; A=new double[9]; b=new double[3]; A[0]=0.5; A[1]=1.0/3; A[2]=1; A[3]=1; A[4]=5.0/3; A[5]=3; A[6]=2; A[7]=4.0/3; A[8]=5; b[0]=1; b[1]=3; b[2]=2; cout<<"initial data:"<<endl; for(int i=0;i<3;i++) cout<<A[i*3]<<" "<<A[i*3+1]<<" "<<A[i*3+2]<<" "<<b[i]<<endl; gaussj(A, b, 3); cout<<"final data:"<<endl; for(int i=0;i<3;i++) cout<<A[i*3]<<" "<<A[i*3+1]<<" "<<A[i*3+2]<<" "<<b[i]<<endl; delete[] A; delete[] b; } */