[問題]反矩陣

看板Fortran作者shockbon (bonbon)時間15年前 (2009/05/13 22:19)推噓1(1推 0噓 0→)

留言1則, 1人參與討論串1/1

我是用彭國倫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)