Inverse Matrix and Plane Equation

Zadkiel · June 22, 2001, 2:03pm

Hi, I need help solving these two maths problems.

First, Inverse 4*4 Matrix

I need code or equations or just a explanation of how to get the inverse of a 44 matrix. I do know that not all 44 matrices have an inverse, but if someone can give me pointers on how to find the inverse for the general case I’ll make sure I never use a matrix which wont work.

Second, How do I get the plane equation form 3 points.

I think theres a method whereby you calculate the normal of the plane (could someone post the code for that, I lost it)and then using that along with one of the points you make a plane equation.

Can anyone help me with either that or another method.

Thanks in advance
Zadkiel

imported_john · June 22, 2001, 5:40pm

Howdy,

answering your second question first, then;

two ways of finding the equation of a plane from three points. The way I am most confident with answering is to find the determinant of a 4x4 matrix:

|  i  |
| pt1 |
| pt2 |
| pt3 |

where i is the ‘parameter’ vector [ x y z w ] and pt1…3 are the three points in homogeneous coordinates. The determinant of this matrix is a 1x4 vector representing the plane in homogeneous coordinates, which just so happens to also be the co-efficients to

Ax+By+Cz+D = 0

that’s method one.

The second method:

suppose you have the normal to your plane (n) and some point that lies on this plane §. Well, the equation of the plane

Ax+By+Cz+D = 0

has the normal encoded as <A B C>, so since we know the normal, then we immediatly have A=n[0], B=n[1], C=n[2]. To find D, all you need to do is subsititute x, y, z for the point that lies on the plane and solve for D.

The question, then, is how to find the normal. Well, that’s the cross product of two vectors that lie on the plane, which you can take to be (p1…p2) and (p1…p3).

In answer to yoru FIRST question… well… hmm. I’m not sure of a good computer way of doing this… I remember computing inverses from 2x2 by finding the determinant… hmm. The more general way is to use gauss-jordan elimination on the matrix to get the identity matrix, and apply the same ops to the identity matrix. When your orig matrix is the identiy, your former identity matrix is the inverse… but how to do this in code is tricky, i’d imagine =)

you could try looking at some maths packages. There are tonnes of them around; I use COOL/IUE Numerics which is bundled with target junior. It’s free.

As an aside, tho’, don’t forget that the inverse for the projection matricies are given to you in the red book =)

cheers,
John

ffish · June 22, 2001, 6:23pm

Here’s something I always wanted to try but haven’t gotten around to yet - Numerical Recipes in C has a method for matrix multiplication and inversion that is O(N^logbase2(7)) instead of O(N^3). Not much difference but mmmmmm, maths! The link is: http://www.ulib.org/webRoot/Books/Numerical_Recipes/bookcpdf/c2-11.pdf
and the parent node of that link (bookcpdf.html) has all you will need for matrix inversion. At a guess, look at chapters 2.1, 2.2, 2.3 and maybe 2.10. I’m sure inversion is mentioned in there once or twice.

There are several methods of calculating matrix inverses (I know of four) but (apart from the link) they are all O(N^3) so you can use any one, especially for 4x4 matrices (one in particular is elegant and easy to code but recursive - it’s as fast as the others but if you try it on large matrices it takes forever and runs out of memory)

Zadkiel · June 23, 2001, 12:38am

Thanks John and ffish.

John I going to use the first method you suggested for the plane equation, it will suit my program perfectly. However, how do I calculate the detriment of a 44 matrix. I only know 22.

ffish, I’m looking up the pdf to suggested, however you also say you know of 3 other ways, could you post a link to them just in cae they may fit my prograsm better than this one.

I just found a good inverse calculationg routine, says its written specificaslly of OpenGL matrices; http://www.cs.unc.edu/~gotz/code/affinverse.html

But thanks both of you for the help, I was stuck at this bit of the program.

[This message has been edited by Zadkiel (edited 06-23-2001).]

ffish · June 23, 2001, 2:50am

(i) The first is Gauss-Jordan elimination. I presume you know how to do that. It’s fairly easy to code.

(ii) The second uses the formula A^-1(ij) = C(ji)/det(A) or A^-1 = C’/det(A) where A^-1 is the inverse matrix of A, i and j are subscripts and det(A) is the determinant of A. C is the “cofactor matrix” and C(ij) is the ith row and jth column of C (note the formula above has (ji) as the subscript).

The easiest way to explain how to calculate the inverse according to this formula is like this: To compute A(ij), cross out row j and column i of A. Multiply the determinant of the resulting matrix by -1^(i+j) to get the cofactor, then divide by det(A) for the result. It’s easiest to do examples with the formula with 3x3 matrices, but I can verify it works with all nxn matrices. I find this method the easiest and quickest to do on paper for small matrices.

