module newton_mod

   use dual_problem_mod
   use dwr_mod
   use element_mod
   use elemental_mod
   use estimates
   use euler_problem
   use matrix_oper
   use matrix_oper_int
   use main_data
   use mesh_mod
   use paramets
   use solution_mod
   use terms_mod

implicit none

   public :: EstimateNewtonADwr
   public :: EstimateNewtonADwrFixed
   public :: EstimateNewtonARes
   public :: localNewtonEstimate
   public :: Newton_output
   public :: NewtonSolve
   public :: performOneNewtonStdgmStep
   public :: setAndersonInNewton
   public :: setDampingInNewton
   !public :: setDampingInNewtonLocal
   public :: setDampingInNewton_DWR
   public :: setLocalDampingDWR


contains


   ! seeking of optimal damping parameter in the Newton method,
   ! this should help the convergence of Newton when the approximation is still rough
   subroutine setDampingInNewton( newton, grid, deg_plus, eta, loc_implicitly, &
                                  res_max_val, max_l, l, ratio, DWR)
      class( Newton_t ), intent ( inout ) :: newton
      class( mesh ), intent(inout) :: grid
      logical, intent(in) :: deg_plus
      real, intent(in) :: eta
      logical, intent(in) :: loc_implicitly
      real, intent(in) :: res_max_val
      integer, intent(in) :: max_l ! max number of damping iterations
      integer, intent(out) :: l ! number of damping iterations
      real, intent(out) :: ratio ! number of damping iterations
      type( DWR_t), intent(inout), optional :: DWR
      class(element), pointer :: elem
      integer :: kk, elemDof, i, k, j, p_mod, q_mod, nsize
      real :: lambda_exit
      logical :: plus, BiCG


      ! the use of BiCG method for DWR 
      BiCG = .false.
      if(present(DWR) ) then
         if( DWR%BI) BiCG = .true.
      endif

      lambda_exit = ratio**(max_l - 1)

      newton%lambda = min( newton%lambda1, 1.0)   ! initialization of lambda (= damping factor)
      !newton%lambda = 1.0   ! initialization of lambda (= damping factor)
      !newton%lambda1 = 1.0   ! initialization of lambda (= damping factor)
       !max_l = 10

      !print*,'                            ___  _ _ -'
       do l=1,max_l   ! iterations, seeking the optimal damping factor
       !   write(*,'(a10, 3i5, 3es12.4)') &
       !        'ITERS:',state%time%iter, Newton%iter, l, newton%lambda, newton%lambda1

          ! update of the solution
          kk = 0
          do i = 1, grid%nelem
             !!lambda_loc = newton%lambda
             elem => grid%elem(i)
             if (elem%ncv /= kk + 1) then
                print*, 'Problem while copying update to wST (compute.f90)'
                stop
             endif
             elemDof = elem%dof * ndim * elem%Tdof
             elem%vec( res_vec, 1:elemDof) =  newton%rr(kk+1:kk+ elemDof)

             do k = 1, elem%Tdof
                do j = 1, ndim

                   elem%wST(j,1:elem%dof,k) = elem%wST(j,1:elem%dof,k) &
                        + (newton%lambda - newton%lambda_old) &
                        * newton%x(kk + 1 : kk + elem%dof)

                   kk = kk + elem%dof
                   !if(elem%i <=4) &
                   !     write(*,'(a10, 2i5, 300es14.6)') 'b_sol:',elem%i, elemDof,  elem%wST(j,1:elem%dof,k)
                enddo !j
             enddo !k

             if(state%linSolver%lin_solver_not_conv > 0) then
                if(sqrt( dot_product(elem%vec( res_vec, 1:elemDof), elem%vec( res_vec, 1:elemDof))) > &
                     0.5 * res_max_val ) then
                   write(53,*), elem%xc(:), elem%i, &
                        sqrt(dot_product(elem%vec(res_vec, 1:elemDof), elem%vec( res_vec, 1:elemDof))),&
                        res_max_val, elem%wST(1:ndim,1:elem%dof,:), &
                        '  state%linSolver%lin_solver_not_conv te3tg73'
                endif
             endif

          enddo !i

          if(state%linSolver%lin_solver_not_conv > 0 ) close(53)

          newton%lambda_old = newton%lambda

          !call ComputeSTDGM_Terms( deg_plus )
          !print*, 'ComputeSTDGM_Terms was replaced by ComputeST_Terms in setDampingInNewton, deg_plus=', &
          !  deg_plus, "impl = ", newton%implicitly
          call ComputeST_Terms( deg_plus )

          !print*,' PM test (3)'
          !call ComputeSTDGM_Terms(deg_plus ) ! fill elem%rhST

          call FillVectorST( newton%b1(1:state%nsize), eta )

          newton%res0 = VectorPrecondNorm( newton%b1 )
          !print*,'EEEEEEEEEEEEEE!!!!',norm2(Newton%b1)



!          print*, 'new Newton res :' , l,  newton%res0
!          print*, 'old Newton res :' , l, newton%res
          newton%theta  =  newton%res0  / max(1E-15, newton%res)

          !print*,'#########################################################'
          !if(l == 1) print*
          !write(*,'(a25, 3i5, 30es12.4)') '####### new vector:', &
          !     l, max_l, Newton%iter, &
          !     newton%lambda,lambda_exit, newton%res0,newton%res, newton%theta
          !print*,'#########################################################'
          !newton%b1(1:5)

          if(newton%iter > 1) call Newton_output(grid, 55, l, loc_implicitly )

          if(loc_implicitly) call Newton_output( grid, 54, l, loc_implicitly )

          newton%newton_count = newton%newton_count + 1


          newton%lambda1 = newton%lambda/ratio/ratio   ! used for the new Newton iteration
          if(newton%lambda <= lambda_exit ) goto 15  ! too small lambda, NOT CONVERGE
          !if(newton%lambda < 2E-1) goto 15  ! too small lambda, NOT CONVERGE


          ! residuum is not decreasing, reduce lambda
          !if(l < max_l .and. newton%theta >= 1 .and. newton%res0 > 1E-8 ) then  problematic for shock
          if(l < max_l .and. newton%theta >= 1  ) then
          !if(l < max_l .and. newton%theta >= 0.9  ) then
          !if(l < max_l .and. newton%theta >= 0.9 .and. newton%res0 > 1E-8 ) then
             !if(newton%theta > 1 - newton%lambda/4  .and. newton%res0 > 1E-8 )then
             !newton%lambda = min(mu1, 0.5* newton%lambda)
             !newton%lambda = 0.85* newton%lambda
             !newton%lambda = 0.75* newton%lambda
             !newton%lambda = 0.5* newton%lambda
             newton%lambda = ratio * newton%lambda

          else

             ! residuum decreases, quit iterations
             goto 15
          endif


       enddo   !l=1,10    ! iterations, seeking the optimal damping factor

15     continue


       ! update of the dual solution computed by BiCG
       if( BiCG ) then

          p_mod = state%getP_mod()
          q_mod = state%getQ_mod()

          nsize = state%bigNSize( p_mod, q_mod )
          if( nsize /= size( DWR%x, 1 ) ) then
             print*,'Troubles in arrays sizes in Newton.f90 73y7u3',  nsize, size( DWR%x,1 )
             stop
          endif


          if( p_mod > 0) then
             print*,'state%getP_mod() > 0 = ', state%getP_mod()
             print*,'this subroutine has to be re-written, also the elem%wST'
             print*, 'setting of plus below'
             stop
          endif


          plus = .false.

          ! copy the solution vector back to elements zST
          if (plus) then
             call CopyZSTplus_fromLongVector( grid, nsize, p_mod, q_mod, DWR%x(1:nsize) )
          else
             ! optional parameter .true. means "update", zST = zST + DWR%x
             ! otherwise zST = DWR%x
             call CopyZST_fromLongVector( grid, nsize, DWR%x(1:nsize), .true. )
             !call CopyZST_fromLongVector( grid, nsize, DWR%x(1:nsize) )
          end if

          ! Fill residual vector - the residuum is always of the BIGGEST SIZE
          ! CHANGE P_MOD, Q_MOD
          call state%setP_mod( state%p_mod_max)
          call state%setQ_mod( state%q_mod_max)

          ! fills the Dual residuals
          call computeDualResidualPlus(grid, DWR, plus)

          ! setting of the RHS = primal residuals
          call ComputeST_Terms( .true. )

          !print*,' PM test (4)'
          !call ComputeSTDGM_Terms(.true. ) ! fill elem%rhST


          !write(*,'(a10, i5, 200es12.4)') 'damp WX:', grid%elem(1)%i,  grid%elem(1)%rhsST

          call state%setP_mod( p_mod )
          call state%setQ_mod( q_mod )

          DWR%dualProblem_computed = .true.

          ! ONLY OUTPUT
          ! do i = 1, grid%nelem
          !    elem => grid%elem(i)
          !    do k = 1, elem%Tdof
          !       do j = 1, ndim
          !          kk = kk + elem%dof
          !          if(elem%i <=4) &
          !               write(*,'(a10, 2i5, 300es14.6)') 'D_sol:',elem%i, elemDof,  elem%zST(j,1:elem%dof,k)
          !       enddo !j
          !    enddo !k
          ! enddo

       endif

   end subroutine setDampingInNewton



   !> seeking of optimal damping parameter in the Newton method,
   !> using Anderson acceleration
   !> Not tested, copied from the old subroutine
   subroutine setAndersonInNewton( newton, grid, deg_plus, eta, loc_implicitly)
      class( Newton_t ), intent ( inout ) :: newton
      class( mesh ), intent(inout) :: grid
      logical, intent(in) :: deg_plus
      real, intent(in) :: eta
      logical, intent(in) :: loc_implicitly
      class(element), pointer :: elem
      integer :: l, i, j, nsize
      integer :: mk, AA_F, AA_m, AA_w, AA_dF, AA_max
      real, dimension(:,:), allocatable :: dd1
      real, dimension(:,:), pointer :: dd
      real, dimension(:), pointer :: alpha
      real :: vall

      !write(*,'(a8, 3i5, 1es10.2, a1, 30es12.4)') 'newton:',state%time%iter, Newton%iter, l, &
      !     Newton%lambda,'|', &
      !     Newton%res0, Newton%res, Newton%theta, state%linSolver%residuum

      nsize = state%nsize

      AA_F = newton%andersonAcc%AA_F
      AA_m = newton%andersonAcc%AA_m
      AA_dF = newton%andersonAcc%AA_dF
      AA_w = newton%andersonAcc%AA_w
      AA_max = newton%andersonAcc%AA_max


      !! ANDERSON ACCELERATION
      ! maximal number of used vectors
      !mk = min( newton%iter, AA_max) - 1
      mk = min( newton%iter -1, AA_max)

      allocate( dd1(1:mk, mk) )  ! temporary arrat

      ! (AA_w, i, :) contains \Delta g_i = g_{k-1}- g_{k-i+1},  g_i:= G(u^i) - u^i
      ! (AA_F, i, :) contains \Delta G_i = G_{k-1}- G_{k-i+1},

      alpha => newton%andersonAcc%alpha(1:mk)
      dd => newton%andersonAcc%dd(1:mk, 1:mk)

      ! shifting the vectors from the previous iterations
      do l =  mk, 2, -1
         newton%andersonAcc%AA_x(1:AA_m, l, 1:nsize) = newton%andersonAcc%AA_x(1:AA_m, l-1, 1:nsize)
      enddo

      ! shifting the matrix
      dd1(2:mk, 2:mk) = dd(1:mk-1, 1:mk-1)
      dd(2:mk, 2:mk) = dd1(2:mk, 2:mk)

      !! new values for the  Anderson acceleration:
      ! newton%andersonAcc%AA_update(1:nsize) is now a temporary field, now it contains g_k

      ! var (i) : difference  G(w^k) - w^k := d^{k+1} = w^{k+1} - w^{k}
      newton%andersonAcc%AA_update(1:nsize) =  newton%x(1:nsize) * newton%lambda

      ! var (ii) : difference  G(w^k) - w^k := F^{k+1}
      !newton%andersonAcc%AA_update(1:nsize) =  newton%b1(1:nsize) !* Newton%lambda

      ! var (iii) : difference  G(w^k) - w^k :=  P F^{k+1}
      !call VectorPrecond(newton%b1(1:nsize), newton%andersonAcc%AA_update(1:nsize) )


      ! the new difference  \Delta g_1,  (AA_w, 0, :) contains g_{k-1} at this moment
      newton%andersonAcc%AA_x(AA_w, 1, 1:nsize) = &
           newton%andersonAcc%AA_update(1:nsize) -  newton%andersonAcc%AA_x(AA_w, 0, 1:nsize)

      ! storing of the new g_{k-1} = g_k
      newton%andersonAcc%AA_x(AA_w, 0, 1:nsize) = newton%andersonAcc%AA_update(1:nsize)


      !  newton%andersonAcc%AA_update(1:nsize) is now a temporary field, now it contains G(u^k)
      call CopyWST_toLongVector( grid, nsize, newton%andersonAcc%AA_update(1:nsize) )

      ! var (i) already OK

      ! var (ii)
      !newton%andersonAcc%AA_update(1:nsize) = newton%andersonAcc%AA_update(1:nsize) &
      !     -  newton%x(1:nsize) * newton%lambda


      ! the new difference  \Delta G_1,  (AA_F, 0, :) contains G_{k-1}
      newton%andersonAcc%AA_x(AA_F, 1, 1:nsize) = &
           newton%andersonAcc%AA_update(1:nsize) -  newton%andersonAcc%AA_x(AA_F, 0, 1:nsize)

      ! storing of the new G_{k-1} = G_k
      newton%andersonAcc%AA_x(AA_F, 0, 1:nsize) = newton%andersonAcc%AA_update(1:nsize)


      !do i=1,max(mk,0)
      !   write(*,'(a8, i5, 30es12.4)') ' Delta g::',i, newton%andersonAcc%AA_x(AA_w, i,1:6)
      !enddo
      !print*
      !do i=1,max(mk,0)
      !   write(*,'(a8, i5, 30es12.4)') ' Delta G::',i, newton%andersonAcc%AA_x(AA_F, i,1:6)
      !enddo
      !print*
      !write(*,'(a8, i5, 30es12.4)') ' G(u^k)::',0, newton%andersonAcc%AA_update(1:6)

      ! Anderson update
      if(mk == 0) then
         ! NO ACTION

       else

          ! computing of the scalar products
          do i=1, mk
             ! RHS
             alpha(i) = -dot_product(newton%andersonAcc%AA_x(AA_w, 0, :),  &
                  newton%andersonAcc%AA_x(AA_w, i,: ))

             ! matrix
             dd(1,i) = dot_product(newton%andersonAcc%AA_x(AA_w, 1, :), &
                     newton%andersonAcc%AA_x(AA_w, i, : )  )
             dd(i,1) = dd(1,i)

          enddo

          !print*,'.,.'
          !do i=1,mk
          !   write(*,'(a8, i5, 30es12.4)') 'dd:',i, dd(i, :)
          !enddo
          !print*,' ___'
          !write(*,'(a8, i5, 30es12.4)') 'RHS:',mk, alpha(:)
          !print*

          ! assembling of the linear algebraic problem
          dd1(1:mk, 1:mk) = dd(1:mk, 1:mk)
          !newton%andersonAcc%alpha(1:mk) =  -newton%andersonAcc%dd(1:mk, 0)
          call SolveLocalMatrixProblem(mk, dd1(1:mk, 1:mk), 1, alpha(1:mk) )

          !write(*,'(a8, i5, 30es12.4)') 'gamma:',mk, alpha(:)

          ! setting of the new update,   newton%andersonAcc%AA_update  already contains G(u^k)
          do i=1,mk
             newton%andersonAcc%AA_update(:) = newton%andersonAcc%AA_update(:) &
                  + newton%andersonAcc%alpha(i) *  newton%andersonAcc%AA_x(AA_F, i, :)
             !write(*,'(a10, 300es12.4)') 'update:', Newton%AA_update(1:10)
          enddo

          ! setting the new solution to wST
          !call Fill_wST_from_VectorST( newton%andersonAcc%AA_update(:) )
          ! in solution_mod
          call CopyWST_fromLongVector( grid, nsize, newton%andersonAcc%AA_update(:) )

          ! update rhs of the new system
          !call ComputeSTDGM_Terms( deg_plus )
