Использование Lapack с 128-битной точностью

Question

Я пытаюсь использовать Lapack для вычисления точности 128-битной матричной декомпозиции по сингулярным значениям (SVD), и я обнаружил, что для этого есть некоторая магия черного компилятора. Компилятор Intel Fortran (ifort) поддерживает опцию -r16 который инструктирует компилятор принимать все переменные, объявленные как DOUBLE PRECISION быть 128-битными реалами. Поэтому я скомпилировал Lapack и BLAS, используя:

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

Чтобы включить это в мою программу (которая является C++), я могу использовать компилятор Intel C++ (icc) с опцией -Qoption,cpp,--extended_float_type который создает тип данных _Quad это 128-битная переменная с плавающей точкой. Мой пример SVD выглядит так:

#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;
}

составлено с

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

Пока все отлично работает. Но вывод:

Размер скаляра: 16
Входная матрица: 
1       2
2       3

Необходимый размер рабочего пространства: 134

Единственные значения: 
инф
инф

U: 
-0,525731       -0,850651
-0,850651       0,525731
VT: 
-0,525731       0,850651
-0,850651 -0,525731

Я проверил, что U и VT правильны, но единственные значения, очевидно, нет. Кто-нибудь понял, почему это происходит или как можно обойти это?
Спасибо за вашу помощь.

Answer 1

При использовании внешних библиотек с повышенной точностью также проверьте, используют ли они старый стиль d1mach.f, r1mach.f, i1mach.f для получения информации об арифметике машины. Здесь могут быть некоторые значения для настройки.

Это не может быть проблемой с Lapack, который использует dlamch.f (здесь http://www.netlib.org/lapack/util/dlamch.f), который использует встроенные функции Fortran 90 для получения этих машинных констант.

Но это может стать проблематичным, например, с BLAS или SLATEC, если вы их используете.