为什么mkl_sparse_d_ev稀疏特征求解器这么慢？

Question

我试图在intel-mkl中使用mkl_sparse_d_ev以获取最小的特征向量。但是，我发现mkl似乎非常慢。我的test_FtF.txt文件使用matlab eigs的时间只有0.02s，但是mkl需要5s。有人知道吗？测试文件已上传到https://github.com/qiqiSipangyiyou/mkl_test/tree/master。这是我的代码：

#include<iostream>
#include<fstream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include<vector>
#include "mkl.h"
#include "mkl_solvers_ee.h"
#include<chrono>

#define max(a, b) (a) < (b) ? (b): (a)

using namespace std;
int main()
{
    /* read Matrix A from file in COO format */
    //fstream FtF_file("/home/qicai/ap_build/mex_cmake/build-src-unknown-Default/mex/mexAdd/test.txt",std::ios::in);
    fstream FtF_file("/home/qicai/ap_build/mex_cmake/test_FtF.txt",std::ios::in);

    MKL_INT num_pts;//matrix size
    FtF_file>>num_pts;

    MKL_INT M = num_pts;               /* number of columns in matrix A */
    MKL_INT N = num_pts;               /* number of columns in matrix A */
    MKL_INT nnz = 0;             /* number of non-zeros in matrix */


    double tmp;
    cout<<num_pts<<endl;
    for(int i=0;i<num_pts;i++){
        for(int j=0;j<num_pts;j++){
            FtF_file>>tmp;
            if(j<i)continue;
        //    cout<<"data:"<<tmp<<endl;
            if(abs(tmp)>1e-10){
                nnz++;
            }
        }
    }
    cout<<"nonzeros:"<<nnz<<endl;
    FtF_file.clear();
    FtF_file.seekg(0);

    MKL_INT *ia=new MKL_INT[nnz];
    MKL_INT *ja=new MKL_INT[nnz];
    double *a=new double[nnz];
    MKL_INT id=0;
    FtF_file>>tmp;
    for(int i=0;i<num_pts;i++){
        for(int j=0;j<num_pts;++j){
            FtF_file>>tmp;
            if(j<i)continue;
            if(abs(tmp)>1e-10){
      //          cout<<"data:"<<tmp<<endl;

                ia[id]=i;
                ja[id]=j;
                a[id]=tmp;

                id++;
            }
        }
    }

    cout<<endl<<"ia:";
    for(int i=0;i<nnz;i++){
        cout<<ia[i]<<',';
    }
    cout<<endl<<"ja:";
    for(int i=0;i<nnz;i++){
        cout<<ja[i]<<',';
    }
    cout<<endl<<"a:";
    for(int i=0;i<nnz;i++){
        cout<<a[i]<<',';
    }

    /* mkl_sparse_d_ev input parameters */
    char         which = 'S'; /* Which eigenvalues to calculate. ('L' - largest (algebraic) eigenvalues, 'S' - smallest (algebraic) eigenvalues) */
    MKL_INT      pm[128];     /* This array is used to pass various parameters to Extended Eigensolver Extensions routines. */
    MKL_INT      k0  = 1;     /* Desired number of max/min eigenvalues */

    /* mkl_sparse_d_ev output parameters */
    MKL_INT      k;           /* Number of eigenvalues found (might be less than k0). */
    double *E=new double[k0];
    double *X=new double[num_pts];
    double *res=new double[num_pts]; /* Residual */
    /* Local variables */
    MKL_INT      info;               /* Errors */
    MKL_INT      compute_vectors = 1;/* Flag to compute eigenvecors */
    MKL_INT      tol = 7;            /* Tolerance */

    /* Sparse BLAS IE variables */
    sparse_status_t status;
    sparse_matrix_t A = NULL; /* Handle containing sparse matrix in internal data structure */
    sparse_matrix_t B = NULL; /* Handle containing sparse matrix in internal data structure */

    struct matrix_descr descr; /* Structure specifying sparse matrix properties */

    /* Create handle for matrix A stored in CSR format */
    descr.type =SPARSE_MATRIX_TYPE_SYMMETRIC;
    descr.mode =SPARSE_FILL_MODE_UPPER;
    descr.diag =SPARSE_DIAG_NON_UNIT;

    status = mkl_sparse_d_create_coo ( &A, SPARSE_INDEX_BASE_ZERO,N, N, nnz, ia, ja, a );
    mkl_sparse_convert_csr(A,SPARSE_OPERATION_NON_TRANSPOSE,&B);
    /* Step 2. Call mkl_sparse_ee_init to define default input values */
    mkl_sparse_ee_init(pm);

    pm[1] = tol; /* Set tolerance */
    pm[6] = compute_vectors;
    pm[2]=1;
    pm[8]=1;

    /* Step 3. Solve the standard Ax = ex eigenvalue problem. */
    auto start_time=chrono::steady_clock::now();
    info = mkl_sparse_d_ev(&which, pm, B, descr, k0, &k, E, X, res);
    auto end_time=chrono::steady_clock::now();
    auto dr_ms=std::chrono::duration<double>(end_time-start_time).count();
    cout<<"time cost for mkl_sparse_d_ev:"<<endl;
    cout << dr_ms <<endl;

    printf("mkl_sparse_d_ev output info %d \n",info);
    if ( info != 0 )
    {
        printf("Routine mkl_sparse_d_ev returns code of ERROR: %i", (int)info);
        return 1;
    }

    printf("*************************************************\n");
    printf("************** REPORT ***************************\n");
    printf("*************************************************\n");
    printf("#mode found/subspace %d %d \n", k, k0);
    printf("Index/Exact Eigenvalues/Estimated Eigenvalues/Residuals\n");

    mkl_sparse_destroy(A);
    return 0;
}

Answer 1

您运行的mkl的问题大小和版本是什么？是32位还是64位代码？