Hi, I am having a little problem with this matrix class, it just drives me crazy.
I use this matrix class in my natural cubic spline implementation, so I instantiate some instance in a function to unit test the cubic splines
Instantiation like this here works fine
cross::matrix ma,ma2;
cross::natural_cubic_spline ncs;
Now when I change ma to ma(5) to create a 5x5 matrix the program crashes somewhen, after the test function terminated.
instatiating the matrix ma only works fine either as does, ma2 only or ncs only.
But as soon as I instatiate 2 cross::matrix objects the program crashes somewhen, not during usage of the matrix classes.
This bug really drives me crazy since I have no clue whats wrong with the code underneath.
What makes the whole thing even stranger is that I get an access violation, but not inside this code, however instantiation of the matrix class causes the access violation in the first place.
Any help is appreciated, please have a look at the constructors and destructors primarily maybe I am doing something wrong there that I am not aware of that corrupts the stack or something.
thx in advance
/*
* You may use this code for nonprofit projects.
* You have to include this header in any copy of that code.
*
* Copyright @ Eduard Aumüller 2006
* basiror@gmx.de
*/
#ifndef CROSS_MATRIX
#define CROSS_MATRIX
#include <vector>
#include <algorithm>
#include <iostream>
using std::cout;
using std::endl;
namespace cross
{
class matrix
{
public:
matrix()
{
m_uiDimension = 0;
m_puiRowIndices = NULL,
m_puiColumnIndices = NULL;
m_pflElements = NULL;
m_pflInternalStorage = NULL;
};
matrix(const matrix& m)
{
if(m.dimension() == 0)
return;
init(m.dimension());
std::copy(&m.m_puiRowIndices[0], &m.m_puiRowIndices[m.dimension()], &m_puiRowIndices[0]);
std::copy(&m.m_puiColumnIndices[0], &m.m_puiColumnIndices[m.dimension()], &m_puiColumnIndices[0]);
std::copy(&m.m_pflElements[0], &m.m_pflElements[m.dimension()*m.dimension()], &m_pflElements[0]);
std::copy(&m.m_pflInternalStorage[0], &m.m_pflInternalStorage[m.dimension()], &m_pflInternalStorage[0]);
};
matrix(uint32 dim)
{
init(dim);
};
~matrix()
{
cout << "deconst"<<endl;
/* if(m_puiRowIndices != NULL)
delete [] m_puiRowIndices;
if(m_puiColumnIndices != NULL)
delete [] m_puiColumnIndices;
if(m_pflElements != NULL)
delete [] m_pflElements;
if(m_pflInternalStorage != NULL)
delete [] m_pflInternalStorage;*/
};
void init(uint32 dim)
{
uint32 i=0;
if(dim != m_uiDimension)
{
m_uiDimension = dim;
if(m_puiRowIndices != NULL)
delete [] m_puiRowIndices;
if(m_puiColumnIndices != NULL)
delete [] m_puiColumnIndices;
if(m_pflElements != NULL)
delete [] m_pflElements;
if(m_pflInternalStorage != NULL)
delete [] m_pflInternalStorage;
m_puiRowIndices = new uint32[dimension()];
m_puiColumnIndices = new uint32[dimension()];
m_pflElements = new float32[dimension()*dimension()];
m_pflInternalStorage = new float32[dimension()];
}
std::fill<float32*,float32>(&m_pflElements[0],&m_pflElements[dimension()*dimension()],0.0f);
std::fill<float32*,float32>(&m_pflInternalStorage[0],&m_pflInternalStorage[dimension()],0.0f);
//set row & column indices
for(i=0;i<dimension();i++)
{
m_puiRowIndices=i;
m_puiColumnIndices=i;
}
};
bool invert()
{
uint32 i=0;
uint32 j=0;
float32 flSign = 1.0f;
float32 flInvDet;
flInvDet = det();
if(fabs(flInvDet) < EPSILON)
return false;
flInvDet = 1.0f / flInvDet;
float32 *pflInv = new float32[dimension()*dimension()];
for(i=0;i<dimension();i++)
{
for(j=0;j<dimension();j++)
{
flSign = pow(-1.0f,static_cast<float32>((i+j)%2));
pflInv[j*dimension()+i] = minor(i,j)*flSign*flInvDet;
}
}
std::copy< float32*,float32* >(&pflInv[0],&pflInv[dimension()*dimension()],&m_pflElements[0]);
delete [] pflInv;
return true;
};
public:
float32 elem(uint32 i, uint32 j) const
{
return m_pflElements[row_index(i)*dimension()+col_index(j)];
};
void set_elem(uint32 i, uint32 j, float32 flValue)
{
m_pflElements[row_index(i)*dimension()+col_index(j)] = flValue;
};
void swap_rows(uint32 i1, uint32 i2)
{
std::swap(m_puiRowIndices[i1],m_puiRowIndices[i2]);
};
void swap_columns(uint32 j1, uint32 j2)
{
std::swap(m_puiColumnIndices[j1],m_puiColumnIndices[j2]);
};
uint32 dimension() const { return m_uiDimension;};
//determinant
float32 det()
{
uint32 i=0;
float32 flResult = 0.0f;
float32 flSign = 1.0f;
for(i=0;i<dimension();i++)
{
flResult += flSign*elem(i,0)*minor(i,0);
flSign *= -1.0f;
}
return flResult;
};
//returns the minor determinant of coordinate i,j
float32 minor(uint32 i, uint32 j)
{
float32 flResult = 0.0f;
float32 flSign = 1.0f;
uint32 k=0;
//map rows & columns to get a dimension()-1 sized matrix
for(k=i;k>0;k--)
swap_rows(k,k-1);
for(k=j;k>0;k--)
swap_columns(k,k-1);
flResult = flSign*minor_mapped(1,dimension()-1);
//unmap rows & columns to reconstruct the original matrix
for(k=0;k<i;k++)
swap_rows(k,k+1);
for(k=0;k<j;k++)
swap_columns(k,k+1);
return flResult;
};
//solves the system for the result input
bool solve_system(const float32 *pflResult, float32 *pflCoefficients)
{
float32 flDet = 0.0f;
uint32 j=0;
flDet = det();
if(fabs(flDet) < EPSILON)
return false;
flDet = 1.0f/flDet;
for(j=0;j<dimension();j++)
{
save_column(j,m_pflInternalStorage);
set_column(j,pflResult);
pflCoefficients[j] = det()*flDet;
set_column(j,m_pflInternalStorage);
}
return true;
};
private:
//stores column j in a buffer
void save_column(uint32 j,float32 *pflInternalStorage)
{
uint32 i=0;
for(i=0;i<dimension();i++)
{
pflInternalStorage = elem(i,j);
}
};
//replaces column j with new values
void set_column(uint32 j,const float32 *pflResult)
{
uint32 i=0;
for(i=0;i<dimension();i++)
{
set_elem(i,j,pflResult);
}
};
//expects a mapped squared matrix from [start][start] to [end][end]
//sub determinants are evaluated with a single swap
float32 minor_mapped(uint32 start, uint32 end)
{
uint32 i=0;
uint32 j=0;
float32 flSign = 1.0f;
float32 flResult = 0.0f;
float32 flElem = 0.0f;
if((end+1) - start == 3)
{
flResult = (
elem(start,start)*elem(start+1,start+1)*elem(start+2,start+2)
+elem(start,start+1)*elem(start+1,start+2)*elem(start+2,start)
+elem(start,start+2)*elem(start+1,start)*elem(start+2,start+1)
-elem(start+2,start)*elem(start+1,start+1)*elem(start,start+2)
-elem(start+2,start+1)*elem(start+1,start+2)*elem(start,start)
-elem(start+2,start+2)*elem(start+1,start)*elem(start,start+1)
);
}
else
{
//loop through the left column and alter the sign each iteration
for(i=start;i<=end;i++)
{
//get multiplication element before swapping
flElem = elem(i,start);
for(j=i;j>start;j--)
swap_rows(j,j-1);
flResult += flSign*flElem*minor_mapped(start+1,end);
flSign *= -1.0f;
for(j=start;j<i;j++)
swap_rows(j,j+1);
}
}
return flResult;
};
uint32 row_index(uint32 i) const { return m_puiRowIndices;};
uint32 col_index(uint32 j) const { return m_puiColumnIndices[j];};
private:
uint32 m_uiDimension;
uint32* m_puiRowIndices;
uint32* m_puiColumnIndices;
float32* m_pflElements;
float32* m_pflInternalStorage;
};
inline std::ostream& operator << ( std::ostream& strm, const matrix& m)
{
uint32 i=0;
uint32 j=0;
for(i=0;i<m.dimension();i++)
{
cout <<"| ";
for(j=0;j<m.dimension();j++)
{
cout<<std::fixed<<m.elem(i,j)<<" ";
}
cout <<"|"<<endl;
}
return strm;
};
};
#endif /*CROSS_MATRIX*/
