Solved

C++ code in Gauss Jordan Matrix reverse

Posted on 2012-03-31
6
1,184 Views
Last Modified: 2016-02-10
I have found a few on the net but i dont want to just copy them.
Here is part of my code.
What am i missing ?
HWK6.txt
0
Comment
Question by:c_hockland
[X]
Welcome to Experts Exchange

Add your voice to the tech community where 5M+ people just like you are talking about what matters.

  • Help others & share knowledge
  • Earn cash & points
  • Learn & ask questions
  • 3
  • 2
6 Comments
 
LVL 53

Expert Comment

by:Infinity08
ID: 37790191
>> What am i missing ?

An explanation of what problem you're having, and what you need our help with ;)
0
 

Author Comment

by:c_hockland
ID: 37790193
i need to use the Gauss Jordan method to find the inverse of a matrix nxn dimensions.
0
 
LVL 22

Expert Comment

by:ambience
ID: 37790709
Your choice of languages for the question is pretty intetersting. Few points

- Inverse is only possible for a square matrix so there is no need to read dimensions as row, col. This will simplify logic.
- I would suggest building a c++ abstraction say Matrix to handle matrices easily. A quick search took me here and looks like a good starting point. Use the Matrix class from here

http://stackoverflow.com/questions/2286991/c-two-dimensional-stdvector-best-practices

>> What am i missing ?

1) The identity matrix to start with.
2) You should use double for algorithms involving mathematical operations like divisions, otherwise chances are you'd get zeroed out results because of truncations.
3) I would suggest to simplify things by writing basic matrix operations needed for guass jordan method. In other words, build a DSL (domain specific language)

Whereever you are performing operations on A you should also be performing identical operation on I (identity matrix).

				 for (i=0 ; i < n  ; i++ )
				 {
						double buffer = A[i][i];
						for  (j=1 ; j < n ; j++ )
						{
							A[i][j] /= buffer ;
							I[i][j] /= buffer ;
						}
				 }

Open in new window


However, there are so many other issues with the code that I think I'll suffice with these for now and hopefully in a few iterations we can fix all.
0
Salesforce Has Never Been Easier

Improve and reinforce salesforce training & adoption using WalkMe's digital adoption platform. Start saving on costly employee training by creating fast intuitive Walk-Thrus for Salesforce. Claim your Free Account Now

 
LVL 22

Accepted Solution

by:
ambience earned 500 total points
ID: 37790815
This was an interesting problem and I rolled out a quick example of what I was saying using the Matrix class from the link. This is minimalistic and only sets up the basics and does not take into account error conditions. You can fatten the solver class. Hope this helps ...

#include "stdafx.h"
#include <assert.h>
#include <algorithm>
#include <memory>
#include <vector>
#include <iostream>

template<class T, class A=std::allocator<T> >
struct Matrix {
	typedef T value_type;
	typedef std::vector<value_type, A> Container;

	Matrix() : _b(0) {}
	Matrix(int rows, int cols, value_type const& initial=value_type())
		: _b(0)
	{
		resize(rows, cols, initial);
	}
	Matrix(Matrix const& other)
		: _data(other._data), _b(other._b)
	{}

	Matrix& operator=(Matrix copy) {
		swap(*this, copy);
		return *this;
	}

	bool empty() const { return _data.empty(); }
	void clear() { _data.clear(); _b = 0; }

	int numRows() const { return _b ? _data.size() / _b : 0; }
	int numCols() const { return _b; }

	bool rowIsNull(int row) {
		for(int i=0; i<numCols(); i++)
			if(_data[row * _b + i] != 0)
				return false;
		return true;
	}

	bool colIsNull(int col) {
		for(int i=0; i<numRows(); i++)
			if(_data[i * _b + col] != 0)
				return false;
		return true;
	}

	value_type& operator()(int i) {
		return _data[i];
	}

	const value_type& operator()(int i) const {
		return _data[i];
	}

	value_type& operator()(int row, int col) {
		return _data[row * _b + col];
	}

	const value_type& operator()(int row, int col) const {
		return _data[row * _b + col];
	}

	void resize(int a, int b, value_type const& initial=value_type()) {
		if (a == 0) {
			b = 0;
		}
		_data.resize(a * b, initial);
		_b = b;
	}

	friend void swap(Matrix& a, Matrix& b) {
		using std::swap;
		swap(a._data, b._data);
		swap(a._b,    b._b);
	}

