Is there a Lapack implementation in extended precision?

1) There are some implementations out there. If you google you will find a few. LAPACK does not support extended precision.

2) For Cholesky you can almost ``do it yourself``. Cholesky is pretty easy. Sometimes compiler with a flag can change the precision at compile time.

You can have a look at DPOTF2+DPOTRS ( or even better DPOTRF2+DPOTRS ) and try to use a multi-precision package below it. Should work. This is not ``too`` hard.

3) You can look at: (

https://github.com/langou/latl )

It uses C++ templates. If you know C++, you can use multiprecision with with it. It is a while I did not work on this, but, at the time, it was working nicely.

4) One of the main problems is that you will not have an optimized BLAS library. So computation will take time. If the matrices are small, this is fine. If they are large, this will take time.

I know there is the XBLAS implementation of BLAS in extended precision?

XBLAS uses extra precision internally. The idea is to use quadruple precision inside the routine to do accurate computations, then returns the result as double. What I am saying is that using XBLAS below Cholesky will not help in your case. (I am not suggesting that this is what you are suggesting.) We would need to write a full on multi-precision LAPACK.

I must deal with a symmetric system of equations, A*x=B with a very badly conditioned matrix (cond ~ 1.0e21). Solving the problem with double precision routines (using GNU Octave, which calls Lapack) I obtain a solution, but as the forward error bound (if I remember correctly) is FE<=cond(A)*epsilon, and epsilon is 2.22e-16, the forward error is bounded by the number 2e5.

This is correct. If the ill-conditioning comes from ``scaling``, maybe you can do something about it. Using DPOSVX or DPOSVXX could help.

I would like to use the DPOTRS subroutine in extended precision in order to solve my problem.

Please note that this is not only the solve (TRS) that needs to be in extended precision, but you also want the factorization (TRF) in extended precision. Unless you want to use iterative refinement, which would be a good idea. So you factor in double, and you use iterative refinement in quadruple to improve the quality of the solution. (This is the idea behind DPOSVXX.)

I must deal with a symmetric system of equations, A*x=B with a very badly conditioned matrix (cond ~ 1.0e21).

Please note that what often happens, unfortunately, is that, if the condition number of A in double precision is computed in double precision as 1e16, then, unfortunately, the condition number of A in quadruple precision is computed in quadruple precision as 1e34. I.e. what is really happening is that the condition number of A is "infinite" and increasing the precision does not really help. The more precision, the more ``you see`` the ill-conditioning. The ill-conditioning is intrinsic to the problem. Let us hope that this is not the case for your problems.

Cheers,

Julien.