I am now using some c++ code to compute matrix inverse and the matrix is stored in a two-dimensional array.
However, I met some possible memory issues that I had no idea how to fix it. Sometimes, I can get good return while sometimes error occurs. Really appreciate if you can give advice.
I got error message "malloc: *** error for object 0x163faaa38: incorrect checksum for freed object - object was probably modified after being freed. ".
Also, please see this picture and I use it to indicate which part of the code might go wrong.
Thank you very much.
long double Determinant(long double **a, int n)
{
int i, j, j1, j2 ; // general loop and matrix subscripts
long double det = 0 ; // init determinant
// pointer to pointers to implement 2d
// square array
long double **m;
if (n == 1) { // should not get here
det = a[0][0] ;
}
else if (n == 2) { // basic 2X2 sub-matrix determinate
// definition. When n==2, this ends the
det = a[0][0] * a[1][1] - a[1][0] * a[0][1] ;// the recursion series
}
// recursion continues, solve next sub-matrix
else { // solve the next minor by building a
// sub matrix
det = 0 ; // initialize determinant of sub-matrix
// for each column in sub-matrix
for (j1 = 0 ; j1 < n ; j1++) {
// get space for the pointer list
m = new long double*[n-1];
//m = (double **) malloc((n-1)* sizeof(double *)) ;
for (i = 0 ; i < n-1 ; i++)
m[i] = new long double[n-1];
// i[0][1][2][3] first malloc
// m -> + + + + space for 4 pointers
// | | | | j second malloc
// | | | +-> _ _ _ [0] pointers to
// | | +----> _ _ _ [1] and memory for
// | +-------> _ a _ [2] 4 doubles
// +----------> _ _ _ [3]
//
// a[1][2]
// build sub-matrix with minor elements excluded
for (i = 1 ; i < n ; i++) {
j2 = 0 ; // start at first sum-matrix column position
// loop to copy source matrix less one column
for (j = 0 ; j < n ; j++) {
if (j == j1) continue ; // don't copy the minor column element
m[i-1][j2] = a[i][j] ; // copy source element into new sub-matrix
// i-1 because new sub-matrix is one row
// (and column) smaller with excluded minors
j2++ ; // move to next sub-matrix column position
}
}
det += pow(-1.0, 1.0 + j1 + 1.0) * a[0][j1] * Determinant(m, n-1);
// sum x raised to y power
// recursively get determinant of next
// sub-matrix which is now one
// row & column smaller
for (i = 0 ; i < n-1 ; i++) delete m[i];// free the storage allocated to
// to this minor's set of pointers
// free the storage for the original
// pointer to pointer
delete [] m;
}
}
return(det);
}
// calculate the cofactor of element (row,col)
void GetMinor(long double **src, long double **dest, int row, int col, int order)
{
// indicate which col and row is being copied to dest
int colCount=0,rowCount=0;
for(int i = 0; i < order; i++ )
{
if( i != row )
{
colCount = 0;
for(int j = 0; j < order; j++ )
{
// when j is not the element
if( j != col )
{
dest[rowCount][colCount] = src[i][j];
colCount++;
}
}
rowCount++;
}
}
}
// Calculate the determinant recursively.
long double CalcDeterminant(long double **mat, int order)
{
// order must be >= 0
// stop the recursion when matrix is a single element
if( order == 1 )
return mat[0][0];
// the determinant value
long double det = 0;
// allocate the cofactor matrix
long double **minor;
minor = new long double*[order-1];
for(int i=0;i<order-1;i++)
minor[i] = new long double[order-1];
for(int i = 0; i < order; i++)
{
// get minor of element (0,i)
GetMinor(mat, minor, 0, i, order);
// the recusion is here!
det += (i%2==1?-1.0:1.0) * mat[0][i] * CalcDeterminant(minor,order-1);
//det += pow( -1.0, i ) * mat[0][i] * CalcDeterminant( minor,order-1 );
}
// release memory
for(int i=0;i<order-1;i++)
delete [] minor[i];
delete [] minor;
return det;
}
// matrix inversioon
void MatrixInversion(long double **InputM, int order, long double **InvM)
{
// get the determinant of a
long double det = 1.0/CalcDeterminant(InputM, order);
// memory allocation
long double *temp = new long double[(order-1)*(order-1)];
long double **minor = new long 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++)
{
// get the co-factor (matrix) of A(j,i)
GetMinor(InputM,minor,j,i,order);
InvM[i][j] = det*CalcDeterminant(minor,order-1);
if((i+j)%2 == 1)
InvM[i][j] = -InvM[i][j];
}
}
// release memory
//delete [] minor[0];
delete [] temp;
for(int i=0; i<order-1; i++)
delete [] minor[i];
delete [] minor;
}
