Matrix Inversion help :'(

[xb211]

Hi All,

I'm currently working on my project (which try to perform degree reduction for a Bezier curve).

I've spent a long time trying to find library or useful class to help me perform matrix inversion.

For example degree reduction in Bezier's curve (let's say from degree 3 to degree 2)

B2: bezier curve for degree of 2 (which is unknown in this case)
B3: bezier curve for degree of 3 (which is a given in this case)
M: A matrix involving bezier curves (I already wrote a function to get this matrix)

The fomula to find B2 is as followed:

B2 = inverse(Mtranspose * M) * M * B3

I wrote function to perform matrices multiplication and matrix transpose. I'm having a great deal of difficulty integrating TNT::JAMA_LU and other classes I've found. I just want to do an inverse of a matrix (dimension NxN)

I've spent too much time trying to figure out how to do matrix inversion and that's not even the main point of my project. I chose to use C++ because I'm using OpenGL to plot my bezier curve... and I can't seem to integrate TNT::JAMA_LU to my program. *sob*

Also... i've tried the information from (http://www.dreamincode.net/forums/showtopic43006.htm) .. and for some reason i just have NAN as my result :S Or maybe im not integrating it correctly. But.. I just copy and pasted the code. :-S

Can anyone provide me with good resources? I really appreciate any help in advance.

** additional information
My current class Matrix has the following variables:
- double* data <-- which contain the list of double values in the matrix
- rows <-- number of row in the matrix
- cols <-- number of column in matrix

