/* 
 * File:   matrixClass.h
 * Author: bkwiatek
 * Oryginal project: http://www.dreamincode.net/forums/topic/55772-determinant-of-matrix-in-c/
 */

#ifndef MATRIXCLASS_H
#define	MATRIXCLASS_H



#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <cstdio>
#include <math.h>
#include <omp.h>
using namespace std;
#include "exception.h"
#include "pararell.h"
#include "mt19937.h"

// Declarations
template <typename MxType> class Matrix;
template <typename MxType> MxType Det(const MxType & a);
template <typename MxType> Matrix<MxType> Diag(const unsigned long n);
template <typename MxType> Matrix<MxType> Diag(const MxType& v);
template <typename MxType> Matrix<MxType> Inv(const MxType& a);
template <typename MxType> Matrix<MxType> Ones(const unsigned long rows, const unsigned long cols);
template <typename MxType> unsigned long Size(const MxType& a, const unsigned long i);
template <typename MxType> Matrix<MxType> Zeros(const unsigned long rows, const unsigned long cols);


// Matrix class

template <typename MxType> class Matrix {
public:
    /***************************************************************************
     * Constructors and destructors
     */

    /*
     * Create a Matrix object without content
     */
    Matrix() {
	this->p = NULL;
	this->rows = 0;
	this->cols = 0;
    }

    /*
     * Create a Matrix object with given number of rows and columns and initially fill in zeros for all values in the matrix
     */
    Matrix(const unsigned long row_count, const unsigned long column_count) {
	if (para.hasThreads()) {
	    this->p = NULL;
	    if (row_count > 0 && column_count > 0) {
		/*matrix config*/
		this->rows = row_count;
		this->cols = column_count;
		this->p = new double*[this->rows];
		/*pararell config*/
		unsigned long chunk = para.calcChunk(this->rows + this->cols);
		unsigned long r;
		unsigned long c;
		double timer = 0;
		/*timer start*/
		timer -= omp_get_wtime();
		/*pararell compute*/
#pragma omp parallel shared(timer,chunk) private(r,c)
		{
		    /*compute*/
#pragma omp for schedule (static,chunk)
		    for (r = 0; r < this->rows; r++) {
			this->p[r] = new double[this->cols];
			for (c = 0; c < this->cols; c++) {
			    this->p[r][c] = 0;
			}
		    }
		}
		/*timer stop*/
		timer += omp_get_wtime();
		para.printTimer(timer);
	    }
	}
    }

    /*
     * Create a Matrix object as copy of a Matrix A
     */
    Matrix(const Matrix<MxType>& A) {
	if (para.hasThreads()) {
	    this->rows = A.rows;
	    this->cols = A.cols;
	    this->p = new double*[A.rows];
	    /*pararell config*/
	    unsigned long chunk = para.calcChunk(this->rows + this->cols);
	    unsigned long r;
	    unsigned long c;
	    double timer = 0;
	    /*timer start*/
	    timer -= omp_get_wtime();
	    /*pararell compute*/
#pragma omp parallel shared(timer) private(r,c)
	    {
		/*compute*/
#pragma omp for schedule (static,chunk)
		for (r = 0; r < A.rows; r++) {
		    this->p[r] = new double[A.cols];
		    /*copy the values from the matrix A*/
		    for (c = 0; c < A.cols; c++) {
			this->p[r][c] = A.p[r][c];
		    }
		}
	    }
	    /*timer stop*/
	    timer += omp_get_wtime();
	    para.printTimer(timer);
	}
    }

    /*
     * Randomize contents of the matrix
     */
    void Randomize() {
	if (para.hasThreads()) {
	    /*pararell config*/
	    unsigned long chunk = para.calcChunk(this->rows + this->cols);
	    unsigned long r;
	    unsigned long c;
	    double timer = 0;
	    /*timer starts*/
	    timer -= omp_get_wtime();
	    /*pararell compute*/
#pragma omp parallel shared(timer) private(r,c)
	    {
		/*initialize mt19937*/
		init_genrand(omp_get_thread_num() * timer);
		/*compute*/
#pragma omp for schedule (static,chunk)
		for (r = 0; r < this->rows; r++) {
		    /*copy the values from the matrix A*/
		    for (c = 0; c < this->cols; c++) {
			this->p[r][c] = genrand_real1();
		    }
		}
	    }
	    /*timer stop*/
	    timer += omp_get_wtime();
	    para.printTimer(timer);
	}
    }

    /*
     * Destructor - clean up allocated memory
     */
    ~Matrix() {
	for (unsigned long r = 0; r < this->rows; r++) {
	    delete[] this->p[r];
	}
	delete this->p;
	this->p = NULL;
    }

    /***************************************************************************
     * Operations
     */

    /*
     * Add a double value (elements wise)
     */
    Matrix<MxType>& Add(const double no) {
	if (para.hasThreads()) {
	    /*pararell config*/
	    unsigned long chunk = para.calcChunk(this->rows + this->cols);
	    unsigned long r;
	    unsigned long c;
	    double timer = 0;
	    /*timer start*/
	    timer -= omp_get_wtime();
	    /*pararell compute*/
#pragma omp parallel shared(timer) private(r,c)
	    {
		/*compute*/
#pragma omp for schedule (static,chunk)
		for (r = 0; r < this->rows; r++) {
		    for (c = 0; c < this->cols; c++) {
			this->p[r][c] += no;
		    }
		}
	    }
	    /*timer stop*/
	    timer += omp_get_wtime();
	    para.printTimer(timer);
	    return *this;
	}
    }

    /*
     * Subtract a double value (elements wise)
     */
    Matrix<MxType> & Subtract(const double no) {
	return Add(-no);
    }

    /*
     * Multiply a double value (elements wise)
     */
    Matrix<MxType> & Multiply(const double no) {
	if (para.hasThreads()) {
	    /*pararell config*/
	    unsigned long chunk = para.calcChunk(this->rows + this->cols);
	    unsigned long r;
	    unsigned long c;
	    double timer = 0;
	    /*timer start*/
	    timer -= omp_get_wtime();
	    /*pararell compute*/
#pragma omp parallel shared(timer) private(r,c)
	    {
		/*compute*/
#pragma omp for schedule (static,chunk)
		for (r = 0; r < this->rows; r++) {
		    for (c = 0; c < this->cols; c++) {
			this->p[r][c] *= no;
		    }
		}
	    }
	    /*timer stop*/
	    timer += omp_get_wtime();
	    para.printTimer(timer);
	    return *this;
	}
    }

    /*
     * Divide a double value (elements wise)
     */
    Matrix<MxType> & Divide(const double no) {
	return Multiply(1 / no);
    }

    /***************************************************************************
     * Operators
     */

    /*
     * Index operator. You can use this class like myMatrix(col, row)
     * The indexes are one-based, not zero based!
     */
    double& operator()(const unsigned long r, const unsigned long c) {
	if (this->p != NULL && r > 0 && r <= this->rows && c > 0 && c <= this->cols) {
	    return this->p[r - 1][c - 1];
	} else {
	    throw Exception("Subscript out of range");
	}
    }

    /*
     * Index operator. You can use this class like myMatrix.get(col, row)
     * The indexes are one-based, not zero based!
     * Use this function get if you want to read from a const Matrix.
     */
    double get(const unsigned long r, const unsigned long c) const {
	if (this->p != NULL && r > 0 && r <= this->rows && c > 0 && c <= this->cols) {
	    return this->p[r - 1][c - 1];
	} else {
	    throw Exception("Subscript out of range");
	}
    }

    /*
     * Assignment operator
     */

    Matrix<MxType>& operator=(const Matrix<MxType> & A) {
	if (para.hasThreads()) {
	    this->rows = A.rows;
	    this->cols = A.cols;
	    this->p = new double*[A.rows];
	    /*pararell config*/
	    unsigned long chunk = para.calcChunk(this->rows + this->cols);
	    unsigned long r;
	    unsigned long c;
	    double timer = 0;
	    /*timer start*/
	    timer -= omp_get_wtime();
	    /*pararell compute*/
#pragma omp parallel shared(timer) private(r,c)
	    {
		/*compute*/
#pragma omp for schedule (static,chunk)
		for (r = 0; r < A.rows; r++) {
		    this->p[r] = new double[A.cols];
		    // copy the values from the matrix a
		    for (c = 0; c < A.cols; c++) {
			this->p[r][c] = A.p[r][c];
		    }
		}
	    }
	    /*timer stop*/
	    timer += omp_get_wtime();
	    para.printTimer(timer);
	    return *this;
	}
    }


    /***************************************************************************
     * Operators: Addition
     */

    /*
     * Addition of Matrix with Matrix
     */
    friend Matrix<MxType> operator+(const Matrix<MxType>& A, const Matrix<MxType> & B) {
	if (para.hasThreads()) {
	    /*check if the dimensions match*/
	    if (A.rows == B.rows && A.cols == B.cols) {
		Matrix<MxType> res(A.rows, A.cols);
		/*pararell config*/
		unsigned long chunk = para.calcChunk(res.rows + res.cols);
		unsigned long r;
		unsigned long c;
		double timer = 0;
		/*timer start*/
		timer -= omp_get_wtime();
		/*pararell compute*/
#pragma omp parallel shared(timer,res) private(r,c)
		{
		    /*compute*/
#pragma omp for schedule (static,chunk)
		    for (r = 0; r < A.rows; r++) {
			for (c = 0; c < A.cols; c++) {
			    res.p[r][c] = A.p[r][c] + B.p[r][c];
			}
		    }
		}
		/*timer stop*/
		timer += omp_get_wtime();
		para.printTimer(timer);
		return res;
	    } else {
		throw Exception("Dimensions does not match");
	    }
	}
	return Matrix(); // return an empty matrix (this never happens but just for safety)
    }

    /*
     * Addition of Matrix with number
     */
    friend Matrix<MxType> operator+(const Matrix<MxType>& A, const double b) {
	Matrix<MxType> res = A;
	res.Add(b);
	return res;
    }

    /*
     * Addition of number with Matrix
     */
    friend Matrix<MxType> operator+(const double b, const Matrix<MxType> & A) {
	Matrix<MxType> res = A;
	res.Add(b);
	return res;
    }

    /***************************************************************************
     * Operators: Substraction
     */

    /*
     * Substraction of Matrix with Matrix
     */
    friend Matrix<MxType> operator-(const Matrix<MxType>& A, const Matrix<MxType> & B) {
	if (para.hasThreads()) {
	    // check if the dimensions match
	    if (A.rows == B.rows && A.cols == B.cols) {
		Matrix<MxType> res(A.rows, A.cols);
		/*pararell config*/
		unsigned long chunk = para.calcChunk(res.rows + res.cols);
		unsigned long r;
		unsigned long c;
		double timer = 0;
		/*timer start*/
		timer -= omp_get_wtime();
		/*pararell compute*/
#pragma omp parallel shared(timer,res) private(r,c)
		{
		    /*compute*/
#pragma omp for schedule (static,chunk)
		    for (r = 0; r < A.rows; r++) {
			for (c = 0; c < A.cols; c++) {
			    res.p[r][c] = A.p[r][c] - B.p[r][c];
			}
		    }
		}
		/*timer stop*/
		timer += omp_get_wtime();
		para.printTimer(timer);
		return res;
	    } else {
		throw Exception("Dimensions does not match");
	    }
	}
	return Matrix(); // return an empty matrix (this never happens but just for safety)
    }

    /*
     * Substraction of Matrix with number
     */
    friend Matrix<MxType> operator-(const Matrix<MxType>& A, const double b) {
	Matrix<MxType> res = A;
	res.Subtract(b);
	return res;
    }

    /*
     * Substraction of number with Matrix
     */
    friend Matrix<MxType> operator-(const double b, const Matrix<MxType> & A) {
	Matrix<MxType> res = -A;
	res.Add(b);
	return res;
    }

    /*
     * Operator unary minus
     */
    friend Matrix<MxType> operator-(const Matrix<MxType> & A) {
	if (para.hasThreads()) {
	    Matrix<MxType> res(A.rows, A.cols);
	    /*pararell config*/
	    unsigned long chunk = para.calcChunk(res.rows + res.cols);
	    unsigned long r;
	    unsigned long c;
	    double timer = 0;
	    /*timer start*/
	    timer -= omp_get_wtime();
	    /*pararell compute*/
#pragma omp parallel shared(timer,res) private(r,c)
	    {
		/*compute*/
#pragma omp for schedule (static,chunk)
		for (r = 0; r < A.rows; r++) {
		    for (c = 0; c < A.cols; c++) {
			res.p[r][c] = -A.p[r][c];
		    }
		}
	    }
	    /*timer stop*/
	    timer += omp_get_wtime();
	    para.printTimer(timer);
	    return res;
	}
    }


    /***************************************************************************
     * Operators: Multiplication
     */

    /*
     * Multiplication of Matrix with Matrix
     */
    friend Matrix<MxType> operator*(const Matrix<MxType>& A, const Matrix<MxType> & B) {
	if (para.hasThreads()) {
	    // check if the dimensions match
	    if (A.cols == B.rows) {
		Matrix<MxType> res(A.rows, B.cols);
		/*pararell config*/
		unsigned long chunk = para.calcChunk(A.rows + B.cols + A.cols);
		unsigned long r;
		unsigned long c;
		unsigned long c_res;
		double timer = 0;
		/*timer start*/
		timer -= omp_get_wtime();
		/*pararell compute*/
#pragma omp parallel shared(timer,res) private(r,c,c_res)
		{
		    /*compute*/
#pragma omp for schedule (static,chunk)
		    for (r = 0; r < A.rows; r++) {
			for (c_res = 0; c_res < B.cols; c_res++) {
			    for (c = 0; c < A.cols; c++) {
				res.p[r][c_res] += A.p[r][c] * B.p[c][c_res];
			    }
			}
		    }
		}
		/*timer stop*/
		timer += omp_get_wtime();
		para.printTimer(timer);
		return res;
	    } else {
		throw Exception("Dimensions does not match");
	    }
	}
	return Matrix(); // return an empty matrix (this never happens but just for safety)
    }

    /*
     * Multiplication of Matrix with number
     */
    friend Matrix<MxType> operator*(const Matrix<MxType>& A, const double b) {
	Matrix<MxType> res = A;
	res.Multiply(b);
	return res;
    }

    /*
     * Multiplication of number with Matrix
     */
    friend Matrix<MxType> operator*(const double b, const Matrix<MxType> & A) {
	Matrix<MxType> res = A;
	res.Multiply(b);
	return res;
    }

    /***************************************************************************
     * Operators: Division
     */

    /*
     * Division of Matrix with Matrix
     */
    friend Matrix<MxType> operator/(const Matrix<MxType>& A, const Matrix<MxType> & B) {
	if (para.hasThreads()) {
	    /*check if the dimensions match: must be square and equal sizes*/
	    if (A.rows == A.cols && A.rows == A.cols && B.rows == B.cols) {
		Matrix<MxType> res(A.rows, A.cols);
		res = A * Inv(B);
		return res;
	    } else {
		throw Exception("Dimensions does not match");
	    }
	}
	return Matrix(); // return an empty matrix (this never happens but just for safety)
    }

    /*
     * Division of Matrix with number
     */
    friend Matrix<MxType> operator/(const Matrix<MxType>& A, const double b) {
	Matrix<MxType> res = A;
	res.Divide(b);
	return res;
    }

    /*
     * Division of number with Matrix
     */
    friend Matrix<MxType> operator/(const double b, const Matrix<MxType> & A) {
	Matrix<MxType> b_matrix(1, 1);
	b_matrix(1, 1) = b;

	Matrix<MxType> res = b_matrix / A;
	return res;
    }

    /***************************************************************************
     * Operations
     */

    /*
     * Returns the minor from the given matrix where the selected row and column are removed
     */
    Matrix<MxType> Minor(const unsigned long row, const unsigned long col) const {
	if (para.hasThreads()) {
	    Matrix<MxType> res;
	    if (row > 0 && row <= this->rows && col > 0 && col <= this->cols) {
		res = Matrix<MxType > (this->rows - 1, this->cols - 1);
		/*pararell config*/
		unsigned long chunk = para.calcChunk((this->rows - (row >= this->rows)) + (this->cols - (col >= this->cols)));
		unsigned long r;
		unsigned long c;
		double timer = 0;
		/*timer start*/
		timer -= omp_get_wtime();
		/*pararell compute*/
#pragma omp parallel shared(timer,res) private(r,c)
		{
		    /*compute*/
#pragma omp for schedule (static,chunk)
		    /*copy the content of the matrix to the minor, except the selected*/
		    for (r = 1; r <= (this->rows - (row >= this->rows)); r++) {
			for (c = 1; c <= (this->cols - (col >= this->cols)); c++) {
			    res(r - (r > row), c - (c > col)) = this->p[r - 1][c - 1];
			}
		    }
		}
		/*timer stop*/
		timer += omp_get_wtime();
		para.printTimer(timer);
	    } else {
		throw Exception("Index for minor out of range");
	    }
	    return res;
	}
    }

    /*
     * Returns the size of the i-th dimension of the matrix.
     * for i=1 the function Returns the number of rows
     * for i=2 the function Returns the number of columns
     * for else the function Returns 0
     */
    unsigned long Size(const unsigned long i) const {
	if (i == 1) {
	    return this->rows;
	} else if (i == 2) {
	    return this->cols;
	}
	return 0;
    }

    /*
     * Returns the number of rows
     */
    unsigned long GetRows() const {
	return this->rows;
    }

    /*
     * Returns the number of cols
     */
    unsigned long GetCols() const {
	return this->cols;
    }

    /*
     * Prints the contents of the matrix
     */
    void Print() const {
	if (this->p != NULL) {
	    unsigned long r;
	    unsigned long c;
	    cout << "[";
	    for (r = 0; r < this->rows; r++) {
		if (r > 0) {
		    cout << " ";
		}
		for (c = 0; c < this->cols; c++) {
		    printf("%.2f", this->p[r][c]);
		    if (c < this->cols - 1) {
			cout << ", ";
		    }
		}
		if (r < this->rows - 1) {
		    cout << "; " << endl;
		}
	    }
	    cout << "]" << endl;
	} else {
	    cout << "[ ]" << endl;
	}
    }

private:
    /*matrix data*/
    unsigned long rows;
    unsigned long cols;
    MxType** p; // pointer to a matrix with doubles
};

