Back in 1975, I wrote a paper which was never published on backwards heat flow problems. Ill-posed by nature, one tries to exploit additional information to create well-posed extensions. The problem ultimately inverted with linear equation solvers is well-posed - but certainly still numerically challenging. When solved originally, I usedthe SLEroutine from what was then the "Michigan Terminal System" library. Results obtained were consistent with my theoretical limits.
It is ironic that now - when it has become desirable to exhume that old work - I am having trouble duplicating what I got back then. My current suspicion is embodied in the following piece from the HELP function for the GESV routine in MKL.
The dsgesv and zcgesv are mixed precision iterative refinement subroutines for exploiting fast single precision hardware. They first attempt to factorize the matrix in single precision (dsgesv) or single complex precision (zcgesv) and use this factorization within an iterative refinement procedure to produce a solution with double precision (dsgesv) / double complex precision (zcgesv) normwise backward error quality (see below). If the approach fails, the method switches to a double precision or double complex precision factorization respectively and computes the solution.
The iterative refinement is not going to be a winning strategy if the ratio single precision performance over double precision performance is too small. A reasonable strategy should take the number of right-hand sides and the size of the matrix into account. This might be done with a call to ilaenv in the future. At present, iterative refinement is implemented.
My question is whether I can force the double precision solution. My fear is that the system will deem the solution via "mixed precision" with "iterative refinement" to have succeeded when, in fact, in this challenging problem, it is bringing in non-trivial round-off. The "saving" in this instance, from using any "fast single precision hardware" may be counter-productive. Any ideas?
If you have previously called MKL to do the requisite L-U decomposition, and you are interested only in obtaining a solution, the call to make is to the generic routine gesv:
call gesv( a, b [,ipiv] [,info] )
If the required arguments are double precision, the call will be directed to the routine dgesv.
If the exhumed code contains specific assumptions regarding the floating point word size, precision, etc., those need to be changed from values applicable to the fossil system to values appropriate to the present platform.