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