module dual_problem_mod
!  use compute_oper
  use ama_L2interpol
  use data_mod
  use dwr_mod
  use element_mod
  use eval_sol
  use euler_problem
  use lin_solvers
  use main_data
  use matrix_oper_int
  use mesh_mod
  use solution_mod
  use terms_mod
  use target_functional_mod

implicit none

   public :: InitDualLinSolver
   public :: DualLinSolver
   public :: PrimalLinSolver
   public :: prepareDualProblem
   public :: SolveDualProblem_stat
   !public :: updateProblemPlus
   public :: WriteAdwrErrorsScreen
   public :: WriteOutputScreenDP
   public :: WriteDWRErrorsScreen


contains

  !> solve \f$ C(u(T))^T z = J(\phi) \f$ for nonlinear stationary problem
  !> if state%p_mod == 0: saves to elem%zST, else: saves to zST_plus
  subroutine SolveDualProblem_stat(DWR, grid, convfile)
    type( DWR_t ), intent(inout) :: DWR
    class( mesh ), intent(inout) :: grid
    character(len=*), intent(in) :: convfile
    class( element ), pointer :: elem
    integer :: i, j, kk, k, elemDof, dof, degP, dofP, R_type, Qnum, nsize
    real :: res_max_val
!    logical :: precond_update
    character(len=50) :: plot_file =  'plot_sol.gnu' !'../DWR/plot_sol.gnu'
    integer :: iPlot = 39
    !real, allocatable, dimension(:,:,:,:) :: Dw ! for the least square reconstruction
    real, allocatable, dimension(:) :: bb, bc_z_error
    real :: l_norm
    real :: residuum
    integer :: lin_solver_not_conv, restart, nloops, p_mod, q_mod
    logical :: plus, loc_implicitly

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

    if ( p_mod > 0) then
        plus = .true.
    else
        plus = .false.
    end if

!    print*, 'SolveDualProblem_stat called !!!'

    call state%cpuTime%startEstimTime()

    if ( state%time_dependent ) &
      stop 'DWR method is not implemented for time-dependent problems!'

    ! TODO - control allocation of b,x,rr
    if (allocated(DWR%b) .or. allocated(DWR%x) .or. allocated(DWR%rr) ) then
        print*, 'DWR%b,x or rr already allocated in SolveDualProblem_stat'
        if (allocated(DWR%b)) deallocate(DWR%b)
        if (allocated(DWR%x)) deallocate(DWR%x)
        if (allocated(DWR%rr)) deallocate(DWR%rr)
    endif

    ! change in future to dual_size add control for nonlinear problems size can be ok

    nsize = state%bigNSize( p_mod, q_mod ) !%nsize
    allocate(DWR%b(1:nsize), DWR%x(1:nsize), DWR%rr(1:nsize) )

    ! do i=1, grid%nelem
    !    elem => grid%elem(i)
    !    !write(*,'(a8, 40es12.4)') "%zST A1:",elem%zST(1:ndim, 1:elem%dof, 1:elem%Tdof)
    !    DWR%x(1:elem%dof*ndim) =  &
    !         Transfer_funST_to_fun( elem%zST(1:ndim, 1:elem%dof, 1:elem%Tdof), &
    !         elem%dof, elem%Tdof, 0, elem%TQnum)

    !    call PlotElemFunction3D(800+state%space%adapt%adapt_level, elem, &
    !         elem%dof, DWR%x(1:elem%dof*ndim)  )
    ! enddo


    ! ZST -> DWR%x - INIT guess
    !call copyZST_toLongVector( grid,  nsize, )
    if (plus) then
        call copyZSTplus_toLongVectorMod( grid,  p_mod, q_mod, DWR%x )
    else
        call copyZST_toLongVectorMod( grid,  p_mod, q_mod, DWR%x )
    end if

    eta = 0.0

    !filling elem%rhsST into global rhs
    ! DWR%rhs -> DWR%b
    call DWR%fillRHSvector( grid, plus ) !used to be: .false.


    !Nullify the number of lin. iterations before dual problem.
    DWR%linSolver_iter = 0
    !maybe globaly ?
    DWR%residuum = 0.0 !state%linSolver%residuum
    lin_solver_not_conv = state%linSolver%lin_solver_not_conv

    ! FR  restart and nloops should be set globally
    !set DWR_PARAMETERS
    restart = 45
    nloops =  40 ! 20 ! 55
    !if(state%time%maxiter == 43) nloops = 5 !2  ! SIMULATION of  aDWR stopping criterion

    write(*,'(a50, i8, a10, i8, a10, es12.4)') &
         "DWR - params for dual lin problem: Restart = ", restart, &
         " Nloops = ", nloops, " tol = ", DWR%linSolver_tol

