128 비트 정밀도의 Lapack 사용

Question

(SVD)의 128 비트 정밀도 계산에 Lapack을 사용하려고하는데이를 수행하기 위해 일부 검은 색 컴파일러 마술이 있음을 알게되었습니다. 인텔 포트란 컴파일러 (ifort)는 컴파일러에 DOUBLE PRECISION으로 선언 된 모든 변수를 128 비트 실수로 취하도록 지시하는 옵션 -r16을 지원합니다.

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

나는 128이다 데이터 유형 _Quad를 생성 옵션 -Qoption,cpp,--extended_float_type와 인텔 C++ 컴파일러 (ICC)를 사용할 수 있습니다 (C++ 인) 내 프로그램이 통합하기 : 그래서 LAPACK 및 BLAS를 사용하여 컴파일 비트 부동 소수점 변수. 내 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 

모든 컴파일 잘 여기까지 작동합니다. 하지만 출력은 다음과 같습니다 :

 
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 

나는 U와 VT는 정확하지만 특이 값은 분명하지 않습니다. 왜 이런 일이 일어나는가 또는 어떻게 우회 할 수 있는지에 대한 아이디어를 얻은 사람이 있습니까?
도움 주셔서 감사합니다.

Answer 1

외부 라이브러리를 확장 된 정밀도로 사용하는 경우 기계 형식 산술에 대한 정보를 얻으려면 구형 d1mach.f, r1mach.f, i1mach.f, i1mach.f을 사용하는지 확인하십시오. 여기 몇 가지 값을 조정할 수 있습니다.

Fortran 90 내장 함수를 사용하여 이러한 기계 상수를 얻는 dlamch.f (여기 http://www.netlib.org/lapack/util/dlamch.f)를 사용하는 Lapack에서는 문제가되지 않습니다.

그러나 BLAS 또는 SLATEC을 사용하는 경우 문제가 될 수 있습니다.