Matrix inverse

[sfotovati]

Friends,

I am a beginner programmer. I am trying to write a matrix inverse code. The code is as following:

#include "stdafx.h"
#include<iostream>
#include<fstream>
#include <math.h>
#include <string>
using namespace std;

int GetMinor(double **src, double **dest, int row, int col, int order)
{
int colCount=0,rowCount=0;
for(int i = 0; i < order; i++ )
{
if( i != row )
{
colCount = 0;
for(int j = 0; j < order; j++ )
{
if( j != col )
{
dest[rowCount][colCount] = src[i][j];
colCount++;
}
}
rowCount++;
}
}
return 1;
}

double Determinant( double **mat, int order)
{
if( order == 1 )
return mat[0][0];
double det = 0;
double **minor;
minor = new double*[order-1];
for(int i=0;i<order-1;i++)
minor[i] = new double[order-1];

for(int i = 0; i < order; i++ )
{
GetMinor( mat, minor, 0, i , order);
det += (i%2==1?-1.0:1.0) * mat[0][i] * Determinant(minor,order-1);
}

for(int i=0;i<order-1;i++)
delete [] minor[i];
delete [] minor;
return det;
}

void MatrixInv(double **A, int order, double **Y)
{
double det = 1.0/Determinant(A,order);
double *temp = new double[(order-1)*(order-1)];
double **minor = new double*[order-1];
for(int i=0;i<order-1;i++)
minor[i] = temp+(i*(order-1));
for(int j=0;j<order;j++)
{
for(int i=0;i<order;i++)
{
GetMinor(A,minor,j,i,order);
Y[i][j] = det*Determinant(minor,order-1);
if( (i+j)%2 == 1)
Y[i][j] = -Y[i][j];
}
}
delete [] temp;
delete [] minor;
}

int main()
{
double S[3][3]={0.};
double C[3][3]={0.};

S[0][0]=1;
S[0][1]=1;
S[0][2]=2;

S[1][0]=S[0][1];
S[1][1]=3;
S[1][2]=1;

S[2][0]=S[0][2];
S[2][1]=S[1][2];
S[2][2]=2;

MatrixInv(S,3,C);

cout << "determinant of S = " << Determinant(S,3) << endl;
cin.get();

return 0;
}

The User defined functions are what I have found from some forums. But it does not work.

I get two errorw :
- error C2664: 'Determinant' : cannot convert parameter 1 from 'double [3][3]' to 'double **'
- error C2664: 'MatrixInv' : cannot convert parameter 1 from 'double [3][3]' to 'double **'

I have no idea what goes wrong. I appreciate if you could help me with it.

Regards,

Shawn

[vlad from moscow]

You defined arrays as

double S[3][3]={0.};
double C[3][3]={0.};


When they are passed by value to a function they are implicitly converted to pointers to their first elements. The element of a two dimensional array is one-dimensional array. So when you pass your arrays to a function they are converted to

double ( * )[3]

It is not the same as double **

[sfotovati]

Thanks for the reply,

But I still don't get it. How should I declare Matrices S, and C in the main program? Note that as they are passed to functions, their size will be different at each time. So I know they must be in a dynamic form. But I have no idea how to specify them at the main.

[vlad from moscow]

For example you could declare function MatrixInv as

void MatrixInv(double ( *A )[3], int order, double ( * Y )[3] );

[sfotovati]

I did what you proposed. But i got the same error. the error is related to the time that the program wants to pass matrix to find its minor (adjugate) matrices. So, whenever I have "GetMinor(A,minor,j,i,order)" the error shows up.

I think I should not use matrix declaration in that way. I intend to calculate matrix inverse. Therefore, for a [3x3] matrix, it must also consider [2x2] matrices to bcalculate "adjugate" matrix for the initial one. Is there any other method I can specify the matrices S and C?

[keskiverto]

Many.

You actually create and destroy one here:

double *temp = new double[(order-1)*(order-1)];
double **minor = new double*[order-1];
for(int i=0;i<order-1;i++)
  minor[i] = temp+(i*(order-1));
..
delete [] temp;
delete [] minor;

And here, with a twist:

double **minor;
minor = new double*[order-1];
for(int i=0;i<order-1;i++)
  minor[i] = new double[order-1];
..
for(int i=0;i<order-1;i++)
  delete [] minor[i];
delete [] minor;

Or, you coud define a Matrix-class that offers neat access methods.

Or, you coud use std::vector< std::vector<double> >

[vlad from moscow]

You should use the declaration I showed everywhere where you are passing the original arrays to functions.

The other way is to use casting. For example

void MatrixInv(double **A, int order, double **Y] );
//...
MatrixInv( ( double **)S, 3, (double **)C);

or

MatrixInv( reinterpret_cast<double **>( S ), 3, reinterpret_cast<double **>( C ));




.

[sfotovati]

I got it keskiverto,

This is how I had to define them:

double **S;
double **C;

S=new double *[3];
C=new double *[3];
for(int i=0;i<3;i++)
{
S[i] = new double[3];
C[i] = new double[3];
}

Thanks everyone for helping beginners :)

[keskiverto]

The order is known at compile-time, and small. Therefore, templates could work too:

template <unsigned int N>
double Determinant( double mat[N][N] )
{
  if( N == 1 )
    return mat[0][0];
  double det = 0;
  double minor[N-1][N-1];

  for(int i = 0; i < N; ++i )
    {
      GetMinor( mat, minor, 0, i, N );
      det += (i%2==1?-1.0:1.0) * mat[0][i] * Determinant( minor );
    }
  return det;
}