Using Lapack with 128 bit precision

Question

I am trying to use Lapack for a 128 bit precision calculation of a matrix singular value decomposition (SVD) and I found out that there is some black compiler magic to accomplish this. The Intel Fortran compiler (ifort) supports the option -r16 which instructs the compiler to take all variables declared as DOUBLE PRECISION to be 128 bit reals. So I compiled Lapack and BLAS using:

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

To incorporate this in my program (which is C++) I can use the Intel C++ compiler (icc) with the option -Qoption,cpp,--extended_float_type which creates a data type _Quad that is a 128 bit floating point variable. My SVD example looks like this:

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

compiled with

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

Everything works fine this far. But the output is:

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

I checked that U and VT are correct, but the singular values are obviously not. Has anyone got an idea why this happens or how one could circumvent it?
Thanks for your help.

Does this example work correctly with ordinary double precision arithmetics? — ev-br, Apr 25 '12 at 10:17
@Zhenya Yes it does. It calculates the correct singular values when computed with ordinary double precision. (4.23607, 0.236068) — Maxwell, Apr 25 '12 at 12:12
In that case, I'd check the `DBDSQR` routine: as far as I can see from the source of the reference implementation (http://www.netlib.org/lapack/double/dgesvd.f), it's computing the singular values given the `U` and `VT` matrices. — ev-br, Apr 25 '12 at 14:23
What kind of computation are you doing that requires 128bits float precision? — Adrien, Aug 06 '12 at 14:17
Actually I resolved the issue by now. It was, to my shame, due to an error I introduced to the make.inc of lapack, which lead to not all files being compiled with the -r16 option. Specifically the files which needed to be compiled without optimization. They have got an extra "compile-options" variable which I overlooked. Never the less, thanks for your help. — Maxwell, Aug 08 '12 at 09:25

score 2 · Answer 1 · answered Oct 03 '12 at 15:19

When using external libraries with extended precision, also check whether they use old-style d1mach.f, r1mach.f, i1mach.f to get information on machine arithmetic. There may be some values to tweak here.

It can't be a problem with Lapack, which uses dlamch.f (here http://www.netlib.org/lapack/util/dlamch.f), which uses Fortran 90 intrinsics to get these machine constants.

But it may become problematic for example with BLAS or SLATEC, if you use them.

Using Lapack with 128 bit precision

1 Answers1