!          print*, 'ComputeSTDGM_Terms was replaced by ComputeST_Terms in setAndersonInNewton, impl = ', &
!            state%nlSolver%implicitly
          call ComputeST_Terms( deg_plus )

          !print*,' PM test (5) '
          !call ComputeSTDGM_Terms(deg_plus ) ! fill elem%rhST


          ! put the rhs into newton%b1
          call FillVectorST( newton%b1(1:nsize), eta )
          ! the actual residuum

          newton%res0 = VectorPrecondNorm( newton%b1 )


          !write(*,'(a8, 3i5, 1es10.2, a1, 30es12.4)') 'anderson:',state%time%iter, Newton%iter, -1, &
          !     Newton%lambda,'|', &
          !     Newton%res0, Newton%res, Newton%theta, state%linSolver%residuum
          !
          !print*,'___________________________________________________________________________'
          !stop "i30id03ik3ew"

       endif

       !print*,'_______________________________________________________', mk


       deallocate(dd1)

   end subroutine setAndersonInNewton


!   !> seeking of optimal damping parameter in the Newton method,
!   !> using Anderson acceleration
!   !> Not tested, copied from the old subroutine
!   subroutine setAndersonInNewton_OLD( newton, grid, deg_plus, eta, loc_implicitly)
!      class( Newton_t ), intent ( inout ) :: newton
!      class( mesh ), intent(inout) :: grid
!      logical, intent(in) :: deg_plus
!      real, intent(in) :: eta
!      logical, intent(in) :: loc_implicitly
!      class(element), pointer :: elem
!      integer :: l, i, j, nsize
!      integer :: mk, AA_F, AA_m, AA_w, AA_dF, AA_max
!      real, dimension(:,:), allocatable :: dd1
!      real, dimension(:,:), pointer :: dd
!      real, dimension(:), pointer :: alpha
!      real :: vall
!
!      write(*,'(a8, 3i5, 1es10.2, a1, 30es12.4)') 'newton:',state%time%iter, Newton%iter, l, &
!           Newton%lambda,'|', &
!           Newton%res0, Newton%res, Newton%theta, state%linSolver%residuum
!
!      nsize = state%nsize
!
!      AA_F = newton%andersonAcc%AA_F
!      AA_m = newton%andersonAcc%AA_m
!      AA_dF = newton%andersonAcc%AA_dF
!      AA_w = newton%andersonAcc%AA_w
!      AA_max = newton%andersonAcc%AA_max
!
!      !! ANDERSON ACCELERATION
!      ! maximal number of used vectors
!      !mk = min( newton%iter, AA_max) - 1
!      mk = min( newton%iter -1, AA_max)
!
!      alpha => newton%andersonAcc%alpha(1:mk)
!      dd => newton%andersonAcc%dd(1:mk, 1:mk)
!
!      ! storing for the vectors for the Anderson acceleration from the previous iterations
!      do l =  mk, 1, -1
!         newton%andersonAcc%AA_x(1:AA_m, l, 1:nsize) = newton%andersonAcc%AA_x(1:AA_m, l-1, 1:nsize)
!      enddo
!
!      !! new values for the  Anderson acceleration:
!      ! var (i) : difference  G(w^k) - w^k := d^{k+1} = w^{k+1} - w^{k}
!      newton%andersonAcc%AA_x(AA_F, 0, 1:nsize) =  newton%x(1:nsize) * newton%lambda
!
!      ! var (ii) : difference  G(w^k) - w^k := F^{k+1}
!      !newton%andersonAcc%AA_x(AA_F, 0, 1:nsize) =  newton%b1(1:nsize) !* Newton%lambda
!
!      ! var (iii) : difference  G(w^k) - w^k :=  P F^{k+1}
!      !call VectorPrecond(newton%b1(1:nsize), newton%andersonAcc%AA_x(AA_F, 0, 1:nsize) )
!
!      ! difference of the old updates and the new one
!      do l =  mk, 1, -1
!         newton%andersonAcc%AA_x(AA_dF, l, 1:nsize) = &
!              newton%andersonAcc%AA_x(AA_F, l, 1:nsize) - newton%andersonAcc%AA_x(AA_F, 0, 1:nsize)
!      enddo
!
!      ! result of updated w^{k+1}= G(w^k)
!      !call FillVectorST_from_wST( newton%andersonAcc%AA_x(AA_w, 0, 1:nsize) )
!      call CopyWST_toLongVector( grid, nsize, newton%andersonAcc%AA_x(AA_w, 0, 1:nsize) )
!
!
!      ! Anderson update
!      if(mk == 0) then
!         ! only computation of the scalar product
!         ! NO ACTION
!         !newton%andersonAcc%dd(0,0) = dot_product(newton%andersonAcc%AA_x(AA_F, 0, :), &
!         !      newton%andersonAcc%AA_x(AA_F, 0, : )  )
!          !write(*,'(a8,i5, 20es12.4)') 'dd:',0, dd(0, :)
!
!       else
!
!          ! computing of the scalar products
!          do i=1, mk
!             ! RHS
!             alpha(i) = -dot_product(newton%andersonAcc%AA_x(AA_F, 0, :),  &
!                  newton%andersonAcc%AA_x(AA_dF, i,: ))
!
!             ! matrix
!             do j=i, mk
!                dd(i,j) = dot_product(newton%andersonAcc%AA_x(AA_dF, i, :), &
!                     newton%andersonAcc%AA_x(AA_dF, j, : )  )
!                dd(j,i) = dd(i,j)
!             enddo
!
!          enddo
!
!          ! assembling of the linear algebraic problem
!          allocate( dd1(1:mk, mk) )
!
!          !dd1(1:mk, 1:mk) = dd(1:mk, 1:mk)
!          !newton%andersonAcc%alpha(1:mk) =  -newton%andersonAcc%dd(1:mk, 0)
!          call SolveLocalMatrixProblem(mk, dd(1:mk, 1:mk), 1, newton%andersonAcc%alpha(1:mk) )
!
!
!          newton%andersonAcc%alpha(0) = 1. - sum (newton%andersonAcc%alpha(1:mk))
!          ! alpha(:) contains the correct values of the minimum
!          !deallocate(dd1)
!
!
!          ! setting of the better update
!          newton%andersonAcc%AA_update = 0.
!          do i=0,mk
!             newton%andersonAcc%AA_update(:) = newton%andersonAcc%AA_update(:) &
!                  + newton%andersonAcc%alpha(i) *  newton%andersonAcc%AA_x(AA_w, i, :)
!             !write(*,'(a10, 300es12.4)') 'update:', Newton%AA_update(1:10)
!          enddo
!
!          ! setting the new solution to wST
!          !call Fill_wST_from_VectorST( newton%andersonAcc%AA_update(:) )
!          ! in solution_mod
!          call CopyWST_fromLongVector( grid, nsize, newton%andersonAcc%AA_update(:) )
!
!          ! update rhs of the new system
!          !call ComputeSTDGM_Terms( deg_plus )
!          print*, 'ComputeSTDGM_Terms was replaced by ComputeST_Terms in setAndersonInNewton'
!          call ComputeST_Terms( deg_plus )
!
!          ! put the rhs into newton%b1
!          call FillVectorST( newton%b1(1:nsize), eta )
!          ! the actual residuum
!
!          newton%res0 = VectorPrecondNorm( newton%b1 )
!
!          write(*,'(a8, 3i5, 1es10.2, a1, 30es12.4)') 'anderson:',state%time%iter, Newton%iter, -1, &
!               Newton%lambda,'|', &
!               Newton%res0, Newton%res, Newton%theta, state%linSolver%residuum
!
!          print*,'___________________________________________________________________________'
!          !stop "i30id03ik3ew"
!
!       endif
!
!     end subroutine setAndersonInNewton_OLD

