Rcpp 中的英特尔 MKL SVD:结果不一致
Intel MKL SVD in Rcpp: inconsistent results
我写了一个共享的 C++ 库,我现在正尝试使用 Rcpp 的 .Call
函数来做一个 R 扩展来调用这个库中的一个函数。该函数使用英特尔 MKL 函数 LAPACKE_dgesdd,它执行矩阵的 SVD。
当我有一个大小为 m=5000 和 n=8 的矩阵时,R 扩展和本机 C++ 版本的输出(U、S、V' 矩阵)完全相同(小数点后第 15 位)点),正如人们所期望的那样。然而,当 m=5000 和 n=12 时,R 版本给出不一致的输出,它与原生 C++ 版本略有不同。此外,每次我 运行 它时,R 版本都会给我一个稍微不同的输出(不像当 n=8 时,它是一致的)。
我真的不知道如何解释这种奇怪的行为。有人有想法吗?
代码如下:
function.cpp:
#include <iostream>
#include "mkl_lapacke.h"
#include "mkl.h"
#include <chrono>
int testSvd(){
int m=5000;
int l=12;
//Allocate matrices
double * s2 = (double *)mkl_malloc( l*sizeof( double ), 64 );
double * Rt2 = (double *)mkl_malloc( l*l*sizeof( double ), 64 );
double * U_l2 = (double *)mkl_malloc( m*l*sizeof( double ), 64 );
double * AQ2 = (double *)mkl_malloc( m*l*sizeof( double ), 64 );
populate_matrix_random (m,l,AQ2); //Note that this populates AQ2 with the same matrix every time it is run
//Perform the svd
LAPACKE_dgesdd( LAPACK_ROW_MAJOR, 'S', m, l, AQ2, l, s2, U_l2, l, Rt2, l );
for (int i=0; i<50; i++ ){
printf( " %6.12f", U_l2[i] );
}
printf("\n\s:");
for (int i=0; i<50; i++ ){
printf( " %6.12f", s2[i] );
}
printf("\n\Rt\n:");
for (int i=0; i<50; i++ ){
printf( " %6.12f", Rt2[i] );
}
}
mainR.cpp:
#include "function.cpp"
#include <Rcpp.h>
RcppExport SEXP testSvdR() {
testSvd();
}
然后是R代码:
library(Rcpp);
rm(list = ls())
dyn.load("/opt/intel/composer_xe_2015/mkl/lib/intel64/libmkl_rt.so", local=FALSE);
dyn.load("mainR.so", local=FALSE);
result = .Call('testSvdR');
我在编译 function.cpp 时通过静态 linking MKL 解决了这个问题(我之前是动态 linking)。
我发现这个 link 也很有帮助:https://software.intel.com/en-us/articles/intel-mkl-custom-static-linkage
我不清楚为什么我在做动态 linkage 时会得到如此奇怪的结果。我所能想到的是,也许 MKL 的某些依赖项是从 R 自动加载的 OpenBLAS 或 LAPACKE 库中提取的。我真的不知道。
我写了一个共享的 C++ 库,我现在正尝试使用 Rcpp 的 .Call
函数来做一个 R 扩展来调用这个库中的一个函数。该函数使用英特尔 MKL 函数 LAPACKE_dgesdd,它执行矩阵的 SVD。
当我有一个大小为 m=5000 和 n=8 的矩阵时,R 扩展和本机 C++ 版本的输出(U、S、V' 矩阵)完全相同(小数点后第 15 位)点),正如人们所期望的那样。然而,当 m=5000 和 n=12 时,R 版本给出不一致的输出,它与原生 C++ 版本略有不同。此外,每次我 运行 它时,R 版本都会给我一个稍微不同的输出(不像当 n=8 时,它是一致的)。
我真的不知道如何解释这种奇怪的行为。有人有想法吗?
代码如下:
function.cpp:
#include <iostream>
#include "mkl_lapacke.h"
#include "mkl.h"
#include <chrono>
int testSvd(){
int m=5000;
int l=12;
//Allocate matrices
double * s2 = (double *)mkl_malloc( l*sizeof( double ), 64 );
double * Rt2 = (double *)mkl_malloc( l*l*sizeof( double ), 64 );
double * U_l2 = (double *)mkl_malloc( m*l*sizeof( double ), 64 );
double * AQ2 = (double *)mkl_malloc( m*l*sizeof( double ), 64 );
populate_matrix_random (m,l,AQ2); //Note that this populates AQ2 with the same matrix every time it is run
//Perform the svd
LAPACKE_dgesdd( LAPACK_ROW_MAJOR, 'S', m, l, AQ2, l, s2, U_l2, l, Rt2, l );
for (int i=0; i<50; i++ ){
printf( " %6.12f", U_l2[i] );
}
printf("\n\s:");
for (int i=0; i<50; i++ ){
printf( " %6.12f", s2[i] );
}
printf("\n\Rt\n:");
for (int i=0; i<50; i++ ){
printf( " %6.12f", Rt2[i] );
}
}
mainR.cpp:
#include "function.cpp"
#include <Rcpp.h>
RcppExport SEXP testSvdR() {
testSvd();
}
然后是R代码:
library(Rcpp);
rm(list = ls())
dyn.load("/opt/intel/composer_xe_2015/mkl/lib/intel64/libmkl_rt.so", local=FALSE);
dyn.load("mainR.so", local=FALSE);
result = .Call('testSvdR');
我在编译 function.cpp 时通过静态 linking MKL 解决了这个问题(我之前是动态 linking)。
我发现这个 link 也很有帮助:https://software.intel.com/en-us/articles/intel-mkl-custom-static-linkage
我不清楚为什么我在做动态 linkage 时会得到如此奇怪的结果。我所能想到的是,也许 MKL 的某些依赖项是从 R 自动加载的 OpenBLAS 或 LAPACKE 库中提取的。我真的不知道。