!    call SolveBlockLinearSTDGMProblem_Dual( nsize, eta, DWR%b, DWR%x, &
!            DWR%residuum, DWR%linSolver_iter, &
!            lin_solver_not_conv, DWR%linSolver_tol, restart, nloops)

    call SolveBlockLinearBigProblem_Dual( nsize, eta, DWR%b, DWR%x, &
            DWR%residuum, DWR%linSolver_iter, &
            lin_solver_not_conv, DWR%linSolver_tol, restart, nloops)

    call bMVprodBIG_Dual(DWR%rr, DWR%x, nsize)    ! rr = Ax


    !print*, "---xxx"
    ! ? is the sign correct?
    DWR%rr(:) = DWR%rr(:) - DWR%b(:)
    !print*, 'norm od dual res 3 = ' , norm2(DWR%rr)

    if( lin_solver_not_conv > 0 ) then
      print*, 'lin solver for DUAL problem has not converged yet '
      print*
          ! does not do anything !
          open(53, file='GMRES_DP_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) =  DWR%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
          close(53)
    endif

    ! REZ VD MOD
    !DWR%x(:) = DWR%x(:)  + DWR%b1(:)

    ! 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
        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)
    call computeDualResidualPlus(grid, DWR, plus)
    call state%setP_mod( p_mod )
    call state%setQ_mod( q_mod )

    ! FR_aDWR What is this for???
!     do i=1, grid%nelem
!        elem => grid%elem(i)
!        !write(*,'(a8, 40es12.4)') "%zST A1:",elem%zST(1:ndim, 1:elem%dof, 1:elem%Tdof)
!        DWR%x(1:elem%dof*ndim) =  &
!             Transfer_funST_to_fun( elem%zST(1:ndim, 1:elem%dof, 1:elem%Tdof), &
!             elem%dof, elem%Tdof, 0, elem%TQnum)
!
!        call PlotElemFunction3D(900+state%space%adapt%adapt_level, elem, &
!             elem%dof, DWR%x(1:elem%dof*ndim)  )
!     enddo

    ! PLOT solution
    ! splot 'plot_sol.gnu' with lines notitle
!    open( iPlot, file = plot_file, action="write", status="replace" )
!    do i = 1, grid%nelem
!      elem => grid%elem(i)
!      ! from elem%w( 0, ...)
!      call PlotElemSolution3D(iPlot, elem)
!    end do
!
!    print*, " Dual solution saved to plot_sol.gnu .."
!    close(iPlot)

    call PlotSolDual(0)
    !print*, " Dual solution saved to plot_sol_dual_0.gnu .."

      ! number of iter = 1, used only with aDWR
    call WriteOutputScreenDP(DWR, 1)

    call state%cpuTime%addEstimTime()

    DWR%dualProblem_computed = .true.

    associate( JJ => DWR%J )
        select type ( JJ )
            ! integrate over edges
            class is ( DragAndLift_t)
                print*, "Test in SolveDualProblem_stat - z*n on Gamma_W"