!      ! seeking of optimal damping parameter in the Newton method,
!   ! this should help the convergence of Newton when the approximation is still rough
!   subroutine setDampingInNewtonLocal( newton, grid, deg_plus, eta, loc_implicitly, res_max_val, max_l, l)
!      class( Newton_t ), intent ( inout ) :: newton
!      class( mesh ), intent(inout) :: grid
!      logical, intent(in) :: deg_plus
!      real, intent(in) :: eta
!      logical, intent(in) :: loc_implicitly
!      real, intent(in) :: res_max_val
!      integer, intent(in) :: max_l ! max number of damping iterations
!      integer, intent(out) :: l ! number of damping iterations
!      class(element), pointer :: elem
!      integer :: kk, elemDof, i, k, j
!      real, dimension(1:ndim) :: temp
!      real, dimension(:), allocatable :: estimOld, estim
!      logical, dimension(:), allocatable :: damping
!      real :: theta
!      integer :: numberOfUpdates, nelem, dof, ncv, nsize, nDamp
!
!      print*, 'setLocalDamping called! '
!
!      nelem = grid%nelem
!
!      newton%lambda = 1.0   ! initialization of lambda (= damping factor)
!
!      allocate( estimOld(1:nelem), source = 0.0 )
!      allocate( estim(1:nelem), source = 0.0 )
!      ! continue with damping on this element
!      allocate( damping(1:nelem) , source = .true. )
!      ! TODO
!      ! there should be some flag for the case when all damping == .false.
!
!      print*, 'newton%res 1 ' , newton%res
!      print*, 'why ere these residuals different??? '
!
!      !!!!!!! compute locally the OLD estimate
!      call localNewtonEstimate( newton, grid, eta, estimOld(1:nelem) )
!      ! global estimate
!      newton%res = VectorPrecondNorm( newton%b1 )
!      print*, 'newton%res 2 ' , newton%res
!      ! should be same as newton%res
!      !try = dot_product( estimOld(:), estimOld(:) )**0.5
!
!      do l=1,max_l   ! iterations, seeking the optimal damping factor
!          ! update of the solution
!          kk = 0
!          do i = 1, grid%nelem
!             !!lambda_loc = newton%lambda
!             elem => grid%elem(i)
!             if (elem%ncv /= kk + 1) then
!                print*, 'Problem while copying update to wST (compute.f90)'
!                stop
!             endif
!             elemDof = elem%dof * ndim * elem%Tdof
!             elem%vec( res_vec, 1:elemDof) =  newton%rr(kk+1:kk+ elemDof)
!
!             do k = 1, elem%Tdof
!                do j = 1, ndim
!                  if ( damping(i) ) then
!                     elem%wST(j,1:elem%dof,k) = elem%wST(j,1:elem%dof,k) &
!                        + (newton%lambda - newton%lambda_old) &
!                        * newton%x(kk + 1 : kk + elem%dof)
!                  endif ! damping
!                   kk = kk + elem%dof
!                   !write(*,'(a10, 2i5, 300es12.4)') 'b_sol:',elem%i, elemDof,  elem%wST(j,1:elem%dof,k)
!                enddo !j
!             enddo !k
!
!             if(state%linSolver%lin_solver_not_conv > 0) then
!                if(sqrt( dot_product(elem%vec( res_vec, 1:elemDof), elem%vec( res_vec, 1:elemDof))) > &
!                     0.5 * res_max_val ) then
!                   write(53,*), elem%xc(:), elem%i, &
!                        sqrt(dot_product(elem%vec(res_vec, 1:elemDof), elem%vec( res_vec, 1:elemDof))),&
!                        res_max_val, elem%wST(1:ndim,1:elem%dof,:)
!                endif
!             endif
!
!          enddo !i
!
!          if(state%linSolver%lin_solver_not_conv > 0 ) close(53)
!          newton%lambda_old = newton%lambda
!
!          ! compute the error locally
!          !call ComputeSTDGM_Terms( deg_plus )
!          print*, 'ComputeSTDGM_Terms replaced by ComputeSTDGM_Terms in '
!          call ComputeST_Terms( deg_plus )
!
!          ! compute the local version of the residual estimate for each element
!          call localNewtonEstimate( newton, grid, eta, estim(1:nelem) )
!          ! global estimate
!
!          newton%res0 = VectorPrecondNorm( newton%b1 )
!
!!          print*, 'new Newton res :' , l,  newton%res0
!!          print*, 'old Newton res :' , l, newton%res
!
!          newton%theta  =  newton%res0  / max(1E-15, newton%res)
!
!          if(newton%iter > 1) call Newton_output(grid, 55, l, loc_implicitly )
!          if(loc_implicitly) call Newton_output( grid, 54, l, loc_implicitly )
!          newton%newton_count = newton%newton_count + 1
!
!          ! write # of damped elements
!          nDamp = 0
!          do i =1,nelem
!            if (damping(i) ) then
!               nDamp = nDamp + 1
!            endif
!          enddo
!          print*, '#elem for which damping in the', l, &
!                  ' cycle was used', nDamp
!
!          ! residuum is not decreasing, reduce lambda
!          if(newton%theta >= 1 .and. newton%res0 > 1E-8 ) then
!             newton%lambda = 0.5* newton%lambda
!
!             ! local dumping
!             !!!!!!!!!!!!!!!!!! IF WE COMMENT this it should work as standard estimates !!
!             do i=1,nelem
!                theta = estim(i) /  max(1E-15, estimOld(i) )
!                if (theta < 1) then
!                   damping(i) = .false.
!                 ! we cannot start damping again - we would need special lambda for each element  !
!!                else
!!                  damping(i) = .true.
!                endif
!             end do
!
!          ! the norm has decreased globally -> no dumping
!          else
!             newton%lambda1 = min(1., 1.5*newton%lambda )
!             ! residuum decreases, quit iterations
!             goto 19
!          endif
!
!          if(newton%lambda < 1E-2) goto 19  ! too small lambda, NOT CONVERGE
!
!
!       enddo   !l=1,10    ! iterations, seeking the optimal damping factor
!
!19     continue
!
!       print*, 'number of used dampings:', l
!       ! number of damped elements in the last cycle
!
!       deallocate( estim, estimOld, damping )
!
!   end subroutine setDampingInNewtonLocal

   ! seeking of optimal damping parameter in the Newton method,
   ! this should help the convergence of Newton when the approximation is still rough
   subroutine setDampingInNewton_DWR( newton, grid, DWR, deg_plus, eta, &
                                      loc_implicitly, res_max_val, max_l, l)
      class( Newton_t ), intent ( inout ) :: newton
      class( mesh ), intent(inout) :: grid
      class( DWR_t ), intent(inout) :: DWR
      logical, intent(in) :: deg_plus
      real, intent(in) :: eta   ! not used
      logical, intent(in) :: loc_implicitly
      real, intent(in) :: res_max_val
      integer, intent(in) :: max_l ! max number of damping iterations
      integer, intent(out) :: l ! number of damping iterations
      class(element), pointer :: elem
      integer :: kk, elemDof, i, k, j

      newton%lambda = 1.0   ! initialization of lambda (= damping factor)

      ! may not be needed
      ! compute the estimate for the old wST
      call computeNonlinDWRestimates( DWR, grid )
      newton%res = DWR%estimNL

      do l=1,max_l   ! iterations, seeking the optimal damping factor
         ! update of the solution
         !print*, 'newton lambda old:', newton%lambda_old

          kk = 0
          do i = 1, grid%nelem

             elem => grid%elem(i)
             if (elem%ncv /= kk + 1) then
                print*, 'Problem while copying update to wST (compute.f90)'
                stop
             endif
             elemDof = elem%dof * ndim * elem%Tdof
             elem%vec( res_vec, 1:elemDof) =  newton%rr(kk+1:kk+ elemDof)

             do k = 1, elem%Tdof
                do j = 1, ndim
                   elem%wST(j,1:elem%dof,k) = elem%wST(j,1:elem%dof,k) &
                        + (newton%lambda - newton%lambda_old) &
                        * newton%x(kk + 1 : kk + elem%dof)
                   kk = kk + elem%dof
                enddo !j
             enddo !k

             if(state%linSolver%lin_solver_not_conv > 0) then
                if(sqrt( dot_product(elem%vec( res_vec, 1:elemDof), elem%vec( res_vec, 1:elemDof))) > &
                     0.5 * res_max_val ) then
                   write(53,*), elem%xc(:), elem%i, &
                        sqrt(dot_product(elem%vec(res_vec, 1:elemDof), elem%vec( res_vec, 1:elemDof))),&
                        res_max_val, elem%wST(1:ndim,1:elem%dof,:)
                endif
             endif

          enddo !i

          if(state%linSolver%lin_solver_not_conv > 0 ) close(53)
          newton%lambda_old = newton%lambda

          !call ComputeSTDGM_Terms( deg_plus )
          ! rhsST -> newton%b1
          !call FillVectorST( newton%b1(1:state%nsize), eta )

          call computeNonlinDWRestimates( DWR, grid )
          newton%res0 = DWR%estimNL
!          print*, 'new Newton res :' , l,  newton%res0
!          print*, 'old Newton res :' , l, newton%res

          newton%theta  =  newton%res0  / max(1E-15, newton%res)

          if(newton%iter > 1) call Newton_output(grid, 55, l, loc_implicitly )
          if(loc_implicitly) call Newton_output( grid, 54, l, loc_implicitly )

          newton%newton_count = newton%newton_count + 1

          ! residuum is not decreasing, reduce lambda
          !
          if(newton%theta >= 1 .and. newton%res0 > 1E-8 ) then
             newton%lambda = 0.5* newton%lambda
          else
             newton%lambda1 = min(1., 1.5*newton%lambda )
             goto 16
          endif

          if(newton%lambda < 1E-2) goto 16  ! too small lambda, NOT CONVERGE

      enddo   !l=1,10    ! iterations, seeking the optimal damping factor

16    continue

   end subroutine setDampingInNewton_DWR

      ! seeking of optimal damping parameter in the Newton method,
   ! this should help the convergence of Newton when the approximation is still rough
   subroutine setLocalDampingDWR( newton, grid, DWR, deg_plus, eta, &
                                      loc_implicitly, res_max_val, max_l, l)
      class( Newton_t ), intent ( inout ) :: newton
      class( mesh ), intent(inout) :: grid
      class( DWR_t ), intent(inout) :: DWR
      logical, intent(in) :: deg_plus
      real, intent(in) :: eta   ! not used
      logical, intent(in) :: loc_implicitly
      real, intent(in) :: res_max_val
      integer, intent(in) :: max_l ! max number of damping iterations
      integer, intent(out) :: l ! number of damping iterations
      class(element), pointer :: elem
      integer :: kk, elemDof, i, k, j
      real, dimension(1:ndim) :: temp
      real, dimension(:), allocatable :: estimOld
      logical, dimension(:), allocatable :: damp
      real :: theta
      logical :: flag
      integer :: numberOfUpdates, nelem

      print*, 'setLocalDampingDWR called! '

      nelem = grid%nelem
      newton%lambda = 1.0   ! initialization of lambda (= damping factor)

      allocate( estimOld(1:nelem), source = 0.0 )
      allocate( damp(1:nelem) , source = .true. )


      ! compute the estimate for the old wST
      ! not needed -- was already called before
      call computeNonlinDWRestimates( DWR, grid )

      do i = 1,nelem
         elem => grid%elem(i)
         estimOld(i) = elem%eta(dwrA,1)
      end do !i
      newton%res = DWR%estimNL


      do l=1,max_l   ! iterations, seeking the optimal damping factor
          ! update of the solution
          !print*, 'newton lmbda old:', newton%lambda_old

          kk = 0
          do i = 1, grid%nelem

             elem => grid%elem(i)
             if (elem%ncv /= kk + 1) then
                print*, 'Problem while copying update to wST (compute.f90)'
                stop
             endif
             elemDof = elem%dof * ndim * elem%Tdof
             elem%vec( res_vec, 1:elemDof) =  newton%rr(kk+1:kk+ elemDof)

             do k = 1, elem%Tdof
                do j = 1, ndim
                  if ( damp(i) ) then
                     elem%wST(j,1:elem%dof,k) = elem%wST(j,1:elem%dof,k) &
                        + (newton%lambda - newton%lambda_old) &
                        * newton%x(kk + 1 : kk + elem%dof)
                   end if
                   kk = kk + elem%dof
                enddo !j
             enddo !k

             if(state%linSolver%lin_solver_not_conv > 0) then
                if(sqrt( dot_product(elem%vec( res_vec, 1:elemDof), elem%vec( res_vec, 1:elemDof))) > &
                     0.5 * res_max_val ) then
                   write(53,*), elem%xc(:), elem%i, &
                        sqrt(dot_product(elem%vec(res_vec, 1:elemDof), elem%vec( res_vec, 1:elemDof))),&
                        res_max_val, elem%wST(1:ndim,1:elem%dof,:)
                endif
             endif

          enddo !i

          if(state%linSolver%lin_solver_not_conv > 0 ) close(53)
          newton%lambda_old = newton%lambda

          !call ComputeSTDGM_Terms( deg_plus )
          ! rhsST -> newton%b1
          !call FillVectorST( newton%b1(1:state%nsize), eta )

          call computeNonlinDWRestimates( DWR, grid )
          newton%res0 = DWR%estimNL
