稀疏矩阵高性能求解器PARDISO的使用

2023年8月3日 浏览：1784 评论：4 收藏：2

PARDISO是一个商用的高速的稀疏矩阵线性方程组直接法求解器，其可以用于大型稀疏对称或非对称线性方程组的求解，商用版支持的并行方式包括共享内存，分布式内存和NVIDA的GPU并行。多次实践已经表明，PARDISO求解器具有优秀的性能，在各种稀疏矩阵直接求解器的比较中通常处于第一梯队，其他常见的直接法求解器包括SuiteSparse，MUMPS，SPOOLES等，商用有限元软件COMSOL的官方文档甚至直接指出在COMSOL中PARDISO求解速度快于MUMPS。。

对于线性方程组求解来说，一般常见的求解方法可分为直接法和迭代法。直接法常见的方法是高斯消去，LU分解法等，在求解过程中不需要迭代，但是在求解中可能需要较大内存，商业软件中经常使用的一种直接法叫做多波前法。迭代法又可分为两个大类，分裂迭代和Krylov迭代。前者通过将系数矩阵进行拆分形成迭代列式，常见的是雅克比迭代，高斯赛德尔迭代等，分裂迭代这一类型的迭代法在商业软件中很少直接使用，但是经常作为Krylov迭代的预处理基础；后者以Krylov子空间为基础上，常见的包括共轭梯度法，广义最小残差法等，Krylov迭代法在商业软件中较为常见。

在稀疏矩阵直接法求解中，一个常见的技术是矩阵重排，通过矩阵重排，可以使得矩阵的带宽减小（即非0元素位置相对更趋近于对角线），从而使得直接求解效率更高，一种常见的矩阵重排方法是Cuthill_McKee，常见的可以实现矩阵重排的库是METIS。以下是采用Cuthill_McKee重排前和重排后的非0元素分布。

稀疏矩阵高性能求解器PARDISO的使用的图1

重排前非0元素分布

稀疏矩阵高性能求解器PARDISO的使用的图2

重排后非0元素分布

商用版本的PARDISO个人一般较少直接使用（主要是要买）。实际中较为广泛使用的其中一个版本是配置在MKL库的PARDISO版本。PARDISO求解需要将矩阵采用CSR格式存储，并且在求解中需要正确的设置较多的参数才能正确运行。

以下提供MKL_PARDISO的一个官方例子：

Fortran版本：

INCLUDE 'mkl_PARDISO.f90'
PROGRAM PARDISO_sym_f90
USE mkl_PARDISO
IMPLICIT NONE
INTEGER, PARAMETER :: dp = KIND(1.0D0)
!.. Internal solver memory pointer 
TYPE(MKL_PARDISO_HANDLE), ALLOCATABLE  :: pt(:)
!.. All other variables
INTEGER maxfct, mnum, mtype, phase, n, nrhs, error, msglvl, nnz
INTEGER error1
INTEGER, ALLOCATABLE :: iparm( : )
INTEGER, ALLOCATABLE :: ia( : )
INTEGER, ALLOCATABLE :: ja( : )
REAL(KIND=DP), ALLOCATABLE :: a( : )
REAL(KIND=DP), ALLOCATABLE :: b( : )
REAL(KIND=DP), ALLOCATABLE :: x( : )
INTEGER i, idum(1)
REAL(KIND=DP) ddum(1)
!.. Fill all arrays containing matrix data.
n = 8 
nnz = 18
nrhs = 1 
maxfct = 1 
mnum = 1
ALLOCATE(ia(n + 1))
ia = (/ 1, 5, 8, 10, 12, 15, 17, 18, 19 /)
ALLOCATE(ja(nnz))
ja = (/ 1,    3,       6, 7,    &
           2, 3,    5,          &
              3,             8, &
                 4,       7,    &
                    5, 6, 7,    &
                       6,    8, &
                          7,    &
                             8 /)
