preconditioner_apply Subroutine

  subroutine preconditioner_apply(this, jacobian, vector, estimate)
    !* Author: Chris MacMackin
    !  Date: December 2016
    !
    ! Use Picard iteration to approximately multiply the state vector
    ! by the inverse Jacobian. The details for this procedure are in
    ! the documentation of the [[preconditioner(type)]] type.
    !
    class(preconditioner), intent(in)                    :: this
    class(jacobian_block), dimension(:,:), intent(inout) :: jacobian
      !! An \(n\times n\) matrix approximating the Jacobian for which
      !! the preconditioner is used.
    class(scalar_field), dimension(:), intent(in)        :: vector
      !! A vector of size \(n\) which is to be preconditioned.
    class(scalar_field), dimension(:), intent(inout)     :: estimate
      !! On entry, an initial guess for the preconditioned vector. On
      !! exit, the iteratively determined value of the preconditioned
      !! vector.

    character(len=77), parameter :: success_format = '("Picard solver '// &
         'converged with error of ",es12.5," after ",i2," iterations.")'
    character(len=77), parameter :: failure_format = '("Picard solver '// &
         'reached maximum (",i3,") iterations, with error ",es12.5,".")'
    integer, parameter                             :: msg_len = 72
    class(scalar_field), dimension(:), allocatable :: prev_estimate
    integer                                        :: i, j, k, n
    real(r8)                                       :: max_err, old_max_err
    class(scalar_field), allocatable               :: tmp_field
    logical                                        :: first
    character(len=msg_len)                         :: msg
    
#ifdef DEBUG
    call logger%debug('preconditioner_apply','Entering function `precondition`.')
#endif
    call vector%guard_temp(); call estimate%guard_temp()
    n = size(vector)
    allocate(prev_estimate(n), mold=estimate)
    allocate(tmp_field, mold=vector(1))
#ifdef DEBUG
    if (size(jacobian,1) /= size(jacobian,2)) then
      error stop('Jacobian is not a square matrix.')
    end if
    if (size(jacobian,1) /= size(vector)) then
      error stop('Vector is of different size than Jacobian.')
    end if
    if (size(estimate) /= size(vector)) then
      error stop('Estimate is of different size than vector.')
    end if
#endif
    ! Until reached maximum number of iterations...
    do i = 1, this%max_iterations
      ! For each row of the Jacobian, solve for the diagonal block,
      ! with off-diagonal blocks applied to the previous guess and
      ! subtracted from the right-hand-side.
      max_err = 0._r8
      prev_estimate = estimate
      do j = 1, n
        first = .true.
        do k = 1, n
          if (k == j) cycle
          if (first) then
            first = .false.
            tmp_field = jacobian(j,k) * estimate(k)
          else
            tmp_field = tmp_field + jacobian(j,k) * estimate(k)
          end if
        end do
        estimate(j) = jacobian(j,j)%solve_for(vector(j) - tmp_field)
        max_err = max(max_err, &
                      maxval(abs( (estimate(j)-prev_estimate(j))/(prev_estimate(j)+1e-10_r8) )))
      end do
      ! If difference between result and previous guess is less than
      ! the tolerance, stop iterations
      if (max_err < this%tolerance) then
        write(msg,success_format) max_err, i
#ifdef DEBUG
        call logger%debug('preconditioner%apply',msg)
        call logger%debug('preconditioner%apply','Exiting function `precondition`.')
#endif
        call vector%clean_temp(); call estimate%clean_temp()
        return
      end if
      if (i > 1 .and. old_max_err <= max_err) then
        call logger%trivia('preconditioner%apply','Iterations diverging. Exiting '// &
                           'and returning previous iterate, with maximum error '//    &
                           str(old_max_err)//'.')
#ifdef DEBUG
        call logger%debug('preconditioner%apply','Exiting function `precondition`.')
#endif
        estimate = prev_estimate
        call vector%clean_temp(); call estimate%clean_temp()
        return
      end if
      old_max_err = max_err
    end do
    write(msg,failure_format) this%max_iterations, max_err
    call logger%warning('preconditioner%apply',msg)
    call vector%clean_temp(); call estimate%clean_temp()
    call logger%debug('preconditioner%apply','Exiting function `precondition`.')
  end subroutine preconditioner_apply