	std::ostream& print(std::ostream& s) 
	{
		s << "<Matrix" << numRows() << 'x' << numCols();
		if (!empty()) 
		{
			bool first = true;
			typedef typename Container::const_iterator Iter;
			Iter i = _data.begin(), end = _data.end();
			while (i != end) 
			{
				s << (first ? " [[" : "], [");
				first = false;
				s << *i;
				++i;
				for (int b = _b - 1; b; --b) {
					s << ", " << *i;
					++i;
				}
			}
			s << "]]";
		}
		s << '>' << std::endl;
		return s;
	}

private:
	Container _data;
	int _b;
};


template <class M>
struct MatrixBuilder
{
	typedef typename M::value_type value_type;
	typedef M return_type;

	MatrixBuilder(int a) :_m(a, a) {
	}

	MatrixBuilder(int a, int b) :_m(a, b) {
	}

	MatrixBuilder(const M& m) :_m(m) {
	}

	MatrixBuilder& zero() {
		_m = M(_m.numRows(), _m.numCols(), value_type(0));
		return *this;
	}


	MatrixBuilder& random() {
		_m = M(_m.numRows(), _m.numCols(), value_type(0));
		for(int i=0;i<_m.numRows();i++)
		for(int j=0;j<_m.numCols();j++)
			_m(i,j) = rand();
		return *this;
	}

	MatrixBuilder& identity() {
		assert(_m.numRows() == _m.numCols());
		_m = M(_m.numRows(), _m.numCols(), value_type(0));
		for(int i=0;i<_m.numRows();i++) {
			_m(i,i) = 1;
		}
		return *this;
	}

	MatrixBuilder& addToRow(int row, value_type v) {
		assert(!_m.empty());
		assert(row <= _m.numRows());
		for(int i=0;i<_m.numCols();i++) {
			_m(row, i) += v;
		}
		return *this;
	}

	MatrixBuilder& addToColumn(int col, value_type v) {
		assert(!_m.empty());
		assert(col <= _m.numCols());
		for(int i=0;i<_m.numRows();i++) {
			_m(i, col) += v;
		}
		return *this;
	}

	MatrixBuilder& multiplyRowBy(int row, value_type v) {
		assert(!_m.empty());
		assert(row <= _m.numRows());
		for(int i=0;i<_m.numCols();i++) {
			_m(row, i) *= v;
		}
		return *this;
	}

	MatrixBuilder& multiplyColumnBy(int col, value_type v) {
		assert(!_m.empty());
		assert(col <= _m.numCols());
		for(int i=0;i<_m.numRows();i++) {
			_m(i, col) *= v;
		}
		return *this;
	}

	MatrixBuilder& divideRowBy(int row, value_type v) {
		assert(!_m.empty());
		assert(row <= _m.numRows());
		for(int i=0;i<_m.numCols();i++) {
			_m(row, i) /= v;
		}
		return *this;
	}

	MatrixBuilder& divideColumnBy(int col, value_type v) {
		assert(!_m.empty());
		assert(col <= _m.numCols());
		for(int i=0;i<_m.numRows();i++) {
			_m(i, col) /= v;
		}
		return *this;
	}

	MatrixBuilder& swapRows(int row1, int row2) {
		assert(!_m.empty());
		assert(row1 <= _m.numRows());
		assert(row2 <= _m.numRows());

		for(int i=0;i<_m.numCols();i++) {
			std::swap( _m(row1, i), _m(row2, i) );
		}
		return *this;
	}

	MatrixBuilder& swapColumns(int col1, int col2) {
		assert(!_m.empty());
		assert(col1 <= _m.numCols());
		assert(col2 <= _m.numCols());

		for(int i=0;i<_m.numRows();i++) {
			std::swap( _m(i, col1), _m(i, col2) );
		}
		return *this;
	}

	MatrixBuilder& addRowElements(int row1, int row2, value_type multiplier) {
		assert(!_m.empty());
		assert(row1 <= _m.numRows());
		assert(row2 <= _m.numRows());

		for(int i=0;i<_m.numCols();i++) {
			_m(row1, i) += _m(row2, i) * multiplier;
		}
		return *this;
	}

