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]
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 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
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 )
    	return det;

    // 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;

    // 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 *