(iii) The third method is one that I got from a parallel programming unit. Looking at it now, it may be the method from Numerical Recipes in C. Anyway, I’ll post all the code from my implementation here for everyone’s benefit to look at if they want.

(iv) Numerical Recipes method - same as (iii)?

I’d do things quite differently if I was to implement this file again (in fact I have, and I have ). You’ll have to get rid of some of the lines at the top regarding MPI (parallel programming API) but I’ll leave them in for completeness. I don’t even know if this code will compile as is, but I think it should from a quick glance:

/* Must be defined at the top for documentation purposes. */
#ifndef _MATRIX_H
#define _MATRIX_H

/** @name pmatrix.h

@type Interface file

@memo pmatrix.h - template file for a generic Matrix class,

including lots of legacy functions that may not be used for this

assignment, and lots of extra functions that I may have used at

some stage but not in the end.

@version 0.2 (27/7/99) Modified: 20/4/2000
*/

//@{

extern “C” {
#include <mpi.h>
}

#include
#include
#include
#include
#include
#include
#include
#include
#include “pmatrix.h”

/** Constant decision parameter for Matrix constructors.

@type const int /
const int IDENTITY=0;
/* Constant decision parameter for Matrix constructors.

@type const int /
const int FULL=1;
/* Constant printing parameter.

@type const int /
const int PRECISION=3;
/* Constant default Matrix size.

@type const int */
const int SIZE=5;

// Namespace for various MPI related variables.
namespace globals {
int myid;
int numprocs;
MPI_Comm workers, world(MPI_COMM_WORLD);
MPI_Group world_group, worker_group;
MPI_Status msgStatus;
}

// Class template forward declaration.
template
class Matrix;

// Declare external vars from test3.cpp.
//test3::myid;
//test3::numprocs;

// Function forward declarations for friend functions.
template
const Matrix operator+(const Matrix& m1, const Matrix& m2);

template
const Matrix operator-(const Matrix& m1, const Matrix& m2);

template
const Matrix operator*(const Matrix& m1, const Matrix& m2);

template
const vector operator*(const vector& v, const Matrix& m);

template
const vector operator*(const Matrix& m, const vector& v);

template
const Matrix operator*(const T& s, const Matrix& m);

template
const Matrix operator*(const Matrix& m, const T& s);

template
const Matrix operator*(const int& s, const Matrix& m);

template
const Matrix operator*(const Matrix& m, const int& s);

template
ostream& operator<<(ostream& out, const Matrix& m);

template
const Matrix operator-(const Matrix& m);

template
void swap_rows(Matrix *m, int row1, int row2);

template
void sort_rows(Matrix *m, int col);

template
const vector init_vector(const Matrix& m);

template
const vector row(int r, const Matrix& m);

template
const vector col(int c, const Matrix& m);

template
T dot_product(const vector& v1, const vector& v2);

template
const Matrix Gauss_Jordan_inverse(const Matrix& m);

template
const Matrix transpose(const Matrix& m);

template
const Matrix subMatrix(const Matrix& m1, int r, int c);

template
const Matrix subMatrix(const Matrix& m1, int r1, int c1, int r2,
int c2);

template
const T det(const Matrix& m);

template
const T cofactor(const Matrix& m, int r, int c);

template
const Matrix cofactor_Matrix(const Matrix& m1);

template
const Matrix cofactor_inverse(const Matrix& m1);

template
const Matrix operator/(const Matrix& m, const T& s);

template
const Matrix MPI_Notes_inverse(const Matrix& m1);

/** Declaration of the Matrix class. This Matrix class is by no means
complete. It serves one purpose only, to implement a floating

point benchmarking test. As such it only contains minimal
functionality.
/
template
class Matrix {
public:
/* Matrix constructor.
* @param init is a decision parameter. Valid values are 1 for a
* random Matrix of floats, otherwise the identity. Pass
* one of the const ints defined as FULL=1 or IDENTITY=0.
* The default is IDENTITY.
* @param s is the row and column size of the Matrix. The default
* is defined above and can be changed.
* @return a Matrix object.
*/
Matrix(const int init=IDENTITY, const int s=SIZE);

/** Matrix destructor.
*/
~Matrix();

/** Matrix copy constructor.
* @param m2 is the Matrix to copy.
* @return a Matrix object.
*/
Matrix(const Matrix& m2);

/** Overloaded Matrix constructor.
* @param a is an array of T’s.
* @return a Matrix object.
/
Matrix(const T a, const int& s);

/** Return the Matrix size.
* @return the size of the Matrix. i.e if the Matrix is n * n,
 	returns n.
*/
inline const int get_size() { return size; };

