[petsc-users] Lapack with Quadruple Precision in PETSc and SLEPc

[Jose E. Roman]

Yes. To complement Barry’s answer:

The matrix exponential is a particular case, since it is not directly available in LAPACK. First of all, I would suggest to upgrade to slepc-3.8 that has a new implementation of Higham’s method (Padé up to order 13). This might be more accurate than the basic Padé of order 6, but still it relies on some parameters that assume double precision (theta_m in Higham’s paper SIMAX 2005). So in this case the problem is in the SLEPc implementation, and has nothing to do with LAPACK. I have to think what needs to be changed for quad precision. 

Jose


> El 13 nov 2017, a las 19:50, Smith, Barry F. <bsmith at mcs.anl.gov> escribió:
> 
> 
> 
>  Tobias,
> 
>   When you use PETSc in quad precision you need to ./configure with --download-f2cblaslapack this uses a version of BLAS/LAPACK obtained by running f2c on the reference version of BLAS/LAPACK (that is, fortran code from netlib) and then massages the source code for quad precision. 
> 
> The tool we use to do the conversion and massaging is $PETSC_DIR/bin/maint/toclapack.sh 
> 
> It is possible that we do not do all the conversions for everything in LAPACK properly for quad precision. What we have tested with PETSc all seems to handle the full precision but the truth is we do not necessary test all the parts of LAPACK that may depend on such constants. 
> 
> If you find things we missed we would eagerly accept your corrections.
> 
>  Barry
> 
>> On Nov 13, 2017, at 12:22 PM, Tobias Jawecki <tobias.jawecki at tuwien.ac.at> wrote:
>> 
>> Dear all,
>> 
>> I am interested in computations with higher precision. The application is mainly error analysis of high order Magnus integrators. In some cases the asymptotic behavior of the error can only be observed when the error is already on double precision and round-off errors of the subroutines can matter. Computation with higher precision then helps to compute a reference solution.
>> 
>> Using PETSc/SLEPc with the float128 setting works for me, but I do have some questions on how it works with Lapack/Blas.
>> Does PETSc use the standard Lapack package with double precision changed to quadruple precision?
>> I did not find a lot about using the standard Lapack with quadruple precision online, I just saw that all the routines in Lapack expect double precision input (according to netlib.org).
>> 
>> For example the Lapack routine dstevr for the eigenvalue computation of small problems which appears in the Lanczos code of SLEPc.
>> 
>> I assume some Lapack methods use parameters which are adjusted on double precision?
>> 
>> To give an example what kind of parameters would need to be considered:
>> Using the Matrix Function structure of SLEPc to compute the Matrix Exponential Function with Krylov methods seems to use a Padé approximation for the small problem. The used Padé approximation is of order 6 with scaling and squaring s.t. |H 2^(-s)|<0.5. This choice of parameter leads to an error smaller than 1e-16 (N. Higham book on Matrix Functions) for approximating the Matrix Exponential of the small problem. The solution is then correct only up to double precision even if higher precision is used for the computations.
>> 
>> My question is how PETSc and SLEPc handle the Lapack methods which are needed on quadruple precision. My main concern is that some Lapack methods use parameters which are chosen so that the results are correct up to double precision.
>> 
>> Thanks for your efforts and greetings,
>> Tobias Jawecki
>