ALLOCATE(a(nnz))
a = (/ 7.d0,        1.d0,             2.d0, 7.d0,        &
             -4.d0, 8.d0,       2.d0,                    &
                    1.d0,                         5.d0,  &
                          7.d0,             9.d0,        &
                                5.d0, 1.d0, 5.d0,        &
                                     -1.d0,       5.d0,  &
                                           11.d0,        &
                                                  5.d0 /)
ALLOCATE(b(n))
ALLOCATE(x(n))
!..
!.. Set up PARDISO control parameter
!..
ALLOCATE(iparm(64))
DO i = 1, 64
   iparm(i) = 0
END DO
iparm(1) = 1 ! no solver default
iparm(2) = 2 ! fill-in reordering from METIS
iparm(4) = 0 ! no iterative-direct algorithm
iparm(5) = 0 ! no user fill-in reducing permutation
iparm(6) = 0 ! =0 solution on the first n components of x
iparm(8) = 2 ! numbers of iterative refinement steps
iparm(10) = 13 ! perturb the pivot elements with 1E-13
iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
iparm(13) = 0 ! maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm(13) = 1 in case of inappropriate accuracy
iparm(14) = 0 ! Output: number of perturbed pivots
iparm(18) = -1 ! Output: number of nonzeros in the factor LU
iparm(19) = -1 ! Output: Mflops for LU factorization
iparm(20) = 0 ! Output: Numbers of CG Iterations
error  = 0 ! initialize error flag
msglvl = 1 ! print statistical information
mtype  = -2 ! symmetric, indefinite
!.. Initialize the internal solver memory pointer. This is only
! necessary for the FIRST call of the PARDISO solver.
ALLOCATE (pt(64))
DO i = 1, 64
   pt(i)%DUMMY =  0 
END DO
!.. Reordering and Symbolic Factorization, This step also allocates
! all memory that is necessary for the factorization
phase = 11 ! only reordering and symbolic factorization
CALL PARDISO (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
              idum, nrhs, iparm, msglvl, ddum, ddum, error)
    
WRITE(*,*) 'Reordering completed ... '
IF (error /= 0) THEN
   WRITE(*,*) 'The following ERROR was detected: ', error
   GOTO 1000
END IF
WRITE(*,*) 'Number of nonzeros in factors = ',iparm(18)
WRITE(*,*) 'Number of factorization MFLOPS = ',iparm(19)
!.. Factorization.
phase = 22 ! only factorization
CALL PARDISO (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
              idum, nrhs, iparm, msglvl, ddum, ddum, error)
WRITE(*,*) 'Factorization completed ... '
IF (error /= 0) THEN
   WRITE(*,*) 'The following ERROR was detected: ', error
   GOTO 1000
ENDIF
!.. Back substitution and iterative refinement
iparm(8) = 2 ! max numbers of iterative refinement steps
phase = 33 ! only solving
DO i = 1, n
   b(i) = 1.d0
END DO
CALL PARDISO (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
              idum, nrhs, iparm, msglvl, b, x, error)
WRITE(*,*) 'Solve completed ... '
IF (error /= 0) THEN
   WRITE(*,*) 'The following ERROR was detected: ', error
   GOTO 1000
ENDIF
WRITE(*,*) 'The solution of the system is '
DO i = 1, n
   WRITE(*,*) ' x(',i,') = ', x(i)
END DO
      
1000 CONTINUE
!.. Termination and release of memory
phase = -1 ! release internal memory
CALL PARDISO (pt, maxfct, mnum, mtype, phase, n, ddum, idum, idum, &
              idum, nrhs, iparm, msglvl, ddum, ddum, error1)
IF (ALLOCATED(ia))      DEALLOCATE(ia)
IF (ALLOCATED(ja))      DEALLOCATE(ja)
IF (ALLOCATED(a))       DEALLOCATE(a)
IF (ALLOCATED(b))       DEALLOCATE(b)
IF (ALLOCATED(x))       DEALLOCATE(x)
IF (ALLOCATED(iparm))   DEALLOCATE(iparm)
IF (error1 /= 0) THEN
   WRITE(*,*) 'The following ERROR on release stage was detected: ', error1
   STOP 1
