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

your code is not working

Could you post your error message here?