8 Replies - 30629 Views - Last Post: 29 September 2007 - 08:43 AM Rate Topic: -----

#1 Max_Power82  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 6
  • Joined: 03-December 06

Gauss-Jordan Elimination matrix

Posted 03 December 2006 - 11:53 AM

Hi guys, I'm new to c++, i'm writing code to do Gauss-Jordan elimination on a matrix.

It is doing it almost correctly, except that every other entry in my inverse is coming up as zero! There must be a problem with my loops, but the numbers that ARE there are correct, so it must be something minor.

Any help appreciated. Very sorry for the messy and unproffessional but I am a new student! YOU CAN IGNORE MOST OF THE CODE, The problem function in question is gaussJordan(.,.) called through inv(.,.), it's near the bottom.

Try running the program and ask it to do inverse for N= 5 or 6, you'll see a 'diamond' pattern of zeros in the inverse matrix.

Thanks very much




#include <iostream>
#include <cmath>
#include <ctime>
#include <iomanip>
using namespace std;

double gausstoosmall = 0.000000001;  //global declaration of tolerance for pivot size 

class vector{
	private:
		int size;
		double *vec;  // my vector is going to be a pointer to an array.
		int range(int); //this is used by the indexing operator.  It makes the class more robust, by making it impossible to return a value outside the array of the objet in question.  same thing in the matrix class. 
		
	public:
		vector(int);  //vector constructor for empty vector of given length
		vector(const double, int); //vector constructor for initialised vector
		vector(const vector&); 
		~vector();		//destructor
		
		double A(int i, int j);
		friend class matrix;

		inline double& vector::operator[](int i){ return vec[range(i)];} // declare an indexing operator 
		inline int vector::getsize(){ return size;}
		vector& operator=(const vector &v);

		void swap(int i, int j);
		void vprint(vector& vc); // this is my home-made vector output function. 
		void permute(vector &v, vector &p);
		vector index(int n);
};

class matrix{
	private:
		int numrows;
		int numcols;
		vector **mat;  //i.e. my matrices are pointers to vectors (so pointers to arrays of pointers to arrays. 
		int range(int);
	public:
		matrix(int, int);	//matrix constructor for empty matrix of given size
		
		matrix(const matrix&);
		//~matrix();	 // destructor
		inline vector& matrix::operator [](int i){return *mat[range(i)];}
		inline int matrix::getrows(){ return numrows;}
		inline int matrix::getcols(){ return numcols;}
		matrix& operator=(const matrix &m); 
		double A(int i, int j); //used to initailise the vectors in parts a) and iv)
		void A_Init(matrix &ma);

		void mprint(matrix& ma);
		void triangulate(matrix &a, vector &b, vector &x); 
		void swap(int i, int j);
		void pivot(matrix &a, vector &b, vector &x, int h);
		int pivotSign(matrix &m, vector &p, int i);
		void backSubs(matrix &a, vector &b, vector &x); 
		matrix inverse(matrix &m); 
		void gaussJordan(matrix &a, matrix &b);
		void inv(matrix &ma, matrix &inv);
};