!          print*, 'new Newton res :' , l,  newton%res0
!          print*, 'old Newton res :' , l, newton%res

          newton%theta  =  newton%res0  / max(1E-15, newton%res)

          ! residuum is not decreasing, reduce lambda
          do i = 1, nelem
            elem => grid%elem(i)
            theta = elem%eta(dwrA, 1) / estimOld(i)
            if (theta < 1) then
               !damp(i) = .false.
            endif
          end do

          if(newton%iter > 1) call Newton_output(grid, 55, l, loc_implicitly )
          if(loc_implicitly) call Newton_output( grid, 54, l, loc_implicitly )

          newton%newton_count = newton%newton_count + 1

          newton%lambda = 0.5* newton%lambda
          if(newton%lambda < 1E-2) goto 16  ! too small lambda, NOT CONVERGE

          damp(:) = .false.
          damp(1) = .true.

      enddo   !l=1,10    ! iterations, seeking the optimal damping factor

16    continue

      print*, 'lambda is equal:' , newton%lambda
      stop 'stopped in setLocalDampingDWR'

      deallocate( estimOld, damp )

   end subroutine setLocalDampingDWR


   !> perform one step of the Newton method, i.e.
   !> solve the linear alg. system given by C(u)d^k = -F(u)
   !> solution of the problem is done by GMRES and its parameters: nloops and restart are set here !
   !> but other methods may be implemented in SolveBlockLinearSTDGMProblem later
   !> the tolerance on the linear algebraic error is saved in state%linSolver%tol
   !> When Newton%non_alg_stop == 'aDWR' then the DWR_t structure is needed to compute the estimates of the alg. error
   subroutine performOneNewtonStdgmStep( newton, grid, imp, iter, &
      deg_plus, newtonDone, res_max_val, loc_implicitly, DWR )
      class( NonlinearSol_t ), intent ( inout ) :: newton ! this should the Newton type
      class( mesh ), intent(inout) :: grid
      integer, intent(in) :: imp
      integer, intent(in) :: iter
      logical, intent(in) :: deg_plus
      logical, intent(out) :: newtonDone
      real, intent(out) :: res_max_val
      logical, intent(out) :: loc_implicitly
      type( DWR_t), intent(inout), optional :: DWR

      logical :: update, BiCG
      class(element), pointer :: elem
      integer :: max_l, l, kk, elemDof, i, k, j, p_mod, q_mod
      integer :: restart, nloops, nsize, bigNsize
      logical :: vector_update
      real, allocatable, dimension(:) :: vecTest, outTest, outTest2, vecTestBig
      real, allocatable, dimension(:) :: newtonB, NewtonX
      real, allocatable, dimension(:) :: x0, y0, b0, c0, d0, e0
      real :: xi_initP, xi_initD, J1, J2, J1a, J2a, xi_P, xi_D
      !real :: t1, t2, time_prepare, tt, tt1

      ! the use of BiCG method for DWR 
      BiCG = .false.
      if(present(DWR) ) then
         if( DWR%BI) BiCG = .true.
      endif

      xi_initP = 0.
      xi_initD = 0.

      newtonDone = .false.

      call state%cpuTime%startPrepareTime()
      if(imp == 0)  newton%implicitly = .true.

      if ( state%nlSolver%non_alg_stop == 'aDWR') then
         write(debug,*) 'There may be problems with implicitly!'
         !print*, 'Probably in the matrix there are values with deg_plus!'
         write(debug,*) 'Set implicitly = .true.'
         newton%implicitly = .true.
      endif

      loc_implicitly = newton%implicitly
      update = .false.

      !print*
      !print*,'newton%implicitlyDE#G=',newton%implicitly, state%nlSolver%iter

      if( newton%implicitly ) then
          update = .true.
          !deg_plus = .true.


          !call ComputeST_Terms( .false. )
          if (state%space%estim_space == "DWR" .and. state%model%linear) then ! VD Newton
             !print*, "linProblem: ComputeST_Terms in performOneNewtonStdgmStep called with impl true, deg_plus = True"
             call ComputeST_Terms( .true. )
          else
             !print*, "nonlinProblem: ComputeST_Terms in performOneNewtonStdgmStep called with impl true, deg_plus = False"
             call ComputeST_Terms( .false. )
             !print*,' PM test (A)'
             
          end if
          !call ComputeST_Terms( .true. )

          !print*, 'PM test ComputeSTDGM_Terms commented with impl true! ',' deg_plus=', deg_plus
          !call ComputeSTDGM_Terms( .false. )
          !!call ComputeSTDGM_Terms( deg_plus )

          !call WriteMatrixST_Blocks(0.)

          newton%newton_count = 0
          newton%updates = Newton%updates + 1
          newton%implicitly = .false.

          !print*, 'ComputeST_Terms in performOneSTDGMStep with impl = false degPlus= ', deg_plus
          !write(31,'(a10, i5, 200es12.4)') 'RHS ..:', grid%elem(1)%i,  grid%elem(1)%rhsST

          call ComputeST_Terms( deg_plus ) ! fill elem%rhST
          !print*,' PM test (Q)'
          !write(31,'(a10, i5, 200es12.4)') 'RHS !!:', grid%elem(1)%i,  grid%elem(1)%rhsST

          !print*,' PM test '
          !call ComputeSTDGM_Terms(deg_plus ) ! fill elem%rhST



          !do i=1,grid%nelem
          !   elem => grid%elem(i)
          !   write(*,'(a10, i5, 300es12.4)') '!! elm:', i, elem%rhsST(:,:, :)
          !enddo

           ! FR test only
!           allocate( vecTest(1:state%nsize) , source = 0.0 )
!           allocate( outTest(1:state%nsize) , source = 0.0 )
!
!           eta = 1./state%time%tau(1)
!           call FillVectorST(vecTest(1:state%nsize), eta ) ! value of eta

