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

Just2Moment

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

查看代码

打印?

001	module bicgstabj

002	!----------------------------------------module coment

003	! Version : V1.0

004	! Coded by : syz

005	! Date : 2016-4-1

006	!-----------------------------------------------------

007	! Description : 稳定双共轭梯度法

008

!

009	!-----------------------------------------------------

010	! Parameters :

011	! 1. IMAX--最大允许迭代次数

012	! 2. tol--误差容限

013

!

014	! Contains :

015	! 1. solve 迭代法方法函数

016

! 2.

017	! 3. dr(r,N) 计算向量长度平方函数

018	! 4. Ar(A,r,N) 计算矩阵乘以向量函数，返回向量

019	! 5. rAr(A,r,N) 计算（Ar,r）函数

020	!-----------------------------------------------------

021	! Post Script :

022	! 1. 里面的三个函数，可以简化程序，同时可以用在其他地方

023

! 2.

024	!-----------------------------------------------------

025	implicit real*8(a-z)

026	integer::IMAX=200

027	real*8::tol=1d-7

028

contains

029	subroutine solve(A,b,x,x0,N)

030	!---------------------------------subroutine comment

031	! Version : V1.0

032	! Coded by : syz

033	! Date : 2016-4-1

034	!-----------------------------------------------------

035	! Purpose : 稳定双共轭梯度法

036	! 用于计算方程AX=b

037	!-----------------------------------------------------

038	! Input parameters :

039	! 1. A,b 意义即 AX=b

040	! 2. x0 迭代初值

041	! 3. N 方程的维数

042	! Output parameters :

043	! 1. x 方程的解

044

! 2.

045	! Common parameters :

046

!

047	!----------------------------------------------------

048	implicit real*8(a-z)

049	integer::N

050	integer::i,j,k

051	real*8::A(N,N),b(N),x(N),x0(N)

052	real*8::r0(N),r1(N),r(N),r00(N),p0(N),p1(N)

053	real*8::x1(N),x2(N)

054

!写入标题

055	write(102,501)

