fortranで行列関係の関数をいくつか作ったので、まとめておく。といっても、MATLABの良く使いそうなやつをいくつか真似して作ってみただけ。
fortranは数値計算は何も考えなくても早いけども、ちょっとしたことが不便すぎる。いろいろツールを作っておかないとそのちょっとしたことで時間を取られてしまうし、コードも汚くなってしまう。まだfortran77のソースも良く見るけど、fortran90, fortran95の機能を知りながら使う気には到底なれないなあ。
それにしてもMATLABで便利なのって、結局、転置行列X'とか行列積X*Yが表現できることだったり、変数がそのまま残ることだったりするような気がする（速度を出すにはテクニックが必要で、まだ全然使いこなせない）。
あとは、MATLABのsave関数的なものがあると便利そうではある。
! ---------------------------------------------------------
module matrix
! ---------------------------------------------------------
! inv_lapack
! Ai = inv_lapack(A) : return inverse of A(m,n) using lapack.
!
! whos
! call whos(X) : shows size, bytes, and type of X.
! only real(8) and integer are supported.
!
! linspace
! X = linspace(x0,x1,n) : return linearly spaced vector just same as MATLAB.
!
! eye, ones, zeros (same as MATLAB, Octave)
! X = eye(n,[m]) : get identity matrix
! X = ones(n,[m]) : get matrix of all ones
! X = zeros(n,[m]) : get matrix of all zeros
!
! show_dmat
! call show_dmat(X,[fmt]) : show real(8) X(n,m).
!
! ---------------------------------------------------------
! Note:
! build: gfortran matrix.f90 main.f90 -llapack
! ---------------------------------------------------------
! Created by sako 2013/05/14
! ---------------------------------------------------------
implicit none
interface whos
module procedure whos_d1, whos_d2, whos_i1, whos_i2
end interface
contains
! --------------------------------
! estimate inverse of array A(m,n).
! Ai = inv_lapack(A)
! --------------------------------
function inv_lapack(A) result (AI)
real(8), intent(in) :: A(:,:)
real(8) :: AI(size(A,1),size(A,2))
integer, parameter :: NB = 64
integer :: info
integer :: m, n
! for lapack
integer :: LWORK
integer, allocatable :: IPIV(:) ! dimension N
real(8), allocatable :: WORK(:) ! dimension LWORK
! get size of A(m,n)
m = size(A,1)
n = size(A,2)
LWORK = n*NB
allocate(IPIV(n))
allocate(WORK(LWORK))
AI = A
call DGETRF(m, n, AI, m, IPIV(1:min(m,n)), info) ! factorize
call DGETRI(n, AI, m, IPIV, WORK, LWORK, info) ! inverse
deallocate(IPIV)
deallocate(WORK)
if (INFO==0) then
!write(*,*) "succeeded."
endif
if (INFO/=0) then
print *, 'error: failed to estimate inverse at "inv_lapack".'
endif
return
end function
! --------------------------------
! show variable size, bytes, type
! like whos of MATLAB. (call whos(X))
! Size Bytes Kind
! ==== ===== ====
! 1x100 800 double
! --------------------------------
subroutine whos_d1(x)
real(8), intent(in) :: x(:)
call show_whos(1, size(x), kind(x), 'double')
return
end subroutine whos_d1
subroutine whos_d2(x)
real(8), intent(in) :: x(:,:)
call show_whos(size(x,1), size(x,2), kind(x), 'double')
return
end subroutine whos_d2
subroutine whos_i1(x)
integer, intent(in) :: x(:)
call show_whos(1, size(x), kind(x), 'integer')
return
end subroutine whos_i1
subroutine whos_i2(x)
integer, intent(in) :: x(:,:)
call show_whos(size(x,1), size(x,2), kind(x), 'integer')
return
end subroutine whos_i2
subroutine show_whos(n, m, b, valtype)
integer, intent(in) :: n, m, b
character(len=*), intent(in) :: valtype
character(len=11) :: strn, strm
character(len=16) :: strb
write(strn,'(i11)') n
write(strm,'(i11)') m
write(strb,'(i16)') b*n*m ! bytes
write(*,'(a48)') " Size Bytes Kind"
write(*,'(a48)') " ==== ===== ===="
write(*,'(a11,"x",a11,1x,a16,1x,a7)') adjustr(strn), adjustl(strm), adjustr(strb), adjustr(valtype)
return
end subroutine
! --------------------------------
! return linearly spaced matrix.
! just same as linspace of MATLAB.
! --------------------------------
function linspace(x0,x1,nd) result(X)
real(8), intent(in) :: x0, x1
integer, intent(in), optional :: nd
real(8), allocatable :: X(:)
real(8) :: dx
integer :: n, i, ierr
if (present(nd)) then
n = nd
else
n = 100
endif
! error check
if (n < 1) n = 1 ! force to be 1
allocate(X(n), stat=ierr)
if (ierr /= 0) then
print *, 'error: failed to allocate memory at "linspace".'
stop
endif
! special case for n==1
if (n == 1) then
X(1) = x1
return
endif
dx = (x1-x0)/dble(n-1)
do i = 1,n
X(i) = x0 + dx*dble(i-1)
enddo
return
end function linspace
! --------------------------------
! return idendity matrix E(n[,m])
! --------------------------------
function eye(n,m) result(E)
!implicit none
integer, intent(in) :: n
integer, intent(in), optional :: m
real(8), allocatable :: E(:,:)
integer :: i, n2, ierr
if (present(m)) then
n2 = m
else
n2 = n
endif
allocate(E(n,n2), stat=ierr)
if (ierr /= 0) then
print *, 'error: failed to allocate memory at "matrix_eye".'
stop
endif
! initialize
E(1:n,1:n2) = 0.d0
do i = 1,n
E(i,i) = 1.d0
enddo
return
end function eye
! --------------------------------
! return matrix X(n[,m]) of all ones
! --------------------------------
function ones(n,m) result(X)
!implicit none
integer, intent(in) :: n
integer, intent(in), optional :: m
real(8), allocatable :: X(:,:)
integer :: i, n2, ierr
if (present(m)) then
n2 = m
else
n2 = n
endif
allocate(X(n,n2), stat=ierr)
if (ierr /= 0) then
print *, 'error: failed to allocate memory at "matrix_ones".'
stop
endif
! initialize
X(1:n,1:n2) = 1.d0
return
end function ones
! --------------------------------
! return matrix X(n[,m]) of all zeros
! --------------------------------
function zeros(n,m) result(X)
!implicit none
integer, intent(in) :: n
integer, intent(in), optional :: m
real(8), allocatable :: X(:,:)
integer :: i, n2, ierr
if (present(m)) then
n2 = m
else
n2 = n
endif
allocate(X(n,n2), stat=ierr)
if (ierr /= 0) then
print *, 'error: failed to allocate memory at "matrix_zeros".'
stop
endif
! initialize
X(1:n,1:n2) = 0.d0
return
end function zeros
! --------------------------------
! show matrix A(m,n)
! --------------------------------
subroutine show_dmat(A,fmtf)
!implicit none
real(8), intent(in) :: A(:,:)
character(len=*), intent(in), optional :: fmtf
character(len=25) :: fmtc
integer :: m, n, i, j
m = size(A,1); n = size(A,2)
if (present(fmtf)) then
fmtc = fmtf
else
fmtc = '(f13.4," ")'
endif
do i=1,m
do j=1,n
write(*,fmtc, advance='no') A(i,j)
enddo
write(*,*) ""
enddo
return
end subroutine show_dmat
end module matrix