/*
 * Returns the size of the i-th dimension of the matrix.
 * for i=1 the function Returns the number of rows
 * for i=2 the function Returns the number of columns
 * for else the function Returns 0
 */
template <typename MxType> unsigned long Size(const MxType& A, const unsigned long i) {
    return A.Size(i);
}

/*
 * Returns a matrix with size cols x rows with ones as values
 * @param   rows of matrix
 * @param   cols of matrix
 * @return  a ones matrix
 */
template <typename MxType> Matrix<MxType> Ones(const unsigned long rows, const unsigned long cols) {
    Matrix<MxType> res = Matrix<MxType > (rows, cols);
    /*pararell config*/
    unsigned long chunk = para.calcChunk(rows + cols);
    unsigned long r;
    unsigned long c;
    double timer = 0;
    /*timer start*/
    timer -= omp_get_wtime();
    /*pararell compute*/
#pragma omp parallel shared(timer,res) private(r,c)
    {
	/*compute*/
#pragma omp for schedule (static,chunk)
	for (r = 1; r <= rows; r++) {
	    for (c = 1; c <= cols; c++) {
		res(r, c) = 1;
	    }
	}
    }
    /*timer stop*/
    timer += omp_get_wtime();
    para.printTimer(timer);
    return res;
}

/**
 * Returns a matrix with size cols x rows with zeros as values
 * @param   rows of matrix
 * @param   cols of matrix
 * @return  a zeros matrix
 */