!           state%p_mod = state%p_mod_max
!           state%q_mod = state%q_mod_max
!           bigNsize = state%bigNSize(state%p_mod, state%q_mod)
!           allocate( vecTestBig(1:bigNsize) , source = 0.0 )
!           vecTestBig(1:state%nsize) = vecTest(1:state%nsize)
!           allocate( outTest2(  1: bigNsize) , source = 0.0 )
!
!           call bMVprodST(outTest,vecTest, state%nsize)
!           print*, 'Norm of the result = ' , norm2( outTest )
!
!           call bMVprodBIG(outTest2, vecTestBig, bigNsize )
!           print*, 'Norm of sec the result = ' , norm2( outTest2 )
!           state%p_mod = 0
!           call bMVprodBIG(outTest2(1:state%nsize), vecTestBig(1:state%nsize), state%nsize )
!           print*, 'Diff in the result = ' , norm2( outTest2(1:state%nsize) - outTest )

           ! allocate long vector from wST
           !deallocate(vecTest, vecTestBig, outTest, outTest2 )
           !stop 'test - not done !'

          newton%implicitly = .true.

          vector_update = .true.
          newton%implicitly = .false.

      ! computing with the matrices from previous time step ?
      else

          ! for iter > 1, array b(:) = F(x^k)
          if(iter == 1) then
             !print*, 'ComputeST_Terms in performOneSTDGMStep, degPlus=', deg_plus, &
             !  "impl = ", newton%implicitly
             call ComputeST_Terms(deg_plus )

             !print*,' PM test (2)'
             !call ComputeSTDGM_Terms(deg_plus )

             !write(*,'(a30, 2es12.4)') 'CT_ends   3 :', tt - state%start_time,  tt - tt1
            vector_update = .true.
          else
             vector_update = .false.
             Newton%b(1:state%nsize) = Newton%b1(1:state%nsize)
             Newton%res = Newton%res0
          endif
      endif ! end of if(state%nlSolver%implicitly)


       eta = 1./state%time%tau(1)

       if(vector_update) then  !filling elem%rhsST into global rhs
          !state%nlSolver%pseudo_time%PTS_rho = 1.
          !print*,' PM test (X)'
          call FillVectorST(newton%b(1:state%nsize), eta ) ! value of eta is not used, only its sign
          Newton%res = VectorPrecondNorm(Newton%b)
          !print*,'EEEEEEEEEEEEEE????',norm2(Newton%b)  ! BICG, residuum can increase after one NEWTON STEP,   WHY???????????????????
       endif

       if(iter == 1) newton%res_ref = Newton%res

       ! first exit based on the absolute value of the residuum vector || F(w) ||
       !if(iter > 1 .and. newton%res < newton%tol &  DOES NOT WORKS !!
       !if(iter > 1 .and. newton%res/state%nsize < newton%tol &
       if(iter > 1 ) then
          if( newton%res < newton%tol .and. state%nlSolver%non_alg_stop /= 'aDWR' ) then
             print*,' #E#E#E#  NEWTON DONE STOP:' , newton%res/state%nsize , newton%tol
             newtonDone = .true.
          endif
       endif

       newton%iter  = newton%iter + 1
       newton%Aiter = newton%Aiter  + 1

       newton%lambda_old = 0.
       newton%x(1:state%nsize) = 0.

       if (.not. state%linSolver%tol_fixed) then
          if(iter == 1) then
             state%linSolver%tol = 0.25
          else
             state%linSolver%tol = min(0.25, abs(newton%res - state%linSolver%residuum) / newton%res1)
             state%linSolver%tol = max(state%linSolver%tol, 1E-6)
          endif
       endif

       call state%cpuTime%addPrepareTime()

       !call WriteMatrixST_Blocks(eta)

       !call test_bMVprodST()
       !call Write_rhsST()
       !stop "u834398dj3"

       ! set number of iterations and number of nloops of lin alg. solver
       if (newton%non_alg_stop == 'aDWR') then
         ! set number of iter after which restart is done
         if ( present(DWR) ) then
            restart = DWR%aDWR%restart_primal
            nloops = 1
            DWR%aDWR%iter = DWR%aDWR%iter + 1

            ! set the right tolerance
            state%linSolver%tol = DWR%aDWR%linTol * DWR%aDWR%C_Safe
            !print*, 'nlTol set to : ' , DWR%aDWR%nlTol
            !print*, 'linTol set to : ' , state%linSolver%tol
            !if ( DWR%aDWR%nlTol < DWR%aDWR%linTol ) &
            !   stop 'nlTol is under linTol in DWR nonlinear solve'

            ! SOLVE THE SYSTEM
            !print*, 'SOLVE PRIMAL problem !'
            call PrimalLinSolver( DWR, grid, newton, state%linSolver%tol )

         else
            stop 'performOneNewtonStdgmStep must be called with DWR, if Newton%non_alg_stop == aDWR '
         endif
       else
         ! parameters for the linear (GMRES) solver
         restart = 50 !30
         nloops =  10 !5 ! 10  ! HERE TO CHANGE 1
         !if(state%time%maxiter == 43) nloops = 5 ! 2 ! SIMULATION of  aDWR stopping criterion

         !if( state%modelName == 'porous') nloops = 1


         if(state%model%varying_time_term) eta = 0. ! time deriv terms already included
         if(state%nlSolver%pseudoT) eta =  state%nlSolver%pseudo_time%PTS_rho


         !print*,'dwr iter state%space%estim_space == ', state%space%estim_space, iter
         ! preparation of the right hand-side for dual problem
         ! if ( .not. present(DWR) ) stop 'DWR is not presented'

         ! BiCG technique according [Strakos, Tichy, SISC 2011]
         !   A x = b_0, initial approximation x_0   <==>    A x~ = b_0 -   A x_0 = b
         ! A^T y = c_0, initial approximation y_0   <==>  A^T y~ = c_0 - A^T y_0 = c
         !  x~ = Newton%x      y~ = DWR%x
         !  b  = Newton%b      c  = DWR%b
         if(BiCG) then
            if(iter == 1 .and. state%time%iter_loc == 1) then
               !print*,'DWR cleaning and preparing '
               call DWR%clean()

               !!call PrepareDualProblem( DWR )  - NOT NECESSARY
               !print*,'call DWR %J%findSupp( grid )'
               call DWR%J%findSupp( grid )

            endif

            p_mod = state%getP_mod()
            q_mod = state%getQ_mod()

            ! setting of the RHS for the dual problem
            ! RHS must be computed after ComputeSTDGM_Terms( )
            !print*,"call DWR%setRHS_new( grid )"
            call DWR%setRHS_new( grid )

            !print*,'call PrepareDualProblem_RHS( DWR, grid )'
            ! assembles dwr%rhs  -> DWR%b and
            !           elem%zST -> DWR%x
            call PrepareDualProblem_RHS( DWR, grid )

            bigNsize = state%getBigNSize()

            if(p_mod > 0) then
               print*,'state%getP_mod() > 0 = ', state%getP_mod()
               print*,'this subroutine has to be re-written, also the elem%wST GESEJ'
               stop
            endif

            allocate(x0(1:bigNsize), y0(1:bigNsize), b0(1:bigNsize), c0(1:bigNsize))
            allocate(d0(1:bigNsize), e0(1:bigNsize))

            !k = 6
            !do i=0,bigNsize, k
            !   write(31,'(a6, i5, 300es14.6)') 'DWRx:',i+1, DWR%x(i+1: min(bigNsize, i+k) )
            !enddo
            !DWR%x = 0.


            ! update of the RHS: A^T  z_0 (z_0 = zST)   DWR%rr is temporary array now
            call bMVprodBIG_Dual( DWR%rr(1: bigNsize), DWR%x(1:bigNsize) , bigNsize)

            ! new RHS for the dual problem b = b - A^T z_0  ( a la Newton update)
            ! is set several lines below

            ! assemple elem%w_ST -> Newton%x
            call FillVectorST_from_wST( Newton%x(1: bigNsize ) )

            ! A x_0   !  NOT NECESSARY for REAL computation
            !XXXxxx
            call bMVprodBIG( b0(1: bigNsize), Newton%x(1:bigNsize) , bigNsize)
            ! b0 = A x_0 + b   <==>  b = b0 - A x_0
            !XXXxxx
            b0(1:bigNsize) = b0(1:bigNsize) + Newton%b(1:bigNsize)
            c0(1:bigNsize) = DWR%b(1:bigNsize)
            !XXX

            ! init the value of the functional (see [Strakos, Tichy, SISC 2011]
            J1 = dot_product (Newton%x(1: bigNsize), DWR%b(1: bigNsize) )
            J2 = dot_product (DWR%x(1: bigNsize), Newton%b(1: bigNsize) )
            xi_initP = J1 + J2
            xi_P = J1 + J2
            !write(*, '(a30, 3es22.14)') &
            !     'INIT J(u) 1 (xi_i, c*x, y*b~ =', xi_initP,  J1, J2

            ! new RHS for the dual problem c = c - A^T z_0  ( a la Newton update)
            ! MUST be COMPATIBLE with the update of the dual solution ! VD DUAL BICG
            DWR%b(1:bigNsize) = DWR%b(1:bigNsize) - DWR%rr(1:bigNsize)

            !XXXxxx
            J1 = dot_product (Newton%x(1: bigNsize), DWR%b(1: bigNsize) )
            J2 = dot_product( DWR%x(1: bigNsize), b0(1: bigNsize) )
            xi_initD = J1 + J2
            xi_D = J1 + J2

            !write(*, '(a30, 3es22.14)') &
            !    'INIT J(u) 2 (xi_i, c~*x, y*b =', J1 + J2, J1, J2
            !XXXxxx

            !k = 6
            !do i=0,bigNsize-1, k
            !!do i=300,370-1, k
            !!do i=300, 306, k
            !   write(*,'(a6, i5, 300es14.6)') 'RHsu:',i+1, DWR%b(i+1: min(bigNsize, i+k) )
            !enddo

            !XXX Newton%rr(1: bigNsize ) = Newton%x(1: bigNsize )
            !XXX NEEDED FOR new estimates
            x0(1: bigNsize ) = Newton%x(1: bigNsize )  ! used in estP,estD in bicg.f90
            y0(1: bigNsize ) = DWR%x(1: bigNsize )

            ! k = 6
            ! do i=0,3*k, k
            !    !write(*,'(a10, 2i5, 300es14.6)') &
            !    write(*,'(a10, 2i5, 300es24.16)') &
            !         'init solP:',i+1, k, Newton%rr(i+1: min(bigNsize, i+k) )
            ! enddo
            ! print*
            !k = 6
            !do i=0,bigNsize, k
            !    !write(*,'(a10, 2i5, 300es14.6)') 'init solD:', &
            !   write(20+state%time%iter,'(a10, 2i5, 300es24.16)') 'init solD:', &
            !        i+1, k, x0(i+1: min(bigNsize, i+k) )
            !enddo



            ! initialization of updates
            Newton%x(1: bigNsize ) = 0.
            DWR%x(1:bigNsize) = 0.

         endif ! if(BiCG)



         ! SOLVE THE SYSTEM
         call state%cpuTime%startSolveTime()

         ! for Ben Southworth
         ! writing of the sparse block matrix
         !call WriteBlockLinearSTDGMsystem(state%nsize, eta, newton%b, newton%x, &
         !     newton%iter, state%time%iter)

         !call WriteMatrixA_ST(eta)
         !print*,'eta = ', eta, state%time%tau(1)
         !call WriteMblock_Screene(grid%elem(1)%blockST(0) )


         call state%setP_mod( 0 )
         call state%setQ_mod( 0 )
         bigNsize = state%getBigNSize()

         !write(*,'(a10, 200es12.4)') 'SOL b:', newton%b
         !stop

         !  simultaneous solution of the primal and dual problems  using of BiCG method
         !if (state%space%estim_space == 'DWR' .and.  DWR%BI) then
         if (BiCG ) then
            restart = 10 !30   ! DEFINES ALSO  update_iter = \NU
            !if(iter > 1) restart = 3
            !nloops = 100 !5 ! 10
            nloops = 80 !40 !1000  !40 ! 50 !5 ! 10
            !if(state%time%maxiter == 43) nloops = 10  ! SIMULATION of  aDWR stopping criterion
            
            
            call SolveBlockLinearBigPrimalDualProblem(bigNsize, &
                 eta, newton%b,  DWR%b, newton%x, DWR%x, &
                 state%linSolver%residuum, state%linSolver%iter, &
                 state%linSolver%lin_solver_not_conv, restart, nloops, xi_initP, xi_initD, &
                 x0(1:bigNsize), y0(1:bigNsize), b0(1:bigNsize), c0(1:bigNsize))
            
            !print*,'after  sx:',norm2( newton%x), norm2( DWR%x)
            !k = 6
            !do i=0,bigNsize, k
            !!    !write(*,'(a10, 2i5, 300es14.6)') 'init solD:', &
            !   write(30+state%time%iter,'(a10, 2i5, 300es24.16)') 'init solD:', &
            !        i+1, k, newton%x(i+1: min(bigNsize, i+k) )
            !enddo
            !k = 6
            !do i=0,bigNsize, k
            !!    !write(*,'(a10, 2i5, 300es14.6)') 'init solD:', &
            !   write(40+state%time%iter,'(a10, 2i5, 300es24.16)') 'init solD:', &
            !        i+1, k, x0(i+1: min(bigNsize, i+k) )
            !enddo
            
            ! writting of the matrices and vectors
            !if(mod(state%space%adapt%adapt_level, 4) == 0) &
            !     call WriteBiCG_Matlab_All(bigNsize, state%space%adapt%adapt_level, &
            !     newton%b,  DWR%b, newton%x,  DWR%x )
            
            !!call WriteBIGmatrixMatlab_Naive(bigNsize)
            !!call WriteBIGmatrixMatlab_Simple(bigNsize)
            !call WriteBIGmatrixMatlab(bigNsize)  ! faster but not preconditioner
            !call WriteBIGVectorMatlab(bigNsize, newton%b,  DWR%b, "rhs_b", "rhs_c")
            !call WriteBIGVectorMatlab(bigNsize, newton%x,  DWR%x, "sol_b", "sol_c")
            
            
            !print*, '(2) Size and norm of the solution with bigBlock x = ', size(newton%x) , norm2(newton%x)
            ! output quantity resulting from BiCG solver
            DWR%J%Ju_BiCG = xi_initP
            
            !x0(1:bigNsize) = x0(1:bigNsize) + Newton%x(1:bigNsize)
            !y0(1:bigNsize) = y0(1:bigNsize) + DWR%x(1:bigNsize)
            
            ! ! k = 6
            ! ! do i=0,3*k, k
            ! !    !write(*,'(a10, 2i5, 300es14.6)') &
            ! !    write(*,'(a10, 2i5, 300es24.16)') &
            ! !         'solP:',i+1, k, x0(i+1: min(bigNsize, i+k) )
            ! ! enddo
            ! ! print*
            ! ! do i=0,3*k, k
            ! !    !write(*,'(a10, 2i5, 300es14.6)') &
            ! !    write(*,'(a10, 2i5, 300es24.16)') &
            ! !         'solD:', i+1, k, y0(i+1: min(bigNsize, i+k) )
            ! ! enddo
            
            !print*
            !J1 = dot_product(c0(1:bigNsize), x0(1:bigNsize) )
            !J2 = dot_product(b0(1:bigNsize), y0(1:bigNsize) )
            !print*,'output quantities, J1, J2, aver = ', J1, J2, (J1+J2)/2
            
            !J1a = xi_P + dot_product(DWR%b(1:bigNsize), Newton%x(1:bigNsize) )
            !J2a = xi_D + dot_product(Newton%b(1:bigNsize), DWR%x(1:bigNsize) )
            !print*,'output quant loc , J1, J2, aver = ', J1a, J2a
            !print*,'____________________________________________________________'
            
            ! call bMVprodBIG( d0(1: bigNsize), Newton%x(1:bigNsize) , bigNsize)
            ! d0(1: bigNsize) = d0(1: bigNsize) - Newton%b(1:bigNsize)
            
            ! call bMVprodBIG_Dual( e0(1: bigNsize), DWR%x(1:bigNsize) , bigNsize)
            ! e0(1: bigNsize) = e0(1: bigNsize) - DWR%b(1:bigNsize)
            
            ! print*,'difference , J1, J2, aver = ', &
            !      dot_product( d0(1: bigNsize) , y0(1:bigNsize) ), &
            !      dot_product( e0(1: bigNsize) , x0(1:bigNsize) ), &
            !      J1 - J1a - dot_product( d0(1: bigNsize) , y0(1:bigNsize) ), &
            !      J2 - J2a - dot_product( e0(1: bigNsize) , x0(1:bigNsize) )
            
            
            ! print*,'____________________________________________________________'
            
            deallocate(x0, y0, b0, c0, d0, e0)

         else

            call SolveBlockLinearBigProblem(bigNsize, &
                 eta, newton%b, newton%x, &
                 state%linSolver%residuum, & ! state%linSolver%iter, &
                 state%linSolver%lin_solver_not_conv, restart, nloops)
            !print*, '(1) Size and norm of the solution with bigBlock x = ', size(newton%x) , norm2(newton%x)
         endif


         !k = 6
         !do i=0,bigNsize-1, k
         !   write(97,'(a6, i5, 300es14.6)') 'rhss:',i+1, newton%b(i+1: min(bigNsize, i+k) )
         !   write(98,'(a6, i5, 300es14.6)') 'sols:',i+1, newton%x(i+1: min(bigNsize, i+k) )
         !enddo


         !write(*,'(a10, 200es12.4)') 'SOL x:', newton%x

!         state%p_mod = 1
!         bigNsize = state%getBigNSize()
!
!         allocate( newtonX(1:bigNsize), source = 0.0 )
!         allocate( newtonB(1:bigNsize), source = 0.0 )
!         call FillBigVector( newtonB(1:bigNsize), eta ) ! fill RHS
!         !print*, 'newton RHS = ', newtonB
!
!         call SolveBlockLinearBigProblem( bigNSize, eta, newtonB, newtonX, &
!            state%linSolver%residuum, state%linSolver%iter, &
!            state%linSolver%lin_solver_not_conv, restart, nloops)
!
!         print*, 'Size and norm of the solution with bigBlock x = ', size(newtonX) , norm2(newtonX)
!         deallocate( newtonB, newtonX )

         !call bMVprodST(newton%rr, newton%x, state%nsize)    ! rr = Ax
         !print*, 'norm of the residual = ', norm2(newton%rr)
         call bMVprodBIG( newton%rr, newton%x, state%nsize)    ! rr = Ax
         !print*, 'norm of the residual1 = ', norm2(newton%rr)

         newton%rr(:) = newton%rr(:) - newton%b(:)
         !print*, 'norm of the residual2 = ', norm2(newton%rr)

         call state%cpuTime%addSolveTime()

       endif

       call state%cpuTime%startPrepareTime()

       ! some manipulation for the fail of the computation
       if(state%linSolver%lin_solver_not_conv > 0) then
          open(53, file='GMRES_failed', status='UNKNOWN', position='append')
          res_max_val = 0.
          do i = 1, grid%nelem
             elem => grid%elem(i)

             kk = elem%ncv - 1
             elemDof = elem%dof * ndim * elem%Tdof
             elem%vec( res_vec, 1:elemDof) =  newton%rr(kk+1:kk+ elemDof)
             res_max_val = max( res_max_val, &
                  sqrt( dot_product(elem%vec( res_vec, 1:elemDof), elem%vec( res_vec, 1:elemDof) )))
          enddo
          ! FR added
          close(53)
       endif

       call state%cpuTime%addPrepareTime()

   end subroutine performOneNewtonStdgmStep

   !> write the outputs of the Newton method into the file 'ifile' -- 'criter'
   subroutine Newton_output(grid, ifile, l, loc_implicitly )
    class( mesh ), intent(in) :: grid
    integer, intent(in) :: ifile
    integer, intent(in) :: l
    logical, intent(in) :: loc_implicitly
    real :: t2

    associate ( Newton => state%nlSolver )

      call cpu_time(t2)


      write(ifile,'(3i4,14es13.5, i10,2l3,F10.1,i9,i11, es12.4, 4i8, 4es12.4, i8)') &
           state%time%iter,Newton%iter,l, &                             ! 1..3
           state%time%iter + 1.*(Newton%iter-1)/(Newton%max_iter+2), &      ! 4
           state%estim(resA_ST,1), &                                       ! 5
           Newton%res/state%nsize,  Newton%res0/state%nsize, &             ! 6..7
           Newton%res,  Newton%res0, &                                     ! 8..9
           state%linSolver%tol,  Newton%lambda, Newton%theta, &            !10 .. 12
           state%L_estim(resA:resST), state%time%tau(1), &                 !13 .. 17
           state%linSolver%iter, loc_implicitly, &                         !18 .. 19
           state%linSolver%precond_update, t2- state%start_time, &         !20 .. 21
           Newton%Aiter, state%linSolver%iter_tot, &                       !22 .. 23
           state%time%ttime+ state%time%tau(1),  &   !24
           state%num_call_C, state%num_call_F, grid%nelem, state%nsize, &  !25 .. 28
           state%L_estim(resA) / max(1E-15, state%L_estim(resS) ),   &     ! 29
           state%L_estim(resA) / max(1E-15, state%L_estim(resT) ),   &     ! 30
           state%L_estim(resT) / max(1E-15, state%L_estim(resS) ), &       ! 31
           state%linSolver%residuum, state%linSolver%iter                  ! 32..33
           end associate

   end subroutine Newton_output

   !> Estimate the nonlinear error in Newton for aRES
   !> it should be called only when Newton%non_alg_stop == 'aRES'
   subroutine EstimateNewtonARes( Newton, time_prepare, time_estim, done, estS)
      class( Newton_t ), intent ( in ) :: Newton
      real, intent(inout) :: time_prepare
      real, intent(inout) :: time_estim
      logical, intent(inout) :: done
      real, intent(inout) :: estS
      class(element), pointer :: elem
      integer :: i
      real :: t1, t2, estim

      call cpu_time(t2)
      time_prepare = time_prepare + t2 - t1

      !print*, 'PerformOneSTDGMstep calling RezidErrorEstimate(true)'
      !call RezidErrorEstimates( .true., .false. )  !onlyAS =.true.  -- almost the same
      call RezidErrorEstimates( .false., .false. )

      call cpu_time(t1)
      time_estim = time_estim + t1 - t2
      !write(*,'(a6, 3es12.4)') 'timesD',time_prepare, time_solve , time_estim

      ! if(Newton%iter == 1) &
      !write(*,'(a4,a4,5a12,a7,2a12)') '  ', 'iter', '|F(x^k)| ','|F(x^k+1)|', &
      !      'eta_A', 'eta_S', &
      !      'rel. err. ', '  LA it',&
      !      'LA resid', '   CPU(s) '

      ! write(*,'(a4,2i4,5es12.4,a2,i5,2es12.4,l2)') 'NEW:', Newton%iter, l-1, &
      !       Newton%res,Newton%res0, &
      !       state%L_estim(resA), state%L_estim(resST), &
      !       state%L_estim(resA)/ state%L_estim(resST),' |', &
      !       state%linSolver%iter, state%linSolver%residuum, t1 -  state%start_time,&
      !       loc_implicitly

      if(.not. state%time_dependent) then
         if(state%estim(resA_ST,1) < Newton%tol2 ) then
            !if(state%estim(resA_ST_loc,1) < Newton%tol2 ) then  ! HERE TO CHANGE
            done = .true.

            !goto 200  ! Newton method has achieved the prescribed tolerance
         endif

      else !  state%time_dependent      ! HERE TO CHANGE 

         !write(*,'(a8, 40es12.4)') &
         !     'STG:',state%L_estim(resA),state%L_estim(resS),state%L_estim(resT),&
         !     state%estim(resA_ST,1) , Newton%tol2, &
         !      state%L_estim(resA) /  max(1E-15, state%L_estim(resT))
         
         ! OLD VARIANT
         !if(state%estim(resA_ST,1) < Newton%tol2  ) then 
         !.and. state%L_estim(resA) /  max(1E-15, state%L_estim(resT)) < 0.1) then

         ! NEW VARIANT
         estim = state%L_estim(resA) /max(1E-15, minval(state%L_estim(resS:resT)) )
         if( estim < Newton%tol2  ) then 
            !if(state%estim(resA_ST_loc,1) < Newton%tol2 ) then
            done = .true.
            if(state%time%iter <=1) print*,'NEW aRES criterion'
            !goto 200  ! Newton method has achieved the prescribed tolerance
         endif
      endif

      estS = state%L_estim(resS)

      ! associate ( Newton => state%nlSolver )
      !   if(Newton%iter == 1) then
      !      open(103, file = 'criter_loc', status='replace')
      !   else
      !      open(103, file = 'criter_loc', status='UNKNOWN', position='append')
      !   endif

      !   do i=1,grid%nelem
      !      elem => grid%elem(i)
      !      if( elem%eta(resA, 1)/ elem%eta(resS, 1) > 1E-3 .and. elem%eta(resA, 1)> 1E-3 ) then
      !           write(103, *) elem%i + 1.*(Newton%iter-1)/(Newton%max_iter+2), &
      !           elem%eta(resA, 1), elem%eta(resS, 1), elem%eta(resT, 1)
      !           write(1000+elem%i, *) elem%i + 1.*(Newton%iter-1)/(Newton%max_iter+2), &
      !           elem%eta(resA, 1), elem%eta(resS, 1), elem%eta(resT, 1), '76hhju73'
      !        endif
      !   enddo

      !   close(103)

      ! end associate


   end subroutine EstimateNewtonARes


   !> Estimate the nonlinear error in Newton for aDWR
   !> AND compute the dual problem, then compute nl estimate again
   !> it should be called only when Newton%non_alg_stop == 'aDWR'
   subroutine EstimateNewtonADwr( grid, DWR, iter, newtonDone)
      class( mesh ), intent(inout) :: grid
      type( DWR_t), intent(inout) :: DWR
      integer, intent(in) :: iter
      logical, intent(out) :: newtonDone
      real :: etaD, some_criter, t2, t1, dualTol

      ! a) estimate nlError for the 1st time
      ! if we are the 1st mesh then no ZST is computed
      if ( grid%adapt_level > 0 .or. iter > 1 ) then
         call computeNonlinDWRestimates( DWR, grid )

         if ( DWR%estimNL < state%nlSolver%tol ) then
            newtonDone = .true.
         else
            newtonDone = .false.
         endif

         print*, 'After ',iter,'th nonlin est! estNL <? TOL' , DWR%estimNL , state%nlSolver%tol
      else
         newtonDone = .true.
      endif

      if (iter >= state%nlSolver%max_iter) then
         newtonDone = .true.
         print*, 'Maximal number of Newton iters achieved -> compute dual problem!'
      endif

      ! update the nonlinear estimate to see if we really can end Newton iterations !
      if (newtonDone) then
            ! b) Dual problem
            call UpdateDualProblem( DWR )
            ! solve the dual problem with tol = C_A * estimNL
            dualTol = DWR%aDWR%linTol * DWR%aDWR%C_Safe
            call DualLinSolver( DWR, grid, iter, dualTol)

            !DWR%aDWR%iter_lin_dual =  -10  !DWR%linSolver_iter !DWR%aDWR%iter_lin_dual  +
            !call PlotSolDual( DWR%aDWR%iter_lin_dual )
            call PlotSolDual( 0 )
            ! too much memory
            !call PlotSolDual( iter )

            ! c) estimate nlError for the 2nd time
            call computeNonlinDWRestimates( DWR, grid )

            if ( DWR%estimNL < state%nlSolver%tol ) then
               newtonDone = .true.
            else
               newtonDone = .false.
            endif
            !print*, 'After second nonlin est! estNL <? TOL' ,  DWR%estimNL , state%nlSolver%tol
      end if

   end subroutine EstimateNewtonADwr

   !> Estimate the nonlinear error in Newton for aDWR
   !> AND compute the dual problem, then compute nl estimate again
   !> it should be called only when Newton%non_alg_stop == 'aDWR' and
   !> fixedIter - fixed number of lin iters -> decrease graphs !!!
   !> TESTED only for linear problems
   subroutine EstimateNewtonADwrFixed( grid, DWR, iter, newtonDone)
      class( mesh ), intent(inout) :: grid
      type( DWR_t), intent(inout) :: DWR
      integer, intent(in) :: iter
      logical, intent(out) :: newtonDone
      real :: etaD, some_criter, t2, t1, dualTol

      print*, 'EstimateNewtonADwrFixed called'


      ! b) Dual problem
      call UpdateDualProblem( DWR )
            ! solve the dual problem with tol = C_A * estimNL
      dualTol = 1.E-12 ! TOL IS NOT USED DWR%aDWR%linTol * DWR%aDWR%C_Safe
      call DualLinSolver( DWR, grid, iter, dualTol)

            !DWR%aDWR%iter_lin_dual =  -10  !DWR%linSolver_iter !DWR%aDWR%iter_lin_dual  +
            !call PlotSolDual( DWR%aDWR%iter_lin_dual )
            !call PlotSolDual( 0 )
            ! too much memory
            !call PlotSolDual( iter )

            ! c) estimate nlError for the 2nd time
      call computeNonlinDWRestimates( DWR, grid )

      print*, 'TOlerance: ', state%nlSolver%tol
      print*, 'Primal :' ,DWR%estimNL , DWR%estimLP
      print*, 'Dual:' , DWR%estimLD
      print*

      newtonDone = .true.
      print*, 'Always NEWTON DONE !!!'

   end subroutine EstimateNewtonADwrFixed



   !> solve the Newton problem F(w) = 0
   !> with tolerance given in state%nlSolver%tol with respect to one of the
   !> estimating methods: aRES, aDWR, rezL2
   subroutine NewtonSolve( grid, imp, iter, deg_plus, DWR)
      class( mesh ), intent(inout) :: grid
      integer, intent(inout) :: imp ! integer which is used to define implicitly - Should be changed
      integer, intent(inout) :: iter ! in - previous iterations , out - total iterations
      logical, intent(in) :: deg_plus
      type( DWR_t), intent(inout), optional :: DWR
      class(element) , pointer :: elem
      logical :: newtonDone, PTSdone  ! Newton and pseudo-time stepping finished
      real :: time_prepare, time_estim, t1, t2
      integer :: i, nDamping, maxNumberOfDamping
      real :: ratio_Damping
      integer :: iPTS, nPTS, iter_total, nsize
      real :: res_max_val, time, lost
      real :: estS, estS_ini, fac
      logical :: loc_implicitly, bicg_used

      bicg_used = .false.    ! use of the bicg method for the primal and dual solutions

