From: Ian Scott <ian.scott@st...>  20050112 09:16:26

Riccardo Lattanzi wrote: > Thanks for your answer, I've search the web a little and I came out with > an idea to solve the problem (see below). I'd really appreciate if you > could give me your opinion on that. > > >>I take it using an SVD based inverse with double precision is not >>sufficient? >> > > > Unfortunately it's not sufficient. I have to deal with 750 by 750 > matrices of complex elements, which are badly conditioned. > I've implemented using Mathematica with 50 digits precision, but I need > the performance of C++ for my purposes. > > >>......... >> >>vnl should then allow you to create a matrix with a complex version of >>this new type: >> >>vnl_matrix<vcl_complex<vnl_arbitrary_precision> > >> >>The biggest problem however is that most (and probably all) of the >>matrix decomposition and inversion classes in vnl_algo are actually just >>thin wrappers around the netlib library. This Fortran library can only >>cope with matrices of float, double, complex<float> and complex<double> >>types. Writing a basic LU decomposition based class would be possible. >> >>Ian. >> > > > I found a package that allows to create a type with arbitrary precision > ( ARPREC, http://crd.lbl.gov/~dhbailey/mpdist/ ). The library also > includes the redefinition of the basic arithmetic operators to smartly > handle the new type. > I could write my own "matrixinverter" code using this package, but I'm > not sure I would end up with something as efficient as the VNL matrix > inversion. > > Do you think it's feasible to create a new type with ARPREC, use it inside > VNL and edit the VNL source codes for matrix inversion and matrix > multiplication forcing the program to use the typespecific operators > instead of the standard +  * / ? > > Or do you suggest to write my own LU decomposition using the ARPREC > library? Having had a quick look at ARPREC, it appears to provide the standard operator+, operator*, etc for the main numeric classes mpreal and mpcomplex. So vnl would should be able to create matrices of mpcomplex fairly easily. You might need to add some simple definitions of vnl_numeric_traits<mpreal>, vnl_complex_traits<mpcomplex>, etc. The basic matrix operations  addition, multiplication  should work fine on vnl_matrix<mpcomplex>. You would still need to write a decomposition class in order to invert your matrix (as I mentioned above, vnl_algo's decomposition and inversion classes only work on a strictly limited set of types) but I guess Numerical Recipes would help there. If you have to use mpcomplex to do the numbers, and write your own decomposition, then it is hard to say if the advantage of using vnl's matrix handling, is worth your effort in setting up vnl. On the plus side, VXL is pretty easy to setup, especially if you switch off everything but the numerics in the build options. On the plus side for us, we would be interested in your experience of using ARPREC+vnl, and we would be happy to put your decomposition code in the library. If you do want to write a LU decomposition class for vnl_matrix<mpcomplex> I suggest you write it for vnl_matrix<vcl_complex<double> > first and then turn it into a templated class. Ian. 