我正在 Fortran 中研究双共轭梯度算法，并按照 Saad, Y.“稀疏线性系统的迭代方法”(普通 BiCG 方法)中的算法编写了完整的代码。但是，它没有收敛到所需的迭代次数，也没有返回正确的结果。
该算法在维基百科的“无预处理版本”中给出( http://en.wikipedia.org/wiki/Biconjugate_gradient_method#Unpreconditioned_version_of_the_algorithm )
我对 Fortran 还比较陌生，并且不明白为什么这没有按预期运行，因为据我所知，它的编码完全按照指定的方式进行。如果有人看到任何非正统的代码，或者算法中的错误，我将非常感激!
为了简单起见，我添加了一个测试矩阵:
!
!////////////////////////////////////////////////////////////////////////
!
! BiCG_main.f90
! Created: 19 February 2013 12:01
! By: Robin Fox
!
!////////////////////////////////////////////////////////////////////////
!
PROGRAM bicg_main
!
IMPLICIT NONE
!-------------------------------------------------------------------
! Program to implement the Bi-Conjugate Gradient method
! follows algorithm in Saad
!-------------------------------------------------------------------
!
COMPLEX(KIND(0.0d0)), DIMENSION(:,:), ALLOCATABLE ::A
COMPLEX(KIND(0.0d0)), DIMENSION(:), ALLOCATABLE ::b
COMPLEX(KIND(0.0d0)), DIMENSION(:), ALLOCATABLE ::x0, x0s
COMPLEX(KIND(0.0d0)), DIMENSION(:), ALLOCATABLE ::x, xs
COMPLEX(KIND(0.0d0)), DIMENSION(:), ALLOCATABLE ::p, ps
COMPLEX(KIND(0.0d0)) ::alpha, rho0, rho1, r_rs
COMPLEX(KIND(0.0d0)), DIMENSION(:), ALLOCATABLE ::r,rs, res_vec
COMPLEX(KIND(0.0d0)), DIMENSION(:), ALLOCATABLE ::Ax, ATx
COMPLEX(KIND(0.0d0)), DIMENSION(:), ALLOCATABLE ::Ap, Aps
COMPLEX(KIND(0.0d0)) ::beta
!
REAL(KIND(0.0d0)) ::tol,res, n2b, n2r0, rel_res
!
INTEGER ::n,i,j,k, maxit
!////////////////////////////////////////////////////////////////////////
!----------------------------------------------------------
n=2
ALLOCATE(A(n,n))
ALLOCATE(b(n))
A(1,1)=CMPLX(-0.73492,7.11486)
A(1,2)=CMPLX(0.024839,4.12154)
A(2,1)=CMPLX(0.274492957,3.7885537)
A(2,2)=CMPLX(-0.632557864,1.95397735)
b(1)=CMPLX(0.289619736,0.895562183)
b(2)=CMPLX(-0.28475616,-0.892163111)
!----------------------------------------------------------
ALLOCATE(x0(n))
ALLOCATE(x0s(n))
!Use all zeros initial guess
x0(:)=CMPLX(0.0d0,0.0d0)
DO i=1,n
x0s(i)=CONJG(x0(i))
END DO
ALLOCATE(Ax(n))
ALLOCATE(ATx(n))
ALLOCATE(x(n))
ALLOCATE(xs(n))
! Multiply matrix A with vector x0
DO i=1,n
Ax(i)=CMPLX(0.0,0.0)
DO j=1,n
Ax(i)=Ax(i)+A(i,j)*x0(j) !==Ax=A*x0
END DO
END DO
! Multiply matrix A^T with vector x0
DO i=1,n
ATx(i)=CMPLX(0.0,0.0)
DO j=1,n
ATx(i)=ATx(i)+CONJG(A(j,i))*x0s(j) !==A^Tx=A^T*x0
END DO
END DO
res=0.0d0
n2b=0.0d0
x=x0
ALLOCATE(r(n))
ALLOCATE(rs(n))
ALLOCATE(p(n))
ALLOCATE(ps(n))
!Initialise
DO i=1,n
r(i)=b(i)-Ax(i)
rs(i)=CONJG(b(i))-ATx(i)
p(i)=r(i) !p0=r0
ps(i)=rs(i) !p0s=r0s
END DO
DO i=1,n
n2b=n2b+(b(i)*CONJG(b(i)))
res=res+(r(i)*CONJG(r(i))) !== inner prod(r,r)
END DO
n2b=SQRT(n2b)
res=SQRT(res)/n2b
!Check that inner prod(r,rs) =/= 0
n2r0=0.0d0
DO i=1,n
n2r0=n2r0+r(i)*CONJG(rs(i))
END DO
IF (n2r0==0) THEN
res=1d-20 !set tol so that loop doesn't run (i.e. already smaller than tol)
PRINT*, "Inner product of r, rs == 0"
END IF
WRITE(*,*) "n2r0=", n2r0
!----------------------------------------------------------
ALLOCATE(Ap(n))
ALLOCATE(Aps(n))
ALLOCATE(res_vec(n))
tol=1d-6
maxit=50 !for n=720
k=0
!Main loop:
main: DO WHILE ((res>tol).AND.(k<maxit))
k=k+1
! Multiply matrix A with vector p
DO i=1,n
Ap(i)=CMPLX(0.0,0.0)
DO j=1,n
Ap(i)=Ap(i)+A(i,j)*p(j)
END DO
END DO
! Multiply matrix A^T with vector p
! N.B. transpose is also conjg.
DO i=1,n
Aps(i)=CMPLX(0.0,0.0)
DO j=1,n
Aps(i)=Aps(i)+CONJG(A(j,i))*ps(j)
END DO
END DO
rho0=CMPLX(0.0d0,0.0d0)
DO i=1,n
rho0=rho0+(r(i)*CONJG(rs(i)))
END DO
WRITE(*,*) "rho0=", rho0
rho1=CMPLX(0.0d0,0.0d0)
DO i=1,n
rho1=rho1+(Ap(i)*CONJG(ps(i)))
END DO
WRITE(*,*) "rho1=", rho1
!Calculate alpha:
alpha=rho0/rho1
WRITE(*,*) "alpha=", alpha
!Update solution
DO i=1,n
x(i)=x(i)+alpha*p(i)
END DO
!Update residual:
DO i=1,n
r(i)=r(i)-alpha*Ap(i)
END DO
!Update second residual:
DO i=1,n
rs(i)=rs(i)-alpha*Aps(i)
END DO
!Calculate beta:
r_rs=CMPLX(0.0d0,0.0d0)
DO i=1,n
r_rs=r_rs+(r(i)*CONJG(rs(i)))
END DO
beta=r_rs/rho0
!Update direction vectors:
DO i=1,n
p(i)=r(i)+beta*p(i)
END DO
DO i=1,n
ps(i)=rs(i)+beta*ps(i)
END DO
!Calculate residual for convergence check
! res=0.0d0
! DO i=1,n
! res=res+(r(i)*CONJG(r(i))) !== inner prod(r,r)
! END DO
!----------------------------------------------------------
!Calculate updated residual "res_vec=b-A*x" relative to current x
DO i=1,n
Ax(i)=CMPLX(0.0d0, 0.0d0)
DO j=1,n
Ax(i)=Ax(i)+A(i,j)*x(j)
END DO
END DO
DO i=1,n
res_vec(i)=b(i)-Ax(i)
END DO
DO i=1,n
rel_res=rel_res+(res_vec(i)*CONJG(res_vec(i)))
END DO
res=SQRT(res)/REAL(n2b)
WRITE(*,*) "res=",res
WRITE(*,*) " "
END DO main
!----------------------------------------------------------
!Output message
IF (k<maxit) THEN
WRITE(*,*) "Converged in",k,"iterations"
ELSE
WRITE(*,*) "STOPPED after",k, "iterations because max no. of iterations was reached"
END IF
!Output solution vector:
WRITE(*,*) "x_sol="
DO i=1,n
WRITE(*,*) x(i)
END DO
!----------------------------------------------------------
DEALLOCATE(x0,x0s, Ax, ATx, x, xs, p, ps ,r, rs, Ap, Aps, res_vec)
DEALLOCATE(A,b)
!
END PROGRAM
!
!////////////////////////////////////////////////////////////////////////
结果:我的脚本的结果如下:
STOPPED after 50 iterations because max no. of iterations was reached
x_sol=
(-2.88435711452590705E-002,-0.43229898544084933 )
( 0.11755325208241280 , 0.73895038053993978 )
而实际结果是使用 MATLAB 内置的 bicg.m 函数给出的:
-0.3700 - 0.6702i
0.7295 + 1.1571i
最佳答案
以下是您的程序中的一些缺陷。它们是否是错误有点主观，是否修改代码完全取决于您。
在这一行
IF (n2r0==0) THEN您测试(可能长时间运行)循环的结果是否相加 精确到 0。这对于浮点来说总是一个坏主意 数字。如果你不知道这一点，请回顾很多很多问题 这里带有标签
floating-point，它源于广泛的 对合理期望的理解不准确 f-p 算术。我认为您在比较的左侧使用实数并在右侧使用整数不会使情况变得更糟，但它不会使它们变得更好。在代码中至少有两个地方可以计算矩阵向量乘积。您可以用对内在
matmul例程的调用来替换这些循环(我想，我没有像您那样仔细检查您的代码)。这实际上可能会减慢您的代码速度，但这不是现阶段的问题。调用经过良好测试的库例程而不是自行编写将 (a) 减少您必须维护/测试/修复的代码量，并且 (b) 更有可能提供一次正确的解决方案。一旦代码可以运行，如果必须的话，请担心性能。您声明了许多具有 double 的实数和复数变量，但是 使用以下语句初始化它们:
A(1,1)=CMPLX(-0.73492,7.11486)double 变量大约有 15 位小数位可用， 但在这里您只提供前 6 个值。你 不能依赖编译器将其他数字设置为任何 特定的值(value)观。相反，像这样初始化:
A(1,1)=CMPLX(-0.73492_dp,7.11486_dp)这将导致这些值被初始化为 double 最接近
-0.73492和7.11486的精度数字。当然，您之前必须编写过类似dp = kind(0d0)的内容，并且还有其他方法可以强制文字常量的精度，但这是我通常采用的方法。如果您有一个提供内在iso_fortran_env模块的现代 Fortran 编译器，您可以将_dp替换为现在标准的_real64。这段代码
x0(:)=CMPLX(0.0d0,0.0d0) DO i=1,n x0s(i)=CONJG(x0(i)) END DO可以替换为
x0 = CMPLX(0.0d0,0.0d0) x0s = x0使用数组语法将第一个归零似乎有点奇怪 数组，然后循环将第二个归零；看起来更加奇特 当
CONJG(0,0)==(0,0)时重复调用CONJG。这段代码
DO i=1,n n2b=n2b+(b(i)*CONJG(b(i))) res=res+(r(i)*CONJG(r(i))) !== inner prod(r,r) END DO n2b=SQRT(n2b) res=SQRT(res)/n2b如果我理解正确的话，可以替换为
n2b = sqrt(dot_product(b,b)) res = sqrt(dot_product(r,r))/n2b我实际上没有发现您的代码有任何问题，但是使用内在函数可以减少您需要编写和维护的行数，就像上面的
matmul的情况一样。
可能还有其他不太明显的瑕疵，但这些应该可以帮助您开始。
关于fortran90 - Fortran 中的 BiCG 算法无法正常工作？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/15065908/