/** Overloaded assignment operator for the Matrix class.
* @type Matrix&
* @param m2 is the Matrix to copy.
* @return a Matrix that is identical to m2.
*/
Matrix& operator=(const Matrix& m2);

/** Overloaded addition operator for the Matrix class.
* @type friend const Matrix
* @param m1 is the lhs Matrix to add with.
* @param m2 is the rhs Matrix to add with.
* @return a Matrix which is the sum of m1 and m2.
*/
friend const Matrix operator+ <>(const Matrix& m1,
const Matrix& m2);

/** Overloaded subtraction operator for the Matrix class.
* @type friend const Matrix
* @param m1 is the lhs Matrix to subtract from.
* @param m2 is the rhs Matrix to subtract with.
* @return a Matrix which equals m1 - m2.
*/
friend const Matrix operator- <>(const Matrix& m1,
const Matrix& m2);

/** Overloaded multiplication operator for the Matrix class.
* @type friend const Matrix
* @param m1 is the lhs Matrix to multiply with.
* @param m2 is the rhs Matrix to be multiplied against.
* @return a Matrix which is the product of m1 and m2.
/
friend const Matrix operator (const Matrix& m1,
const Matrix& m2);

/** Overloaded multiplication operator for the Matrix class.
* @type friend const vector
* @param v is the vector to multiply with.
* @param m is the Matrix to be multiplied against.
* @return a vector which is the product of v and m.
/
friend const vector operator (const vector& v,
const Matrix& m);

/** Overloaded multiplication operator for the Matrix class.
* @type friend const vector
* @param m is the Matrix to multiply with.
* @param v is the vector to be multiplied against.
* @return a vector which is the product of m and v.
/
friend const vector operator (const Matrix& m,
const vector& v);

/** Overloaded scalar multiplication operator for the Matrix
* class.
* @type friend const Matrix
* @param s is the scalar to multiply with.
* @param m1 is the Matrix to be multiplied against.
* @return a Matrix which is the product of s and m.
/
friend const Matrix operator (const int& s,
const Matrix& m);

/** Overloaded scalar multiplication operator for the Matrix
* class.
* @type friend const Matrix
* @param m is the Matrix to be multiplied against.
* @param s is the scalar to multiply with.
* @return a Matrix which is the product of s and m.
/
friend const Matrix operator (const Matrix& m,
const int& s);

/** Overloaded scalar multiplication operator for the Matrix
* class.
* @type friend const Matrix
* @param s is the scalar to multiply with.
* @param m1 is the Matrix to be multiplied against.
* @return a Matrix which is the product of s and m.
/
friend const Matrix operator (const T& s,
const Matrix& m);

/** Overloaded scalar multiplication operator for the Matrix
* class.
* @type friend const Matrix
* @param m is the Matrix to be multiplied against.
* @param s is the scalar to multiply with.
* @return a Matrix which is the product of s and m.
/
friend const Matrix operator (const Matrix& m,
const T& s);

/** Overloaded scalar division operator for the Matrix
* class.
* @type friend const Matrix
* @param m1 is the Matrix to be divided against.
* @param s is the scalar to divide with.
* @return a Matrix which is the elementwise division of m1 by s.
*/
friend const Matrix operator/ (const Matrix& m,
const T& s);

/** Overloaded output operator for the Matrix class. Allows
* printing of a Matrix object.
* @type friend ostream&
* @param out is the stream to write onto.
* @param m is the Matrix to write.
* @return an ostream reference.
*/
friend ostream& operator<< (ostream& out, const Matrix& m);

/** Overloaded unary negation operator for the Matrix class.
* @type friend const Matrix
* @param m1 is the Matrix to negate.
* @return a Matrix which is the negative of m1.
*/
// friend const Matrix operator- (const Matrix& m1);

/** Swap row1 of the Matrix m into row2 and swap row2 into row1.
* @type friend void
* @param m is a pointer to a Matrix.
* @param row1 is the index of the first row.
* @param row2 is the index of the second row.
* @return void
*/
friend void swap_rows(Matrix *m, int row1, int row2);

/** A slow, inefficient sort algorithm.
* @type friend void
* @param m is a pointer to the Matrix to sort rows on.
* @param col is the column to sort on.
*/
friend void sort_rows(Matrix *m, int col);

/** Initializes a vector that has the same number of rows, cols
* as the Matrix m.
* @type friend const vector
* @param m is a reference to the Matrix used to initialize the
* vector.
* @return a vector with the same size as m
*/
friend const vector init_vector(const Matrix& m);