int main(){
//part a) 
	int N; //user selects the value of N (i.e. size of square matrix A.
	cout << "\n\nPart a)\n\nEnter the dimension you require for the NxN matrix [A].  N= ";
	cin >> N;
	
	matrix Amatrix(N,N); //create the empty vector [A] 
		Amatrix.A_Init(Amatrix);
		cout << "\n\nMatrix [A] has been initialised as follows:\n\n";
		Amatrix.mprint(Amatrix);  //this function prints the matrix to the monitor
	vector Bvector(1, N); //create the vector B, using the second constructor to initialise to all ones. 
		//for (int v=0; v<Bvector.getsize(); ++v){	Bvector[v] = 1;}; //initialise B vector as vector of ones.
		cout << "\n\nThe vector B is initialised as follows: \n\n"; 
		Bvector.vprint (Bvector); //custom vector print function
	vector X(N); //the purpose of this vector is to keep a track of which x_n value is which after row swapping has taken place.  It will be initialised to {1,2,3,4,etc...}
		for(int i=0; i<N; ++i) X[i] = i + 1;
		cout << "\n\nThe vector of variables X is initialised as follows: \n\n"; 
		X.vprint(X); 
		cout << "\n(Where the number indicates the subscript of the x variable within the X vector.)\n\n"; 

	cout << "\n\nPress enter to perform Gaussian elimination on this system of linear equations.\n\n";
		cin.get(); cin.get();
		cout  << "After diagonalisation, Matrix [A] is transformed as follows:\n\n"; 
		Amatrix.triangulate(Amatrix, Bvector, X);
		Amatrix.backSubs(Amatrix, Bvector, X);
	
//Part i)
	//cout << "\n\nPart i)\n\nNow we set N = 3.  The matrix [A] is now defined and initialised as:\n\n"; 
	cout << "Enter value of  N for size of matrix, matrix will be initialised for you, a[i][j] = (i*j )/ (i+j-1): \n";
cin >> N;
	//N = 3;	
		matrix E(N,N);
		E.A_Init(E);
		E.mprint(E);
		matrix Einv(N,N);  //here we create an empty matrix that will become the inverse. 
	cout << "\n\nTo calculate to Inverse of matrix [A], press enter...\n";
		cin.get();
		E.inv(E, Einv);
		



		cout << "Press enter to continue..."; 
//part iv)
		cin.get();
	cout << "\n\nPart iv)\n\nFor the special case where N = 10 , the NxN matrix C is initialised to: [c] =\n\n";
		matrix C(10, 10);
		C.A_Init(C); //initialise the 10x10 matrix. 
		C.mprint(C);
		matrix Cinv(10, 10);  //create a new matrix object (empty).  This will become inverse of Cmatrix.
	cout << "\n\nTo calculate to Inverse of matrix [c], press enter...\n"; 
		cin.get();	
		C.inv(C, Cinv); //
	cout << "Press enter to continue...";
	cin.get();

//part ii)
	cout << "\n\nPart ii)\n\n";
	



	return 0;
}
vector::vector(int n){
	size = n;
	vec = new double [size]; 
	if (!vec)
		cout << "ERROR\tMemory allocation failed in: vector::vector(int).\n\n"; //this makes the program more robust.
	for (int i = 0; i < size; i++) vec[i] = 0;
}
vector::vector(const double d, int n){ 
	size = n;
	vec = new double [size];
	if (!vec)
		cout << "ERROR\tMemory allocation failed in: vector::vector(double*, int).";
	for (int i = 0; i < size; i++) vec[i] = d; 
}
vector::~vector() {delete vec;}

int vector::range(int i)
{if(i <0|| i>= size) return -1; else return i;}

double vector::A(int i, int j){
	return (i * j)/(i + j - 1.0);
}
void vector::vprint(vector& vc){ 
	cout << "\n\n";
	for (int s=0; s < vc.getsize(); ++s) 
				{cout << "[" << vc[s] << "]\n";};
}

void vector::swap(int i, int j){ 
	double temp = vec[range(i)];
	vec[i] = vec[range(j)];
	vec[j] = temp;

}
vector& vector::operator=(const vector &v){
	if (size != v.size)
		{cout << "ERROR: Size mismatch."; return *this;} 
	for(int i = 0; i <size; ++i) vec[i] = v.vec[i];
	return *this;
}

void permute(vector &v, vector &p){
	int n = v.getsize();
	vector u(n);
	for (int i = 0; i < n; ++i) 
		u[i] = v[int(p[i])];
	v = u;
}

matrix::matrix(int nrows, int ncols){		//constructor for NxM matrix
	numrows = nrows; numcols = ncols;
	mat = new vector* [numrows];
	if (!mat) cout << "ERROR:\tRow allocation failed in matrix::matrix(int, int).\n"; 
	for (int i = 0; i < numrows; ++i)
	{	mat[i] = new vector(numcols);
		if (!mat[i]) cout << "ERROR:\tColumn allocation failed in matrix::matrix(int, int).\n";
	};
}
int matrix::range(int i) 
{if (i < 0 || i >= numrows) return -1; else return i;}

double matrix::A(int i, int j){
	return (i * j)/(i + j - 1.0);
}
void matrix::mprint(matrix& ma){
	for (int r =0; r<ma.getrows (); ++r)
	{cout << "["; for (int c=0; c < ma.getcols(); ++c) 
				{cout << ma[r][c] << "\t";}
	cout << "]\n";}
}



