Uso de Lapack con precisión de 128 bits

Question

Estoy tratando de usar Lapack para un cálculo de precisión de 128 bits de una matriz singular value decomposition (SVD) y descubrí que hay algo de magia en el compilador negro para lograr esto. El compilador Intel Fortran (ifort) admite la opción -r16 que indica al compilador que tome todas las variables declaradas como DOUBLE PRECISION como reales de 128 bits. Así que he recopilado Lapack y BLAS usando:

ifort -O3 -r16 -c isamax.f -o isamax.o 
ifort -O3 -r16 -c sasum.f -o sasum.o 
... 

Incorporar esto en mi programa (que es C++) Puedo utilizar el compilador Intel C++ (ICC) con la opción -Qoption,cpp,--extended_float_type que crea un tipo de datos _Quad que es un 128 variable de coma flotante de bit Mi ejemplo SVD se ve así:

#include "stdio.h" 
#include "iostream" 
#include "vector" 

using namespace std; 
typedef _Quad scalar; 

//FORTRAN BINDING 
extern "C" void dgesvd_(char *JOBU, char *JOBVT, int *M, int *N, 
    scalar *A, int *LDA, 
    scalar *S, 
    scalar *U, int *LDU, 
    scalar *VT, int *LDVT, 
    scalar *WORK, int *LWORK, int *INFO); 

int main() { 
    cout << "Size of scalar: " << sizeof(scalar) << endl; 
    int N=2; 
    vector<scalar> A(N*N); 
    vector<scalar> S(N); 
    vector<scalar> U(N*N); 
    vector<scalar> VT(N*N); 

    // dummy input matrix 
    A[0] = 1.q; 
    A[1] = 2.q; 
    A[2] = 2.q; 
    A[3] = 3.q; 
    cout << "Input matrix: " << endl; 
    for(int i = 0; i < N; i++) { 
     for(int j = 0;j < N; j++) 
      cout << double(A[i*N+j]) << "\t"; 
     cout << endl; 
    } 
    cout << endl; 

    char JOBU='A'; 
    char JOBVT='A'; 
    int LWORK=-1; 
    scalar test; 
    int INFO; 

    // allocate memory 
    dgesvd_(&JOBU, &JOBVT, &N, &N, 
     &A[0], &N, 
     &S[0], 
     &U[0], &N, 
     &VT[0], &N, 
     &test, &LWORK, &INFO); 
    LWORK=test; 
    int size=int(test); 
    cout<<"Needed workspace size: "<<int(test)<<endl<<endl; 
    vector<scalar> WORK(size); 

    // run... 
    dgesvd_(&JOBU, &JOBVT, &N, &N, 
     &A[0], &N, 
     &S[0], 
     &U[0], &N, 
     &VT[0], &N, 
     &WORK[0], &LWORK, &INFO); 
    // output as doubles 
    cout << "Singular values: " << endl; 
    for(int i = 0;i < N; i++) 
     cout << double(S[i]) << endl; 
    cout << endl; 
    cout << "U: " << endl; 
    for(int i = 0;i < N; i++) { 
    for(int j = 0;j < N; j++) 
     cout << double(U[N*i+j]) << "\t"; 
    cout << endl; 
    } 
    cout << "VT: " << endl; 
    for(int i = 0;i < N; i++) { 
    for(int j = 0;j < N; j++) 
     cout << double(VT[N*i+j]) << "\t"; 
    cout << endl; 
    } 
    return 0; 
} 

compilado con

icc test.cpp -g -Qoption,cpp,--extended_float_type -lifcore ../lapack-3.4.0/liblapack.a ../BLAS/blas_LINUX.a 

todo funciona bien hasta aquí. Pero la salida es:

 
Size of scalar: 16 
Input matrix: 
1  2 
2  3 

Needed workspace size: 134 

Singular values: 
inf 
inf 

U: 
-0.525731  -0.850651 
-0.850651  0.525731 
VT: 
-0.525731  0.850651 
-0.850651  -0.525731 

Comprobé que U y VT son correctos, pero los valores singulares obviamente no lo son. ¿Alguien tiene una idea de por qué sucede esto o cómo se podría eludir?
Gracias por su ayuda.

Answer 1

Al utilizar bibliotecas externas con precisión extendida, también verifique si usan el estilo antiguo d1mach.f, r1mach.f, i1mach.f para obtener información sobre la aritmética de la máquina. Puede haber algunos valores para modificar aquí.

No puede ser un problema con Lapack, que usa dlamch.f (aquí http://www.netlib.org/lapack/util/dlamch.f), que usa Fortran 90 intrinsics para obtener estas constantes de máquina.

Pero puede ser problemático, por ejemplo, con BLAS o SLATEC, si los usa.