HomerReid / scuff-em

A comprehensive and full-featured computational physics suite for boundary-element analysis of electromagnetic scattering, fluctuation-induced phenomena (Casimir forces and radiative heat transfer), nanophotonics, RF device engineering, electrostatics, and more. Includes a core library with C++ and python APIs as well as many command-line applications.
http://www.homerreid.com/scuff-em
GNU General Public License v2.0
126 stars 50 forks source link

LUSolve() for an already factorized Hmatrix #47

Open eunnieverse opened 9 years ago

eunnieverse commented 9 years ago

I am importing a Hmatrix from a HDF5 file that has already been LUFactorized. LUSolve() currently requires that LUFactorize() be run again, which is double work. It would help if this step can be skipped.

HomerReid commented 9 years ago

The difficulty is that LUFactorize() creates some internal state that gets stored inside the HMatrix class but is not written to disk by the Export...() routines. In principle this could be rectified by modifying the Export() routines to write the extra state, although this would clutter the .hdf5 files.

Here is a somewhat tortured workaround.

Export:

M = (an HMatrix created somehow)
int N=M->NR;
M->LUFactorize();
ipivVector=new HVector(N);
for(int n=0; n<N; n++) 
 ipivVector->SetEntry(n, (double) (M->ipiv[n])) ;
M->ExportToHDF5(Context,"M");
ipivVector->ExportToHDF5(Context,"ipiv");
delete ipivVector;

Import:

M = new HMatrix("MyFile.hdf5","M");
int N=M->NR;
ipivVector=new HVector("MyFile.hdf5","ipiv");
M->ipiv=(int *)mallocEC(N*sizeof(int));
for(int n=0; n<N; n++) 
 M->ipiv[n]=(int)(ipivVector->GetEntryD(n));
delete ipivVector;