PDA

View Full Version : Inverse Matrix and Plane Equation



Zadkiel
06-22-2001, 02:03 PM
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 4*4 matrix. I do know that not all 4*4 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

john
06-22-2001, 05:40 PM
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 (p). 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
06-22-2001, 06:23 PM
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
06-23-2001, 12:38 AM
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 4*4 matrix. I only know 2*2.

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).]

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

ffish
06-23-2001, 02:50 AM
(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 http://www.opengl.org/discussion_boards/ubb/smile.gif). 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 <cassert>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <iostream>
#include <iomanip>
#include <stack>
#include <vector>
#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 T>
class Matrix;

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

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

template<class T>
const Matrix<T> operator-(const Matrix<T>&amp; m1, const Matrix<T>&amp; m2);

template<class T>
const Matrix<T> operator*(const Matrix<T>&amp; m1, const Matrix<T>&amp; m2);

template<class T>
const vector<T> operator*(const vector<T>&amp; v, const Matrix<T>&amp; m);

template<class T>
const vector<T> operator*(const Matrix<T>&amp; m, const vector<T>&amp; v);

template<class T>
const Matrix<T> operator*(const T&amp; s, const Matrix<T>&amp; m);

template<class T>
const Matrix<T> operator*(const Matrix<T>&amp; m, const T&amp; s);

template<class T>
const Matrix<T> operator*(const int&amp; s, const Matrix<T>&amp; m);

template<class T>
const Matrix<T> operator*(const Matrix<T>&amp; m, const int&amp; s);

template<class T>
ostream&amp; operator<<(ostream&amp; out, const Matrix<T>&amp; m);

template<class T>
const Matrix<T> operator-(const Matrix<T>&amp; m);

template<class T>
void swap_rows(Matrix<T> *m, int row1, int row2);

template<class T>
void sort_rows(Matrix<T> *m, int col);

template<class T>
const vector<T> init_vector(const Matrix<T>&amp; m);

template<class T>
const vector<T> row(int r, const Matrix<T>&amp; m);

template<class T>
const vector<T> col(int c, const Matrix<T>&amp; m);

template<class T>
T dot_product(const vector<T>&amp; v1, const vector<T>&amp; v2);

template<class T>
const Matrix<T> Gauss_Jordan_inverse(const Matrix<T>&amp; m);

template<class T>
const Matrix<T> transpose(const Matrix<T>&amp; m);

template<class T>
const Matrix<T> subMatrix(const Matrix<T>&amp; m1, int r, int c);

template<class T>
const Matrix<T> subMatrix(const Matrix<T>&amp; m1, int r1, int c1, int r2,
int c2);

template<class T>
const T det(const Matrix<T>&amp; m);

template<class T>
const T cofactor(const Matrix<T>&amp; m, int r, int c);

template<class T>
const Matrix<T> cofactor_Matrix(const Matrix<T>&amp; m1);

template<class T>
const Matrix<T> cofactor_inverse(const Matrix<T>&amp; m1);

template<class T>
const Matrix<T> operator/(const Matrix<T>&amp; m, const T&amp; s);

template<class T>
const Matrix<T> MPI_Notes_inverse(const Matrix<T>&amp; 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 T>
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<T>&amp; m2);

/** Overloaded Matrix constructor.
* @param a is an array of T's.
* @return a Matrix object.
*/
Matrix(const T* a, const int&amp; 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<T>&amp;
* @param m2 is the Matrix to copy.
* @return a Matrix that is identical to m2.
*/
Matrix<T>&amp; operator=(const Matrix<T>&amp; 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<T> operator+ <>(const Matrix<T>&amp; m1,
const Matrix<T>&amp; 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<T> operator- <>(const Matrix<T>&amp; m1,
const Matrix<T>&amp; 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<T> operator* <T>(const Matrix<T>&amp; m1,
const Matrix<T>&amp; m2);

/** Overloaded multiplication operator for the Matrix class.
* @type friend const vector<T>
* @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<T> operator* <T>(const vector<T>&amp; v,
const Matrix<T>&amp; m);

/** Overloaded multiplication operator for the Matrix class.
* @type friend const vector<T>
* @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<T> operator* <T>(const Matrix<T>&amp; m,
const vector<T>&amp; v);

/** Overloaded scalar multiplication operator for the Matrix
* class.
* @type friend const Matrix<T>
* @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<T> operator* <T>(const int&amp; s,
const Matrix<T>&amp; m);

/** Overloaded scalar multiplication operator for the Matrix
* class.
* @type friend const Matrix<T>
* @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<T> operator* <T>(const Matrix<T>&amp; m,
const int&amp; s);

/** Overloaded scalar multiplication operator for the Matrix
* class.
* @type friend const Matrix<T>
* @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<T> operator* <T>(const T&amp; s,
const Matrix<T>&amp; m);

/** Overloaded scalar multiplication operator for the Matrix
* class.
* @type friend const Matrix<T>
* @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<T> operator* <T>(const Matrix<T>&amp; m,
const T&amp; s);

/** Overloaded scalar division operator for the Matrix
* class.
* @type friend const Matrix<T>
* @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<T> operator/ <T>(const Matrix<T>&amp; m,
const T&amp; s);

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

/** Overloaded unary negation operator for the Matrix class.
* @type friend const Matrix<T>
* @param m1 is the Matrix to negate.
* @return a Matrix which is the negative of m1.
*/
// friend const Matrix<T> operator- <T>(const Matrix<T>&amp; 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<T>(Matrix<T> *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<T>(Matrix<T> *m, int col);

/** Initializes a vector that has the same number of rows, cols
* as the Matrix m.
* @type friend const vector<T>
* @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<T> init_vector<T>(const Matrix<T>&amp; m);

/** Returns a row from the Matrix m, where the rows go from 0 to
* m.size - 1.
* @type friend const vector<T>
* @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<T> row<T>(int r, const Matrix<T>&amp; m);

/** Returns a col from the Matrix m, where the cols go from 0 to
* m.size - 1.
* @type friend const vector<T>
* @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<T> col<T>(int c, const Matrix<T>&amp; m);

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

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

/** Computes a sub-Matrix from a Matrix.
* @type friend const Matrix<T>
* @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<T> subMatrix<T>(const Matrix<T>&amp; 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<T>
* @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<T> subMatrix<T>(const Matrix<T>&amp; 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<T>(const Matrix<T>&amp; 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<T>(const Matrix<T>&amp; 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<T> cofactor_Matrix<T>(const Matrix<T>&amp; 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<T> cofactor_inverse<T>(const Matrix<T>&amp; 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<T> MPI_Notes_inverse<T>(const Matrix<T>&amp; m1);

// Data in the Matrix. Ended up making this private so main()
// could access it.
vector< vector<T> > 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<class T>
const int Matrix<T>: http://www.opengl.org/discussion_boards/ubb/tongue.gifrecision;

/** 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<class T>
Matrix<T>::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<T>(0));
}
m[i][i] = static_cast<T>(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<T>(rand()) /
static_cast<T>(rand()));
}
}
break;
}
}

/** Matrix destructor. */
template<class T>
Matrix<T>::~Matrix() {}

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

/** Matrix copy constructor.
* @param m2 is the Matrix to copy.
* @return a Matrix object.
*/
template<class T>
Matrix<T>::Matrix(const Matrix<T>&amp; 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<class T>
Matrix<T>::Matrix(const T* a, const int&amp; 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<T>&amp;
* @param m2 is the Matrix to copy.
* @return a Matrix that is identical to m2.
*/
template<class T>
Matrix<T>&amp; Matrix<T>: http://www.opengl.org/discussion_boards/ubb/redface.gifperator=(const Matrix<T>&amp; m2)
{
if (this == &amp;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&amp;
* @param out is the stream to write onto.
* @param m is the Matrix to write.
* @return an ostream reference.
*/
template<class T>
ostream&amp; operator<<(ostream&amp; out, const Matrix<T>&amp; 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] << "\t";
}
out << "]\n";
}
out << endl;
return out;
}

/** Overloaded output operator for a vector.
* @type ostream&amp;
* @param out is the stream to write onto.
* @param m is the Matrix to write.
* @return an ostream reference.
*/
template<class T>
ostream&amp; operator<<(ostream&amp; out, const vector<T>&amp; v)
{
out << std::setprecision(PRECISION);
out << "[";
for (int i(0); i<v.size(); i++) {
out << v[i]<< "\t";
}
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<class T>
void swap_rows(Matrix<T> *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<class T>
const Matrix<T> operator+(const Matrix<T>&amp; m1, const Matrix<T>&amp; m2)
{
assert(m1.size == m2.size);
Matrix<T> 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<class T>
const Matrix<T> operator-(const Matrix<T>&amp; m1, const Matrix<T>&amp; 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<class T>
const Matrix<T> operator*(const Matrix<T>&amp; m1, const Matrix<T>&amp; m2)
{
assert(m1.size == m2.size);
Matrix<T> 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<T>
* @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<class T>
const vector<T> operator*(const vector<T>&amp; v, const Matrix<T>&amp; m)
{
assert(v.size() >= m.size);
vector<T> result;
result.reserve(m.size);
T sum;
for (int i(0); i<m.size; i++) {
sum = static_cast<T>(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<T>
* @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<class T>
const vector<T> operator*(const Matrix<T>&amp; m, const vector<T>&amp; v)
{
assert(v.size() >= m.size);
vector<T> result;
result.reserve(m.size);
T sum;
for (int i(0); i<m.size; i++) {
sum = static_cast<T>(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<T>
* @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<class T>
const Matrix<T> operator*(const T&amp; s, const Matrix<T>&amp; m1)
{
Matrix<T> 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<T>
* @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<class T>
const Matrix<T> operator*(const Matrix<T>&amp; m, const T&amp; s)
{
return s * m;
}

/** Overloaded scalar multiplication operator for the Matrix
* class.
* @type const Matrix<T>
* @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<class T>
const Matrix<T> operator*(const int&amp; s, const Matrix<T>&amp; m)
{
return static_cast<T>(s)*m;
}

/** Overloaded scalar multiplication operator for the Matrix
* class.
* @type const Matrix<T>
* @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<class T>
const Matrix<T> operator*(const Matrix<T>&amp; m, const int&amp; s)
{
return static_cast<T>(s)*m;
}

/** Overloaded unary negation operator for the Matrix class.
* @type const Matrix<T>
* @param m is the Matrix to negate.
* @return a Matrix which is the negative of m.
*/
template<class T>
const Matrix<T> operator-(const Matrix<T>&amp; m)
{
return -1*m;
}

/** Overloaded scalar division operator for the Matrix
* class.
* @type const Matrix<T>
* @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<class T>
const Matrix<T> operator/(const Matrix<T>&amp; m1, const T&amp; s)
{
assert(s); // s != 0
Matrix<T> 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<T>
* @param m is a reference to the Matrix used to initialize the vector.
* @return a vector with the same size as m
*/
template<class T>
const vector<T> init_vector(const Matrix<T>&amp; m)
{
vector<T> v;
v.reserve(m.size);
srand(time(0));
for (int n(0); n<m.size; n++) {
v.push_back(static_cast<T>(rand()) / static_cast<T>(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<class T>
void sort_rows(Matrix<T> *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<T>
* @param m1 is the Matrix to find the inverse of.
* @return the inverse of m.
*/
template<class T>
const Matrix<T> Gauss_Jordan_inverse(const Matrix<T>&amp; m1)
{
Matrix<T> m2(m1);
Matrix<T> 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(&amp;m2, p, i);
// Mirror for the RHS Matrix.
swap_rows(&amp;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<T>
* @param r is the row to choose.
* @param m is the Matrix to choose the row from.
* @return a row vector from m.
*/
template<class T>
const vector<T> row(int r, const Matrix<T>&amp; m)
{
vector<T> 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<T>
* @param c is the col to choose.
* @param m is the Matrix to choose the col from.
* @return a col vector from m.
*/
template<class T>
const vector<T> col(int c, const Matrix<T>&amp; m)
{
vector<T> 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<T>
* @param v1 is a vector<T>.
* @param v2 is a vector<T>.
* @return the dot product of v1 and v2.
*/
template<class T>
T dot_product(const vector<T>&amp; v1, const vector<T>&amp; v2)
{
assert(v1.size() == v2.size());
T result = static_cast<T>(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<T>
* @param m is the Matrix to find the transpose of.
* @return the transpose of m.
*/
template<class T>
const Matrix<T> transpose(const Matrix<T>&amp; m)
{
Matrix<T> 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<class T>
const T det(const Matrix<T>&amp; m)
{
Matrix<T> m1(m);
Matrix<T> m2(m);
T result(static_cast<T>(0));
T temp(static_cast<T>(0));
T node(static_cast<T>(0));
T t2(static_cast<T>(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<T>
* @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<class T>
const Matrix<T> subMatrix(const Matrix<T>&amp; m1, int r, int c)
{
assert((r<m1.size) &amp;&amp; (c<m1.size));
Matrix<T> 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<T>
* @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<class T>
const Matrix<T> subMatrix(const Matrix<T>&amp; m1, int r1, int c1,
int r2, int c2)
{
assert((r2-r1) == (c2-c1));
Matrix<T> 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<class T>
const T cofactor(const Matrix<T>&amp; 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<class T>
const Matrix<T> cofactor_Matrix(const Matrix<T>&amp; m1)
{
Matrix<T> 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<class T>
const Matrix<T> cofactor_inverse(const Matrix<T>&amp; m1)
{
T t(det(m1));
Matrix<T> 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<class T>
const Matrix<T> MPI_Notes_inverse(const Matrix<T>&amp; m1)
{
int size(m1.size);
Matrix<T> inverse(IDENTITY, size);
int half_size(size>>1);
// Create Matrices that are used in the calculations.
Matrix<T> a11(subMatrix(m1, 0, 0, half_size-1, half_size-1));
Matrix<T> a12(subMatrix(m1, 0, half_size, half_size-1, m1.size-1));
Matrix<T> a21(subMatrix(m1, half_size, 0, m1.size-1, half_size-1));
Matrix<T> a22(subMatrix(m1, half_size, half_size,
m1.size-1, m1.size-1));
// Create some more Matrices using those already created.
Matrix<T> a11_inverse(Gauss_Jordan_inverse(a11));
Matrix<T> b1(a11_inverse*a12);
Matrix<T> b2(a21*a11_inverse);
Matrix<T> h(a22-a21*a11_inverse*a12);
Matrix<T> h_inverse(Gauss_Jordan_inverse(h));
// Four parts of the inverse are the upper left, upper right,
// lower left and lower right.
Matrix<T> ul(a11_inverse + b1*h_inverse*b2);
Matrix<T> ur(-b1*h_inverse);
Matrix<T> ll(-h_inverse*b2);
// Matrix<T> 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