template <typename MxType> Matrix<MxType> Zeros(const unsigned long rows, const unsigned long cols) {
    return Matrix<MxType > (rows, cols);
}

/**
 * Returns a diagonal matrix with size n x n with ones at the diagonal
 * @param   n a size of matrix
 * @return  a diagonal matrix with ones on the diagonal
 */
template <typename MxType> Matrix<MxType> Diag(const unsigned long n) {
    Matrix<MxType> res = Matrix<MxType > (n, n);
    /*pararell config*/
    unsigned long chunk = para.calcChunk(n);
    unsigned long i;
    double timer = 0;
    /*timer start*/
    timer -= omp_get_wtime();
    /*pararell compute*/
    //#pragma omp parallel shared(timer) private(i)
    //    {
    /*compute*/
    //#pragma omp for schedule (static,chunk)
    for (i = 1; i <= n; i++) {
	res(i, i) = 1;
    }
    //    }
    /*timer stop*/
    timer += omp_get_wtime();
    para.printTimer(timer);
    return res;
}

/**
 * Returns a diagonal matrix
 * @param   V a vector
 * @return  a diagonal matrix with the given vector v on the diagonal
 */
template <typename MxType> Matrix<MxType> Diag(const Matrix<MxType>& V) {
    Matrix<MxType> res;
    if (V.GetCols() == 1) {
	/*the given matrix is a vector n x 1*/
	unsigned long rows = V.GetRows();
	res = Matrix<MxType > (rows, rows);
	/*pararell config*/
	unsigned long chunk = para.calcChunk(rows);
	unsigned long r;
	double timer = 0;
	/*timer start*/
	timer -= omp_get_wtime();
	/*pararell compute*/
	//#pragma omp parallel shared(timer,res,rows) private(r)
	//	{
	/*compute*/
	//#pragma omp for schedule (static,chunk)
	/*copy the values of the vector to the matrix*/
	for (r = 1; r <= rows; r++) {
	    res(r, r) = V.get(r, 1);
	}
	//	}
	/*timer stop*/
	timer += omp_get_wtime();
	para.printTimer(timer);
    } else if (V.GetRows() == 1) {
	/*the given matrix is a vector 1 x n*/
	unsigned long cols = V.GetCols();
	res = Matrix<MxType > (cols, cols);
	/*pararell config*/
	unsigned long chunk = para.calcChunk(cols);
	unsigned long c;
	double timer = 0;
	/*timer start*/
	timer -= omp_get_wtime();
	/*pararell compute*/
	//#pragma omp parallel shared(timer,res,cols) private(c)
	//	{
	/*compute*/
	//#pragma omp for schedule (static,chunk)
	/*copy the values of the vector to the matrix*/
	for (c = 1; c <= cols; c++) {
	    res(c, c) = V.get(1, c);
	}
	//	}
	/*timer stop*/
	timer += omp_get_wtime();
	para.printTimer(timer);
    } else {
	throw Exception("Parameter for diag must be a vector");
    }
    return res;
}