void matrix::triangulate(matrix &a, vector &b, vector &x){ //this is the centrepiece of my class, it does gauss elimination, and prepares the matrix for backsubstitution by getting ones on the diagonals.  it also deals with vector B and the vector X of x_i variables.  Devised by me entirely. 

for(int j = 0; j < a.getcols()-1; ++j)
{
	for(int p = j + 1; p < a.getrows(); ++p)
	{	a.pivot(a, b, x, j);
		double l = (-1.0*a[p][j])/(a[j][j]);
		b[p]  += l*b[j];						 // this deals with the b vector 
			for(int i=j; i< a.getcols(); i++)
			{	
				a[p][i] += l * a[j][i];
			}
	}		
}
	for(int k = 0; k < a.getrows(); ++k) // here we get ones on the diagonal.  I am assuming that there will be no zeros on the diagonal, which is reasonable since this can not happen in this question. 
	{	double y = a[k][k];
		for(int m = 0; m < a.getcols(); ++m)
		{	
			a[k][m] = (a[k][m])/y;
			b[k] = (b[k])/y;
			if(k>m) {
						a[k][m] = fabs(a[k][m]); //this is purely cosmetic, it gets rid of minus signs on the lower diagonal entries, which are of course now zero.  If you wish to check this, just comment this line out, the program still works in exactly the same way. 
						if((a[k][m])!=0 && (a[k][m])< 0.000001) (a[k][m])=0; };  // likewise, this is cosmetic, it just gets rid or tiny values which remain on the lower diagonal due to rounding errors, we're talking 10exp-16 small.  You can get rid of it if you like. 
		}
	
	}

}
void matrix::swap(int i, int j){
	vector *temp = mat[range(i)];
	mat[i] = mat[range(j)];
	mat[j] = temp;
}
void matrix::pivot(matrix &a, vector &b, vector &x, int h){ 
	int n = a.getrows();
	int g = h;
	double t = 0;
	for (int k = h; k<n; ++k)
		{double elm = fabs(a[k][h]);  //finds the largest value in the column, the "pivot".
			if(elm > t) {t = elm; g=k;} 
		}
	if (g > h)
		{	a.swap(h,g);
			b.swap(h,g);
			x.swap(h,g);  //this is to swap the entries of the X vector, so that we can keep track of the solution properly.  This may not be useful in this example because x2 to xN are all zero, but in general it would be useful.	 
		}
} 
void matrix::backSubs(matrix &a, vector &b, vector &x){
	cout << "The solution to the system of linear equations [A]*X = B, is as follows:\n\n";
	cout << "After Gaussian Elimination with Full Pivoting, the matrix [A] = \n\n"; 
	a.mprint(a);
	cout << "\n\n, the vector B = \n";
	b.vprint(b);
	cout << "\n\n, and the variable vector X = \n";
	x.vprint(x); cout << "(Where the number indicates the subscript of the x variable.)\n\nTherefore:\n\n"; 
	
	for (int i = 0; i < a.getcols(); ++i)
		{	cout << "x_" << x[i] << " = " << (b[i])/a[i][i] << endl;
		};
	cout << "\n(The order in which the x_i variables appear reflects their respective final order in the X vector.)\n\n"; 
}

