D:/demona/2.0/manual/doc/2-31set/processor/krylov/krylov.f90

Go to the documentation of this file.

module krylov


use mF

use myFunctions

use demonaParameters

implicit none


contains

subroutine gmres(b,myX,m,restart,tolerance,residualVector,val,colInd,rowPtr,diag,valsize,totalit,preconditioner)

implicit none
integer i,j,ii,jj,dir(9)
integer m,restart,RestartCounter
integer q,rIndex

integer totalit
real(8) beta,residual,aux,BetaZero,oldresidual
!real(8),parameter :: tol=1e-6
real(8) tolerance
! relative residual check, use beta for ||r0||
!static arrays
real(8) POINTER, myX(:)
real(8) w(size(myX)), b(size(myX))  
real(8) f(size(myX)) !deltax(size(x))
!preconditioner
real(8) p(size(myX)) !,zero(size(x))

real(8) residualVector(size(myX))  ! later avoid this and use one of the vectors on the algorithm

!real(8) tempx(size(x))

!dynamic arrays
real(8), POINTER, DIMENSION(:) :: y !max:y(m)
real(8), POINTER, DIMENSION(:) :: H !max:H(m+1)
real(8), POINTER, DIMENSION(:) :: sint !max:sint(m+1)
real(8), POINTER, DIMENSION(:) :: cost !max:cost(m+1)
real(8), POINTER, DIMENSION(:) :: gamma !max:gamma(m+1)
real(8), POINTER, DIMENSION(:,:) :: V !max:V(nx*ny,m+1)
real(8), POINTER, DIMENSION(:,:) :: R !max:R(m+1,m)
logical NP,LP,RP
real(8) res_test(size(myX)),progressCounter,temp1,temp2
integer happy


! crs jacobian definitions

integer valSize
real(8) val(valSize)
integer colInd(valSize)
integer rowPtr(size(myX)+1)
integer diag(size(myX))

integer preconditioner


if (gsSolver.eqv.no) then

print*,' '
print*,'###################'
print*,'GMRES'
print*,'###################'

end if



RestartCounter=0
Open (1,File='GmresResidual.dat',Status='unknown')




!preconditioner=0  no preconditioner
!preconditioner=1  left preconditioner
!preconditioner=2  right preconditioner







! initialize

totalit=0

10 w=0.d0


allocate(V(size(myX),1),H(2),R(2,1),sint(1),cost(1),gamma(1))

V=0.d0
H=0.d0
R=0.d0
sint=0.d0
cost=1.d0

p=residualVector


20 continue



!v(:,1)=b-matVecMult(myX)
v(:,1)=b-matVec(myX)



if (preconditioner.eq.1) then

V(:,1)=LUyX(val,colInd,rowPtr,size(myX),diag,valSize,V(:,1))

end if



V(:,1)=invscaler(V(:,1),p)  !diagonal scaling


beta=dsqrt(dot_product(v(:,1),v(:,1)))



v(:,1)=v(:,1)/beta   ! v1=ro/beta



gamma(1)=beta ! note that beta="res_zero"
!residual=beta

if (RestartCounter.eq.0) then
        residual=beta
        BetaZero=beta
        !write(1,'(F20.12,1X)')residual/BetaZero
        write(1,'(E20.12,1X)')residual/1.d0
end if


! Arnoldi Orthogonalization Procedure

do j=1,m
        
        totalit=totalit+1


        if (preconditioner.eq.0) then

                w=matVec(V(:,j))
        
        elseif (preconditioner.eq.1) then

                w=LUyX(val,colInd,rowPtr,size(myX),diag,valSize,matVec(V(:,j)))

        elseif (preconditioner.eq.2) then

                w=matVec(LUyX(val,colInd,rowPtr,size(myX),diag,valSize,V(:,j)))

        end if



        w=invscaler(w,p) !diagonal scaling

        do i=1,j   ! reduction openmp
                H(i)=dot_product(w,v(:,i))
                w=w-H(i)*v(:,i)
        end do
        ! normalization


        H(j+1)=dsqrt(dot_product(w,w))  !??? if için belki bu satýr gerekmeyebilir, kontrol et

        ! happy breakdown was here ----- 1

        R=>reallocate_matrix(R,j+1,j)

        !apply previous  transformations
        R(:,j)=previous(H,sint,cost) ! <--- and add newly modified H to R

        if (j.gt.1) R(j+1,j-1)=0.d0  ! infact that value is not needed during the backward substitution

        !calculate last angle
        sint(j)=R(j+1,j)/dsqrt(R(j+1,j)**2+R(j,j)**2)
        cost(j)=R(j,j)/dsqrt(R(j+1,j)**2+R(j,j)**2)

        !if (H(j+1).eq.0.d0) then
        if (dabs(H(j+1)).le.epsilon(0.d0)) then
                if (gsSolver.eqv.no) print*,'                   |Happy Breakdown!'
                happy=1
                exit  
        end if
