Sử dụng Lapack với độ chính xác 128 bit

Question

Tôi đang cố gắng sử dụng Lapack để tính toán chính xác 128 bit của ma trận singular value decomposition (SVD) và tôi phát hiện ra rằng có một số phép thuật trình biên dịch màu đen để thực hiện việc này. Trình biên dịch Intel Fortran (ifort) hỗ trợ tùy chọn -r16 hướng dẫn trình biên dịch nhận tất cả các biến được khai báo là DOUBLE PRECISION là các bit thực 128 bit. Vì vậy, tôi biên soạn LAPACK và BLAS sử dụng:

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

Để kết hợp này trong chương trình của tôi (đó là C++) Tôi có thể sử dụng Intel C++ (ICC) với các tùy chọn -Qoption,cpp,--extended_float_type mà tạo ra một kiểu dữ liệu _Quad đó là một 128 biến số dấu chấm động bit. Ví dụ SVD My trông như thế này:

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

biên soạn với

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

Tất cả mọi thứ hoạt động tốt này đến nay. Nhưng đầu ra là:

 
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 

Tôi đã kiểm tra U và VT là chính xác, nhưng giá trị số ít rõ ràng là không. Có ai có một ý tưởng tại sao điều này xảy ra hoặc làm thế nào người ta có thể phá vỡ nó?
Cảm ơn sự giúp đỡ của bạn.

Answer 1

Khi sử dụng thư viện bên ngoài với độ chính xác mở rộng, cũng kiểm tra xem họ sử dụng kiểu cũ d1mach.f, r1mach.f, i1mach.f để có được thông tin trên máy tính số học. Có thể có một số giá trị để tinh chỉnh ở đây.

Nó không thể là một vấn đề với Lapack, sử dụng dlamch.f (ở đây http://www.netlib.org/lapack/util/dlamch.f), sử dụng Fortran 90 nội tại để có được các hằng số máy này.

Nhưng nó có thể trở thành vấn đề ví dụ với BLAS hoặc SLATEC, nếu bạn sử dụng chúng.