非対称，構造対称マトリクスの解法 MKL PARADISO

備忘録。
水熱連成で，非対称マトリクスになってしまう場合について。
MKLのParadisoで非対称マトリクスを解くためのサブルーチンの一例。

サンプルコード：

INCLUDE 'mkl_pardiso.f90'
subroutine MainMKL()
use MKL_PARDISO
real(8) :: A(5,5)
real(8) :: val(13)
integer :: col(13)
integer :: ptrB(5)
integer :: ptrE(5)
integer :: m
real(8) :: alpha

!MKL PARAMETERS
type(mkl_pardiso_handle) :: pt(64)
integer :: maxfct=1, mnum=1, mtype=11, phase=13, n, nrhs=1, error, msglvl=0
integer :: iparm(64)
integer :: nnz
integer :: ncol
integer :: idum
real(8) :: ddum
real(8), allocatable :: value(:)
real(8), allocatable :: b(:)
real(8), allocatable :: x(:)
integer, allocatable :: columns(:)
integer, allocatable :: rowindex(:)
integer, allocatable :: perm(:)

integer :: i, j, k
integer :: nzidx
integer :: MAX_OF_NODES

pt(:)=mkl_pardiso_handle(0)
iparm(1)=0

A=0.D0
A(1,1)=1.D0
A(1,2)=-1.D0
A(1,4)=-3.D0
A(2,1)=-2.D0
A(2,2)=5.D0
A(3,3)=4.D0
A(3,4)=6.D0
A(3,5)=4.D0
A(4,1)=-4.D0
A(4,3)=2.D0
A(4,4)=7.D0
A(5,2)=8.D0
A(5,5)=-5.D0


!n x nを定義
MAX_OF_NODES=5
n=MAX_OF_NODES


!rowIndexを定義
allocate(rowindex(1:(n+1)))
allocate(b(1:n))
allocate(x(1:n))
allocate(perm(1:n))
rowindex(1)=1
nnz=0
do i = 1, n
ncol=0
do j = 1, n
if(A(i,j)/=0.D0) nnz=nnz+1
if((A(i,j)/=0.D0)) ncol=ncol+1
enddo
rowindex(i+1)=rowindex(i)+ncol
enddo
allocate(value(1:nnz))
allocate(columns(1:nnz))

nzidx=0
do i =1, n
do j=1, n
if(A(i,j)/=0.D0)then
nzidx=nzidx+1
value(nzidx)=A(i,j)
columns(nzidx)=j
endif
enddo
enddo

do i=1, nnz
write(*,*) value(i), columns(i)
enddo
!-------------------------------------------------------------
!Factorization & Solution
!-------------------------------------------------------------
phase = 13 ! only factorization

!ここで，OutFvecを代入
do i = 1, n
b(i) = 1.D0
end do

call pardiso (pt, maxfct, mnum, mtype, phase, n, value, rowindex, columns,&
& perm, nrhs, iparm, msglvl, b, x, error)
write(*,*) 'error', error
do i = 1, n
write(*,*) ' x(',i,') = ', x(i), b(i)
end do


b=matmul(A, x)
do i = 1, n
write(*,*) ' b(',i,') = ', b(i)
end do

!-------------------------------------------------------------
!Termination and release of memory
!-------------------------------------------------------------
phase = -1 ! release internal memory
CALL pardiso (pt, maxfct, mnum, mtype, phase, MAX_OF_NODES, value, rowindex, columns,&
& perm, nrhs, iparm, msglvl, b, x, error)
return
end subroutine MainMKL