[关闭]
@BravoWA 2016-11-25T17:15:37.000000Z 字数 6921 阅读 2490

PCG_MATLAB vs PCG_Eigen_MEX

LinearAlgebraSolver PCG Eigen

Movtivation

When I running my FSI code, I use Matlab bulit-in Preconditioned Conjugate Gradient (pcg) funtion to solve my pressure poisson equation. However, I found this built-in Matlab pcg solver is not multi-threaded. Hence, I want to use c-mex to warp a external c/cpp sparse linear algebra equation solver to speedup calculation. Here, my choice is Eigen (its development branch).

comparison code in Matlab

Platform: Win7 64bit with Intel C/CPP compiler, Matlab 2014b.

  1. ...
  2.     %%
  3.     % +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  4.     %  solve
  5.     % +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  7.     fprintf('MSG: Test of Eigen_pcg begin!!\n');
  8.     f_nnode = size(f_node_coord,1);
  9.     load('rhs_foreigentest.mat');
  10.     b = rhs;
  11.     x0 = zeros(f_nnode,1);
  13.     [f_H,b] = sol_CBS_impose_pbc(f_H,b,f_pbc_dof,f_pbc_val);
  14.     tic
  15.     [x1] = pcg(f_H, b,1e-5,1000 ); 
  16.     toc
  18.     b = rhs;
  19.     tic
  20.     [x] = c_sol_pcg_eigen(f_H_ir, f_H_jc, f_H_val, b, x0, 1000, 1e-5,int32(f_pbc_dof), f_pbc_val);
  21.     toc
  23.     b = rhs;
  24.     tic
  25.     [x2] = c_sol_pcg_eigen_OMP...
  26.     (f_H_ir, f_H_jc, f_H_val, b, x0, 1000, 1e-5,int32(f_pbc_dof), f_pbc_val); 
  27.     toc
  28. end

Sparse matrix H is [29292x29292], nnz of H is 389062.
Serial CPP-mex code for pcg_Eigen(dev)