	MatrixBuilder& subRowElements(int row1, int row2, value_type multiplier) {
		assert(!_m.empty());
		assert(row1 <= _m.numRows());
		assert(row2 <= _m.numRows());

		for(int i=0;i<_m.numCols();i++) {
			_m(row1, i) -= _m(row2, i) * multiplier;
		}
		return *this;
	}

	M result() const  {
		return _m;
	}

private:
	M _m;
};


template <class M>
struct GuassJordanSolver
{
	typedef typename M::value_type value_type;
	typedef M return_type;

	GuassJordanSolver(const M& m) :_m(m), _b(m), _i(m.numRows()) {
		assert(!_m.empty());
		assert(_m.numRows() == _m.numCols());
		_i.identity();
	}

	bool solve()
	{
		for(int i=0; i<_m.numCols(); i++)
		{
			value_type v = _b.result()(i, i);

			// Interchange rows to bring a non-zero element, only look for rows greater than current
			if(v == 0 && !fixNullPivot(i, i))
			{
				return false;
			}

			v = _b.result()(i, i);
			_b.divideRowBy(i, v);
			_i.divideRowBy(i, v);

			_b.result().print(std::cout);

			for(int j=0; j<_m.numRows(); j++)
			{
				if(i == j) continue;

				v = _b.result()(j, i);

				if(v == 0) continue;

				_b.subRowElements(j, i, v);
				_i.subRowElements(j, i, v);

				_b.result().print(std::cout);
			}

		}

		return true;
	}

	M original() const  {
		return _m;
	}
	
	M transformed() const  {
		return _b.result();
	}

	M inverse() const  {
		return _i.result();
	}

private:
	bool fixNullPivot(int r, int c)
	{
		for(int i = r+1; i<_m.numRows(); i++)
		{
			if(_b.result()(i, c) != 0)
			{
				_b.swapRows(r, i);
				_i.swapRows(r, i);
				_b.result().print(std::cout);
				return true;
			}
		}
		return false;
	}

private:
	M _m;
	MatrixBuilder<return_type> _b;
	MatrixBuilder<return_type> _i;
};



Matrix<double> fromData() 
{
	Matrix<double> _m = Matrix<double>(3, 3, 0);
	_m(0, 0) = 0;
	_m(0, 1) = -2;
	_m(0, 2) = 3;

	_m(1, 0) = 2;
	_m(1, 1) = 0;
	_m(1, 2) = 6;

	_m(2, 0) = 3;
	_m(2, 1) = 5;
	_m(2, 2) = 0;
	return _m;
}

int _tmain(int argc, _TCHAR* argv[])
{
	MatrixBuilder<Matrix<double>> b(3);
	GuassJordanSolver<Matrix<double>> s(fromData());

	s.solve();
	s.original().print(std::cout);
	s.transformed().print(std::cout);
	s.inverse().print(std::cout);

	return 0;
}

Open in new window

0
 

Author Comment

by:c_hockland
ID: 37790825
thnks for the quick input. We are not supposed to use classes. Only matrices and functions if needed. Let me try the code you sent me...
0
 
LVL 22

Expert Comment

by:ambience
ID: 37790833
In that case it shouldnt be that difficult to convert c++ code into equivalent c code, just pass matrix as first parameter in all methods and make them functions.
0

Featured Post

[Webinar] Code, Load, and Grow

Managing multiple websites, servers, applications, and security on a daily basis? Join us for a webinar on May 25th to learn how to simplify administration and management of virtual hosts for IT admins, create a secure environment, and deploy code more effectively and frequently.

Question has a verified solution.

If you are experiencing a similar issue, please ask a related question

Suggested Solutions

Windows programmers of the C/C++ variety, how many of you realise that since Window 9x Microsoft has been lying to you about what constitutes Unicode (http://en.wikipedia.org/wiki/Unicode)? They will have you believe that Unicode requires you to use…
Performance in games development is paramount: every microsecond counts to be able to do everything in less than 33ms (aiming at 16ms). C# foreach statement is one of the worst performance killers, and here I explain why.
The goal of the tutorial is to teach the user how to use functions in C++. The video will cover how to define functions, how to call functions and how to create functions prototypes. Microsoft Visual C++ 2010 Express will be used as a text editor an…
The viewer will learn additional member functions of the vector class. Specifically, the capacity and swap member functions will be introduced.

737 members asked questions and received personalized solutions in the past 7 days.

Join the community of 500,000 technology professionals and ask your questions.

Join & Ask a Question