Matrix Inversion for complex number on C++

To implement Discrete Wavenumber Method (DWM) on C++, I selected boost library + FFTW. Because DWM needs support of complex number, matrix inversion and FFT, all of which can be found in the boost, FFTW & STL library. However I later realized boost ublas does not handle matrices with large condition number. For example, A=[6,6]
(((0,0.000109621),(0.0167211,-0.000109621),(0.250979,2.90726e-15),(-0.250979,-38.2832),(-0,-0),(0,0)),
((-9.33491e-05,0),(9.33491e-05,-0.0213859),(2.90726e-15,0.213724),(-48.9633,-0.213724),(0,0),(0,0)),
((0,-0),(0,0),(0,0),(0,0),(-1.26981e-18,0.000109621),(0.0167211,-0.000109621)),
((-0,-4.86612e+08),(5.4689e+10,4.86612e+08),(-1.11411e+12,0.00657373),(1.11411e+12,-1.25212e+14),(0,-0),(0,0)),
((3.96682e+08,0),(-3.96682e+08,-6.05079e+10),(0.0105204,-9.08209e+11),(-1.38534e+14,9.08209e+11),(0,-0),(0,0)),
((0,-0),(0,0),(0,-0),(0,0),(-5.73716e-07,-3.17712e+08),(1.59949e+10,3.17712e+08)))
will raise boost “internal logic error” which stems in this commend “BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, unit_lower> (m), e), cm1), internal_logic ());matrix_type cm2 (e);”. Google says it is due to the ill-condition of the matrix and the boost ublas authors wanted to ensure the matrix is invertible within the machine precision (epsilon > 1.11e-16 for double).  Indeed, adding ” #define BOOST_UBLAS_NDEBUG 1″ to the code to avoid checking produced a wrong inverse matrix. Interestingly, the determinant is not zero ( -2.9762e+003 +2.4557e+016i), and Matlab can invert the above matrix with a warning, which is actually desirable at this stage. Therefore using analytic solution, i.e., the way of directly calculating an inverse matrix may be a way out. I modified the C++ code from http://chi3x10.wordpress.com/2008/05/28/calculate-matrix-inversion-in-c/ and share it here with those who are in the middle of C++ coding using boost library. Caution: do make sure your inverse matrix multiple the original matrix is not distant from the identity matrix.


 /*
  * Analytical Matrix inversion routine
  */
template 
int GetMinor(matrix& src, matrix& dest, int row, int col)
{
    // indicate which col and row is being copied to dest
    int colCount=0,rowCount=0,order=src.size1();

    for(int i = 0; i
T CalcDeterminant( matrix& mat)
{
    // order must be >= 0
	// stop the recursion when matrix is a single element
	int order=mat.size1();
	T det;

    if( order == 1 )
    {
    	det=mat(0,0);
    	return det;
    }

	det=0.0;
    // the determinant value
   // T det = 0;

    // allocate the cofactor matrix
    matrix minor_mat (order-1,order-1);

    for(int i = 0; i
void MatrixInversion(matrix& A, matrix& Y)
{
	int order=A.size1();
    // get the determinant of a
    T detA,det;
    detA=CalcDeterminant(A);
	detA=(RealType)1.0/detA;

    // memory allocation
	 matrix minor_a (order-1, order-1);

    for(int j=0;j

It is good for small square matrices. NOT tested for large matrices. If your matrix is well scaled within the machine precision, use "fancier" inverse algorithm please.

2 thoughts on “Matrix Inversion for complex number on C++

Leave a Reply

Your email address will not be published. Required fields are marked *