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

Go to the documentation of this file.

subroutine bicgstab_no_s(A,b,x,n,tol,mainloop,k)
implicit none
integer n,i,j,k,mainloop
real(8) A(n,n),b(n),x(n),res(n)
real(8) rhat(n),p(n),v(n),t(n) !,s(n)
real(8) alpha,beta,omega,rho,rhonew,norm_value,aux,res_one
real(8) tol
integer,parameter :: display=100
! "mainloop" is necessary to stop the iteration loop is it does not converge after a while
! "display" is the variable which determines the on-screen data print interval...
!  ... meaning that if "display=1" it shows every residual norm in every iteration step...
!  ... or if "display=100" it shows the residual norm in every 100 iteration step...

! this interface is necessary to use vector values functions
! mfmultiplier multiplies the A matrix with any given vector is such a way that
! the A is stored in five different vectors namely;mfD,mfW,mfE,mfS,mfN
! "mf" for matrix free and D,W,E,S,N for Diagonal,West,East,South and North respectively.
!
! so you don't have the use "call", you can use this routine inline.
! check the code to understand how it is used.


!initialize the parameters
v=0.d0
p=0.d0
t=0.d0
!s=0.d0
res=0.d0
rhat=0.d0



!-----------------
!res=b-matmul(A,x)

call MatVecMult(A,x,res,n)
!res=b-MatVecMultFunc(A,x) !needs interface
res=b-res

!-----------------



! normally rhat must be selected in such a way that their dot product is not zero
! generally it is selected as the residual(res) itself...
! an idea may be using the normalized residual(res) as rhat but I don't know if 
! it will came out useful or not
rhat=res

! bicgstab variables
rho=1.d0
alpha=1.d0
omega=1.d0




Open (111,File='linearResidual.dat',Status='unknown')
norm_value=dsqrt(dot_product(res,res)) ! first norm
write(111,'(E30.20,1X)')norm_value

res_one=norm_value

do i=1,mainloop
        
        !print*,i


        !rhonew=dotter(res,rhat,nh*nv)
        rhonew=dot_product(res,rhat)


        if (rhonew.eq.0) then
                print*,'Method fails: rho=0'
                exit
        end if

        beta=rhonew/rho*alpha/omega



        p=res+beta*(p-omega*v)

        

        !v=matmul(A,p)
        call MatVecMult(A,p,v,n)
        !v=backwardsub(U,forwardsub(L,matmul(A,p)))
        !v=mfmultiplier(mfD,mfW,mfE,mfS,mfN,p,nh,nv)

        !alpha=rhonew/dotter(v,rhat,nh*nv)
        alpha=rhonew/dot_product(v,rhat)

        !s=res-alpha*v
        !print*,'----->',res !-alpha*v

        res=res-alpha*v

        ! res=s: s is zero - converged solution (happy breakdown?)
        if (dot_product(res,res).le.epsilon(0.d0)) then

                goto 100
                

        end if

        
        !t=matmul(A,s)
        !t=mfmultiplier(mfD,mfW,mfE,mfS,mfN,s,nh,nv)
        !t=matmul(A,res)
        call MatVecMult(A,res,t,n)
        !t=backwardsub(U,forwardsub(L,matmul(A,res)))
        !t=mfmultiplier(mfD,mfW,mfE,mfS,mfN,res,nh,nv)


        !omega=dotter(t,s,nh*nv)/dotter(t,t,nh*nv)
        !omega=dotter(t,res,nh*nv)/dotter(t,t,nh*nv)
        omega=dot_product(t,res)/dot_product(t,t)

        if (omega.eq.0) then
                print*,'Method cannot continue:w=0'
                exit
        end if

        !norm_value=norm2(x,x+alpha*p+omega*s,n)

        !write(111,'(F20.12,1X)')norm_value

        !if (norm_value.le.tol) then
        !       x=x+alpha*p+omega*s
        !       print*,'Solver converged'
        !       print*,'i=',i ,'norm=',norm_value
        !       exit
        !else

        !x=x+alpha*p+omega*s
100     x=x+alpha*p+omega*res


        !end if

        
        !res=s-omega*t
        res=res-omega*t

        !print*,res


        rho=rhonew


        norm_value=dsqrt(dot_product(res,res))

        !print*,norm_value

        write(111,'(E30.20,1X)')norm_value

        if (norm_value/res_one.le.tol) then
                print*,'Solver converged'
                print*,'i=',i ,'norm=',norm_value
                return
        end if

        if (mod(i,display).eq.0) then
                print*,'i=',i ,'relative res=',norm_value/res_one
        end if

end do


close(111)


end subroutine

---