作者shockbon (bonbon)
看板Fortran
標題[問題]反矩陣
時間Wed May 13 22:19:54 2009
我是用彭國倫FORTRAN95書裡面的"inverse.f90"作修該
module LinearAlgebra
implicit none
contains
! 求反矩陣
subroutine inverse(A,IA)
implicit none
real :: A(:,:), IA(:,:)
real, allocatable :: B(:,:)
integer :: i,j,N
N = size(A,1)
allocate(B(N,N))
! 先把IA設定成單位矩陣
forall(i=1:N,j=1:N,i==j) IA(i,j)=1.0
forall(i=1:N,j=1:N,i/=j) IA(i,j)=0.0
! 保存原先的矩陣A, 使用B來計算
B=A
! 把B化成對角線矩陣(除了對角線外,都為0)
call Upper(B,IA,N) ! 先把B化成上三角矩陣
call Lower(B,IA,N) ! 再把B化成下三角矩陣
! 求解
forall(i=1:N) IA(i,:)=IA(i,:)/B(i,i)
return
end subroutine
! 輸出矩陣的副程式
subroutine output(matrix)
implicit none
real :: matrix(:,:)
integer :: m,n,i
character(len=20) :: for='(??(1x,f6.3))'
m = size(matrix,1)
n = size(matrix,2)
! 用字串來設定輸出格式
write( FOR(2:3), '(I2)' ) N
do i=1,N
write( *, FMT=FOR ) matrix(i,:)
end do
return
end subroutine output
! 求上三角矩陣的副程式
subroutine Upper(M,S,N)
implicit none
integer :: N
real :: M(N,N)
real :: S(N,N)
integer :: I,J
real :: E
do I=1,N-1
do J=I+1,N
E=M(J,I)/M(I,I)
M(J,I:N)=M(J,I:N)-M(I,I:N)*E
S(J,:)=S(J,:)-S(I,:)*E
end do
end do
return
end subroutine Upper
! 求下三角矩陣的副程式
subroutine Lower(M,S,N)
implicit none
integer :: N
real :: M(N,N)
real :: S(N,N)
integer :: I,J
real :: E
do I=N,2,-1
do J=I-1,1,-1
E=M(J,I)/M(I,I)
M(J,1:N)=M(J,1:N)-M(I,1:N)*E
S(J,:)=S(J,:)-S(I,:)*E
end do
end do
return
end subroutine Lower
end module
! 求解聯立式
program main
use LinearAlgebra
implicit none
integer, parameter :: N=3 ! Size of Matrix
real :: A(N,N) = (/1,2,3,4,5,6,7,8,8 /)
real :: IA(N,N)
integer :: i
write(*,*) "原矩陣"
call output(A)
call inverse(A,IA)
write(*,*) "反矩陣"
call output(IA)
stop
end program
我要得到n維反矩陣
小弟程式跑出來的結果反矩陣一直都是零
不知道哪裡出了問題
有大大可以幫忙看一下嗎謝謝
以下程式暫時用3x3作測試
而g.dat 使用以下數值
1
2
3
4
5
6
7
8
8
算出結果應該如下
ans =
-2.6667 2.6667 -1.0000
3.3333 -4.3333 2.0000
-1.0000 2.0000 -1.0000
修改後程式碼
program inv
implicit none
integer, parameter :: N=3 ! Size of Matrix
real :: A(N,N)
real :: IA(N,N)
integer :: i,yy,zz
open(unit=101,file='g.dat',status='old')
!write(101,*) "Matrix:"
DO yy = 1, 3
read(101,*) (A(yy,zz), zz = 1, 3)
END DO
open(30, file="result.dat",status="unknown")
write(30,*) "原矩陣"
DO yy = 1, 36
do zz= 1, 36
write(30,*) yy,zz,A(yy,zz)
END DO
END DO
call output(A)
call output(IA)
call inverse(A,IA)
write(30,*) "反矩陣"
DO yy = 1, 36
do zz= 1, 36
write(30,*) yy,zz,IA(yy,zz)
END DO
END DO
stop
end program
! 輸出矩陣的副程式
subroutine output(matrix)
real :: matrix(:,:)
integer :: m,n,i
character(len=20) :: for='(??(1x,f6.3))'
m = size(matrix,1)
n = size(matrix,2)
! 用字串來設定輸出格式
write( FOR(2:3), '(I2)' ) N
do i=1,N
write( *, FMT=FOR ) matrix(i,:)
end do
return
end subroutine output
! 求反矩陣
subroutine inverse(A,IA)
real :: A(:,:), IA(:,:)
real, allocatable :: B(:,:)
integer :: i,j,N
N = size(A,1)
allocate(B(N,N))
! 先把IA設定成單位矩陣
forall(i=1:N,j=1:N,i==j) IA(i,j)=1.0
forall(i=1:N,j=1:N,i/=j) IA(i,j)=0.0
! 保存原先的矩陣A, 使用B來計算
B=A
! 把B化成對角線矩陣(除了對角線外,都為0)
call Upper(B,IA,N) ! 先把B化成上三角矩陣
call Lower(B,IA,N) ! 再把B化成下三角矩陣
! 求解
forall(i=1:N) IA(i,:)=IA(i,:)/B(i,i)
return
end subroutine
! 求上三角矩陣的副程式
subroutine Upper(M,S,N)
integer :: N
real :: M(N,N)
real :: S(N,N)
integer :: I,J
real :: E
do I=1,N-1
do J=I+1,N
E=M(J,I)/M(I,I)
M(J,I:N)=M(J,I:N)-M(I,I:N)*E
S(J,:)=S(J,:)-S(I,:)*E
end do
end do
return
end subroutine Upper
! 求下三角矩陣的副程式
subroutine Lower(M,S,N)
integer :: N
real :: M(N,N)
real :: S(N,N)
integer :: I,J
real :: E
do I=N,2,-1
do J=I-1,1,-1
E=M(J,I)/M(I,I)
M(J,1:N)=M(J,1:N)-M(I,1:N)*E
S(J,:)=S(J,:)-S(I,:)*E
end do
end do
return
end subroutine Lower
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 140.115.66.121
※ 編輯: shockbon 來自: 140.115.66.121 (05/13 22:27)
1F：推 GP03:直接用IMSL的函式試試看? 05/13 22:31