/*
 * Returns the determinant of Matrix A
 * Compute non pararell due Minor is computed pararell
 */
template <typename MxType> double Det(const Matrix<MxType>& A) {
    double d = 0; // value of the determinant
    unsigned long rows = A.GetRows();
    unsigned long cols = A.GetRows();
    if (rows == cols) {
	/*this is a square matrix*/
	if (rows == 1) {
	    /*this is a 1 x 1 matrix*/
	    d = A.get(1, 1);
	} else if (rows == 2) {
	    /*this is a 2 x 2 matrix*/
	    /*the determinant of [a11,a12;a21,a22] is det = a11*a22-a21*a12*/
	    d = A.get(1, 1) * A.get(2, 2) - A.get(2, 1) * A.get(1, 2);
	} else {
	    /*this is a matrix of 3 x 3 or larger*/
	    for (unsigned long c = 1; c <= cols; c++) {
		Matrix<MxType> M = A.Minor(1, c);
		/*d += pow(-1, 1+c)  * a(1, c) * Det(M);*/
		d += (c % 2 + c % 2 - 1) * A.get(1, c) * Det(M); /*faster than with pow()*/
	    }
	}
    } else {
	throw Exception("Matrix must be square");
    }
    return d;
}

/*
 * Swap two values
 */
void Swap(double& a, double& b) {
    double temp = a;
    a = b;
    b = temp;
}

