Fortran：用csr存储格式并使用pardiso求解稀疏矩阵

使用链表，将大型稀疏对称矩阵存储为CSR格式，并用mkl函数库中的pardiso求解器进行快速求解
Program CSR
    !// author: luk
    !// qq: 735343320
    implicit none
    integer, parameter        :: n = 5  !// depend on your equation
    integer                   :: i, j, mm, tmp, nn, fileid, first, num
    real(kind=8)              :: matrix(n,n), b(n), x(n)
    real(kind=8), allocatable :: aa(:)
    integer, allocatable      :: ia(:), ja(:)
  
    !// =====================the variables of pardiso====================
    integer(kind=8)           :: pt(64)   !// kind=8:x64; kind=4:win32
    integer                   :: maxfct, mnum, mtype, phase, nrhs, error, msglvl
    integer, allocatable      :: perm(:)
    integer                   :: iparm(64)
    integer, external         :: mkl_get_max_threads
    !// =====================the variables of pardiso====================
  
    type                      :: node
        real(kind=8)          :: s
        integer               :: k1, k2
        type (node), pointer  :: next
    end type node 
    type (node), pointer      :: head, tail, p, q  
  
    !// input data
    allocate( perm( n ) )
    open ( newunit = fileid, file = 'SparseMatrix.txt' )  !// The sparse matrix data needs to be given by itself
    Do i = 1, n
        do j = 1, n
            read ( fileid,* ) matrix(i,j)
        end do
    End do
    read( fileid,* ) b
    close ( fileid )

    !// =========================store data by CSR format=======================
    first = 1
    If ( .not.associated(p) ) then
        allocate( p )
        q    => p
        nullify( p%next )
        q%k2 = first
        tmp  = q%k2
    End if
  
    num = 0  !// calculate the number of no-zero
    Do i =  1, n
        mm = 0  !// calculate the position of no-zero in every row
        do j = i, n
            If ( matrix(i,j) /= 0.0 ) then
                if ( .not.associated(head) ) then
                    allocate( head )
                    tail    => head
                    nullify( tail%next )
                    tail%s  = matrix(i,j)
                    tail%k1 = j
                    num     = num + 1
                    mm      = mm + 1
                else
                    allocate( tail%next )
                    tail    => tail%next
                    tail%s  = matrix(i,j)
                    tail%k1 = j
                    num     = num + 1  
                    mm      = mm + 1  
                    nullify( tail%next )
                end if
            End if
        end do
        if ( associated(p) ) then
            allocate( q%next )
            q    => q%next
            q%k2 = tmp + mm
            tmp  = q%k2
            nullify( q%next )
        end if
    End do
  
    allocate( aa(num), ja(num), ia(n+1) )  !// store the no-zero element
    !// output a and ja
    tail => head
    j = 1
    Do 
        if ( .not.associated(tail) ) exit
        aa(j) = tail%s
        ja(j) = tail%k1
        j     = j + 1
        tail  => tail%next
    End do
  
    !// output ia
    q => p
    j = 1
    Do 
        if ( .not.associated(q) ) exit
        ia(j) = q%k2
        j     = j + 1
        q     => q%next
    End do
  
    nullify( head, tail, p, q )
    !// =========================store data by CSR format=======================
    write ( *,'(a)' ) '>> Original data'
    write ( *,'(5f7.2)' ) matrix
    write ( *,'(a)' ) '>> CSR data'
    write ( *,'(*(f7.2))' ) aa
    write ( *,'(a)' ) '>> B data'
    write ( *,'(*(f7.2))' ) b
    write ( *,'(a)' ) '>> Colnum of CSR data'
    write ( *,'(*(i4))' ) ja
    write ( *,'(a)' ) '>> The position of CSR data'
    write ( *,'(*(i4))' ) ia
  
    !// ===========================solve equations=============================
    pt = 0  !// pointer initialization
    maxfct = 1; mnum = 1; mtype = -2  !// mtype = -2: symmetric nonpositive definite matrix, mtype = 2: symmetric positive definite matrix
    perm = 0; nrhs = 1
    x = 0.d0
    iparm(1) = 0  !// iparm use default values
    iparm(3) = mkl_get_max_threads()  !// iparm(3): parallel threads are obtained by MKL. In general, the number of cpus is the same as the number of threads
    error = 0  
    msglvl = 0
  
    phase = 12  !// LU decompose
    call pardiso( pt, maxfct, mnum, mtype, phase, n, aa, ia, ja, perm, nrhs, iparm, msglvl, b, x, error )
    phase = 33  !// solve equations
    call pardiso( pt, maxfct, mnum, mtype, phase, n, aa, ia, ja, perm, nrhs, iparm, msglvl, b, x, error )
    write ( *,'(a)' ) '>> The solution as follows ...... '
    do i = 1, n
        write( *,'(f7.4)' ) x(i)
    end do
  
    phase = -1
    call pardiso (pt, maxfct, mnum, mtype, phase, n, aa, ia, ja, perm, nrhs, iparm, msglvl, b, x, error)
  
    deallocate( aa, ja, ia, perm )
    !// ===========================stop solving equations=============================
  
End program CSR
  

!// remark:
!// windows: set mkl libraries or use command line with /Qmkl  [ifort /Qmkl pardiso.f90]
!// linux: use command line with -mkl   [ifort -mkl pardiso.f90]
!// phase = 12  !// LU decompose
!// call pardiso( pt, maxfct, mnum, mtype, phase, n, aa, ia, ja, perm, nrhs, iparm, msglvl, b, x, error )
!// phase = 33  !// solve equations 
!// call pardiso( pt, maxfct, mnum, mtype, phase, n, aa, ia, ja, perm, nrhs, iparm, msglvl, b, x, error )
!// The top four rows can be replaced by the bottom two rows
!// phase = 13  !// solve equations 
!// call pardiso( pt, maxfct, mnum, mtype, phase, n, aa, ia, ja, perm, nrhs, iparm, msglvl, b, x, error )

这里给出测试矩阵文件SparseMatrix.txt, 元素如下：
   1.00  -1.00   0.00  -3.00   0.00
  -1.00   5.00   0.00   0.00   0.00
   0.00   0.00   4.00   6.00   4.00
  -3.00   0.00   6.00   7.00   0.00
   0.00   0.00   4.00   0.00  -5.00
右端项为：
  -13.00   
   9.00  
   56.00  
   43.00 
  -13.00
得其解为：
1
2
3
4
5