E:/turkuaz/krylov/codes/newton/1.2b/linear/krylov/bicgstab.f90

Go to the documentation of this file.

subroutine bicgstab(n,A,b,x,tol,maxIter,k)
!CALL bicgstab(N,JAC,-F0,deltax,linearTol,maxLinearIter,i) !<--- call within the Newton's Method

implicit none

integer n,maxIter

real(8) :: A(n,n),b(n),x(n)

real(8) :: tol


! bicgstab related definitions

real(8) r(n) ! residual

real(8) rhat(n) ! residual hat

real(8) rhoZero,alpha,omegaZero   ! beware: alpha, omega etc should not conflict with other alphas,omegas etc

real(8) v(n),p(n),s(n),t(n)

real(8) rhoOne

real(8) beta


integer decision

integer i,k


r=0.d0

r=b-matmul(A,x)

rhat=r

if ( dot_product(r,rhat).le.epsilon(0.d0) ) then

        print*,'r and rhat are orthogonal or r is zero'
        !print*,'0 to stop else to select random rhat'
        !read*,decision

        !if (decision.eq.0) then

                print*,'kNaN stop!'
                stop
        !else
                !rhat=1.d0  ! should be random or different than 1?
        !end if



end if


rhoZero=1.d0
alpha=1.d0
omegaZero=1.d0

v=0.d0
p=0.d0
s=0.d0




! as suggested by van der vorst, in his bicgstab paper, if s is small, then we must exit the solver with 
! as done in the ||s||2 part of the code. however, for small sized problems, s might be zero
! in that case, t=0 and we might find that omegaZeo=0/0 so run time error
! to avoid this case, another check is put after the calculation of s, if the norm is less than the
! machine epsilon (i.e. if the norm is exactly zero then it is smaller than mach_eps) then the algorithm 
! jumps to x update section - since t will be zero, t will not contribute to x - t can be totally avoided.
! since s is smaller than mach_eps, it will automatically be smaller than the tolerance of the model
! so return will be invoked - if however residual is used cor convergence rather than the s itself
! then a return could be added since we can get no gain in the foregoing computations
! but it might be interesting to check whether the code will exit even if even if the return is not
! explicity invoked in the residual_correction part (which is absent in this code)
! but still - pre_bicgstab_no_s will be used in kNaN - al least for testing purposes





do i=1,maxIter

!if (k.eq.6) print*,'---------------------'



rhoOne=dot_product(rhat,r)

!if (k.eq.6) print*,'rhoOne=',rhoOne

beta=(rhoOne/rhoZero)*(alpha/omegaZero)

!if (k.eq.6) print*,'beta=',beta

p=r+beta*(p-omegaZero*v)



!if (k.eq.6) print*,'p=',p



v=matmul(A,p)

!if (k.eq.6) print*,'v=',v

alpha=rhoOne/dot_product(rhat,v)

!if (k.eq.6) print*,'alpha=',alpha


s=r-alpha*v

if (dot_product(s,s).le.epsilon(0.d0)) then

        goto 1000

end if


!if (k.eq.6) print*,'s=',s

t=matmul(A,s)

!if (k.eq.6) print*,'t=',t

omegaZero=dot_product(t,s)/dot_product(t,t)

!if (k.eq.6) print*,'omegaZero=',omegaZero

1000 x=x+alpha*p+omegaZero*s

!if (k.eq.6) print*,'x=',x

if (dot_product(s,s).le.tol) then

        print*,'bicgstab converged in step=',i
        return

end if



r=s-omegaZero*t

!if (k.eq.6) print*,'r=',r

rhoZero=rhoOne



end do








end subroutine bicgstab

---