!        happy breakdown: zero residual i.e. exact solution
!        do not calculate the last transformation but apply previous transformations

        ! apply last transformation
        aux=R(j,j)
        R(j,j)=cost(j)*R(j,j)+sint(j)*R(j+1,j)  ! burada R(j,j) kendisini update ettiði için sorun oluyor.

        R(j+1,j)=-sint(j)*aux+cost(j)*R(j+1,j)

    gamma=>reallocate_vector(gamma,j+1)
        gamma(j+1)=0.d0
        aux=gamma(j)
        gamma(j)=cost(j)*gamma(j)+sint(j)*gamma(j+1)   ! gamma(j+1)'nýn deðeri olmadýðý için 1e+66 tarzý þeyler ortaya çýktý...
        gamma(j+1)=-sint(j)*aux+cost(j)*gamma(j+1)

        residual=-sint(j)*residual  !in fact residual =gamma(j+1)
        
        !print*,abs(residual)

        !write(1,'(F20.12,1X)')dabs(residual/BetaZero)
        write(1,'(E20.12,1X)')dabs(residual/1.d0)

        !if (((dabs(residual/BetaZero)).LE.tolerance).OR.(j.EQ.m)) exit
        if (((dabs(residual/1.d0)).LE.tolerance).OR.(j.EQ.m)) exit

        H=>reallocate_vector(H,j+2)
        sint=>reallocate_vector(sint,j+1)
        cost=>reallocate_vector(cost,j+1)
        V=>reallocate_matrix(V,size(myX),j+1)

        v(:,j+1)=w/H(j+1) 

end do

if (gsSolver.eqv.no) print*,'                   |gmres exit @ j=',j,'          |'
!print*,'                       |residual is=',dabs(residual/BetaZero),' |'
if (gsSolver.eqv.no) print*,'                   |residual is=',dabs(residual/1.d0),' |'


allocate(y(size(R,DIM=2)))

y=0.d0
y=backwardsub(R(1:size(R,DIM=2),:),gamma(1:size(R,DIM=2)))


if (preconditioner.ne.2) then

        myX=myX+matmul(V(:,1:size(R,DIM=2)),y)

else

        myX=myX+LUyX(val,colInd,rowPtr,size(myX),diag,valSize,matmul(V(:,1:size(R,DIM=2)),y))

end if



deallocate(y,H,cost,gamma,V,R)




deallocate(sint)


if (happy.eq.1) goto 101

!if ((abs(residual/BetaZero)).GT.tolerance) then
if ((dabs(residual/1.d0)).GT.tolerance) then
        RestartCounter=RestartCounter+1
        if (RestartCounter.le.restart) then
                if (gsSolver.eqv.no) print*,'                   |restart ',RestartCounter
                goto 10
        end if
end if

101  continue


close(1)


!residualVector=0.d0
!residualVector=b-matVec(myX)

!print*,dsqrt(dot_product(residualVector,residualVector))
!stop

end subroutine gmres



subroutine bicgstabRP(b,myx,tol,maxIter,residualVector,val,colInd,rowPtr,diag,valsize,totalit)
! later bicgstab and gmres will be used a krylovSolvers as an module and
! n as size of deltax will not be needed.

!use matrixFree

implicit none

integer i,j

integer n

integer maxIter

real(8) myx(:),b(:)


real(8) tol


!real(8) f(n)  !<--- use this either with or pointer or whatever to avoid the numerous calls in matVecMult