/** Returns a row from the Matrix m, where the rows go from 0 to
* m.size - 1.

@type friend const vector

@param r is the row to choose.

@param m is the Matrix to choose the row from.

@return a row vector from m.
*/
friend const vector row(int r, const Matrix& m);

/** Returns a col from the Matrix m, where the cols go from 0 to
* m.size - 1.

@type friend const vector

@param c is the col to choose.

@param m is the Matrix to choose the col from.

@return a col vector from m.
*/
friend const vector col(int c, const Matrix& m);

/** Computes the inverse of a Matrix using the Gauss-Jordan
* method.
* @type friend const Matrix
* @param m is the Matrix to find the inverse of.
* @return the inverse of m.
*/
friend const Matrix Gauss_Jordan_inverse
(const Matrix& m);

/** Computes the transpose of a Matrix.
* @type friend const Matrix
* @param m is the Matrix to find the transpose of.
* @return the transpose of m.
*/
friend const Matrix transpose(const Matrix& m);

/** Computes a sub-Matrix from a Matrix.
* @type friend const Matrix
* @param m1 is the Matrix to compute the sub-Matrix from.
* @param r is the row to remove from m1.
* @param c is the column to remove from m1.
* @return a sub-Matrix of m1, reduced by one row and one column.
*/
friend const Matrix subMatrix(const Matrix& m1, int r,
int c);

/** Computes a sub-Matrix from a Matrix. Currently I’m only
* supporting square Matrices, so r2-r1 should equal c2-c1.
* @type friend const Matrix
* @param m1 is the Matrix to compute the sub-Matrix from.
* @param r1 is the start row.
* @param c1 is the start column.
* @param r2 is the end row.
* @param c2 is the end column.
* @return a sub-Matrix of m1, from r1 to r2, c1 to c2.
*/
friend const Matrix subMatrix(const Matrix& m1, int r1,
int c1, int r2, int c2);

/** Computes the determinant of a Matrix.
* @type friend const T
* @param m is the Matrix to find the determinant of.
* @return the determinant of m.
*/
friend const T det(const Matrix& m);

/** Computes the cofactor of a Matrix element.
* @type friend const T
* @param m is the Matrix to use in the computation.
* @param r is the row of the element.
* @param c is the column of the element.
* @return the cofactor of m[r][c].
*/
friend const T cofactor(const Matrix& m, int r, int c);

/** Computes the cofactor Matrix of a Matrix.
* @type friend const T
* @param m1 is the Matrix to use in the computation.
* @return the cofactor Matrix of m1.
*/
friend const Matrix cofactor_Matrix(const Matrix& m1);

/** Computes the inverse of a Matrix using its cofactors.
* @type friend const T
* @param m1 is the Matrix to use in the computation.
* @return the inverse Matrix of m1.
*/
friend const Matrix cofactor_inverse(const Matrix& m1);

/** Computes the inverse of a Matrix using the method in the MPI
* course notes, pp79-80. Currently, I’m only implementing
* square Matrices.
* @type friend const T
* @param m1 is the Matrix to use in the computation.
* @return the inverse Matrix of m1.
*/
friend const Matrix MPI_Notes_inverse(const Matrix& m1);

// Data in the Matrix. Ended up making this private so main()
// could access it.
vector< vector > m;
private:
// Precision for printing out fields of the Matrix.
static const int precision = PRECISION;
// Size of the Matrix.
int size;
};
/** Constant for setting the precision of printing out a Matrix.

A precision of n means n significant digits will be printed.

@type const int */
template
const int Matrix: recision;

/** Matrix constructor.
@param init is a decision parameter. Valid values are 1 for a
    random Matrix of floats, otherwise the identity. Pass one of 
   the const ints defined as FULL=1 or IDENTITY=0. The default
    is FULL.
@param s is the row and column size of the Matrix. The default
    is defined above and can be changed.
@return a Matrix object.
*/
template
Matrix::Matrix(const int init, const int s) : size(s)
{
m.resize(size);
m.clear();
switch (init) {
case IDENTITY:
for (int i(0); i<size; i++) {
m[i].resize(size);
m[i].clear();
for (int j(0); j<size; j++) {
m[i].push_back(static_cast(0));
}
m[i][i] = static_cast(1);
}
break;
case FULL:
srand(time(0));
using namespace globals;
for (int i(0); i<size; i++) {
m[i].resize(size);
m[i].clear();
for (int j(0); j<size; j++) {
m[i].push_back(static_cast(rand()) /
static_cast(rand()));
}
}
break;
}
}
/** Matrix destructor. */
template
Matrix::~Matrix() {}

/** Explicit Matrix destructor for Matrices of TYPEs. */
template<>
Matrix::~Matrix()
{
for (int i(0); i<size; i++) {
m[i].vector::~vector();
}
// m.vector< vector >::~vector();
}

