MKL svd 버그

svd의 값이 MKLNUMTHREADS에 따라 춤을 춘다. (11.0 Update 3부터 이런 문제를 일으키는 것 같다.)
예를 들어 다음과 같이 armadillo library와 mkl을 사용하는 경우

/* -*- mode: c++; -*- */
#include <armadillo>
int main(int argc, char *argv[])
{
    using namespace arma;
    mat M = randu(957, 957);
    mat U, V;
    vec s;
    svd(U, s, V, M);
    mat newM = U * diagmat(s) * V.t();
    printf("norm: %f\n", arma::norm(M - newM,2));
    return 0;
}

이런 결과가 나오고

$ MKL_NUM_THREADS=1 ./main
norm: 0.000000
$ MKL_NUM_THREADS=2 ./main
norm: 0.000000
$ MKL_NUM_THREADS=3 ./main
norm: 0.000000
$ MKL_NUM_THREADS=4 ./main
norm: 0.000000
$ MKL_NUM_THREADS=5 ./main
norm: 371.303371
$ MKL_NUM_THREADS=6 ./main
norm: 138.622780
$ MKL_NUM_THREADS=7 ./main
norm: 138.622780

mkl만 사용하는 다음 프로그램을 짜서 돌려보면

/* -*- mode: c++; -*- */
#include <mkl.h>
#include <string.h>
#include <math.h>
#include <stdio.h>

int main(int argc, char *argv[])
{
    lapack_int M = 958;
    double *mat = new double[M * M];
    double *mat_orig = new double[M * M];
    double *s = new double[M];
    double *diags = new double[M*M];
    double *U = new double[M * M];
    double *Vt = new double[M * M];
    double *superb = new double[M-1];

    for (int i=0; i<M*M; ++i)
        mat[i] = static_cast<double>(rand())/RAND_MAX;

    memcpy(mat_orig, mat, sizeof(double) * M * M);
    LAPACKE_dgesvd(LAPACK_COL_MAJOR, 'a', 'a', M, M, mat, M, s, U, M, Vt, M, superb);

    memset(diags, 0, sizeof(double) * M * M);
    for (int i=0; i<M; ++i)
        diags[i * M + i] = s[i];

    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, M, M, 1.0, U, M, diags, M, 0.0, mat, M);
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, M, M, 1.0, mat, M, Vt, M, 0.0, U, M);

    double res = 0.0;
    for (int i=0; i<M*M; ++i) {
        double t = mat_orig[i] - U[i];
        res += (t*t);
    }
    printf("norm: %f\n", sqrt(res));

    delete[] superb;
    delete[] Vt;
    delete[] U;
    delete[] diags;
    delete[] s;
    delete[] mat_orig;
    delete[] mat;
    return 0;
}

결과가 매우 심각하다. ㅡㅡ

$ MKL_NUM_THREADS=1 ./main2
norm: 0.000000
$ MKL_NUM_THREADS=2 ./main2
norm: 0.000000
$ MKL_NUM_THREADS=3 ./main2
norm: 0.000000
$ MKL_NUM_THREADS=4 ./main2
norm: 0.000000
$ MKL_NUM_THREADS=5 ./main2
norm: 457.024091
$ MKL_NUM_THREADS=6 ./main2
Segmentation fault (core dumped)
$ MKL_NUM_THREADS=7 ./main2
Segmentation fault (core dumped)

See http://software.intel.com/en-us/forums/topic/390568