Matlab cpp-mex bulid command：
mex -I".\eigen" c_sol_pcg_eigen.cpp

  1. /*==========================================================
  3. This subroutine is written in C for invoking in Matlab 2009a.
  4. Before using it, compile this C code using proper C compiler.
  6. [x] = c_sol_pcg_eigen(A_ir, A_jc, A_val, b, x0, maxiter, tol, pbc_dof, pbc_val)
  7. Input:
  8. (1) : A_ir, row-id vector for CSR Sparse matrix format.
  9. (2) : A_jc, col-id vector for CSR Sparse matrix format.
  10. (3) : A_val, value vector for CSR Sparse matrix format.
  11. (4) : b, right-hands term.
  12. (5) : x0, initial guess of unknowns variables.
  13. (6) : maxiter, max number of iterations.
  14. (7) : tol, tolerance of residuals to check if converged.
  15. (8) : pbc_dof,
  16. (9) : pbc_val.
  17. Output:
  18. (1) : x
  21. Author: Chen Jiang  2016-04-20-15.42 jiangchen2007@gmail.coom
  22. *========================================================*/
  23. #include "mex.h"
  24. #include "math.h"
  25. #include "matrix.h"
  26. #include <Eigen/Sparse>
  27. #include <Eigen/Dense>
  28. #include <vector>
  29. #include <algorithm>
  30. #include <iostream>
  32. using namespace Eigen;
  33. using namespace std;
  35. typedef Eigen::Triplet<double> T;
  38. // the gateway function
  39. void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
  40. {
  41.     //input vars
  42.     int *A_ir;
  43.     int *A_jc;
  44.     double *A_val;
  45.     double *b;
  46.     double *x0;
  47.     int maxiter;
  48.     double tol;
  49.     int *pbc_dof;
  50.     double *pbc_val;
  52.     //output vars
  53.     double *x;
  56.     //temp vars
  57.     int nrows;
  58.     int nnz;
  59.     std::vector<T> tripletList;
  60.     int npbc_dof;
  62.     //-----------------
  63.     // GetData
  64.     //-----------------
  65.     if(nrhs != 9){
  66.         mexErrMsgIdAndTxt("MEX:c_sol_pcg_eigen:rhs",
  67.                 "This function takes too much input arguments.");
  68.     }
  70.     A_ir = (int*)mxGetData(prhs[0]);
  71.     A_jc = (int*)mxGetData(prhs[1]);
  72.     A_val = mxGetPr(prhs[2]);
  73.     b = mxGetPr(prhs[3]);
  74.     x0 = mxGetPr(prhs[4]);
  75.     maxiter = mxGetScalar(prhs[5]);
  76.     tol = mxGetScalar(prhs[6]);
  77.     pbc_dof = (int*) mxGetData(prhs[7]);
  78.     pbc_val = mxGetPr(prhs[8]);
  80.     nrows = mxGetM(prhs[3]);
  81.     npbc_dof = mxGetM(prhs[7]);
  82.     nnz = mxGetM(prhs[0]);
  84.     plhs[0] = mxCreateDoubleMatrix(nrows,1,mxREAL);
  85.     x = mxGetPr(plhs[0]);
  87.     //-----------------
  88.     //calculations
  89.     //-----------------
  91.     //covert A_ir, A_jc, A_val to Eigen Sparse matrix
  92.     tripletList.reserve(nnz);
  93.     for(int i=0; i<nnz; i++){
  94.         tripletList.push_back(T(A_ir[i]-1, A_jc[i]-1, A_val[i]));
  95.     }
  97.     SparseMatrix<double> A_eigen(nrows, nrows);
  98.     A_eigen.setFromTriplets(tripletList.begin(), tripletList.end());
  100.     Map<VectorXd> b_eigen(b, nrows);
  102.     //apply pbc
  103.     for(int idof = 0; idof < npbc_dof; idof++){
  104.         int dof = pbc_dof[idof]-1;
  105.         A_eigen.coeffRef(dof, dof) = 1.0e8;
  106.         b_eigen(dof) = 1.0e8*pbc_val[idof];
  107.     }
  108.     //CG solver of Eigen
  109.     VectorXd x_eigen;
  110.     Map<VectorXd> x0_eigen(x0, nrows);
  111.     ConjugateGradient<SparseMatrix<double>, Lower|Upper, IncompleteCholesky<double> > cg;
  112.     //ConjugateGradient<SparseMatrix<double>, Lower|Upper > cg;
  113.     //cg.preconditioner().setDroptol(0.001);
  114.     cg.setTolerance(tol);
  115.     cg.setMaxIterations(maxiter);
  116.     cg.compute(A_eigen);
  117.     x_eigen = cg.solveWithGuess(b_eigen, x0_eigen);
  118.     //x_eigen = cg.solve(b_eigen);
  120.     //
  121.     for(int i=0; i<nrows; i++){
  122.         x[i] = x_eigen(i);
  123.     }
  124. }

OpenMP CPP-MEX code for PCG_Eigen(dev)

