初学Fortran，根据算法编译的程序结果显示都是0，烦请指点.

Just2Moment

[Fortran] 纯文本查看 复制代码

module bicgstabj
!----------------------------------------module coment
! Version : V1.0
! Coded by : syz
! Date : 2016-4-1
!-----------------------------------------------------
! Description : 稳定双共轭梯度法
!
!-----------------------------------------------------
! Parameters :
! 1. IMAX--最大允许迭代次数
! 2. tol--误差容限
!
! Contains :
! 1. solve 迭代法方法函数
! 2.
! 3. dr(r,N) 计算向量长度平方函数
! 4. Ar(A,r,N) 计算矩阵乘以向量函数，返回向量
! 5. rAr(A,r,N) 计算（Ar,r）函数
!-----------------------------------------------------
! Post Script :
! 1. 里面的三个函数，可以简化程序，同时可以用在其他地方
! 2.
!-----------------------------------------------------
implicit real*8(a-z)
integer::IMAX=200
real*8::tol=1d-7
contains
subroutine solve(A,b,x,x0,N)
!---------------------------------subroutine comment
! Version : V1.0
! Coded by : syz
! Date : 2016-4-1
!-----------------------------------------------------
! Purpose : 稳定双共轭梯度法
! 用于计算方程AX=b
!-----------------------------------------------------
! Input parameters :
! 1. A,b 意义即 AX=b
! 2. x0 迭代初值
! 3. N 方程的维数
! Output parameters :
! 1. x 方程的解
! 2.
! Common parameters :
!
!----------------------------------------------------
implicit real*8(a-z)
integer::N
integer::i,j,k
real*8::A(N,N),b(N),x(N),x0(N)
real*8::r0(N),r1(N),r(N),r00(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
r00=r0
do k=1,IMAX
tmp1=rAr(r00,r0,N)
tmp2=rAr(r00,Ar(A,p0,N),N)
afa=tmp1/tmp2
x3=x1+afa*p0
r3=r0-afa*Ar(A,p0,N)
oumiga=rAr(A,r3,N)/dr(Ar(A,r3,N),N)
x2=x3+oumiga*r3
!记录迭代中间值
write(102,502)k,x2
502 format(I3,4F12.8)
!如果r0 接近于，则停止迭代
dr_s=dsqrt(dr(r0,N))
if (dr_s<tol) exit
r1=r3-oumiga*Ar(A,r,N)
tmp3=rAr(r1,r00,N)
beta=(afa/oumiga)*(tmp3/tmp1)
p1=r1+beta*(p0-oumiga*Ar(A,p0,N))
!全部更新
r0=r1
p0=p1
x1=x2
end do
x=x2
end subroutine solve
function dr(r,N)
!---------------------------------subroutine comment
! Version : V1.0
! Coded by : syz
! Date : 2016-4-1
!-----------------------------------------------------
! Purpose : 计算向量长度平方 (r,r)
!
!-----------------------------------------------------
! Input parameters :
! 1. r 向量
! 2. N 维数
! Output parameters :
! 1. dr 长度平方
! 2.
! Common parameters :
!
!----------------------------------------------------
implicit real*8(a-z)
integer::N,i
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)
!---------------------------------subroutine comment
! Version : V1.0
! Coded by : syz
! Date : 2016-4-1
!-----------------------------------------------------
! Purpose : !计算 A*r,返回N 维向量
!
!-----------------------------------------------------
! Input parameters :
! 1. r 向量
! 2. N 维数
! 3. A 矩阵
! Output parameters :
! 1. Ar 返回向量
! 2.
! Common parameters :
!
!----------------------------------------------------
implicit real*8(a-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)
!---------------------------------subroutine comment
! Version : V1.0
! Coded by : syz
! Date : 2016-4-1
!-----------------------------------------------------
! Purpose : 向量点乘 v1v2=v1(1)*v2(1)+v1(2)*v(2)+...
!
! Post Script :
! 1.
! 2.
! 3.
!
!-----------------------------------------------------
! Input parameters :
! 1. v1,v2 向量
! 2. N 向量维数
! Output parameters :
! 1. v1,v2 向量点乘值
! 2.
! Common parameters :
! 1.
! 2.
!----------------------------------------------------
implicit real*8(a-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)
!---------------------------------subroutine comment
! Version : V1.0
! Coded by : syz
! Date : 2016-4-1
!-----------------------------------------------------
! Purpose : !计算（Ar,r）,返回标量
!
!-----------------------------------------------------
! Input parameters :
! 1. r 向量
! 2. N 维数
! 3. A 矩阵
! Output parameters :
! 1. Ar 返回标量
! 2.
! Common parameters :
!----------------------------------------------------
implicit real*8(a-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 bicgstabj
program main
!--------------------------------------program comment
! Version : V1.0
! Coded by : syz
! Date : 2016-4-1
!-----------------------------------------------------
! Purpose : 采用稳定双共轭梯度法计算线性方程
!
!-----------------------------------------------------
! In put data files :
! 1.
! 2.
! Output data files :
! 1. Im_result.txt 计算的中间数据
! 2. result.txt 计算结果
!-----------------------------------------------------
use bicgstabj
implicit real*8(a-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=(/1d0,1d0,1d0,1d0/)
!系数
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

fcode

1. 你使用什么编译器？
2. 我使用 gfortran 编译上述程序，会报错。
3. 修改代码后，得到结果
  x(1)=    22.85115620
  x(2)=   160.18391519
  x(3)=  -140.46884071
  x(4)=   209.15349955
不知是否正确？
4. 修改方法是：
53 行 real*8::x1(N),x2(N),x3(N),r3(N)!//要定义x3,r3
133 行 integer::i,N,j!//要定义j
5. 强烈建议，每个程序都写上 Implicit None
请回答我的第一个问题，谢谢