current functions i have for my matrix class
- cells(int row, int col) <- will get the exact data from double* data
- declareMatrix
- getActualSize
- getCols
- getRows
- invert <--- trying to do this... but i'm stuck .. and don't know where else to go
- multiply (which takes in two matrices and return a multiplied matrix)
- transpose (do a transpose on a matrix)
- setData (can take a set of double* data and put it in a matrix of dimension NxM

I hope my problem/request is clear...

THANK YOU very much in advance for any help

[xb211]

Addition things I've tried:

I got java code that does inverse matrix:
http://www.daniweb.com/forums/thread118289.html

And converted it into c++ code: when i execute it... my program will hang... :-s

Matrix inverse(Matrix a){
		int n=a.getRows(a);
		
		Matrix x,b;
		x.declareMatrix(n,n);
		b.declareMatrix(n,n);

		// transform the matrix into an upper triangle
		//a = gaussian(a);

		int* index;
		index = new int[n];
		double* c;
		c = new double[n];
		
		int i,j;
		for(i=0; i<n; ++i){
			index[i] = i;
		}

		//find rescaling factors, one from each row
		for(i=0; i<n; ++i){
			double c1 = 0;
			for(j=0; j<n; ++j){
				double c0 = abs(a.cells(i,j));
				if(c0>c1){
					c1 = c0;
				}
				else{
				}
			}
			c[i] = c1;
		}

		// search the pivoting element from each column
		int k=0;
		for(j=0; j<n-1; ++j){
			double pi1 = 0;
			for(i=j; i<n; ++i){
				double pi0 = abs(a.cells(i,j));
				pi0 /= c[index[i]];
				if(pi0 > pi1){
					pi1 = pi0;
					k=i;
				}
			}
		}

		// interchange rows according to pivoting order
		int itmp = index[j];
		index[j] = index[k];
		index[k] = itmp;

		for(i=j+1; i<n; ++i){
			double pj = a.cells(index[i],j)/a.cells(index[j],j);

			// recording pivoting ratios below the diagonal
			a.cells(index[i],j) = pj;

			// modify other elements accordingly
			for(int l=j+1; l<n; ++l){
				a.cells(index[i],l) -= pj*a.cells(index[j],l);
			}
		}
		// update matrix b
		for(i=0; i<n-1; ++i){
			for(j=i+1; j<n; ++j){
				for(k=0; k<n; ++k){
					b.cells(index[n-1],k) -= a.cells(index[j],i)*b.cells(index[i],k);
				}
			}
		}

		// perform backward substitution
		for(i=0; i<n; ++i){
			x.cells(n-1,i) = b.cells(index[n-1],i) / a.cells(index[n-1],n-1);
			for(i=n-2; j>=0; --j){
				x.cells(i,j) = b.cells(index[j],i);
				for(k=j+1; k<n; ++k){
					x.cells(j,i) -= a.cells(index[j],k)*x.cells(k,i);
				}
				x.cells(j,i) /= a.cells(index[j],j);
			}
		}
		return x;
	}

I'm not sure what's wrong.. :(

[ArkM]

It's suprising that you "have spent a long time trying to find library or useful class to help me perform matrix inversion". I have found then download some tested and ready-to-use freeware libraries in C++ with matrix classes and inverse operation for 10 minutes...

Regrettably at the moment I have no time to post all useful links gotten after my Google search (and more deeper surfing after that). May be you can post your test matricies examples (where your home-made product failed). Later I will try to install and test some of these libraries...

[xb211]

Hi...

Here's the test cpp file i've created to show my problem. I think type double is more precise in java than c++. I when i'm doing some comparison... (ex: a>b)... it might not be as precise to detect the nth decimal place. (ex: from line 126 and 150)

#include <iostream>
#include <math.h>

using namespace std;
#define PI 3.14159265

class Matrix{
public:
	int rows;
	int cols;
	int actualSize;
	double* data;

	Matrix(){
		rows = 1;
		cols = 1;
		actualSize = rows*cols;
		data = new double[rows*cols];
	}

	double &cells(int r, int c){
		return data[r*cols+c];
    }

	void declareMatrix(int r, int c){
		Matrix::rows = r;
		Matrix::cols = c;
		Matrix::actualSize=r*c;
		data = new double[r*c];
	}

		Matrix multiply(Matrix a, Matrix b){
		int r = a.rows;
		int c = b.cols;
		int middle = a.cols;
		double sum;
		Matrix result;
		result.declareMatrix(r,c);
		
		double* tempData;
		tempData = new double[r*c];

		int index = 0;

		for(int i=0; i<r; i++){
			for(int j=0; j<c; j++){
				sum = 0;
				for(int k=0; k<middle; k++){
					
					sum = sum + (a.cells(i,k)*b.cells(k,j));
					result.cells(i,j) = sum;
				}
			}
		}

		return result;	
	}

	Matrix transpose(Matrix a, int r, int c){
		Matrix result;
		result.declareMatrix(c,r);

		for(int i=0; i<r; i++){
			for(int j=0; j<c; j++){
				if(i>c){
				}
				else{
					result.cells(j,i) = a.cells(i,j); 
				}
			}
		}

		return result;
	}

	int getActualSize(){
		return actualSize;
	}

	int getRows(Matrix a){
		return a.rows;
	}

	int getCols(Matrix a){
		return a.cols;
	}

	Matrix inverse(Matrix a){
		int n=a.getRows(a);
		int i,j;

		Matrix x,b;
		x.declareMatrix(n,n);
		b.declareMatrix(n,n);

		for(i=0; i<n; ++i){
			for(j=0; j<n; ++j){
				b.cells(i,j) =1;
			}
		}
		
		int* index;
		index = new int[n];
		
		// initializing index
		for(i=0; i<n; ++i){
			index[i] = i;
		}

		double* c;
		c = new double[n];
		
		//find rescaling factors, one from each row
		for(i=0; i<n; ++i){
			double c1 = 0;
			for(j=0; j<n; ++j){
				double c0;
				if(a.cells(i,j) < 0){
					c0 = (a.cells(i,j))*(-1);
				}
				else{
					c0 = (a.cells(i,j));
				}
				if(c0>c1){ //<<---- i think perhaps the problem lies here.. where double is not as precise as in java
					c1 = c0;
				}
			}
			c[i] = c1;
		}

		// search the pivoting element from each column
		int k=0;
		double pi0 = 0;
		int itmp = 0;

		for(j=0; j<n-1; ++j){
			double pi1 = 0;
			for(i=j; i<n; ++i){
				
				if(a.cells(index[i],j)<0){
					pi0 = a.cells(index[i],j) * (-1);
				}
				else{
					pi0 = a.cells(index[i],j);
				}
				pi0 /= c[index[i]];
				
				if(pi0 > pi1){ //<<---- i think perhaps the problem lies here.. where double is not as precise as in java
					pi1 = pi0;
					k=i;
				}
			}
		
			// interchange rows according to pivoting order
			itmp = index[j];
			index[j] = index[k];
			index[k] = itmp;

			for(i=j+1; i<n; ++i){
				double pj = (a.cells(index[i],j))/(a.cells(index[j],j));

				// recording pivoting ratios below the diagonal
				a.cells(index[i],j) = pj;

				// modify other elements accordingly
				for(int l=j+1; l<n; ++l){
					a.cells(index[i],l) -= pj*a.cells(index[j],l);
				}
			}
		}
	
		// update matrix b
		for(i=0; i<n-1; ++i){
			for(j=i+1; j<n; ++j){
				for(k=0; k<n; ++k){
					b.cells(index[j],k) -= a.cells(index[j],i)*b.cells(index[i],k);
				}
			}
		}

		// perform backward substitution
		for(i=0; i<n; ++i){
			x.cells(n-1,i) = b.cells(index[n-1],i) / a.cells(index[n-1],n-1);
			for(j=n-2; j>=0; --j){
				x.cells(j,i) = b.cells(index[j],i);
				for(k=j+1; k<n; ++k){
					x.cells(j,i) -= a.cells(index[j],k)*x.cells(k,i);
				}
				x.cells(j,i) /= a.cells(index[j],j);
			}
		}
		return x;
	}
};

void dumpInfo(Matrix a){
	for(int i=0; i<a.getRows(a);i++){
		for(int j=0; j<a.getCols(a);j++){
			printf("%f\t", a.cells(i,j));
		}
		printf("\n");
	}
}

Matrix build_M_matrix(int n){
	//n=3
	Matrix result;
	double ratio1, ratio2;
	int r = n+2;
	int c = n+1;

	result.declareMatrix(r, c);
	
	int skip = 0;
	double tempI;
	double tempN;
	int k=0;

	for(int i=0; i<r; i++){
		tempI = (double)i;
		tempN = (double)n;
		ratio1 = (1-(tempI/(tempN+1)));
		ratio2 = (tempI/(tempN+1));

		if(i<=0){
			result.cells(i,i) = ratio1;
			result.cells(i,i+1) = ratio2;
			for(k=i+2; k<c; k++){
				result.cells(i,k) = 0;
			}
		}

		else{

			if(i==(r-1)){
				for(k=0; k<skip; k++){
					result.cells(i,k) = 0;
				}
				result.cells(i,skip) = ratio2;
				result.cells(i,skip+1) = ratio1;
			}
			else{
				for(k=0; k<skip; k++){
					result.cells(i,k) = 0;
				}
				result.cells(i,skip) = ratio1;
				result.cells(i,skip+1) = ratio2;

				for(k=skip+2; k<c; k++){
					result.cells(i,k) = 0;
				}
				skip++;
			}
		}
	}

	return result;
}
void testing(){
	Matrix mtest, mtTest;

	int asdf=2;
	mtest = build_M_matrix(asdf); //M

	printf("--- m ---\n");
	dumpInfo(mtest);

	mtTest = mtest.transpose(mtest,mtest.getRows(mtest), mtest.getCols(mtest)); //M transpose

	printf("--- m transpose ---\n");
	dumpInfo(mtTest);

	mtTest = mtest.multiply(mtest, mtTest);//(mtest,mtest.getRows(mtest), mtest.getCols(mtest)); //M transpose

	printf("--- m* m transpose ---\n");
	dumpInfo(mtTest);

	printf("--- Inverted m* m transpose ---\n");
	mtTest = mtTest.inverse(mtTest);
//	result = mtest.invert(testing,mtTest);
	dumpInfo(mtTest);
}//end display


int main(int argc, char* argv[])
{
	
	testing();
	return 0;
}

[xb211]

Nevermind! it wasn't the precision problem. PROBLEM SOLVED! Wooo... :-)