[Fortran] 纯文本查看 复制代码
module conj
implicit real*8(a-h,o-z)
integer :: IMAX=200
real*8 :: tol=1d-7
contains
subroutine solve(A,b,x,x0,N)
implicit real*8(a-h,o-z)
integer :: i,j,k
real*8 :: A(N,N),b(N),x(N),x0(N)
real*8 :: r0(N),r1(N),p0(N),p1(N)
real*8 :: x1(N),x2(N)
write(102,501)
501 format(//,18x,'共轭梯度法中间结果',//)
x1=x0
r0=b-Ar(A,x1,N)
p0=r0
do k=1,IMAX
tmp1=dr(r0,N)
tmp2=rAr(A,p0,N)
afa=tmp1/tmp2
x2=x1+afa*p0
!记录迭代中间值
write(101,502)k,x2
502 format(I3,4F12.8)
!
!
dr_s=dsqrt(dr(r0,N))
if (dr_s<tol) exit
r1=r0-afa*Ar(A,p0,N)
tmp3=dr(r1,N)
beta=tmp3/tmp1
p1=r1+beta*p0
r0=r1
p0=p1
x1=x2
end do
x=x2
end subroutine solve
function dr(r,N)
implicit real*8(a-h,o-z)
real*8 :: r(N),dr
s=0
do i=1,N
s=s+r(i)**2
end do
dr=s
end function dr
function Ar(A,r,N)
implicit real*8(a-h,o-z)
integer :: i,N
real*8 :: A(N,N),r(N),temp(N),Ar(N)
temp=0
do i=1,N
do j=1,n
temp(i)=temp(i)+A(i,j)*r(j)
end do
end do
Ar=temp
end function ar
function v1v2(v1,v2,N)
implicit real*8(a-h,o-z)
integer :: n
real*8 :: v1(n),v2(n)
integer :: i
v1v2=0
do i=1,n
v1v2=v1v2+v1(i)*v2(i)
end do
end function
function rAr(A,r,N)
implicit real*8(a-h,o-z)
integer :: i,N
real*8 :: A(N,N),r(N),temp(N)
temp=Ar(A,r,N)
rAr=v1v2(r,temp,N)
end function rAr
end module conj
program main
use conj
implicit real*8(a-h,o-z)
integer,parameter :: N=4
real*8 :: A(N,N),b(N),x(N),x0(N)
open(unit=101,file='result.txt')
open(unit=102,file='Im_result.txt')
!迭代初值
x0=(/0d0,0d0,0d0,0d0/)
b=(/62d0,87d0,91d0,90d0/)
A=reshape((/5d0,7d0,6d0,5d0,&
7d0,10d0,8d0,7d0,&
6d0,8d0,10d0,9d0,&
5d0,7d0,9d0,10d0/),(/4,4/))
call solve(A,b,x,x0,N)
write(101,501)x
501 format(/,T10,'共轭梯度法',//,&
2x,'x(1)=',F15.8,/,&
2x,'x(2)=',F15.8,/,&
2x,'x(3)=',F15.8,/,&
2x,'x(4)=',F15.8)
end program main