!                print*, "----"
!                print*, "THE theta was set to (1,1) for testing!!!"
!                print*, "----"
                allocate( bc_z_error(1:grid%nelem), source = 0.0 )
                do i = 1, grid%nelem
                    call JJ%testBC_Euler_W( grid%elem(i), bc_z_error(i) )
                end do

                if ( maxval( bc_z_error(1:grid%nelem)) > 0.01 ) then
                    print*, "On some elements the ADJOINT BC is not satisfied!"
                    print*, "Max value: ", maxval( bc_z_error(1:grid%nelem))  , &
                        " on element:", maxloc( bc_z_error(1:grid%nelem))
                end if

                deallocate( bc_z_error )
            class default
               !stop "cadscdsa"
                !print*, "cadscdsa ?????????????????"
        end select
    end associate



!    l_norm = Solution_L8norm( grid )
!    print*, ' L8 norm of Z: ', l_norm
!    print*,

  end subroutine SolveDualProblem_stat




  !> prepare the right-hand side of the dual problem
  !> \f$ C(u(T))^T z = J(\phi) \f$ for nonlinear stationary problem
  !> see subroutine SolveDualProblem_stat above
  subroutine PrepareDualProblem_RHS(DWR, grid)
    type( DWR_t ), intent(inout) :: DWR
    class( mesh ), intent(inout) :: grid
    class( element ), pointer :: elem
    integer :: i, j, kk, k, elemDof, dof, degP, dofP, R_type, Qnum, nsize
    real :: res_max_val
!    logical :: precond_update
    !real, allocatable, dimension(:,:,:,:) :: Dw ! for the least square reconstruction
    real, allocatable, dimension(:) :: bb
    real :: l_norm
    real :: residuum
    integer :: lin_solver_not_conv, restart, nloops, p_mod, q_mod
    logical :: plus, loc_implicitly

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

    if ( p_mod > 0) then
        plus = .true.
    else
        plus = .false.
    end if

    !print*, 'SolveDualProblem_stat called !!!'

    !call state%cpuTime%startEstimTime()

    if ( state%time_dependent ) &
      stop 'DWR method is not implemented for time-dependent problems YET!'

    nsize = state%bigNSize( p_mod, q_mod ) !%nsize

    ! TODO - control allocation of b,x,rr
    call DWRarrays_allocated_size ( DWR, nsize)

    ! do i=1, grid%nelem
    !    elem => grid%elem(i)
    !    !write(*,'(a8, 40es12.4)') "%zST A1:",elem%zST(1:ndim, 1:elem%dof, 1:elem%Tdof)
    !    DWR%x(1:elem%dof*ndim) =  &
    !         Transfer_funST_to_fun( elem%zST(1:ndim, 1:elem%dof, 1:elem%Tdof), &
    !         elem%dof, elem%Tdof, 0, elem%TQnum)

    !    call PlotElemFunction3D(800+state%space%adapt%adapt_level, elem, &
    !         elem%dof, DWR%x(1:elem%dof*ndim)  )
    ! enddo


    ! ZST -> DWR%x - INIT guess
    !call copyZST_toLongVector( grid,  nsize, )
    !print*,'DWR plus ehue', plus, nsize, p_mod, q_mod
    if (plus) then
        call copyZSTplus_toLongVectorMod( grid,  p_mod, q_mod, DWR%x )
    else
        call copyZST_toLongVectorMod( grid,  p_mod, q_mod, DWR%x )
    end if

    eta = 0.0

    !filling elem%rhsST into global rhs
    ! DWR%rhs -> DWR%b
    call DWR%fillRHSvector( grid, plus ) !used to be: .false.

  end subroutine PrepareDualProblem_RHS



  !> solve and estimate the linear primal problem until
  !> DWR%estimLP < tol
  !> only for aDWR method - uses DWR%aDWR class, which may not be initialized in other cases !
  subroutine PrimalLinSolver(DWR, grid, newton, tol)
    type( DWR_t ), intent(inout) :: DWR
    class( mesh ), intent(inout) :: grid
    class( NonlinearSol_t ), intent(inout) :: newton
