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

Go to the documentation of this file.

! Version 1.2a
! March 30, 2009

subroutine newton(X,n)
use solverParameters
use numerical_libraries

use myFunctions

implicit none
integer i,j !< dummy iteration variable
integer n !< size of the unknown vector
integer dummyBacktrackCount
integer NumberOfCalculatedNorms

real(8)  X(n),deltaX(n),F0(n),F1(n),F2(n) !< definitions of the vectors
real(8) fNorm(maxNewtonIter),lambdaVec(maxNewtonIter)
real(8) JAC(n,n)
real(8) initialNorm
real(8) convergenceTolerance
real(8) fZero,fOne,fTwo
real(8) alpha,mybeta,lambda,lambda1,lambda2,lambdaTest
real(8) slope

! inexact newton settings

real(8) inexactLinearTol

real(8) eta




alpha=1e-4
mybeta=0.9d0
!lambda=1;
eta=0.5d0


newtonSuccess=.false.


! calculate the initial F
call F(X,F0,n)
! calculate the initial L2-norm of F
initialNorm=dot_product(F0,F0)
fnorm(1)=initialNorm
fZero=0.5d0*fnorm(1)   ! f0=0.5*F0'*F0;

i=0

if (fnorm(i+1).le.convergenceTolerance) then
                        
                print*,'---------------------------------------------------------------------------'
                print*,'newton converged @ iter=',i
                print*,'current norm is ',fnorm(i+1)
                print*,'convergence tolerance is',convergenceTolerance
                print*,''
                print*,'-----> already converged solution is given!'
                print*,''

                if (enableNewtonDebug.EQV.ON) print*,'convergence option is', NewtonConvergenceOption
                if (enableNewtonDebugExtra.EQV.ON) print*,'absolute residual=',fnorm(i+1)
                if (enableNewtonDebugExtra.EQV.ON) print*,'relative residual=',fnorm(i+1)/initialNorm

                print*,'---------------------------------------------------------------------------'
                newtonSuccess=.true.
                NumberOfCalculatedNorms=i+1
                return                          

end if


! decide on the convergence tolerance as a convergence criteria
select case(NewtonConvergenceOption)

        case(1)

                convergenceTolerance=nonlinearTolAbs

        case(2)

                convergenceTolerance=nonlinearTolRel*initialNorm

        case(3)

                convergenceTolerance=nonlinearTolRel*initialNorm+nonlinearTolAbs

        case default

                print*,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                print*,'non valid convergence option'
                print*,'!!!!!!!!!!kNaN STOP!!!!!!!!!!'
                stop

end select



!x(1)=0.998162047066876
!x(2)=1.00372230866380