real(8) r(size(b)),p(size(b)),rHat(size(b)),s(size(b)),v(size(b)),t(size(b)),shat(size(b)),phat(size(b))

real(8) residualVector(size(b))  ! later avoid this and use one of the vectors on the algorithm


real(8) alphaB,omegaB,betaB,rhoNewB,rhoOldB,rRhat


real(8) normR,normS,normIni

real(8) rho,rhonew

real(8) res_one,norm_value,bNorm

integer totalit

! crs jacobian definitions

integer valSize
real(8) val(valSize)
integer colInd(valSize)
integer rowPtr(size(myX)+1)
integer diag(size(myX))


print*,' '
print*,'###################'
print*,'BICGSTAB'
print*,'###################'







n=size(b)

totalit=0



912 FORMAT(1(F30.20,1X))

!Open (1,File='pre.dat',Status='unknown')
!read(1,912)pre
!close(1)



!r=0.d0

r=b-matvec(myX)




!LUyX(val,colInd,rowPtr,size(myX),diag,valSize,V(:,j))


normIni=dsqrt(dot_product(r,r))


!print*,r






!print*,normIni

!stop

!normIni=1.d0


!pause

!stop


bNorm=dsqrt(dot_product(b,b))


if (bNorm.eq.0.d0) bNorm=1.d0

!print*,bNorm

!stop


rHat=r


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

rhoOldB=1.d0
alphaB=1.d0
omegaB=1.d0

rho=1.d0





Open (111,File='bicgstab.dat',Status='unknown')

write(111,'(F20.12,1X)')normIni


if (normIni.le.tol) then
        
        print*,'no iter convergence'

        goto 100

end if


alphaB=0.d0
betaB=0.d0
omegaB=1.d0
rhoNewB=0.d0
rhoOldB=0.d0


do i=1,maxIter

        
totalit=totalit+1

        !print*,i


        rhoNewB = dot_product(rHat,r)

        !print*,rhoNewB

        !if (dabs(rhoNewB).le.epsilon(0.d0)) then
        if (rhoNewB.eq.0.d0) then
                print*,rhoNewB
                print*,'Method fails: rho=0'
                exit
        end if


        if (i.eq.1) then

                p=r


        else

                betaB = (rhoNewB / rhoOldB) * (alphaB / omegaB)

        
                p = r + betaB * (p - omegaB * v)

        end if

        !phat=LUyX(val,colInd,rowPtr,size(myX),diag,valSize,p)
        phat=p

        v=matvec(phat)


        alphaB = rhoNewB / dot_product(rHat,v)


        !print*,alphaB
        !stop

        s = r-alphaB * v

        if (i.eq.20) then

                !print*,alphaB
                !pause

        end if


        normS=dsqrt(dot_product(s,s))

        if (normS.le.tol)  then
                myX = myX + alphaB * phat
                write(111,'(F20.12,1X)')normS
                print*,'s exit'
                exit
        end if


        !shat=LUyX(val,colInd,rowPtr,size(myX),diag,valSize,s)
        shat=s

        t=matvec(shat)

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

        !print*,omegaB

        !if ((dabs(omegaB).le.epsilon(0.d0)).or.(dot_product(t,t).le.epsilon(0.d0))) then


        myX = myX + alphaB * phat +omegaB * shat



        r = s - omegaB * t







        normR=dsqrt(dot_product(r,r))
        !normR=dsqrt(dot_product(invscaler(r,pre),invscaler(r,pre)))
        !normR=dsqrt(dot_product(b-matVecMult(myX),b-matVecMult(myX)))

        write(111,'(F20.12,1X)')normR

        if (mod(i,100).eq.0) print*,normR
        
        !if ((normR*normIni).le.tol)  then
        if ((normR).le.tol)  then

                print*,'bicgstab converged in iter=',i
                print*,'abs norm=',normR
                !print*,'rel norm=',normR/bNorm

                exit

        end if


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



        rhoOldB = rhoNewB 




end do


100 continue

close(111)


!print*,myX
!stop


residualVector=0.d0
residualVector=b-matVec(myX)


!print*,dsqrt(dot_product(residualVector,residualVector))
!stop



end subroutine bicgstabRP

























end module krylov

---