!    integer, intent(in) :: iter ! this subroutine is called for the iter-th time
    real, intent(in) :: tol ! given tolerance
    class( element ), pointer :: elem
    real :: tt, res_max_val, eta1
    integer :: elemDof, kk, i, nsize
    integer :: iPlot = 68
    integer :: lin_solver_not_conv
    integer :: restart, nloops
    logical :: done
    integer :: it, maxIt

    eta1 = 0
    nsize = state%nsize

    ! maximal number of iterations
    if ( DWR%aDWR%fixedIter ) then
      maxIt = 1 ! number of repetitions of the cycle
      nloops = 1 ! restart
    else
      maxIt = 20
      nloops = 5
    endif
    it = 1
    done = .false.

    ! only one loop for the alg iterative solver is done
    ! TODO - better working with linear solver ???


    do while (.not. done .and. ( it <= maxIt) )
       !Solve
       call state%cpuTime%startSolveTime()
       stop "new big block structure should be used "
       call SolveBlockLinearSTDGMProblem( nsize, eta1, newton%b, newton%x, &
            state%linSolver%residuum, DWR%aDWR%iter_lin_primal, &
            state%linSolver%lin_solver_not_conv, DWR%aDWR%restart_primal, nloops)

       call state%cpuTime%addSolveTime()
       call state%cpuTime%startPrepareTime()

       !print*, 'iter in PrimalLinSolver:' , DWR%aDWR%iter_lin_primal, state%linSolver%iter
   !    time_solve = time_solve + t1 - t2
       !should be changed

       call bMVprodST(newton%rr, newton%x, state%nsize)    ! rr = Ax
       newton%rr(:) = Newton%rr(:) - Newton%b(:)
       call state%cpuTime%addPrepareTime()

       ! estimate -> DWR%estimLP
       call computePrimalLinearDWRestimates( DWR, grid, newton%rr)

       ! is the error measured with respect to the target quantity
       ! under the given tolerance ?
       if ( DWR%estimLP < tol ) &
         done = .true.

       it = it + 1
    end do !while

    if (.not. done) &
      print*, 'Primal linear problem was not solved with needed tol. '
   print*, 'linPP: it =' , DWR%aDWR%iter_lin_primal , 'tol=' , tol, &
           'and estimate=' , DWR%estimLP

    ! PLOT solution
    ! splot 'plot_sol.gnu' with lines notitle
!    open( iPlot, file = plot_file, action="write", status="replace" )
!    do i = 1, grid%nelem
!      elem => grid%elem(i)
!      ! from elem%w( 0, ...)
!      call PlotElemSolution3D(iPlot, elem)
!    end do
!    print*, 'Primal solution saved to plot_sol.gnu'
!    close(iPlot)

      ! the only use of iter
     !call WriteOutputScreenDP(DWR, iter )

  end subroutine PrimalLinSolver




  !> solve and estimate the linear dual problem until
  !> DWR%estimLD < DWR%aDWR%linTol
  !> for the problem \f$  C(u(T))^T z = J(\phi) \f$
  !> only for aDWR method - uses DWR%aDWR class, which may not be initialized in other cases !
  subroutine DualLinSolver(DWR, grid, iter, tol)
    type( DWR_t ), intent(inout) :: DWR
    class( mesh ), intent(inout) :: grid
    integer, intent(in) :: iter ! this subroutine is called for the iter-th time
    real, intent(in) :: tol
    class( element ), pointer :: elem
    real :: t1, res_max_val, eta1
    integer :: elemDof, kk, i, nsize
    integer :: iPlot = 68
    integer :: lin_solver_not_conv
    integer :: restart, nloops
    logical :: done
    integer :: it, maxIt

    eta1 = 0
    nsize = state%nsize

    ! maximal number of iterations
    if( DWR%aDWR%fixedIter ) then
      maxIt = 1 ! only one restart for tables and decrease graph
    else
      maxIt = 50
    endif

    it = 1
    done = .false.

    ! only one loop for the alg iterative solver is done
    ! TODO - better working with linear solver ???
    nloops = 1


    do while (.not. done .and. ( it <= maxIt) )
      !Solve
      print*, 'SolveBlockLinearSTDGMProblem_Dual should be replaced by SolveBigBlockLinearProblem_Dual'
      call SolveBlockLinearSTDGMProblem_Dual( nsize, eta1, DWR%b, DWR%x, &
              DWR%residuum, DWR%aDWR%iter_lin_dual, &
              lin_solver_not_conv, DWR%linSolver_tol, DWR%aDWR%restart_dual, nloops)

       !should be changed
       call bMVprodST_Dual(DWR%rr, DWR%x, nsize)    ! rr = A^Tx

       DWR%rr(:) = DWR%rr(:) - DWR%b(:)

       ! copy the solution vector back to elements WST, W, and wS for the least squares
       call CopyZST_fromLongVector( grid, nsize, DWR%x(1:nsize) )

      !Estimate
      call computeDualLinearDWRestimates(DWR, grid)

      if ( DWR%estimLD < tol ) &
         done = .true.

      it = it + 1
    end do !while

