珊瑚虫 发表于 2015-1-12 21:41
请把implicit real*8(a-z)
换成implicit real*8(a-h,o-z)
[mw_shl_code=fortran,true]module conj

3Q，完美解决！！！！！

请把implicit real*8(a-z)
换成implicit real*8(a-h,o-z)

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

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