/*
 * Returns the inverse of Matrix A
 */
template <typename MxType> Matrix<MxType> Inv(const Matrix<MxType>& A) {
    Matrix<MxType> res;
    double d = 0; // value of the determinant
    unsigned long rows = A.GetRows();
    unsigned long cols = A.GetRows();
    d = Det(A);
    if (rows == cols && d != 0) {
	// this is a square matrix
	if (rows == 1) {
	    // this is a 1 x 1 matrix
	    res = Matrix<MxType > (rows, cols);
	    res(1, 1) = 1 / A.get(1, 1);
	} else if (rows == 2) {
	    // this is a 2 x 2 matrix
	    res = Matrix<MxType > (rows, cols);
	    res(1, 1) = A.get(2, 2);
	    res(1, 2) = -A.get(1, 2);
	    res(2, 1) = -A.get(2, 1);
	    res(2, 2) = A.get(1, 1);
	    res = (1 / d) * res;
	} else {
	    // this is a matrix of 3 x 3 or larger
	    // calculate inverse using gauss-jordan elimination
	    //      http://mathworld.wolfram.com/MatrixInverse.html
	    //      http://math.uww.edu/~mcfarlat/inverse.htm
	    res = Diag<MxType > (rows); // a diagonal matrix with ones at the diagonal
	    Matrix<MxType> ai = A; // make a copy of Matrix a
	    /*pararell config*/
	    unsigned long chunk = para.calcChunk(cols + rows + cols + rows + cols);
	    unsigned long r;
	    unsigned long c;
	    unsigned long s;
	    double f;
	    double timer = 0;
	    /*timer start*/
	    timer -= omp_get_wtime();
	    /*pararell compute*/
#pragma omp parallel shared(timer,res,ai,rows,cols) private(r,c,s,f)
	    {
		/*compute*/
#pragma omp for schedule (static,chunk)
		for (c = 1; c <= cols; c++) {
		    /*element (c, c) should be non zero. if not, swap content of lower rows*/
		    for (r = c; r <= rows && ai(r, c) == 0; r++) {
		    }
		    if (r != c) {
			/*swap rows*/
			for (s = 1; s <= cols; s++) {
			    Swap(ai(c, s), ai(r, s));
			    Swap(res(c, s), res(r, s));
			}
		    }
		    /*eliminate non-zero values on the other rows at column c*/
		    for (r = 1; r <= rows; r++) {
			if (r != c) {
			    /*eleminate value at column c and row r*/
			    if (ai(r, c) != 0) {
				f = -ai(r, c) / ai(c, c);
				/*add (f * row c) to row r to eleminate the value at column c*/
				for (s = 1; s <= cols; s++) {
				    ai(r, s) += f * ai(c, s);
				    res(r, s) += f * res(c, s);
				}
			    }
			} else {
			    /*make value at (c, c) one, divide each value on row r with the value at ai(c,c)*/
			    f = ai(c, c);
			    for (unsigned long s = 1; s <= cols; s++) {
				ai(r, s) /= f;
				res(r, s) /= f;
			    }
			}
		    }
		}
	    }
	    /*timer stop*/
	    timer += omp_get_wtime();
	    para.printTimer(timer);
	}
    } else {
	if (rows == cols) {
	    throw Exception("Matrix must be square");
	} else {
	    throw Exception("Determinant of matrix is zero");
	}
    }
    return res;
}


#endif	/* MATRIXCLASS_H */