do i=1,maxNewtonIter

        !print*,x
        !print*,fnorm(i)

        
    if (enableNewtonBasic.EQV.ON) print*,'Newton iter=',i

        ! initialize the update 
    deltaX=0.d0

        call JF(X,JAC,n)  !  JAC=JF(X);
        
        !--------------------

        !fZero=0.5d0*dot_product(F0,F0)   ! f0=0.5*F0'*F0;
        
    !fnorm(i)=2.d0*fZero


        ! linear solver will be added here
    ! deltaX=JAC\(-F0);
        !CALL DLSLRG (N, JAC, N, -F0, 1, deltaX)

        inexactLinearTol= min(fnorm(i)*eta,maxLinearTol)        ! implement a method to select eta
        inexactLinearTol= max(inexactLinearTol,LinearTol)   ! this line is needed else linear solver fails if nonlinearTol>linearTol
                                                                                                                ! solver cannot satisfy a solution accurate than the nonlinear accuracy 
                                                                                                                ! ( verify! first comment out this line and work it out!)


        ! krylov solvers
        call bicgstab_no_s(JAC,-F0,deltax,n,inexactLinearTol,maxLinearIter,i)
        !CALL bicgstab(N,JAC,-F0,deltax,linearTol,maxLinearIter,i)


    !  check alpha condition

    lambda=1.d0;
    lambda2=lambda;

    call F(X+lambda*deltaX,F2,n)  !F2=F(X+lambda*deltaX);
        
        
    fTwo=0.5d0*dot_product(F2,F2)

    slope=-dot_product(F0,F0)    ! is it still the slope when solve iteratively?????????
    ! or gradient of f times deltaX

    !slope=GF(X)'*deltaX;

    if (fTwo .LT. fZero + alpha*lambda*slope) then

        X=X+lambda*deltaX
        lambdaVec(i)=lambda

                if (enableNewtonDebug.EQV.ON) print*,'valid alpha condition (direct) in iter',i


                !#######################################################################################
                !# convergence and initialization for the next iteration set starts here
                ! reduce F2 and F0 by overwriting and maybe with pointers!

                F0=F2


                if ((i+1).GT.maxNewtonIter) then
                        print*,'!!!----------------------------------------------------------------------'
                        print*,'Newtons method failed to converge in ',maxNewtonIter,'iterations'
                        print*,'current norm is ',fnorm(i)
                        print*,'convergence tolerance is',convergenceTolerance
                        
                        if (enableNewtonDebug.EQV.ON) print*,'convergence option is', NewtonConvergenceOption
                        if (enableNewtonDebugExtra.EQV.ON) print*,'absolute residual=',fnorm(i+1)
                        if (enableNewtonDebugExtra.EQV.ON) print*,'relative residual=',fnorm(i+1)/initialNorm

                        print*,'!!!------------------------------------------------------------------------'

                        NumberOfCalculatedNorms=i

                        exit
                end if


                fnorm(i+1)=fTwo*2.d0


                if (fnorm(i+1).le.convergenceTolerance) then
                        
                        print*,'---------------------------------------------------------------------------'
                        print*,'newton converged @ iter=',i
                        print*,'current norm is ',fnorm(i+1)
                        print*,'convergence tolerance is',convergenceTolerance

                        if (enableNewtonDebug.EQV.ON) print*,'convergence option is', NewtonConvergenceOption
                        if (enableNewtonDebugExtra.EQV.ON) print*,'absolute residual=',fnorm(i+1)
                        if (enableNewtonDebugExtra.EQV.ON) print*,'relative residual=',fnorm(i+1)/initialNorm

                        print*,'---------------------------------------------------------------------------'
                        newtonSuccess=.true.
                        NumberOfCalculatedNorms=i+1
                        exit                            

                end if

                !fZero=0.5d0*fnorm(i+1)   ! f0=0.5*F0'*F0;
                fZero=fTwo    ! reduced by comparing Newton routine in version 1.2

                if (enableNewtonBasic.EQV.ON) print*,'  |----> current norm=',fnorm(i+1)

                

                !# convergence and initialization for the next iteration set ends here
                !#######################################################################################

        cycle

    else  ! start line search backtracking

                dummyBacktrackCount=1

        ! first try quadratic

        lambda1=quadratic(slope,fZero,fTwo)


        if (lambda1.lt.0.1d0*lambda) then
            
            lambda1=0.1d0*lambda
            
        elseif (lambda1.gt.0.5d0*lambda) then
            
            lambda1=0.5d0*lambda;
            
        end if


        ! check alpha condition

        call F(X+lambda1*deltaX,F1,n) !F1=F(X+lambda1*deltaX);
        
        fOne=0.5d0*dot_product(F1,F1)


        if (fOne .lt. fZero + alpha*lambda*slope) then
            
            lambda=lambda1
            X=X+lambda*deltaX
            lambdaVec(i)=lambda
                        
                        if (enableNewtonDebug.EQV.ON) print*,'valid alpha condition (quadratic) in iter',i


                        !#######################################################################################
                        !# convergence and initialization for the next iteration set starts here
                        ! reduce F2 and F0 by overwriting and maybe with pointers!

                        F0=F1


                        if ((i+1).GT.maxNewtonIter) then
                                print*,'!!!----------------------------------------------------------------------'
                                print*,'Newtons method failed to converge in ',maxNewtonIter,'iterations'
                                print*,'current norm is ',fnorm(i)
                                print*,'convergence tolerance is',convergenceTolerance

                                if (enableNewtonDebug.EQV.ON) print*,'convergence option is', NewtonConvergenceOption
                                if (enableNewtonDebugExtra.EQV.ON) print*,'absolute residual=',fnorm(i+1)
                                if (enableNewtonDebugExtra.EQV.ON) print*,'relative residual=',fnorm(i+1)/initialNorm

                                print*,'!!!------------------------------------------------------------------------'
                                NumberOfCalculatedNorms=i

                                exit
                        end if


                        fnorm(i+1)=fOne*2.d0


                        if (fnorm(i+1).le.convergenceTolerance) then
                                
                                print*,'---------------------------------------------------------------------------'
                                print*,'newton converged @ iter=',i
                                print*,'current norm is ',fnorm(i+1)
                                print*,'convergence tolerance is',convergenceTolerance

                                if (enableNewtonDebug.EQV.ON) print*,'convergence option is', NewtonConvergenceOption
                                if (enableNewtonDebugExtra.EQV.ON) print*,'absolute residual=',fnorm(i+1)
                                if (enableNewtonDebugExtra.EQV.ON) print*,'relative residual=',fnorm(i+1)/initialNorm

                                print*,'---------------------------------------------------------------------------'
                                NumberOfCalculatedNorms=i+1
                                newtonSuccess=.true.
                                exit                            

                        end if



                        !fZero=0.5d0*fnorm(i+1)   ! f0=0.5*F0'*F0;
                        fZero=fOne    ! reduced by comparing Newton routine in version 1.2

                        if (enableNewtonBasic.EQV.ON) print*,'  |----> current norm=',fnorm(i+1)

                        !# convergence and initialization for the next iteration set ends here
                        !#######################################################################################

            
                        cycle

        else
            

            do while (fOne .gt. fZero + alpha*lambda*slope)
                
                if (enableNewtonDebugExtra.EQV.ON) print*,'cubic backtrack number ', dummyBacktrackCount
                

                lambdaTest=cubic(slope,fZero,fOne,fTwo,lambda1,lambda2)
                
                if (lambdaTest.lt.0.1d0*lambda1) then

                    lambdaTest=0.1d0*lambda1

                elseif (lambdaTest.gt.0.5d0*lambda1) then

                    lambdaTest=0.5d0*lambda1

                end if
                
                fTwo=fOne;
                lambda2=lambda1
                
                lambda1=lambdaTest



                                call F(X+lambda1*deltaX,F1,n)

                fOne=0.5d0*dot_product(F1,F1)




                                dummyBacktrackCount=dummyBacktrackCount+1


                                if (dummyBacktrackCount.GT.maxBacktrackCount) then

                                        if (enableNewtonDebugExtra.EQV.ON) print*,'too many backtracks=',maxBacktrackCount
                                        exit

                                end if




            end do


            lambda=lambda1
            
            X=X+lambda*deltaX
            lambdaVec(i)=lambda

                        if (enableNewtonDebug.EQV.ON) print*,'valid alpha condition (cubic) in iter',i


                        !#######################################################################################
                        !# convergence and initialization for the next iteration set starts here
                        ! reduce F2 and F0 by overwriting and maybe with pointers!

                        F0=F1


                        if ((i+1).GT.maxNewtonIter) then
                                print*,'!!!----------------------------------------------------------------------'
                                print*,'Newtons method failed to converge in ',maxNewtonIter,'iterations'
                                print*,'current norm is ',fnorm(i+1)
                                print*,'convergence tolerance is',convergenceTolerance

                                if (enableNewtonDebug.EQV.ON) print*,'convergence option is', NewtonConvergenceOption
                                if (enableNewtonDebugExtra.EQV.ON) print*,'absolute residual=',fnorm(i+1)
                                if (enableNewtonDebugExtra.EQV.ON) print*,'relative residual=',fnorm(i+1)/initialNorm

                                print*,'!!!------------------------------------------------------------------------'
                                NumberOfCalculatedNorms=i
                                exit
                        end if


                        fnorm(i+1)=fOne*2.d0


                        if (fnorm(i+1).le.convergenceTolerance) then
                                
                                print*,'---------------------------------------------------------------------------'
                                print*,'newton converged @ iter=',i
                                print*,'current norm is ',fnorm(i+1)
                                print*,'convergence tolerance is',convergenceTolerance

                                if (enableNewtonDebug.EQV.ON) print*,'convergence option is', NewtonConvergenceOption
                                if (enableNewtonDebugExtra.EQV.ON) print*,'absolute residual=',fnorm(i+1)
                                if (enableNewtonDebugExtra.EQV.ON) print*,'relative residual=',fnorm(i+1)/initialNorm

                                print*,'---------------------------------------------------------------------------'
                                NumberOfCalculatedNorms=i+1
                                newtonSuccess=.true.
                                exit                            

                        end if

                        !fZero=0.5d0*fnorm(i+1)   ! f0=0.5*F0'*F0;
                        fZero=fOne    ! reduced by comparing Newton routine in version 1.2

                        if (enableNewtonBasic.EQV.ON) print*,'  |----> current norm=',fnorm(i+1)

                        !# convergence and initialization for the next iteration set ends here
                        !#######################################################################################


        end if


    end if





end do



!NumberOfCalculatedNorms=7



!%%%% POST PROCESSING START
!simple post processing with data output ! no alteration
910 FORMAT(I4,2(E18.10,1X))


select case(NewtonConvergenceOption)

case(1)


        Open (1,File='NewtonResidual.dat',Status='unknown')
        WRITE(1,*) 'TITLE="results"'
        WRITE(1,*) 'VARIABLES="iter","||f||"'
        WRITE(1,*) 'ZONE I=', NumberOfCalculatedNorms
        DO j= 1, NumberOfCalculatedNorms
                        WRITE(1,910) j, fnorm(j)
        END DO
        close(1)


case(2)

        Open (1,File='NewtonResidual.dat',Status='unknown')
        WRITE(1,*) 'TITLE="results"'
        WRITE(1,*) 'VARIABLES="iter","||f||"'
        WRITE(1,*) 'ZONE I=', NumberOfCalculatedNorms
        DO j= 1, NumberOfCalculatedNorms
                        WRITE(1,910) j, fnorm(j)/fnorm(1)
        END DO
        close(1)

end select
!%%%% POST PROCESSING END

end subroutine newton

---