逆行列の数値計算

Lapackを用いて数値計算で逆行列を求めます。
基本的に数値的にあらわに逆行列を求める意味はそんなにありません。
連立一次方程式を解きたければ、行列のLU分解を通じて解けば良いです。この場合、逆行列を求めることはしません。

そうはいっても使ってみたい時があります。そんな状況を考えましょう。

Lapack, もしくはインテルのマス・カーネル・ライブラリ(MKL)を用います。

用いるルーチンは

?getrf, ?getri

の2つです。

ここに用意したのは配列の大きさNと配列Aを入れると逆行列を出力するプログラムです。
倍精度、複素倍精度のルーチンを書きました。

ifort -mkl main.f90

でコンパイルしてください。

program main
  implicit none
  integer::i,j,N
  complex(kind(0d0)),allocatable::A(:,:),B(:,:)
 
  N=3
  allocate(A(1:N,1:N),B(1:N,1:N))
  A(1,1:3)=(/(2d0,0d0),(3d0,0d0),(0d0,0d0)/)
  A(2,1:3)=(/(2d0,0d0),(9d0,0d0),(1d0,3d0)/)
  A(3,1:3)=(/(3d0,4d0),(5d0,0d0),(3d0,1d0)/)

  B=A
 
  do i=1,N
     do j=1,N
        write(6,'(2f15.10,$)')A(i,j)
     enddo
     write(6,*)
  enddo
  write(6,*)

  call zinvmat(size(A,1),A)
 
  do i=1,N
     do j=1,N
        write(6,'(2f15.10,$)')A(i,j)
     enddo
     write(6,*)
  enddo
  write(6,*)

  A=matmul(B,A)
 
  do i=1,N
     do j=1,N
        write(6,'(2f15.10,$)')A(i,j)
     enddo
     write(6,*)
  enddo
  write(6,*)

  stop
end program main
   
subroutine invmat(N,A)
  implicit none
  integer,intent(in)::N
  double precision,intent(inout)::A(1:N)

  integer::info,lwork,ipiv(1:N)
  double precision::tw(1:1)
  double precision,allocatable::work(:)
 
  info=0
  call dgetrf(N,N,A,N,ipiv,info)
  if(info.ne.0)then
     write(6,*)"Error at inverse matrix ?getrf, info",info
     stop
  endif
 
  !size query
  call dgetri(N,A,N,ipiv,tw,-1,info)
  lwork=nint(tw(1))
  allocate(work(1:lwork+1))
  work=0d0
 
  call dgetri(N,A,N,ipiv,work,lwork,info)
  if(info.ne.0)then
     write(6,*)"Error at inverse matrix ?dgetri, info",info
     stop
  endif
 
  return
end subroutine invmat


subroutine zinvmat(N,A)
  implicit none
  integer,intent(in)::N
  complex(kind(0d0)),intent(inout)::A(1:N)

  integer::info,lwork,ipiv(1:N)
  complex(kind(0d0))::tw(1:1)
  complex(kind(0d0)),allocatable::work(:)
 
  info=0
  call zgetrf(N,N,A,N,ipiv,info)
  if(info.ne.0)then
     write(6,*)"Error at inverse matrix ?getrf, info",info
     stop
  endif
 
  !size query
  call zgetri(N,A,N,ipiv,tw,-1,info)
  lwork=nint(dble(tw(1)))
  allocate(work(1:lwork+1))
  work=dcmplx(0d0,0d0)
 
  call zgetri(N,A,N,ipiv,work,lwork,info)
  if(info.ne.0)then
     write(6,*)"Error at inverse matrix ?dgetri, info",info
     stop
  endif
 
  return
end subroutine zinvmat

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です