12

我正在尝试使用Lapack对矩阵奇异值分解(SVD) 进行 128 位精度计算，我发现有一些黑色的编译器魔法可以实现这一点。英特尔 Fortran 编译器 (ifort) 支持-r16指示编译器将所有声明为DOUBLE PRECISION128 位实数的变量的选项。所以我编译了 Lapack 和 BLAS 使用：

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

要将其合并到我的程序（即 C++）中，我可以使用带有选项的英特尔 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

奇异值：
信息
信息

你：
-0.525731 -0.850651
-0.850651 0.525731
VT：
-0.525731 0.850651
-0.850651 -0.525731

我检查了 U 和 VT 是否正确，但奇异值显然不是。有没有人知道为什么会发生这种情况或如何规避它？
谢谢你的帮助。

4

1 回答 1

2

使用扩展精度的外部库时，还要检查它们是否使用旧式d1mach.f、r1mach.f、i1mach.f来获取有关机器算术的信息。这里可能需要调整一些值。

Lapack 不会有问题，它使用 dlamch.f（此处为http://www.netlib.org/lapack/util/dlamch.f），它使用 Fortran 90 内在函数来获取这些机器常量。

但是，如果您使用 BLAS 或 SLATEC，它可能会出现问题。

于 2012-10-03T15:19:37.013 回答