!      print*, "WORKING ON NewtonSolve!"

      ! factor for divergence, maximal increase of the initial residuum
      fac = 100. ! can be small, decrease of the time step will follow
      if (state%time%tau_fixed) fac = 1E+08

      
      ! initialization parameteres for the setting of divergence
      estS = 1.
      estS_ini = 1.

      !nsize = state%nsize
      ! FR_NEWTON
      nsize = state%getBigNSize()

      iter_total = 0   ! sum of all Newton iterations over all pseudo-time steps
      newtonDone = .false.
      PTSdone = .false.

      associate( nlSolver=> state%nlSolver)

        nlSolver%lambda1 = 1.0   ! initialization of lambda (= damping factor)

        !print*,'nlSolver%anderson ===' , nlSolver%anderson
        ! Anderson acceleration - no standard damping is used, only update the solution
        if ( nlSolver%anderson ) then
           ! DRIV: allocate( AndersonAcceleration_t :: nlSolver%andersonAcc )
           nlSolver%andersonAcc%AA_max = 20

           call nlSolver%andersonAcc%init( nsize )
           maxNumberOfDamping = 1
           ratio_Damping = 0.85
           !maxNumberOfDamping = 3
        else
           !maxNumberOfDamping = 5
           !ratio_Damping = 0.5

           maxNumberOfDamping = 9
           ratio_Damping = 0.75
          !maxNumberOfDamping = 15
           if(nlSolver%max_iter == 1) maxNumberOfDamping = 1 ! typically for linear problem
        end if

        ! pseudo-time stepping
        nPTS = 1 ! ONLY one-repetition of the Newton step
        iPTS = 1 ! index of the loop

        ! pseudo-time stepping
        if( nlSolver%pseudoT) then
           nPTS = 20  ! maximal number of the pseudo-time step repetitions
           call nlSolver%pseudo_time%init(nsize )
           nlSolver%pseudo_time%nPTS = nPTS
           nlSolver%pseudo_time%PTS_rho = 10000.
        endif

        ! MAIN CYCLE for the pseudo-time stepping
        do while ( (.not. PTSdone) .and. (iPTS <= nPTS) )
           iter = 1


           if( nlSolver%pseudoT) then
              call CopyWST_toLongVector( grid, nsize, nlSolver%pseudo_time%PTS_x(1, 1:nsize) )

              !write(*,'(a25, 2i5,l3, 30es12.4)') '## pseudo-time stepping:', iPTS, nPTS, PTSdone, &
              !     sqrt(dot_product(nlSolver%pseudo_time%PTS_x(1,:), nlSolver%pseudo_time%PTS_x(1,:))), &
              !     nlSolver%pseudo_time%PTS_x(1, 1:5)
           endif



           ! MAIN NEWTON CYCLE
           do while ( (.not. newtonDone) .and. (iter <= nlSolver%max_iter) )
              ! Three parts:
              !    1) solve the linear system
              !    2) find the damping
              !    3) estimate the error and decide whether to quit - newtonDone = .true.


              if( nlSolver%pseudoT) then
                 call CopyWST_toLongVector( grid, nsize, nlSolver%pseudo_time%PTS_x(2, 1:nsize) )
                 nlSolver%pseudo_time%PTS_x(2, 1:nsize) = nlSolver%pseudo_time%PTS_x(2, 1:nsize) &
                      - nlSolver%pseudo_time%PTS_x(1, 1:nsize)
              endif


              !print*," 1) SOLVE THE LINEAR SYSTEM"
              if(state%space%estim_space == 'DWR' .and. present(DWR) ) then
                 if(nlSolver%non_alg_stop == 'aDWR' .or. DWR%BI )  then

                    if( DWR%BI ) bicg_used = .true.

                    !if(state%space%estim_space == 'DWR' .and. &
                    !     (nlSolver%non_alg_stop == 'aDWR' .or. DWR%BI ) ) then
                    !if (present(DWR)) then
                    ! inner linear primal problem step
                    ! print*, ' performOneNewtonStdgmStep called '
                    call performOneNewtonStdgmStep( nlSolver, grid, imp,&
                         iter, deg_plus, newtonDone, res_max_val, loc_implicitly, DWR )
                 endif
              elseif(nlSolver%non_alg_stop == 'aDWR') then
                 stop 'performOneNewtonStdgmStep should be called with DWR in its &
                      argument for DWR method'

              else
