私は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