!    if( lin_solver_not_conv > 0 ) then
!      print*, 'lin solver for DUAL problem does not converge'
!      ! does not do anything !
!      open(53, file='GMRES_DP_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) =  DWR%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
!      close(53)
!   endif

    ! PLOT solution
    ! splot 'plot_sol.gnu' with lines notitle
!    open( iPlot, file = plot_file, action="write", status="replace" )
!    do i = 1, grid%nelem
!      elem => grid%elem(i)
!      ! from elem%w( 0, ...)
!      call PlotElemSolution3D(iPlot, elem)
!    end do
!    print*, 'Dual solution saved to plot_sol.gnu'
!    close(iPlot)

     ! only for output
     DWR%linSolver_iter =  DWR%aDWR%iter_lin_dual
     ! we have at least some approximation to zST
     DWR%dualProblem_computed = .true.

      ! the only use of iter
     call WriteOutputScreenDP(DWR, iter )
     print*, 'Dual problem solved!'

  end subroutine DualLinSolver


!  !> updates the size of arrays and matrix blocks when the polynomial degree globally changes
!  !> we have to call computeTerms afterwards to fill the blockST...
!  subroutine updateProblemPlus( grid, plus )
!    class( mesh ), intent( inout ) :: grid
!!    type( DWR_t ), intent( in ) :: DWR
!    logical, intent(in) :: plus ! true - plus state%space%plusDeg, false - minus %plusDeg
!    class( element ), pointer :: elem
!    integer :: i, nelem
!    integer :: p_plus
!
!    print*,
!    print*, 'updateProblemPlus called',  plus
!    print*,
!
!    if (state%time_dependent) then
!      stop 'updateProblemPlus time dependent - not implemented'
!      ! reinit the time stepping - reverse from (0,T) to (T,0)
!      !call state%time%reverseTime()
!    endif
!
!    nelem = grid%nelem
!
!    ! increase/decrease the polynomial degree
!    if ( plus ) then
!         p_plus = state%p_mod_max
!    else
!         p_plus = - state%p_mod_max
!    endif
!
!    !print*,' the size of the system changes, also dofs and quadratures', p_plus
!    if ( p_plus /= 0 ) then
!
!      call grid%setNewElementDofs( p_plus )
!
!      ! update the sizes of mesh arrays and matrices
!      call UpdateMesh_p( grid )
!    endif
!
!  end subroutine updateProblemPlus

  !> prepare the structures for computation of the dual problem
  !> it should be called before every solution process, i.e. after every mesh adaptation
  subroutine PrepareDualProblem( DWR )
    class( DWR_t ), intent(inout) :: DWR
    class(element), pointer :: elem
    integer :: i,j
    integer, pointer :: exp1, exp2
    character(len=20) :: loc_adapt_space
    character(len=50) :: plot_file_primal = 'plot_sol_primal.gnu' !'../DWR/plot_sol_primal.gnu'
    integer :: iPlotPrimal = 41
    logical :: impl_loc

    print*, ' ############## subroutine PrepareDualProblem( DWR )'
    call state%cpuTime%startEstimTime()

    if (state%time%quadType /= 'Radau') &
      stop 'Radau quadrature must be used - to compute C(u(T))'

    ! print primal solution
    open( iPlotPrimal, file = plot_file_primal, action="write", status="replace" )
    do i = 1, grid%nelem
      elem => grid%elem(i)
      call Transfer_wST_to_w_Elem( elem , 0, elem%TQnum)
      call PlotElemSolution3D(iPlotPrimal, elem)
    end do
    print*, 'Primal solution saved to plot_sol_primal.gnu'
    close(iPlotPrimal)


    !print*, 'ComputeST_Terms in prepare dual problem!, deg_plus = false, impl =' , state%nlSolver%implicitly
    !!! for nonlinear problems we need also the matrix at u(T)
    ! for DP computation, dual error estimate (we need bigger matrix for the res computation), for Ritz reconstruction
    ! for linear problems the matrix is already computed since it does not ~ on u_h
    impl_loc = state%nlSolver%implicitly
    if ( .not. state%model%linear ) then
      state%nlSolver%implicitly = .true.
      call ComputeST_Terms( .true. )
    end if
    !!! in all case we need the RES vector with degPlus=True to : r_h(u_h) in error est and Ritz reconstruction
    state%nlSolver%implicitly = .false.
    call ComputeST_Terms( .true. )
    state%nlSolver%implicitly = impl_loc


    call DWR%J%findSupp( grid )
    ! RHS must be computed after ComputeSTDGM_Terms( )
    call DWR%setRHS_new( grid )

    if ( DWR%J%isupp == 0 ) then
        print*, 'epsilon: ', DWR%J%eps, 'grid%h/2', grid%h / 2.0
       stop 'zero support of target func (4)'
    endif

    if (state%time_dependent) &
      stop 'we need to compute IC for dual problem in PrepareDualProblem!'

    call state%cpuTime%addEstimTime()

  end subroutine PrepareDualProblem

  !> write outputs to the screen
  !> finishComputationProcess should be called after this subroutine - they used to be together
  subroutine WriteOutputScreenDP( DWR, iter )
    class( DWR_t ), intent(in) :: DWR
    integer :: ifile = 11
    integer :: it_hours, it_min, it_sec, i, is
    character(len=1) :: ch
    character(len=3) :: yes_conv ! convergence?
    logical :: Aname
    real:: t_end, ratio
    real :: err1, err2, err3
    integer linIter, nlIter, iter

!    iter = 1
    linIter = DWR%linSolver_iter ! DWR%aDWR%lin_iter_dual ???
    ! dual problem is linear
    !nlIter = 0 !state%nlSolver%iter

    ! output of files tri*xxxxx and/or sol*xxxxx
    ch = ' '

    ! no adaptation
    if( state%space%adapt%max_adapt_level == 0 ) then
       if(state%time%ttime >= 0.999999*state%time%FinTime  ) ch = '*'  ! Final time ?
       if( iter >= state%time%maxiter ) ch = '*'      ! Last iteration?
    endif

    ! screen output
    print*,'iter    tau       time     Ax=b&
      & res(iter)'
    print*,'---------------------------------------------------'

    !if(iter == 1 .or. ch == '*') then
    if( state%modelName == 'scalar' .or. state%modelName == '2eqs' .or. &
         state%modelName == 'NSe' .or. state%modelName == 'Eul' ) then
       write(*,'(i5,es10.2,es11.3,a1,es9.2,a1,i4,a1)') &
             iter, state%time%tau(1),&
             state%time%ttime,ch, DWR%residuum,'(',linIter,')'
    else
       stop 'WriteOutputScreenDP not done for this type of problems!'
    endif
    !endif

  end subroutine WriteOutputScreenDP

  !> write on screen the table of dual error estimates
  subroutine WriteDWRErrorsScreen( DWR )
   class( DWR_t ), intent(in) :: DWR
   integer :: i

    print*,'  J(u_h)    J(u)       Error     eta_S    eta_Sabs     eta_A'
    print*,'----------------------------------------------------------------'
    write(*,'( es9.2, a1, es9.2, a2 , es9.2, a2 ,  es9.2, a2, es9.2, a2 , es9.2 )') &
        DWR%J%Ju , ' ', DWR%J%Ju_exact , '  ', DWR%J%Ju - DWR%J%Ju_exact , '  ', &
        sqrt(state%estim(dwrS, 1)) , '  ',sqrt(state%estim(dwrS_abs, 1)), '  ', &
        sqrt( state%estim(dwrA, 1) )
    print*,'----------------------------------------------------------------'
    print*, 'eta_dS    eta_dSabs     eta_dA'
    print*,'----------------------------------------------------------------'
    write(*,'( es9.2, a2, es9.2, a2 , es9.2)') &
            sqrt( state%estim(dwr_dualS, 1)) , '  ', &
            sqrt( state%estim(dwr_dualS_abs, 1) ) , '  ', &
            sqrt( state%estim(dwr_dualA, 1) )
    print*,'----------------------------------------------------------------'
    print*, 'Ju:' , DWR%J%Ju , ' Ju_exact:', DWR%J%Ju_exact
    print*,'----------------------------------------------------------------'

    print*, "etaS = ", sqrt(state%estim(dwrS, 1)) , " /= dual etaS = ", sqrt(state%estim(dwr_dualS, 1))
!    do i= 1,grid%nelem
!        print*, "elem%eta =", grid%elem(i)%eta(dwrS,1), grid%elem(i)%eta(dwr_dualS,1)
!    end do

  end subroutine WriteDWRErrorsScreen

    !> write on screen the table of dual error estimates
  subroutine WriteAdwrErrorsScreen( DWR )
   class( DWR_t ), intent(in) :: DWR

    print*,'  Error     eta_Discr    etaNL     etaA    etaAD'
    print*,'----------------------------------------------------------------'
    write(*,'( es9.2, a1, es9.2, a2 , es9.2, a2 ,  es9.2, a2, es9.2)') &
            DWR%J%Ju - DWR%J%Ju_exact , '  ', &
            DWR%estimDiscr, '  ', DWR%estimNL, '  ', &
            DWR%estimLP, '  ', DWR%estimLD
    print*,'----------------------------------------------------------------'

  end subroutine WriteAdwrErrorsScreen


   !> update the structures for computation of the dual problem
   !> it should be called before lin_problem solve
   !> used only for aDWR method for estimation of the alg errors
  subroutine UpdateDualProblem( DWR )
    class( DWR_t ), intent(inout) :: DWR
    class(element), pointer :: elem
    integer :: i,j
    integer, pointer :: exp1, exp2
    logical :: time_dependent
    character(len=20) :: loc_adapt_space
    character(len=50) :: plot_file_primal = 'plot_sol_primal.gnu' !'../DWR/plot_sol_primal.gnu'
    integer :: iPlotPrimal = 41
    logical :: impl_loc
    logical :: deg_plus

    deg_plus = .false.  !??

    ! we need to compute the matrix at u(T)
    impl_loc = state%nlSolver%implicitly
    state%nlSolver%implicitly = .true.
    ! ComputeSTDGM_Terms( ) call:
    ! print*, 'calling compute terms in UpdateDualProblem. Although expensive, necessary for nonl problems.'
    ! FR: calling compute terms in UpdateDualProblem. Although expensive, necessary for nonl problems.
    !call ComputeSTDGM_Terms( deg_plus )
    print*, 'ComputeSTDGM_Terms replaced by ComputeST_Terms in UpdateDualProblem! NOT DONE',deg_plus, state%nlSolver%implicitly
    call ComputeST_Terms( deg_plus )
    ! put back
    state%nlSolver%implicitly = impl_loc

!    call InitDualLinSolver( DWR, grid )

   end subroutine UpdateDualProblem

end module dual_problem_mod
