行列特異値分解（SVD）の128ビット精度の計算にLapackを使用しようとしていますが、これを実現するための黒いコンパイラーの魔法があることがわかりました。Intel Fortranコンパイラー（ifort）は、128ビット実数として宣言されたすべての変数を取得するようにコンパイラーに指示するオプションをサポートしています。そこで、以下を使用してLapackとBLASをコンパイルしました。-r16DOUBLE PRECISION
ifort -O3 -r16 -c isamax.f -o isamax.o
ifort -O3 -r16 -c sasum.f -o sasum.o
...
これを私のプログラム（C ++）に組み込むために、128ビット浮動小数点変数で-Qoption,cpp,--extended_float_typeあるデータ型を作成するオプションを指定してIntel C ++コンパイラー（icc）を使用できます。_Quad私の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 特異値： inf inf U： -0.525731 -0.850651 -0.850651 0.525731 VT： -0.525731 0.850651 -0.850651 -0.525731
UとVTが正しいことを確認しましたが、特異値は明らかに正しくありません。なぜこれが起こるのか、またはどのようにそれを回避できるのか、誰かが知っていますか？
ご協力いただきありがとうございます。