ENDIF
IF (error /= 0) STOP 1
END PROGRAM PARDISO_sym_f90


具体使用：安装VS和INTEL oneAPI后，在VS中新建Fortran工程，点项目-属性--配置属性-Fortran-Libraries，将Use Intel Math Kernal Library选项改为Parrallel即可运行。

c语言版本：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mkl_PARDISO.h"
#include "mkl_types.h"
// Define the format to printf MKL_INT values
#if !defined(MKL_ILP64)
#define IFORMAT "%i"
#else
#define IFORMAT "%lli"
#endif
int main (void)
{
    /* Matrix data. */
    MKL_INT n = 8;
    MKL_INT ia[9] = { 1, 5, 8, 10, 12, 15, 17, 18, 19};
    MKL_INT ja[18] = 
    { 1,    3,       6, 7,
         2, 3,    5,
            3,             8,
               4,       7,
                  5, 6, 7,
                     6,    8,
                        7,
                           8
    };
    double a[18] = 
    { 7.0,      1.0,           2.0, 7.0,
          -4.0, 8.0,      2.0,
                1.0,                     5.0,
                     7.0,           9.0,
                          5.0, 1.0, 5.0,
                              -1.0,      5.0,
                                   11.0,
                                         5.0
    };
    MKL_INT mtype = -2;       /* Real symmetric matrix */
    /* RHS and solution vectors. */
    double b[8], x[8];
    MKL_INT nrhs = 1;     /* Number of right hand sides. */
    /* Internal solver memory pointer pt, */
    /* 32-bit: int pt[64]; 64-bit: long int pt[64] */
    /* or void *pt[64] should be OK on both architectures */
    void *pt[64];
    /* PARDISO control parameters. */
    MKL_INT iparm[64];
    MKL_INT maxfct, mnum, phase, error, msglvl;
    /* Auxiliary variables. */
    MKL_INT i;
    double ddum;          /* Double dummy */
    MKL_INT idum;         /* Integer dummy. */
/* -------------------------------------------------------------------- */
/* .. Setup PARDISO control parameters. */
/* -------------------------------------------------------------------- */
    for ( i = 0; i < 64; i++ )
    {
        iparm[i] = 0;
    }
    iparm[0] = 1;         /* No solver default */
    iparm[1] = 2;         /* Fill-in reordering from METIS */
    iparm[3] = 0;         /* No iterative-direct algorithm */
    iparm[4] = 0;         /* No user fill-in reducing permutation */
    iparm[5] = 0;         /* Write solution into x */
    iparm[6] = 0;         /* Not in use */
    iparm[7] = 2;         /* Max numbers of iterative refinement steps */
    iparm[8] = 0;         /* Not in use */
    iparm[9] = 13;        /* Perturb the pivot elements with 1E-13 */
    iparm[10] = 1;        /* Use nonsymmetric permutation and scaling MPS */
    iparm[11] = 0;        /* Not in use */
    iparm[12] = 0;        /* Maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm[12] = 1 in case of inappropriate accuracy */
    iparm[13] = 0;        /* Output: Number of perturbed pivots */
    iparm[14] = 0;        /* Not in use */
    iparm[15] = 0;        /* Not in use */
    iparm[16] = 0;        /* Not in use */
    iparm[17] = -1;       /* Output: Number of nonzeros in the factor LU */
    iparm[18] = -1;       /* Output: Mflops for LU factorization */
    iparm[19] = 0;        /* Output: Numbers of CG Iterations */
    maxfct = 1;           /* Maximum number of numerical factorizations. */
    mnum = 1;         /* Which factorization to use. */
    msglvl = 1;           /* Print statistical information in file */
    error = 0;            /* Initialize error flag */
/* -------------------------------------------------------------------- */
/* .. Initialize the internal solver memory pointer. This is only */
/* necessary for the FIRST call of the PARDISO solver. */
/* -------------------------------------------------------------------- */
    for ( i = 0; i < 64; i++ )
    {
        pt[i] = 0;
    }
/* -------------------------------------------------------------------- */
/* .. Reordering and Symbolic Factorization. This step also allocates */
/* all memory that is necessary for the factorization. */
/* -------------------------------------------------------------------- */
    phase = 11;
    PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
             &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
    if ( error != 0 )
    {
        printf ("\nERROR during symbolic factorization: " IFORMAT, error);
        exit (1);
    }
    printf ("\nReordering completed ... ");
    printf ("\nNumber of nonzeros in factors = " IFORMAT, iparm[17]);
    printf ("\nNumber of factorization MFLOPS = " IFORMAT, iparm[18]);
/* -------------------------------------------------------------------- */
/* .. Numerical factorization. */
/* -------------------------------------------------------------------- */
    phase = 22;
    PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
             &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
    if ( error != 0 )
    {
        printf ("\nERROR during numerical factorization: " IFORMAT, error);
        exit (2);
    }
    printf ("\nFactorization completed ... ");
/* -------------------------------------------------------------------- */
/* .. Back substitution and iterative refinement. */
/* -------------------------------------------------------------------- */
    phase = 33;
    iparm[7] = 2;         /* Max numbers of iterative refinement steps. */
    /* Set right hand side to one. */
    for ( i = 0; i < n; i++ )
    {
        b[i] = 1;
    }
    PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
             &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error);
    if ( error != 0 )
    {
        printf ("\nERROR during solution: " IFORMAT, error);
        exit (3);
    }
    printf ("\nSolve completed ... ");
    printf ("\nThe solution of the system is: ");
    for ( i = 0; i < n; i++ )
    {
        printf ("\n x [" IFORMAT "] = % f", i, x[i]);
    }
    printf ("\n");