!                 print*,'call performOneNewtonStdgmStep BB'
                 call performOneNewtonStdgmStep( nlSolver, grid, imp, iter, &
                      deg_plus, newtonDone, res_max_val, loc_implicitly )

                 !print*, "Norm of the res = ", norm2(nlSolver%rr(:))
              endif

              !print*,' bicg_used ============' ,bicg_used

              !print*," 2) seeking of optimal damping parameter",nlSolver%anderson
              call state%cpuTime%startSolveTime()

              select type ( nlSolver )
              type is ( Newton_t )
                 ! aDWR  !
                 if ( nlSolver%non_alg_stop == 'aDWR') then
                    ! AND zST ready
                    if (DWR%dualProblem_computed) then
                       call setDampingInNewton_DWR( nlSolver, grid, DWR, deg_plus, &
                            eta, loc_implicitly, res_max_val, maxNumberOfDamping, nDamping)
                       !   print*, 'temporarily classical damping back!'
                       !   call setDampingInNewton( nlSolver, grid, deg_plus, &
                       !       eta, loc_implicitly, res_max_val, maxNumberOfDamping, nDamping)

                    else
                       print*, 'zST is not ready yet. Classical RES damping is called!'
                       call setDampingInNewton( nlSolver, grid, deg_plus, &
                            eta, loc_implicitly, res_max_val, maxNumberOfDamping, nDamping, ratio_Damping)

                    endif
                 else
                    !print*,'############################################################ 36gd3'
                    !call setDampingInNewtonLocal( nlSolver, grid, deg_plus, &
                    !   eta, loc_implicitly, res_max_val, maxNumberOfDamping, nDamping)
                    if(bicg_used) then
                       !maxNumberOfDamping = 1
                       call setDampingInNewton( nlSolver, grid, deg_plus, &
                            eta, loc_implicitly, res_max_val, maxNumberOfDamping, nDamping,&
                            ratio_Damping, DWR)
                    else
                       call setDampingInNewton( nlSolver, grid, deg_plus, &
                            eta, loc_implicitly, res_max_val, maxNumberOfDamping, nDamping,&
                            ratio_Damping )
                    endif

                 end if

                 ! Anderson Acceleration, 3 variants inside
                 if (nlSolver%anderson ) &
                      call setAndersonInNewton( nlSolver, grid, deg_plus, eta, loc_implicitly)
                 !call setAndersonInNewton_OLD( nlSolver, grid, deg_plus, eta, loc_implicitly)

                 class default
                 stop 'other type of nlSolver'
              end select

              call state%cpuTime%addSolveTime()
              call state%cpuTime%startEstimTime()

              if(nlSolver%newton_count >= 8) nlSolver%implicitly = .true.
              ! !!if(Newton%theta > 0.75) nlSolver%implicitly = .true.
              if(nlSolver%theta > 0.5) nlSolver%implicitly = .true.
              !if(Newton%lambda < 0.9) nlSolver%implicitly = .true. !lambda is small, C(w) update
              if(state%linSolver%lin_solver_not_conv > 0) nlSolver%implicitly = .true.

              !if(nlSolver%implicitly ) write(*,'(a8,i5, a14,i5,a8,es12.4,a8,l2 )')&
              !     'iter = ',state%time%iter, '### Ncount =',nlSolver%newton_count, &
              !     ',  theta = ', nlSolver%theta,',  non conv', state%linSolver%lin_solver_not_conv
              !					state%isol = state%isol + 1
              !					call WriteProgressOutput( 'STE' )
              !call WriteMatrixLU()

              !write(*,'(a8, 2i5, 30es12.4)') 'PTS rez:',iPTS, iter, &
              !     sqrt(dot_product(nlSolver%b1, nlSolver%b1)), nlSolver%res0

              !open(15, file = 'pseudo', status = 'UNKNOWN', position = 'append')
              !if(iter == 1) write(15, '(x)')
              !write(15,'(2i5, 30es12.4)') iPTS, iter, &
              !     1.*( (iPTS-1)* nlSolver%max_iter + iter) / ( nPTS * nlSolver%max_iter + 5), &
              !     sqrt(dot_product(nlSolver%b1, nlSolver%b1))
              !close(15)

              !print*," ! 3 ) ESTIMATE THE ALGEBRACIC (NONLINEAR) ERROR"
              select type ( nlSolver )
              type is ( Newton_t )

                 if( bicg_used ) then   ! BiCG method has its own stopping criterion
                    if( .not. state%model%linear) &
                         stop " BiCG method only for LINEAR problems now!"

                    !call EstimateNewton_DWR( grid, nlSolver, DWR, newtonDone )
                    print*,'Newton-BICG stopping criterion', newtonDone
                    ! one Newton step should be enough
                    !newtonDone = .true.

                 ! nonlinear algebraic residuum criterion
                 elseif(nlSolver%non_alg_stop == 'aRES') then
                    call EstimateNewtonARes( nlSolver, time_prepare, time_estim, newtonDone, estS )
                    if(iter == 1) estS_ini = estS
                    !if(estS > 1E+03 * estS_ini) newtonDone = .true. ! DIVERGENCE

                    if(estS > fac * estS_ini)  newtonDone = .true. ! DIVERGENCE
                    
                        !    if( elem%i == 35) then
                    !       write(100, '(2i5, 300es12.4)' ) i, iter, elem%eta(resA:resST, 1), 99999., &
                    !            elem%wST(1,1:elem%dof,1:elem%Tdof), -99999., &
                    !            elem%rhsST(1,1:elem%dof_plus,1:elem%Tdof_plus)
                    !    endif

                    ! enddo


                 else if (nlSolver%non_alg_stop == 'aDWR') then
                    print*,'aDWR 1'
                    call DWR%J%computeJu(grid)
                    write(debug, *) 'Maybe we should compute Ju_exact - to outputfile'
                    print*,'aDWR 2'
                    call DWR%J%computeJu_exact(grid)
                    print*,'aDWR 3'
                    call DWR%writeNlErrorFile(grid%nelem)
                    print*,'aDWR 4', DWR%aDWR%fixedIter 

                    if ( DWR%aDWR%fixedIter ) then ! fixed number of lin iters - decrease graphs
                       call EstimateNewtonADwrFixed( grid, DWR, iter, newtonDone)
                    else
                       call EstimateNewtonADwr( grid, DWR, iter, newtonDone)
                    endif

                 else if (nlSolver%non_alg_stop == 'rezL2') then
                    write(*,*) 'FR nolinear rezL2 reziduum:' , &
                      nlSolver%res0,  '<' , nlSolver%tol !* nsize

                    !if(nlSolver%res0 < nlSolver%tol * nsize )  then
                    if(nlSolver%res0 < nlSolver%tol )  then
                       !write(*,*) 'Reziduum:' , nlSolver%res0,  '<' , nlSolver%tol * nsize
                       !call cpu_time(t2)
                       !time_prepare = time_prepare + t2 - t1
                       ! Newton method has achieved the prescribed tolerance
                       newtonDone = .true.
                    endif
                    !write(*,*) 'Reziduum:' , nlSolver%res0,  '>' , nlSolver%tol * nsize

                 else
                    stop 'Unknown method for nonlinear algebraic error estimation in NewtonSolve!'
                 endif  ! end of other technique
                 class default
                 stop 'other type of nlSolver'
              end select

              call state%cpuTime%addEstimTime( )

              ! Write Newton output
              open(91, file='criter', status='unknown', position='append')
              if(nlSolver%iter == 1) write(91,'(x)')
              if(nlSolver%Aiter == 1) write(91,'(x)')
              !if(state%time%iter == 1) print*, 'TODO: Repair NewtonOutput for other methods than aRES!'
              call Newton_output(grid, 91, nDamping, loc_implicitly )
              close(91)

              nlSolver%res1 = nlSolver%res

              iter = iter + 1
           end do ! while
           ! end of newton step, save the solution, write down and  go on



           iter_total = iter_total + iter ! sum of all Newton iterations over all pseudo-time steps
           iPTS = iPTS + 1

           if( nlSolver%pseudoT)  nlSolver%pseudo_time%PTS_rho = nlSolver%pseudo_time%PTS_rho / 10.
        end do  ! while
        ! end of the pseudotime stepping

        iter = iter_total

        ! VD no necessary printing, can be deduced from the number of used
        ! NEWTON ITERATIONS printed on the screen
        !if (( iter > nlSolver%max_iter) .and. (.not. newtonDone) ) then
        !   print*, 'Newton solver ended without reaching the tolerance', &
        !           'achieving max. number of iterations:' , nlSolver%max_iter
        !endif

        if ( nlSolver%anderson ) then
           call nlSolver%andersonAcc%dealloc()
        endif

        nlSolver%converged = newtonDone  ! PCQ (23. 8. 2019)
        !nlSolver%converged = .true.

        fac = 10
        if (state%time%tau_fixed) fac = 1E+08
        if(estS > fac * estS_ini) nlSolver%converged = .false. ! DIVERGENCE