056	501 format(//,18x,'稳定双共轭梯度法中间结果',//)

057

x1=x0

058	r0=b-Ar(A,x1,N)

059

p0=r0

060

r00=r0

061	do k=1,IMAX

062	tmp1=rAr(r00,r0,N)

063	tmp2=rAr(r00,Ar(A,p0,N),N)

064	afa=tmp1/tmp2

065	x3=x1+afa*p0

066	r3=r0-afa*Ar(A,p0,N)

067	oumiga=rAr(A,r3,N)/dr(Ar(A,r3,N),N)

068	x2=x3+oumiga*r3

069

!记录迭代中间值

070	write(102,502)k,x2

071	502 format(I3,4F12.8)

072	!如果r0 接近于，则停止迭代

073	dr_s=dsqrt(dr(r0,N))

074	if (dr_s<tol) exit

075	r1=r3-oumiga*Ar(A,r,N)

076	tmp3=rAr(r1,r00,N)

077	beta=(afa/oumiga)*(tmp3/tmp1)

078	p1=r1+beta*(p0-oumiga*Ar(A,p0,N))

079

!全部更新

080

r0=r1

081

p0=p1

082

x1=x2

083

end do

084

x=x2

085	end subroutine solve

086	function dr(r,N)

087	!---------------------------------subroutine comment

088	! Version : V1.0

089	! Coded by : syz

090	! Date : 2016-4-1

091	!-----------------------------------------------------

092	! Purpose : 计算向量长度平方 (r,r)

093

!

094	!-----------------------------------------------------

095	! Input parameters :

096

! 1. r 向量

097

! 2. N 维数

098	! Output parameters :

099	! 1. dr 长度平方

100

! 2.

101	! Common parameters :

102

!

103	!----------------------------------------------------

104	implicit real*8(a-z)

105	integer::N,i

106	real*8::r(N),dr

107

s=0

108

do i=1,N

109	s=s+r(i)**2

110

end do

111

dr=s

112	end function dr

113	function Ar(A,r,N)

114	!---------------------------------subroutine comment

115	! Version : V1.0

116	! Coded by : syz

117	! Date : 2016-4-1

118	!-----------------------------------------------------

119	! Purpose : !计算 A*r,返回N 维向量

120

!

121	!-----------------------------------------------------

122	! Input parameters :

123

! 1. r 向量

124

! 2. N 维数

125

! 3. A 矩阵

126	! Output parameters :

127	! 1. Ar 返回向量

128

! 2.

129	! Common parameters :

130

!

131	!----------------------------------------------------

132	implicit real*8(a-z)

133	integer::i,N

134	real*8::A(N,N),r(N),temp(N),Ar(N)

135

temp=0

136

do i=1,N

137

do j=1,n

138	temp(i)=temp(i)+A(i,j)*r(j)

139

end do

140

end do

141

Ar=temp

142	end function ar

143	function v1v2(v1,v2,N)

144	!---------------------------------subroutine comment

145	! Version : V1.0

146	! Coded by : syz

147	! Date : 2016-4-1

148	!-----------------------------------------------------

149	! Purpose : 向量点乘 v1v2=v1(1)*v2(1)+v1(2)*v(2)+...

150

!

151	! Post Script :

152

! 1.

153

! 2.

154

! 3.

155

!

156	!-----------------------------------------------------

157	! Input parameters :

158	! 1. v1,v2 向量

159	! 2. N 向量维数

160	! Output parameters :

161	! 1. v1,v2 向量点乘值

162

! 2.

163	! Common parameters :

164

! 1.

165

! 2.

166	!----------------------------------------------------

167	implicit real*8(a-z)

168	integer::n

169	real*8::v1(n),v2(n)

170	integer::i

171

v1v2=0

172

do i=1,n

173	v1v2=v1v2+v1(i)*v2(i)

174

end do

175	end function

176	function rAr(A,r,N)

177	!---------------------------------subroutine comment

178	! Version : V1.0

179	! Coded by : syz

180	! Date : 2016-4-1

181	!-----------------------------------------------------

182	! Purpose : !计算（Ar,r）,返回标量

183

!

184	!-----------------------------------------------------

185	! Input parameters :

186

! 1. r 向量

187

! 2. N 维数

188

! 3. A 矩阵

189	! Output parameters :

190	! 1. Ar 返回标量

191

! 2.

192	! Common parameters :

193	!----------------------------------------------------

194	implicit real*8(a-z)

195	integer::i,N

196	real*8::A(N,N),r(N),temp(N)

197	temp=Ar(A,r,N)

198	rAr=v1v2(r,temp,N)

199	end function rAr

200	end module bicgstabj

201	program main

202	!--------------------------------------program comment

203	! Version : V1.0

204	! Coded by : syz

205	! Date : 2016-4-1

206	!-----------------------------------------------------

207	! Purpose : 采用稳定双共轭梯度法计算线性方程

208

!

209	!-----------------------------------------------------

210	! In put data files :

211

! 1.

212

! 2.

213	! Output data files :

214	! 1. Im_result.txt 计算的中间数据

215	! 2. result.txt 计算结果

216	!-----------------------------------------------------

217	use bicgstabj

218	implicit real*8(a-z)

219	integer,parameter::N=4

220	real*8 ::A(N,N),b(N),x(N),x0(N)

221	open(unit=101,file='result.txt')

222	open(unit=102,file='Im_result.txt')

223

!迭代初值

224	x0=(/1d0,1d0,1d0,1d0/)

225

!系数

226	b=(/62d0,87d0,91d0,90d0/)

227	A=reshape((/5d0,7d0,6d0,5d0,&

228	7d0,10d0,8d0,7d0,&

229	6d0,8d0,10d0,9d0,&

230	5d0,7d0,9d0,10d0/),(/4,4/))

231

! 调用函数

232	call solve(A,b,x,x0,N)

233	write(101,501)x

234	501 format(/,T10,'稳定双共轭梯度法',//,&

235	2x,'x(1)=',F15.8,/,&

236	2x,'x(2)=',F15.8,/,&

237	2x,'x(3)=',F15.8,/,&

238	2x,'x(4)=',F15.8)

239	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
请回答我的第一个问题，谢谢