I've made a fully functional Matrix class, to use with OpenGL... but I think you should 1st try to do something youself. Later I can give you some bits...

Here's the header

Code:

#define MATRIX_IDENTITY ((const float**)0x00000004)
#define MATRIX_EMPTY ((const float**)0x00000002)
class Matrix{
private:
float* matrix;
int h,w;//height, width
void buildBuf();
friend bool check1(const Matrix& m);//auxiliaries for invert
friend bool check2(Matrix& m);
Matrix* synched;
bool dosync;
public:
Matrix(int h=1,int w=1,const float** = MATRIX_EMPTY);
Matrix(int,int,const float*);
Matrix(const Matrix&);
~Matrix();
void sync(Matrix&);//all basic operations are made both on this matrix and the one synched
bool unsync();
Matrix* getSynched();
inline int getHeight() const;
inline int getWidth() const;
float* operator[](int) const;
bool invert();
void mkIdentity();
void mkEmpty();
//basic operations
void swapRow(int,int);
void swapCol(int,int);
void addMultRow(int,int,float=1.0);
void addMultCol(int,int,float=1.0);
void multRow(int,float);
void multCol(int,float);
//boolean operators
Matrix& operator=(const Matrix&);//doesn't change synched
bool operator==(const Matrix&) const;
bool operator!=(const Matrix&) const;
//arithmetic operators, aren't aplied on synched
Matrix& operator+=(const Matrix&);
Matrix operator+ (const Matrix&) const;
Matrix& operator-=(const Matrix&);
Matrix operator- (const Matrix&) const;
Matrix& operator*=(const Matrix&);
Matrix operator* (const Matrix&) const;
Matrix& operator*=(float);
Matrix operator* (float) const;
//debuging
void print() const;
};
Matrix operator*(float, const Matrix&);

And an example of it's usage with your problem

[Matrix M]= 1/(AD-BC) * [D -B -C A]

Code:

Matrix A(...),B(...),C(...),D(...);
Matrix temp = A*D - B*C;
if(!temp.invert()){
std::cerr<<"Matrix not invertible.\n"
//break; return; exit();... whatever
}
Matrix M = temp*(D-B-C-A);