!        print*,'#>MNH', fac, estS , fac * estS_ini,  nlSolver%converged 
      end associate ! nlSolver

!      print*,' end subroutine NewtonSolve'
   end subroutine NewtonSolve


!    !> solve the Newton problem F(w) = 0
!   !> with tolerance given in state%nlSolver%tol with respect to one of the
!   !> estimating methods: aRES, aDWR, rezL2
!   subroutine NewtonSolveWSTplus( grid, imp, iter, DWR)
!      class( mesh ), intent(inout) :: grid
!      integer, intent(inout) :: imp ! integer which is used to define implicitly - Should be changed
!      integer, intent(inout) :: iter ! in - previous iterations , out - total iterations
!      type( DWR_t), intent(inout) :: DWR
!      class(element) , pointer :: elem
!      logical :: newtonDone, PTSdone  ! Newton and pseudo-time stepping finished
!      real :: time_prepare, time_estim, t1, t2
!      integer :: i, nDamping, maxNumberOfDamping
!      real :: ratio_Damping
!      integer :: iPTS, nPTS, iter_total, nsize
!      real :: res_max_val, time, lost
!      real :: estS, estS_ini
!      logical :: loc_implicitly
!
!      ! initialization parameteres for the setting of divergence
!      estS = 1.
!      estS_ini = 1.
!
!      !nsize = state%nsize
!
!      iter_total = 0   ! sum of all Newton iterations over all pseudo-time steps
!      newtonDone = .false.
!      PTSdone = .false.
!
!      associate( nlSolver=> state%nlSolver)
!
!        nlSolver%lambda1 = 1.0   ! initialization of lambda (= damping factor)
!        !maxNumberOfDamping = 5
!        !ratio_Damping = 0.5
!        maxNumberOfDamping = 9
!        ratio_Damping = 0.75
!        !maxNumberOfDamping = 15
!
!           ! MAIN NEWTON CYCLE
!           do while ( (.not. newtonDone) .and. (iter <= nlSolver%max_iter) )
!              ! Three parts:
!              !    1) solve the linear system
!              !    2) find the damping
!              !    3) estimate the error and decide whether to quit - newtonDone = .true.
!
!              !1) SOLVE THE LINEAR SYSTEM"
!              if(nlSolver%non_alg_stop == 'aDWR') then
!                stop "aDWR not done for DWR_Ps"
!!                    ! inner linear primal problem step
!!                    ! print*, ' performOneNewtonStdgmStep called '
!!                call performOneNewtonStdgmStep( nlSolver, grid, imp,&
!!                         iter, deg_plus, newtonDone, res_max_val, loc_implicitly, DWR )
!
!              else
!                 ! FR_NEWTON
!                 ! TODO
!                 call performOneNewtonStdgmStep( nlSolver, grid, imp, iter, &
!                      deg_plus, newtonDone, res_max_val, loc_implicitly )
!              endif
!
!              ! 2) seeking of optimal damping parameter"
!              call state%cpuTime%startSolveTime()
!
!              select type ( nlSolver )
!              type is ( Newton_t )
!                 ! aDWR  !
!                 if ( nlSolver%non_alg_stop == 'aDWR') then
!                    stop "aDWR not done for DWR_Ps"
!!                    ! AND zST ready
!!                    if (DWR%dualProblem_computed) then
!!                       call setDampingInNewton_DWR( nlSolver, grid, DWR, deg_plus, &
!!                            eta, loc_implicitly, res_max_val, maxNumberOfDamping, nDamping)
!!                       !   print*, 'temporarily classical damping back!'
!!                       !   call setDampingInNewton( nlSolver, grid, deg_plus, &
!!                       !       eta, loc_implicitly, res_max_val, maxNumberOfDamping, nDamping)
!!
!!                    else
!!                       print*, 'zST is not ready yet. Classical RES damping is called!'
!!                       call setDampingInNewton( nlSolver, grid, deg_plus, &
!!                            eta, loc_implicitly, res_max_val, maxNumberOfDamping, nDamping, ratio_Damping)
!!
!!                    endif
!                 else
!                    ! FR_NEWTON
!                    ! TODO
!                    call setDampingInNewton( nlSolver, grid, deg_plus, &
!                         eta, loc_implicitly, res_max_val, maxNumberOfDamping, nDamping, ratio_Damping)
!                 end if
!
!                 class default
!                    stop 'other type of nlSolver'
!              end select
!
!              call state%cpuTime%addSolveTime()
!              call state%cpuTime%startEstimTime()
!
!              if(nlSolver%newton_count >= 8) nlSolver%implicitly = .true.
!              ! !!if(Newton%theta > 0.75) nlSolver%implicitly = .true.
!              if(nlSolver%theta > 0.5) nlSolver%implicitly = .true.
!              !if(Newton%lambda < 0.9) nlSolver%implicitly = .true. !lambda is small, C(w) update
!              if(state%linSolver%lin_solver_not_conv > 0) nlSolver%implicitly = .true.
!
!              !print*," ! 3 ) ESTIMATE THE ALGEBRACIC (NONLINEAR) ERROR"
!              select type ( nlSolver )
!              type is ( Newton_t )
!
!                 ! nonlinear algebraic residuum criterion
!                 if(nlSolver%non_alg_stop == 'aRES') then
!                    call EstimateNewtonARes( nlSolver, time_prepare, time_estim, newtonDone, estS )
!                    if(iter == 1) estS_ini = estS
!                    if(estS > 1E+03 * estS_ini) newtonDone = .true. ! DIVERGENCE
!                    ! do i=1,grid%nelem
!                    !    elem => grid%elem(i)
!                    !    write(100+iter, '(i5, 300es12.4)' ) i, elem%eta(resA:resST, 1), elem%w(0,:)
!                    !    if( elem%i == 35) then
!                    !       write(100, '(2i5, 300es12.4)' ) i, iter, elem%eta(resA:resST, 1), 99999., &
!                    !            elem%wST(1,1:elem%dof,1:elem%Tdof), -99999., &
!                    !            elem%rhsST(1,1:elem%dof_plus,1:elem%Tdof_plus)
!                    !    endif
!
!                    ! enddo
!
!
!                 else if (nlSolver%non_alg_stop == 'aDWR') then
!                    stop "aDWR not done for DWR_Ps"
!                    call DWR%J%computeJu(grid)
!                    write(debug, *) 'Maybe we should compute Ju_exact - to outputfile'
!                    call DWR%J%computeJu_exact(grid)
!                    call DWR%writeNlErrorFile(grid%nelem)
!
!                    if ( DWR%aDWR%fixedIter ) then ! fixed number of lin iters - decrease graphs
!                       call EstimateNewtonADwrFixed( grid, DWR, iter, newtonDone)
!                    else
!                       call EstimateNewtonADwr( grid, DWR, iter, newtonDone)
!                    endif
!
!                 else if (nlSolver%non_alg_stop == 'rezL2') then
!                    if(nlSolver%res0 < nlSolver%tol * nsize )  then
!                       !write(*,*) 'Reziduum:' , nlSolver%res0,  '<' , nlSolver%tol * nsize
!                       !call cpu_time(t2)
!                       !time_prepare = time_prepare + t2 - t1
!                       ! Newton method has achieved the prescribed tolerance
!                       newtonDone = .true.
!                    endif
!                    !write(*,*) 'Reziduum:' , nlSolver%res0,  '>' , nlSolver%tol * nsize
!                 else
!                    stop 'Unknown method for nonlinear algebraic error estimation in NewtonSolve!'
!                 endif  ! end of other technique
!                 class default
!                  stop 'other type of nlSolver'
!              end select
!
!              call state%cpuTime%addEstimTime( )
!
!              ! Write Newton output
!              open(91, file='criter', status='unknown', position='append')
!              if(nlSolver%iter == 1) write(91,'(x)')
!              if(nlSolver%Aiter == 1) write(91,'(x)')
!              !if(state%time%iter == 1) print*, 'TODO: Repair NewtonOutput for other methods than aRES!'
!              call Newton_output(grid, 91, nDamping, loc_implicitly )
!              close(91)
!
!              nlSolver%res1 = nlSolver%res
!
!              iter = iter + 1
!           end do ! while
!           ! end of newton step, save the solution, write down and  go on
!
!
!        !nlSolver%converged = newtonDone
!        nlSolver%converged = .true.
!        if(estS > 10. * estS_ini) nlSolver%converged = .false. ! DIVERGENCE
!
!      end associate ! nlSolver
!
!      !print*,' end subroutine NewtonSolve'
!   end subroutine NewtonSolveWSTplus

   ! compute the norm of the residual locally for each element in the mesh
   subroutine localNewtonEstimate( newton, grid, eta, estim )
      class( Newton_t ), intent ( inout ) :: newton
      class( mesh ), intent(in) :: grid
      real, intent(in) :: eta ! 1/tau
      real, dimension(1:grid%nelem), intent(inout)  :: estim
      class(element), pointer :: elem
      real, dimension(:), allocatable :: y
      integer :: nsize, dof, ncv, i, nelem

      nelem = grid%nelem

      ! rhsST -> newton%b1
      call FillVectorST( newton%b1(1:state%nsize), eta )
      nsize = size( newton%b1 )
      allocate( y(1:nsize) )

      ! for iLU precond only !!!! ?
      if (state%linSolver%name /= 'GMRES_ILU' ) &
         stop 'localNewtonEstimate only for ILU precond'

      if ( state%time%disc_time == 'STDG' ) then
         call bMViLUprodST( y, newton%b1, nsize )
      else
         stop 'Local damping not tested for BDF'
         call bMViLUprod(y, newton%b1 ,nsize)
      endif

      !print*, 'norm b', norm2( newton%b1 ) , 'y:',norm2(y)

      ! compute locally l2(K) norm of the vector y
      ! newton%res is already computed, but we need it locally
      do i = 1,nelem
         elem => grid%elem(i)
         dof = elem%dof * elem%Tdof *ndim
         ncv = elem%ncv
         ! use sqrt or not???
         estim(i) = dot_product(y(ncv:ncv+dof-1), y(ncv:ncv+dof-1) )**0.5
      end do !i

      deallocate(y)

   endsubroutine localNewtonEstimate

end module newton_mod