/** Matrix copy constructor.

@param m2 is the Matrix to copy.

@return a Matrix object.
*/
template
Matrix::Matrix(const Matrix& m2) : size(m2.size)
{
m.resize(size);
m.clear();
T temp;
for (int i(0); i<size; i++) {
m[i].resize(size);
m[i].clear();
for (int j(0); j<size; j++) {
temp = m2.m[i][j];
m[i].push_back(temp);
}
}
}

/** Overloaded Matrix constructor.

@param a is an array of T’s.

@return a Matrix object.
/
template
Matrix::Matrix(const T a, const int& s) : size(s)
{
m.resize(size);
m.clear();
for (int i(0); i<size; i++) {
m[i].resize(size);
m[i].clear();
int row(i*size);
for (int j(0); j<size; j++) {
m[i][j] = a[row+j];
}
}
}

/** Overloaded assignment operator for the Matrix class.

@type Matrix&

@param m2 is the Matrix to copy.

@return a Matrix that is identical to m2.
*/
template
Matrix& Matrix: perator=(const Matrix& m2)
{
if (this == &m2) return *this;
size = m2.size;
for (int i(0); i<size; i++) {
for (int j(0); j<size; j++) {
m[i][j] = m2.m[i][j];
}
}
return *this;
}

/** Overloaded output operator for the Matrix class. Allows printing

of a Matrix object.

@type ostream&

@param out is the stream to write onto.

@param m is the Matrix to write.

@return an ostream reference.
*/
template
ostream& operator<<(ostream& out, const Matrix& m)
{
out << std::setprecision(m.precision);
for (int i(0); i<m.size; i++) {
out << “[”;
for (int j(0); j<m.size; j++) {
out << m.m[i][j] << " ";
}
out << "]
";
}
out << endl;
return out;
}

/** Overloaded output operator for a vector.

@type ostream&

@param out is the stream to write onto.

@param m is the Matrix to write.

@return an ostream reference.
*/
template
ostream& operator<<(ostream& out, const vector& v)
{
out << std::setprecision(PRECISION);
out << “[”;
for (int i(0); i<v.size(); i++) {
out << v[i]<< " ";
}
out << “]” << endl;
return out;
}

/** Swap row1 of the Matrix m into row2 and swap row2 into row1.

@type void

@param m is a pointer to a Matrix.

@param row1 is the index of the first row.

@param row2 is the index of the second row.

@return void
*/
template
void swap_rows(Matrix *m1, int row1, int row2)
{
const size_t n = m1->size;
T temp[n];
for (int col(0); col<n; col++) {
temp[col] = m1->m[row1][col];
m1->m[row1][col] = m1->m[row2][col];
m1->m[row2][col] = temp[col];
}
}

/** Overloaded addition operator for the Matrix class.

@type const Matrix

@param m1 is the lhs Matrix to add with.

@param m2 is the rhs Matrix to add with.

@return a Matrix which is the sum of m1 and m2.
*/
template
const Matrix operator+(const Matrix& m1, const Matrix& m2)
{
assert(m1.size == m2.size);
Matrix sum(IDENTITY, m1.size);
for (int i(0); i<m1.size; i++) {
for (int j(0); j<m1.size; j++) {
sum.m[i][j] = m1.m[i][j] + m2.m[i][j];
}
}
return sum;
}

/** Overloaded subtraction operator for the Matrix class.

@type const Matrix

@param m1 is the lhs Matrix to subtract from.

@param m2 is the rhs Matrix to subtract with.

@return a Matrix which equals m1 - m2.
*/
template
const Matrix operator-(const Matrix& m1, const Matrix& m2)
{
return m1+(-m2);
}

/** Overloaded multiplication operator for the Matrix class.

@type const Matrix

@param m1 is the lhs Matrix to multiply with.

@param m2 is the rhs Matrix to be multiplied against.

@return a Matrix which is the product of m1 and m2.
/
template
const Matrix operator(const Matrix& m1, const Matrix& m2)
{
assert(m1.size == m2.size);
Matrix result(IDENTITY, m1.size);
for (int i(0); i<m1.size; i++) {
for (int j(0); j<m1.size; j++) {
T sum = dot_product(row(i, m1), col(j, m2));
result.m[i][j] = dot_product(row(i, m1), col(j, m2));
}
}
return result;
}

/** Overloaded multiplication operator for the Matrix class.

@type const vector

@param v is the vector to multiply with.

@param m is the rhs Matrix to be multiplied against.

@return a vector which is the product of v and m.
/
template
const vector operator(const vector& v, const Matrix& m)
{
assert(v.size() >= m.size);
vector result;
result.reserve(m.size);
T sum;
for (int i(0); i<m.size; i++) {
sum = static_cast(0);
for (int j(0); j<m.size; j++) {
sum += v[j] * m.m[j][i];
}
result.push_back(sum);
}
return result;
}

/** Overloaded multiplication operator for the Matrix class.

@type const vector

@param m is the Matrix to multiply with.

@param v is the vector to be multiplied against.

@return a vector which is the product of m and v.
/
template
const vector operator(const Matrix& m, const vector& v)
{
assert(v.size() >= m.size);
vector result;
result.reserve(m.size);
T sum;
for (int i(0); i<m.size; i++) {
sum = static_cast(0);
for (int j(0); j<m.size; j++) {
sum += v[j] * m.m[i][j];
}
result.push_back(sum);
}
return result;
}

/** Overloaded scalar multiplication operator for the Matrix

class.

@type const Matrix

@param s is the scalar to multiply with.

@param m1 is the Matrix to be multiplied against.

@return a Matrix which is the product of s and m1.
/
template
const Matrix operator(const T& s, const Matrix& m1)
{
Matrix m2(m1);
for (int i(0); i<m1.size; i++) {
for (int j(0); j<m1.size; j++) {
m2.m[i][j] *= s;
}
}
return m2;
}

/** Overloaded scalar multiplication operator for the Matrix

class.

@type const Matrix

@param m is the Matrix to be multiplied against.

@param s is the scalar to multiply with.

@return a Matrix which is the product of s and m.
/
template
const Matrix operator(const Matrix& m, const T& s)
{
return s * m;
}

/** Overloaded scalar multiplication operator for the Matrix

class.

@type const Matrix

@param s is the scalar to multiply with.

@param m is the Matrix to be multiplied against.

@return a Matrix which is the product of s and m.
/
template
const Matrix operator(const int& s, const Matrix& m)
{
return static_cast(s)*m;
}

/** Overloaded scalar multiplication operator for the Matrix

class.

@type const Matrix

@param m is the Matrix to be multiplied against.

@param s is the scalar to multiply with.

@return a Matrix which is the product of s and m.
/
template
const Matrix operator(const Matrix& m, const int& s)
{
return static_cast(s)*m;
}

/** Overloaded unary negation operator for the Matrix class.

@type const Matrix

@param m is the Matrix to negate.

@return a Matrix which is the negative of m.
/
template
const Matrix operator-(const Matrix& m)
{
return -1m;
}

/** Overloaded scalar division operator for the Matrix

class.

@type const Matrix

@param m1 is the Matrix to be divided against.

@param s is the scalar to divide with.

@return a Matrix which is the elementwise division of m1 by s.
*/
template
const Matrix operator/(const Matrix& m1, const T& s)
{
assert(s); // s != 0
Matrix m2(m1);
for (int i(0); i<m1.size; i++) {
for (int j(0); j<m1.size; j++) {
m2.m[i][j] /= s;
}
}
return m2;
}

/** Initializes a vector that has the same number of elements as

max(rows, cols) of the Matrix m.

@type vector

@param m is a reference to the Matrix used to initialize the vector.

@return a vector with the same size as m
*/
template
const vector init_vector(const Matrix& m)
{
vector v;
v.reserve(m.size);
srand(time(0));
for (int n(0); n<m.size; n++) {
v.push_back(static_cast(rand()) / static_cast(rand()));
}
return v;
}

/** A slow sorting algorithm based on Knuth’s shell sort.

@type void

@param m is a pointer to the Matrix to sort rows on.

@param col is the column to sort on.

@return void
*/
template
void sort_rows(Matrix *m1, int col)
{
const size_t n = m1->size;
for (int gap(n/2); 0<gap; gap/=2) {
for (int i(gap); i<n; i++) {
for (int j(i-gap); 0<=j; j-=gap) {
if (m1->m[j+gap][col] < m1->m[j][col]) {
swap_rows(m1, j, j+gap);
}
}
}
}
}

/** Computes the inverse of a Matrix using the Gauss-Jordan elimination

method.

@type const Matrix

@param m1 is the Matrix to find the inverse of.

@return the inverse of m.
*/
template
const Matrix Gauss_Jordan_inverse(const Matrix& m1)
{
Matrix m2(m1);
Matrix inv(IDENTITY, m1.size);
const size_t n = m1.size;
// For each col in the Matrix.
for (int i(0); i<n; i++) {
int p=i;
for (; p<n; p++) {
// Search until a non-zero first element is found.
if (m2.m[p][i]) {
break;
}
}
// If there’s a column of all zeros.
if (p==n) {
cout << "No inverse solution exists. Exiting at line: " <<
LINE << endl;
exit(0);
// Swap a row with a non-zero first element.
} else if (p!=i) {
swap_rows(&m2, p, i);
// Mirror for the RHS Matrix.
swap_rows(&inv, p, i);
}
// Divide highest row by its first element.
T divisor = m2.m[i][i];
for (int j(0); j<n; j++) {
m2.m[i][j] /= divisor;
// Mirror for the RHS Matrix.
inv.m[i][j] /= divisor;
}
// Eliminate the rest of the rows.
for (int j(i+1); j<n; j++) {
T mult = m2.m[j][i] / m2.m[i][i];
for (int k(0); k<n; k++) {
m2.m[j][k] -= m2.m[i][k] * mult;
// Mirror for the RHS Matrix.
inv.m[j][k] -= inv.m[i][k] * mult;
}
}
}
// Check that all the diagonals are occupied.
for (int i(0); i<n; i++) {
if (!m2.m[i][i]) {
cout << "No inverse solution exists. Exiting at line: " <<
LINE << endl;
exit(0);
}
}
// Backwards substitution to get the inverse on the LHS.
// For each row.
for (int i(n-1); i>-1; i–) {
// For each subrow.
for (int j(i-1); j>-1; j–) {
// For each column.
for (int k(0); k<n; k++) {
inv.m[j][k] -= m2.m[j][i]*inv.m[i][k];
}
}
}
return inv;
}

/** Returns a row from the Matrix m, where the rows go from 0 to

m.size - 1.

@type const vector

@param r is the row to choose.

@param m is the Matrix to choose the row from.

@return a row vector from m.
*/
template
const vector row(int r, const Matrix& m)
{
vector result;
result.reserve(m.size);
for (int i(0); i<m.size; i++) {
result.push_back(m.m[r][i]);
}
return result;
}

/** Returns a col from the Matrix m, where the cols go from 0 to

m.size - 1.

@type const vector

@param c is the col to choose.

@param m is the Matrix to choose the col from.

@return a col vector from m.
*/
template
const vector col(int c, const Matrix& m)
{
vector result;
result.reserve(m.size);
for (int i(0); i<m.size; i++) {
result.push_back(m.m[i][c]);
}
return result;
}

/** Computes the dot product of two vectors.

@type vector

@param v1 is a vector.

@param v2 is a vector.

@return the dot product of v1 and v2.
*/
template
T dot_product(const vector& v1, const vector& v2)
{
assert(v1.size() == v2.size());
T result = static_cast(0);
for (int i(0); i<v1.size(); i++) {
result += v1[i]*v2[i];
}
return result;
}

/** Computes the transpose of a Matrix.

@type const Matrix

@param m is the Matrix to find the transpose of.

@return the transpose of m.
*/
template
const Matrix transpose(const Matrix& m)
{
Matrix transpose(m);
for (int i(0); i<m.size; i++) {
for (int j(0); j<m.size; j++) {
transpose.m[i][j] = m.m[j][i];
}
}
return transpose;
}

/** Computes the determinant of a Matrix. Nice, easy to understand

recursive function - doesn’t work though for Matrices around size

10x10 or larger due to memory required. This is the most expensive

recursive function I’ve ever written!

@type const T

@param m is the Matrix to find the determinant of.

@return the determinant of m.
*/
template
const T det(const Matrix& m)
{
Matrix m1(m);
Matrix m2(m);
T result(static_cast(0));
T temp(static_cast(0));
T node(static_cast(0));
T t2(static_cast(0));

if (m1.size == 2) {
return m1.m[0][0]*m1.m[1][1] - m1.m[0][1]*m1.m[1][0];
}
for (int i(0); i<m1.size; i++) {
node = m1.m[0][i] * pow(-1, i);
m2 = subMatrix(m1, 0, i);
if (m2.size==2) {
t2 = m2.m[0][0]*m2.m[1][1] - m2.m[0][1]*m2.m[1][0];
} else {
t2 = det(m2);
}
result += node * t2;
}
return result;
}

/** Computes a sub-Matrix from a Matrix.

@type const Matrix

@param m1 is the Matrix to compute the sub-Matrix from.

@param r is the row to remove from m1.

@param c is the column to remove from m1.

@return a sub-Matrix of m1, reduced by one row and one column.
*/
template
const Matrix subMatrix(const Matrix& m1, int r, int c)
{
assert((r<m1.size) && (c<m1.size));
Matrix m2(IDENTITY, m1.size-1);
for (int i1(0), i2(0); i1<m1.size; i1++, i2++) {
if (i1==r) if (++i1>=m1.size) continue;
for (int j1(0), j2(0); j1<m1.size; j1++, j2++) {
if (j1==c) if (++j1>=m1.size) continue;
m2.m[i2][j2] = m1.m[i1][j1];
}
}
return m2;
}

/** Computes a sub-Matrix from a Matrix. Currently I’m only supporting

square Matrices, so r2-r1 should equal c2-c1.

@type const Matrix

@param m1 is the Matrix to compute the sub-Matrix from.

@param r1 is the start row.

@param c1 is the start column.

@param r2 is the end row.

@param c2 is the end column.

@return a sub-Matrix of m1, from r1 to r2, c1 to c2.
*/
template
const Matrix subMatrix(const Matrix& m1, int r1, int c1,
int r2, int c2)
{
assert((r2-r1) == (c2-c1));
Matrix m2(IDENTITY, r2-r1+1);
for (int i1(0), i2(r1); i2<=r2; i1++, i2++) {
for (int j1(0), j2(c1); j2<=c2; j1++, j2++) {
m2.m[i1][j1] = m1.m[i2][j2];
}
}
return m2;
}

/** Computes the cofactor of a Matrix element.

@type const T

@param m is the Matrix to use in the computation.

@param r is the row of the element.

@param c is the column of the element.

@return the cofactor of m[r][c].
*/
template
const T cofactor(const Matrix& m, int r, int c)
{
return det(subMatrix(m, r, c));
}

/** Computes the cofactor Matrix of a Matrix.

@type const T

@param m1 is the Matrix to use in the computation.

@return the cofactor Matrix of m1.
*/
template
const Matrix cofactor_Matrix(const Matrix& m1)
{
Matrix m2(IDENTITY, m1.size);
if (m1.size==2) { // Handle 2x2 Matrices.
m2.m[0][0] = m1.m[1][1]; m2.m[0][1] = -m1.m[1][0];
m2.m[1][0] = -m1.m[0][1]; m2.m[1][1] = m1.m[0][0];
return m2;
}
for (int i(0); i<m1.size; i++) {
for (int j(0); j<m1.size; j++) {
m2.m[i][j] = pow(-1, i+j) * det(subMatrix(m1, i, j));
}
}
return m2;
}

/** Computes the inverse of a Matrix using its cofactors.

@type const T

@param m1 is the Matrix to use in the computation.

@return the inverse Matrix of m1.
*/
template
const Matrix cofactor_inverse(const Matrix& m1)
{
T t(det(m1));
Matrix m2(cofactor_Matrix(m1));
m2 = transpose(m2)/t;
return m2;
}

/** Computes the inverse of a Matrix using the method in the MPI course

notes, pp79-80. Currently, I’m only implementing square Matrices

and only evenly sized Matrices.

@type const T

@param m1 is the Matrix to use in the computation.

@return the inverse Matrix of m1.
/
template
const Matrix MPI_Notes_inverse(const Matrix& m1)
{
int size(m1.size);
Matrix inverse(IDENTITY, size);
int half_size(size>>1);
// Create Matrices that are used in the calculations.
Matrix a11(subMatrix(m1, 0, 0, half_size-1, half_size-1));
Matrix a12(subMatrix(m1, 0, half_size, half_size-1, m1.size-1));
Matrix a21(subMatrix(m1, half_size, 0, m1.size-1, half_size-1));
Matrix a22(subMatrix(m1, half_size, half_size,
m1.size-1, m1.size-1));
// Create some more Matrices using those already created.
Matrix a11_inverse(Gauss_Jordan_inverse(a11));
Matrix b1(a11_inversea12);
Matrix b2(a21a11_inverse);
Matrix h(a22-a21a11_inversea12);
Matrix h_inverse(Gauss_Jordan_inverse(h));
// Four parts of the inverse are the upper left, upper right,
// lower left and lower right.
Matrix ul(a11_inverse + b1h_inverseb2);
Matrix ur(-b1h_inverse);
Matrix ll(-h_inverse*b2);
// Matrix lr(h_inverse); // Redundant.
// Fill the result Matrix with its four subMatrices.
for (int i(0); i<ul.size; i++) {
for (int j(0); j<ul.size; j++) {
inverse.m[i][j] = ul.m[i][j];
}
}
for (int i(0); i<ur.size; i++) {
for (int j(0); j<ur.size; j++) {
inverse.m[i][j+half_size] = ur.m[i][j];
}
}
for (int i(0); i<ll.size; i++) {
for (int j(0); j<ll.size; j++) {
inverse.m[i+half_size][j] = ll.m[i][j];
}
}
for (int i(0); i<h_inverse.size; i++) {
for (int j(0); j<h_inverse.size; j++) {
inverse.m[i+half_size][j+half_size] = h_inverse.m[i][j];
}
}
return inverse;
}

//@}

#endif //_MATRIX_H