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
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