void matrix::A_Init(matrix &ma){
	for (int r =0; r< ma.getrows(); ++r) //this loop initialises the matrix [A] ("Amatrix".) using the memeber function A. 
	{	for (int c=0; c < ma.getcols(); ++c) 
				ma[r][c] = ma.A(r+1,c+1); 
		};
}



  
void matrix::gaussJordan(matrix &a, matrix &b){ 
	matrix GJ(a.getrows(), (a.getcols()*2));
		for(int r = 0; r<a.getrows(); ++r)
		{	for(int c = 0; c<a.getcols(); ++c)
			{	GJ[r][c] = a[r][c];
			
			} 
			GJ[r][r + a.getcols()] = 1;
		}	
		vector u(GJ.getrows()); //dummy vector for pivot

		cout << "\nThe Augmented formed by placing an NxN identity matrix to the right of our matrix is:\n"; 
		GJ.mprint(GJ);	cout << "\n\n";

for(int j = 0; j < GJ.getcols()-1; ++j) //triangulate leftside of composite matrix. this gets zero in every row column, below diagonal.
{
	for(int r = j + 1; r < GJ.getrows(); ++r) //this for each individual column.
	{	GJ.pivot(GJ, u, u, j);
		double l = (-1.0*GJ[r][j])/(GJ[j][j]);
								 
			for(int c=j; c< GJ.getcols(); c++)  //this to multiply each element in the row. 
			{	
				GJ[r][c] += l * GJ[j][c];

			}
	}		
} //GJ.mprint(GJ);
	for(int r = 0; r < GJ.getrows(); ++r) // here we get ones on the leading diagonal.  I am assuming that there will be no zeros on the diagonal, which is reasonable since this can not happen in this question. 
	{	double y = GJ[r][r];
		for(int c = 0; c < GJ.getcols(); ++c)
		{	
			GJ[r][c] = (GJ[r][c])/y;

		}
	}
	for (int k= ((GJ.getcols())/2-1); k > 0; --k) //the outer column loop, this puts zeros in successive upper diag of left side of augmented matrix GJ 
	{	
		for(int r = GJ.getrows()-2 + k -((GJ.getcols())/2-1); r > -1; --r)
		{double y = (GJ[r][k]);	
			for(int c = 0; c < GJ.getcols(); ++c)
			
			{	//cout << "k r c GJ[r][k]   " << k << "\t"<< r << "\t"<< c<< "\t" << GJ[r][k] << "\n\n\n"; 
				GJ[r][c] = (GJ[r][c])-y*(GJ[k][c]);
				//GJ.mprint(GJ);
				//cout << "\n\n";
			}
		}
	}
	for(int r = 0; r < GJ.getrows (); r++)  //this is a clean-up routine, purely aesthetic, after the algorithmhas run, to remove near zero values form the matrix.  this makes it more legible.  small values arise as rounding errors.
	{	for(int c = 0; c < GJ.getcols(); c++)
		{	if(GJ[r][c]<0.000001) GJ[r][c] = 0;
		}
	}
	cout << "After Gauss-Jordan elimination, the Reduced Row-Echelon Form of the augmented matrix is:\n\n"; 
		GJ.mprint(GJ); 
		
	
	for(int r = 0; r < a.getrows(); r++)  //here we write the inverse to the second argument matrix.
	{	for(int c = 0; c < a.getcols(); c++)
		{	b[r][c] = GJ[r][c + a.getcols()];	
		}
	}

	
}	

void matrix::inv(matrix &ma, matrix &inv){
	if (ma.getcols() != inv.getcols() || ma.getrows() != inv.getrows()) //this 'if' prevents size mismatch.  This makes the program more robust. 
		cout << "***ERROR***: Size mismatch.  Please supply an inverse with same dimensions as matrix to be inverted.";
	
	gaussJordan(ma, inv);
	cout << "Press Enter to proceed..."; cin.get();
	cout << "\n\nTherefore the inverse of the supplied matrix is: \n\n";
	inv.mprint(inv);  cout << endl;

}


Is This A Good Question/Topic? 0
  • +

Replies To: Gauss-Jordan Elimination matrix

#2 Max_Power82  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 6
  • Joined: 03-December 06

Re: Gauss-Jordan Elimination matrix

Posted 03 December 2006 - 07:34 PM

