   Matrix inversion and precision problems
  Hi all,I''ve written me a nice little Matrix class (C++), whose elements are objects of my class ''Complex'' (consisting of x and y, both doubles) I need the matrix stuff for a college software project in small signal analysis. I''ve run into some trouble with my matrix inversion routine recently, because even doubles seem to be not quite big enough to guarantee the precision I need to get correct results. And in the course of many subsequent computations, what starts out as a tiny rounding error can turn into something not even remotely resembling the correct result. =( Is there anything I can do to help that at all? I mean PC''s are used for scientific computing, right, so there has to be way to work around the limited precision of floating point numbers. Anyone had a similar problem and found a solution?   
  You can use expression templates to help reduce temporaries and the amount of operations that need to take place.There are lots of articles abut them and MTL and boost.uBLAS use them for their matrices. ;  There's always the long double type, but that's not the best long-term solution.  After any multiplication or inversion, check through each element of the matrix and see if it's less than some small value called epsilon ( ε ), typically 1e-6 or smaller, depending on how precise you want to be.  If an element is less than ε, make it zero.  Also, check the determinant when performing an inversion (since you'll need to calculate it anyways) and see if it's less than ε.  If it is, then your inverted matrix is going to come out very funky with extremely large values.   [edited by - aprosenf on February 11, 2004 3:50:13 PM];  Thanks Aprosenf, I''m gonna try that. ;  are you doing this by brute force or are you using the usual techniques? ;  Note: when I said check if an entry is less than epsilon, I meant check if its absolute value is less than epsilon.  Otherwise, you''d be zeroing out all of the negative entries. ;  When working with real numbers you always want to have the largest possible magnitude in the position where your going to get a leading one.  This is because dividing by a large positive or negative number will give you better accuracy when working with numbers stored in the IEEE floating point format.  You can search the matrix for the largest real number and swap the colum and row to get the largest possible number in the position where you want the leading one.  This causes a problem because you will no longer be reducing down to the identity matrix if you swap columns but you will be reducing to some known permutation matrix.For complex numbers you have a problem with "largest" magnitude but you can use euclidian distance with (x,y) or some other metric space that would make sense and be faster.  This all increases the run time of the algorithm but will give better results.All assuming your using gauss jordan elimination =p ;  I think I was wrong when I said switching the columns would cause you to have to reduce to a permutation matrix.  I think you would still have to reduce to the identity matrix but you would have to swap rows in your resultant inverse matrix.  I'm not even sure anymore I'm obviously not that bright =p   [edited by - O_o on February 12, 2004 8:17:51 PM]