Abstract: Paralution is an open source library for sparse iterative methods with special focus on multi-core and accelerator technology such as GPUs. It has a simple fortran interface and not well designed for multiple right-hand-sides. Here we defined a user friendly interface based on ISO_C_BINDING for Fortran FEA programmes to call Paralution solver.
keywords: Paralution, Fortran, sparse iterative method, iso_c_binding, FEA
Paralution系由瑞典科研工作者开发的基于多种稀疏矩阵迭代求解方法的C++并行计算程序包。它支持包括OpenMP、GPU/CUDA、OpenCL、OpenMP/MIC等多种后台并行库，也允许通过MPI进行多节点并行加速计算。它是多许可证的开源软件，方便易用，且不依赖特定的硬件平台和软件库，支持主流的Linux、Windows和MacOS操作平台。具体内容可以参考其官方主页：www.paralution.com
Paralution附带的算例中有Fortran调用的例子，只是对于有限元程序而言通常求解线性方程组AX=B的时候，右端项B不是向量形式而是矩阵形式。其次，在动力计算中总刚的稀疏结构通常而言并不发生改变，所以迭代求解器和预处理只需进行一次，这也是软件包附带算例没有考虑的。
在有限元分析中调用Paralution计算线性稀疏方程组AX=B的一般步骤是：
1、初始化paralution环境；
2、设置求解器并导入总刚；
3、分析求解X；
4、退出并清空paralution环境
具体的介绍详见代码和算例
CPP code:
1 #include <paralution.hpp> 2 3 namespace fortran_module{ 4 // Fortran interface via iso_c_binding 5 /* 6 Author: T.Q 7 Date: 2014.11 8 Version: 1.1 9 License: GPL v3 10 */ 11 extern "C" 12 { 13 // Initialize paralution 14 void paralution_f_init(); 15 // Print info of paralution 16 void paralution_f_info(); 17 // Build solver and preconditioner of paralution 18 void paralution_f_build(const int[], const int[], const double[], int ,int); 19 // Solve Ax=b 20 void paralution_f_solve(const double[], double[], int , int&, double&, int&); 21 // Stop paralution 22 void paralution_f_stop(); 23 // Subspace method for compute general eigenpairs 24 //void paralution_f_subspace(); 25 } 26 27 int *row = NULL;// Row index 28 int *col = NULL;// Column index 29 int *row_offset = NULL;// Row offset index for CSR 30 double *val = NULL;// Value of sparse matrix 31 32 double *in_x = NULL;// Internal x vector 33 double *in_rhs = NULL;// Internal rhs vector 34 35 bool var_clear = false;// Flag of variables clearance 36 bool ls_build = false;// Flag of solver exist 37 38 using namespace paralution; 39 40 41 LocalMatrix<double> *mat_A = NULL;// Matrix A i.e. Stiffness matrix in FEM 42 LocalMatrix<double> *mat_B = NULL;// Matrix B i.e. Mass matrix in FEM 43 44 // Iterative Linear Solver and Preconditioner 45 IterativeLinearSolver<LocalMatrix<double>, LocalVector<double>, double> *ls = NULL; 46 Preconditioner<LocalMatrix<double>, LocalVector<double>, double> *pr = NULL; 47 48 void paralution_f_clear() 49 { 50 // Clear variables 51 52 if(ls!=NULL){ ls->Clear(); delete ls;} 53 if(pr!=NULL){ pr->Clear(); delete pr;} 54 55 if(row!=NULL)free_host(&row); 56 if(col!=NULL)free_host(&col); 57 if(row_offset!=NULL)free_host(&row_offset); 58 if(val!=NULL)free_host(&val); 59 60 if(in_x!=NULL)free_host(&in_x); 61 if(in_rhs!=NULL)free_host(&in_rhs); 62 63 if(mat_A!=NULL){ mat_A->Clear(); delete mat_A;} 64 if(mat_B!=NULL){ mat_B->Clear(); delete mat_B;} 65 66 var_clear = true; 67 ls_build = false; 68 } 69 70 void paralution_f_init() 71 { 72 paralution_f_clear(); 73 init_paralution(); 74 } 75 76 void paralution_f_info() 77 { 78 info_paralution(); 79 } 80 81 void paralution_f_stop() 82 { 83 paralution_f_clear(); 84 stop_paralution(); 85 } 86 87 88 void paralution_f_build(const int for_row[], const int for_col[], 89 const double for_val[], int n, int nnz) 90 { 91 int i; 92 93 if(var_clear==false)paralution_f_clear(); 94 95 // Allocate arrays 96 allocate_host(nnz,&row); 97 allocate_host(nnz,&col); 98 allocate_host(nnz,&val); 99 100 // Copy sparse matrix data F2C 101 for(i=0;i<nnz;i++){ 102 row[i] = for_row[i] - 1;// Fortran starts with 1 default 103 col[i] = for_col[i] - 1; 104 val[i] = for_val[i];} 105 106 // Load data 107 mat_A = new LocalMatrix<double>; 108 mat_A->SetDataPtrCOO(&row,&col,&val,"Imported Fortran COO Matrix", nnz, n, n); 109 mat_A->ConvertToCSR();// Manual suggest CSR format 110 mat_A->info(); 111 112 // Creat Solver and Preconditioner 113 ls = new CG<LocalMatrix<double>, LocalVector<double>, double>; 114 pr = new Jacobi<LocalMatrix<double>, LocalVector<double>, double>; 115 116 ls->Clear(); 117 ls->SetOperator(*mat_A); 118 ls->SetPreconditioner(*pr); 119 ls->Build(); 120 121 ls_build = true; 122 123 } 124 125 void paralution_f_solve(const double for_rhs[], double for_x[], int n, 126 int &iter, double &resnorm, int &err) 127 { 128 int i; 129 LocalVector<double> x;// Solution vector 130 LocalVector<double> rhs;// Right-hand-side vector 131 132 assert(ls_build==true); 133 134 if(in_rhs!=NULL)free_host(&in_rhs); 135 if(in_x!=NULL)free_host(&in_x); 136 137 allocate_host(n,&in_rhs); 138 allocate_host(n,&in_x); 139 // Copy and set rhs vector 140 for(i=0;i<n;i++){ in_rhs[i] = for_rhs[i];} 141 rhs.SetDataPtr(&in_rhs,"vector",n); 142 // Set solution to zero 143 x.Allocate("vector",n); 144 x.Zeros(); 145 // Solve 146 ls->Solve(rhs,&x); 147 // Get solution 148 x.LeaveDataPtr(&in_x); 149 // Copy to array 150 for(i=0;i<n;i++){ for_x[i] = in_x[i];} 151 // Clear 152 x.Clear(); 153 rhs.Clear(); 154 // Get information 155 iter = ls->GetIterationCount();// iteration times 156 resnorm = ls->GetCurrentResidual();// residual 157 err = ls->GetSolverStatus();// error code 158 } 159 // TODO 160 // Subspace method for solve general eigenpairs for symmetric real positive 161 // defined matrix 162 // A*x=lambda*M*x 163 // A: matrix N*N i.e. stiffness matrix 164 // B: matrix N*N i.e. mass matrix 165 // 166 void paralution_f_subspace(const int for_row[], const int for_col[], 167 const int option[], const double for_stif[], const double for_mass[], 168 double eigval[], double eigvec[], int n, int nnz) 169 { 170 } 171 }