/* -------------------------------------------------------------------- */
/* .. Termination and release of memory. */
/* -------------------------------------------------------------------- */
    phase = -1;           /* Release internal memory. */
    PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
             &n, &ddum, ia, ja, &idum, &nrhs,
             iparm, &msglvl, &ddum, &ddum, &error);
    return 0;
}


以下是计算结果：

稀疏矩阵高性能求解器PARDISO的使用的图3

以上，即是MKL_PARDISO的具体使用例子。在实际中，MKL版本的PARDISO性能是非常高的，在某些稀疏矩阵的条件下，求解速度可以远超matlab。

以下，是某个条件数比较大但是方程组规模不大（61349个未知数）的matlab和PARDISO求解速度对比：

稀疏矩阵高性能求解器PARDISO的使用的图4

稀疏矩阵高性能求解器PARDISO的使用的图5

稀疏矩阵高性能求解器PARDISO的使用的图6

在该问题中，采用AMD3900X的处理器，MKL_PARDISO求解时间为43.6s，求解速度是matlab的两倍多。该稀疏矩阵的下载链接如下：

链接：https://pan.baidu.com/s/1lw-SgNlnSUAB4o9dwJa14Q

提取码：xmim

以上，即是高性能稀疏矩阵线性方程组求解器的简单介绍，感谢您的阅读，欢迎关注公众号  有限元术

二维码1.png
康国政老师的书《非线性本构关系在ABAQUS中的实现》第7章的几处疑似印刷错误 维得WD1117数据手册,引脚图,封装规格
(4条)
默认 最新
寒老师，感谢分享。正好一个问题请教您，之前在找文献的时候看过一些uel，umat源代码，就是用的这个PARDISO求解器求解一个扩充自由度。然后download下来后运行失败，提示是这个求解器（包？）无法使用。请教您如何在子程序里调用这个求解器包呢。
评论 点赞
学习了
评论 点赞

查看更多评论 >

点赞 10 评论 4 收藏 2
关注