Mex build command in Windows7:
mex -I".\eigen_dev" c_sol_pcg_eigen_OMP.cpp COMPFLAGS="$COMPFLAGS /openmp"
Mex build command in Linux (opensuse):
mex c_sol_pcg_eigen_OMP.cpp -I".\eigen_dev" CXXPLFAGS="\$CXXFLAGS -fopenmp" LDFLAGS="\$LDFLAGS -fopenmp"

  1. /*==========================================================
  3. This subroutine is written in C for invoking in Matlab 2009a.
  4. Before using it, compile this C code using proper C compiler.
  6. [x] = c_sol_pcg_eigen(A_ir, A_jc, A_val, b, x0, maxiter, tol, pbc_dof, pbc_val)
  7. Input:
  8. (1) : A_ir, row-id vector for CSR Sparse matrix format.
  9. (2) : A_jc, col-id vector for CSR Sparse matrix format.
  10. (3) : A_val, value vector for CSR Sparse matrix format.
  11. (4) : b, right-hands term.
  12. (5) : x0, initial guess of unknowns variables.
  13. (6) : maxiter, max number of iterations.
  14. (7) : tol, tolerance of residuals to check if converged.
  15. (8) : pbc_dof,
  16. (9) : pbc_val.
  17. Output:
  18. (1) : x
  21. Author: Chen Jiang 2016-04-24-10.59 jiangchen2007@gmail.com
  22. *========================================================*/
  23. #include "mex.h"
  24. #include "math.h"
  25. #include "matrix.h"
  26. #include <Eigen/Sparse>
  27. #include <Eigen/Dense>
  28. #include <vector>
  29. #include <algorithm>
  30. #include <iostream>
  31. #include <omp.h>
  33. using namespace Eigen;
  34. using namespace std;
  36. typedef Eigen::Triplet<double> T;
  39. // the gateway function
  40. void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
  41. {
  42.     //
  43.     Eigen::initParallel();
  44.     Eigen::setNbThreads(2);
  45.     //input vars
  46.     int *A_ir;
  47.     int *A_jc;
  48.     double *A_val;
  49.     double *b;
  50.     double *x0;
  51.     int maxiter;
  52.     double tol;
  53.     int *pbc_dof;
  54.     double *pbc_val;
  56.     //output vars
  57.     double *x;
  60.     //temp vars
  61.     int nrows;
  62.     int nnz;
  63.     std::vector<T> tripletList;
  64.     int npbc_dof;
  66.     //-----------------
  67.     // GetData
  68.     //-----------------
  69.     if(nrhs != 9){
  70.         mexErrMsgIdAndTxt("MEX:c_sol_pcg_eigen:rhs",
  71.                 "This function takes too much input arguments.");
  72.     }
  74.     A_ir = (int*)mxGetData(prhs[0]);
  75.     A_jc = (int*)mxGetData(prhs[1]);
  76.     A_val = mxGetPr(prhs[2]);
  77.     b = mxGetPr(prhs[3]);
  78.     x0 = mxGetPr(prhs[4]);
  79.     maxiter = mxGetScalar(prhs[5]);
  80.     tol = mxGetScalar(prhs[6]);
  81.     pbc_dof = (int*) mxGetData(prhs[7]);
  82.     pbc_val = mxGetPr(prhs[8]);
  84.     nrows = mxGetM(prhs[3]);
  85.     npbc_dof = mxGetM(prhs[7]);
  86.     nnz = mxGetM(prhs[0]);
  88.     plhs[0] = mxCreateDoubleMatrix(nrows,1,mxREAL);
  89.     x = mxGetPr(plhs[0]);
  91.     //-----------------
  92.     //calculations
  93.     //-----------------
  95.     //covert A_ir, A_jc, A_val to Eigen Sparse matrix
  96.     tripletList.reserve(nnz);
  97.     for(int i=0; i<nnz; i++){
  98.         tripletList.push_back(T(A_ir[i]-1, A_jc[i]-1, A_val[i]));
  99.     }
  101.     SparseMatrix<double> A_eigen(nrows, nrows);
  102.     A_eigen.setFromTriplets(tripletList.begin(), tripletList.end());
  104.     Map<VectorXd> b_eigen(b, nrows);
  106.     //apply pbc
  107.     for(int idof = 0; idof < npbc_dof; idof++){
  108.         int dof = pbc_dof[idof]-1;
  109.         A_eigen.coeffRef(dof, dof) = 1.0e8;
  110.         b_eigen(dof) = 1.0e8*pbc_val[idof];
  111.     }
  112.     //CG solver of Eigen
  113.     VectorXd x_eigen;
  114.     Map<VectorXd> x0_eigen(x0, nrows);
  115.     ConjugateGradient<SparseMatrix<double>, Lower|Upper, IncompleteCholesky<double> > cg;
  116.     //ConjugateGradient<SparseMatrix<double>, Lower|Upper > cg;
  117.     //cg.preconditioner().setDroptol(0.001);
  118.     cg.setTolerance(tol);
  119.     cg.setMaxIterations(maxiter);
  120.     cg.compute(A_eigen);
  121.     x_eigen = cg.solveWithGuess(b_eigen, x0_eigen);
  122.     //x_eigen = cg.solve(b_eigen);
  124.     //
  125.     for(int i=0; i<nrows; i++){
  126.         x[i] = x_eigen(i);
  127.     }
  128. }

Result

norm(x1-x)/norm(x1) = 2.4111e-07
Elapsed time:
Matlab built-in pcg function: 2.39255 sec
Serial Cpp-mex pcg of Eigen: 0.56422 sec
OpenMP Cpp-mex pcg of Eigen: 0.518612 sec (2 threads only)

Conclusion

Even CPP-mex PCG of Eigen contains not only the PCG calculation codes but also the data structure transfering codes, it is still faster than its rival.

添加新批注
在作者公开此批注前，只有你和作者可见。
回复批注