128ビット精度のLapackを使用する

Question

（SVD）の128ビット精度の計算にLapackを使用しようとしていますが、これを達成するために黒いコンパイラの魔法があることがわかりました。インテルFortranコンパイラー（ifo​​rt）は、-r16オプションをサポートしています。このオプションは、DOUBLE PRECISIONと宣言されたすべての変数を128ビットの実数にするようコンパイラーに指示します。私は128であるデータ型_Quadを作成するオプション-Qoption,cpp,--extended_float_typeでインテルC++コンパイラー（ICC）を使用することができます

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

を（C++である）私のプログラムでこれを組み込むには：だから私は使用してLAPACKおよびBLASをコンパイルビット浮動小数点変数。

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

すべてはここまで正常に動作してコンパイル

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

：私のSVDの例は次のようになります。しかし出力は：

 
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を使用するかどうかを確認また、拡張精度で外部ライブラリを使用して。ここで調整する価値があるかもしれません。

Fortran 90組み込み関数を使用してこれらのマシン定数を取得するdlamch.f（ここではhttp://www.netlib.org/lapack/util/dlamch.f）を使用するLapackでは、問題はありません。

しかし、BLASやSLATECを使用すると問題が発生することがあります。