Gauss_Jordan Method


kutta
11-21-2008, 02:25 AM
Hello comrade,
Can you please indicate a better way to initialize this coded stuf to enable get output


#include <math.h>

#include "nr3.h"
using namespace std;

int a,b;

class gaussj{

int a,b;
public :
gaussj(MatDoub_IO& a ,MatDoub_IO& b);
gaussj(MatDoub_IO& a);
~gaussj();


};

void gaussj(MatDoub_IO & a, MatDoub_IO& b)
{
Int i,icol,irow,j,k,l,ll,n=a.nrows(),m=b.ncols();
Doub big,dum,pivinv;
VecInt indxc(n),indxr(n),ipiv(n);
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 (abs(a[j][k]) >= big) {
big=abs(a[j][k]);
irow=j;
icol=k;
}
}
}
++(ipiv[icol]);
if (irow != icol) {
for (l=0;l<n;l++) SWAP(a[irow][l],a[icol][l]);
for (l=0;l<m;l++) SWAP(b[irow][l],b[icol][l]);
}
indxr[i]=irow;
indxc[i]=icol;
if (a[icol][icol] == 0.0) throw("gaussj: Singular Matrix");
pivinv=1.0/a[icol][icol];
a[icol][icol]=1.0;
for (l=0;l<n;l++) a[icol][l] *= pivinv;
for (l=0;l<m;l++) b[icol][l] *= pivinv;
for (ll=0;ll<n;ll++)
if (ll != icol) {
dum=a[ll][icol];
a[ll][icol]=0.0;
for (l=0;l<n;l++) a[ll][l] -= a[icol][l]*dum;
for (l=0;l<m;l++) b[ll][l] -= b[icol][l]*dum;
}
}
for (l=n-1;l>=0;l--) {
if (indxr[l] != indxc[l])
for (k=0;k<n;k++)
SWAP(a[k][indxr[l]],a[k][indxc[l]]);

}

}

void gaussj(MatDoub_IO &a)
{
MatDoub b(a.nrows(),0);
gaussj(a,b);
}

int main(){


return 0;
}


Else
A universal method Of initialization that wil solve all the routines of NR.
ie)
The input of any routine requires MatVecDoub_I, MatVecDoub_IO, MatVecInt_I, MatVecInt_IO etc
ie)Master dependancy tool like only for feeding & getting output for C++.
Advance
Thanks
As
Kutta

davekw7x
11-21-2008, 09:43 AM
...indicate a better way
One possibility for finding the inverse: Read the size of the matrix and the matrix elements from a file:


//
// Read matrix values from a file as follows:
//
// The first line in the file is a comment line.
// It is printed out exactly as given.
//
// The second line contains an integer: the size of the
// matrix.
//
// This is followed by the elements of the matrix
//
// The name of the file can be given on the command
// line that invokes this program.
//
// If no command line argument is given, the name
// of the file is "matrxn.dat"
//
//
// davekw7x
//
#include "../code/nr3.h"
#include "../code/gaussj.h"

int main(int argc, char **argv)
{
/*
* Find the inverse of an "nxn" matrix
*/
int n;
int i, j, k;
string str;
char *inname = "matrxn.dat";
if (argc > 1) {
inname = argv[1];
}

ifstream fp(inname);

if (!fp) {
cout << "Can't open file " << inname << " for reading." << endl;
exit(EXIT_FAILURE);
}
cout << fixed << setprecision(6);
getline(fp,str); // The comment line
cout << str << endl << endl;
fp >> n;
if (!fp) {
cerr << "Error reading the size of the matrix" << endl;
exit(EXIT_FAILURE);
}
cout << "Matrix size is " << n << " x " << n << endl;
getline(fp, str); // eat the rest of the matrix size line

MatDoub a(n, n);
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
fp >> a[i][j];
if (!fp) {
cerr << "Error reading a[" << i << "][" << j << "]"
<< endl;
exit(EXIT_FAILURE);
}
}
}

cout << scientific << setprecision(6);
cout << endl << "The matrix [a]: " << endl;
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
cout << setw(14) << a[i][j];
}
cout << endl;
}
cout << endl;

MatDoub ai(a); // Copy of a will hold the inverse

// Now find the inverse
gaussj(ai);

cout << "The calculated inverse of the matrix [a]:" << endl;
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
cout << setw(14) << ai[i][j];
}
cout << endl;
}
cout << endl;

//Multiply [a] * [a-inverse]
MatDoub check(n, n);
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
check[i][j] = 0.0;
for (k = 0; k < n; k++) {
check[i][j] += (a[i][k] * ai[k][j]);
}
}
}
cout << "The matrix [a] times the calculated inverse:"
<< endl;
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
cout << setw(14) << check[i][j];
}
cout << endl;
}
cout << endl;

fp.close();

return EXIT_SUCCESS;
}


Here's matrxn.dat:

Sample problem number 1
3
1.0 0.0 0.0
0.0 2.0 0.0
0.0 0.0 3.0


And here's the output:

Sample problem number 1

Matrix size is 3 x 3

The matrix [a]:
1.000000e+00 0.000000e+00 0.000000e+00
0.000000e+00 2.000000e+00 0.000000e+00
0.000000e+00 0.000000e+00 3.000000e+00

The calculated inverse of the matrix [a]:
1.000000e+00 0.000000e+00 0.000000e+00
0.000000e+00 5.000000e-01 0.000000e+00
0.000000e+00 0.000000e+00 3.333333e-01

The matrix [a] times the calculated inverse:
1.000000e+00 0.000000e+00 0.000000e+00
0.000000e+00 1.000000e+00 0.000000e+00
0.000000e+00 0.000000e+00 1.000000e+00


Else
A universal method Of initialization that wil solve all the routines of NR.
What the heck do you mean by that? Each routine that requires some input must have a program that, somehow, declares variables (matrices, vectors, etc.), and that, somehow, assigns values to the variables, and then calls the routine.

You can put the assignment statements in the program itself and recompile the program for each problem, or you can read values from a file, as I showed in the example.

Or you can get the values in any other way that suits your sense of programming style.

Regards,

Dave