I figured it out, i was missing a call to fabs(.) in the 'if' of my cleanup subroutine, meaning that any negative values were getting wiped out! The interesting thing is that no matter what size matrix you choose, the negatives (and therefore the incrrect zeros) accur in a diamond pattern! I haven't time to work out the maths behind that, it's probably linked to the formula used to initalise the matrix. Or is it something to do eigenvalues? (Don't bother answering that...)
Was This Post Helpful? 0
  • +
  • -

#3 slyth  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 6
  • Joined: 24-September 07

Re: Gauss-Jordan Elimination matrix

Posted 24 September 2007 - 06:42 PM

hey!!! can anybody help me write the program in turbo c??? i'm having a hard time on it...pointers should be used. any help would be appreciated. my homework is due on tuesday next week. thanks!!!
Was This Post Helpful? 0
  • +
  • -

#4 slyth  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 6
  • Joined: 24-September 07

Re: Gauss-Jordan Elimination matrix

Posted 24 September 2007 - 06:52 PM

please help me with my homework!!! i'm having a hard time making a turbo c program about gauss-jordan elimination matrix. our teacher said that we should use pointers. can anybody help me with this??? any help would be appreciated. it's due tuesday next week and i'm still clueless...
again, turbo c..not c++
thanks

i forgot to say this...it has to have arrays...;D ty
Was This Post Helpful? 0
  • +
  • -

#5 PennyBoki  Icon User is offline

  • system("revolution");
  • member icon

Reputation: 53
  • View blog
  • Posts: 2,334
  • Joined: 11-December 06

Re: Gauss-Jordan Elimination matrix

Posted 24 September 2007 - 06:52 PM

View Postslyth, on 24 Sep, 2007 - 06:42 PM, said:

hey!!! can anybody help me write the program in turbo c??? i'm having a hard time on it...pointers should be used. any help would be appreciated. my homework is due on tuesday next week. thanks!!!

No one will do the homework for you, you should give it a try yourself. Then if you run into some trouble try asking for an advice here. Please read the rules before posting and welcome!
Was This Post Helpful? 0
  • +
  • -

#6 PennyBoki  Icon User is offline

  • system("revolution");
  • member icon

Reputation: 53
  • View blog
  • Posts: 2,334
  • Joined: 11-December 06

Re: Gauss-Jordan Elimination matrix

Posted 24 September 2007 - 06:53 PM

As I replied in the other thread:

No one will do the homework for you, you should give it a try yourself. Then if you run into some trouble try asking for an advice here. Please read the rules before posting and welcome!
Was This Post Helpful? 0
  • +
  • -

#7 born2c0de  Icon User is offline

  • printf("I'm a %XR",195936478);
  • member icon

Reputation: 180
  • View blog
  • Posts: 4,667
  • Joined: 26-November 04

Re: Gauss-Jordan Elimination matrix

Posted 25 September 2007 - 03:24 AM

Topics Merged.
Was This Post Helpful? 0
  • +
  • -

#8 kinit  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 1
  • Joined: 29-September 07

Re: Gauss-Jordan Elimination matrix

Posted 29 September 2007 - 07:50 AM

aw! sa cpu ka diba?! hehe! diba C na ang sa inyo? tip lang! maka bulig gid ni.. ang tanan nga cinn>> islan mo lang scanf kag tanan nga cou<< islan mo lang printf.. di ko sigurado kung 100% gid ni nga tsk2! kay lain iban ya nga datatype... hehe! bulig lang ah! chck lang! :ph34r:



[font=Comic Sans Ms][size=7]/*Max_Power82*/
if you dont mind..in java how can u get the data being outputed in the parallelports.. and put it on the textfile? can you convert that to bin digits?...TY...

This post has been edited by kinit: 29 September 2007 - 08:12 AM

Was This Post Helpful? 0
  • +
  • -

#9 PennyBoki  Icon User is offline

  • system("revolution");
  • member icon

Reputation: 53
  • View blog
  • Posts: 2,334
  • Joined: 11-December 06

Re: Gauss-Jordan Elimination matrix

Posted 29 September 2007 - 08:43 AM

View Postkinit, on 29 Sep, 2007 - 07:50 AM, said:

aw! sa cpu ka diba?! hehe! diba C na ang sa inyo? tip lang! maka bulig gid ni.. ang tanan nga cinn>> islan mo lang scanf kag tanan nga cou<< islan mo lang printf.. di ko sigurado kung 100% gid ni nga tsk2! kay lain iban ya nga datatype... hehe! bulig lang ah! chck lang! :ph34r:



[font=Comic Sans Ms][size=7]/*Max_Power82*/
if you dont mind..in java how can u get the data being outputed in the parallelports.. and put it on the textfile? can you convert that to bin digits?...TY...

Hi kinit, welcome! Please post in english.

And if you have any java questions, please post them in the java forum, but I recommend you read the forum rules first :)
Was This Post Helpful? 0
  • +
  • -

Page 1 of 1