solve_for Subroutine

  subroutine solve_for(this, velocity, velocity_dx)
    !* Author: Chris MacMackin
    !  Date: January 2018
    !
    ! Inverts the linear portions of the plume momentum equation with
    ! the provided data. This is done by solving the linear ODE
    ! described in the documentation for the [[coriolis_block]]
    ! type. The block object must first have been initialised using
    ! the constructor.
    !
    ! @Warning Currently this is only implemented for a 1-D field.
    !
    class(coriolis_block), intent(inout) :: this
    class(vector_field), intent(inout)   :: velocity
      !! On input, the velocity value being preconditioned. On output,
      !! the preconditioned velocity.
    class(vector_field), intent(inout)   :: velocity_dx
      !! On input, the velocity derivative being preconditioned. On
      !! output, the preconditioned velocity derivative.

    type(cheb1d_scalar_field), dimension(4) :: F, eDxC_r, eDxC_i, E_r, E_i
    integer :: i
    real(r8), dimension(1) :: rtmp, ctmp
    real(r8), dimension(4) :: bound_vals
    complex(r8), dimension(4) :: rhs, B, C_bounds
    type(cheb1d_scalar_field) :: tmp

    integer :: info
    real(r8) :: rcond

    call velocity%guard_temp(); call velocity_dx%guard_temp()

    ! Construct array of scalar matrices for easier manipulation
    F(1) = velocity%component(1)
    F(2) = velocity%component(2)
    F(3) = velocity_dx%component(1)
    F(4) = velocity_dx%component(2)
    ! Get boundary values
    tmp = F(1)%get_boundary(this%integrate_bound, 1)
    rtmp = tmp%raw()
    bound_vals(1) = rtmp(1)
    tmp = F(2)%get_boundary(this%integrate_bound, 1)
    rtmp = tmp%raw()
    bound_vals(2) = rtmp(1)
    tmp = F(3)%get_boundary(this%integrate_bound, 1)
    rtmp = tmp%raw()
    bound_vals(3) = rtmp(1)
    tmp = F(4)%get_boundary(this%integrate_bound, 1)
    rtmp = tmp%raw()
    bound_vals(4) = rtmp(1)
    ! If have boundary conditions at both boundaries, correct the
    ! input fields
!    if (this%vel_bound_loc /= this%integrate_bound) then
!      F(1) = this%integrator%solve_for(F(1), this%vel_bound_loc, zero)
!      F(1) = F(1)%d_dx(1)
!      F(2) = this%integrator%solve_for(F(2), this%vel_bound_loc, zero)
!      F(2) = F(2)%d_dx(1)
!    end if
!    if (this%dvel_bound_loc /= this%integrate_bound) then
!      F(3) = this%integrator%solve_for(F(3), this%dvel_bound_loc, zero)
!      F(3) = F(3)%d_dx(1)
!      F(4) = this%integrator%solve_for(F(4), this%dvel_bound_loc, zero)
!      F(4) = F(4)%d_dx(1)
!    end if
    ! Calculate \(e^{-\bm{D}x}\bm{V}^{-1}\bm{F}\), aliasing variables
    ! which are not needed yet
    associate (emDxVinvF_r => eDxC_r, emDxVinvF_i => eDxC_i, &
               M_r => this%emDxVinv_r, M_i => this%emDxVinv_i, &
               C_r => E_r, C_i => E_i)
      emDxVinvF_r = [(M_r(i,1)*F(1) + M_r(i,2)*F(2) + M_r(i,3)*F(3) + &
                      M_r(i,4)*F(4), i=1,4)]
      emDxVinvF_i = [(M_i(i,1)*F(1) + M_i(i,2)*F(2) + M_i(i,3)*F(3) + &
                      M_i(i,4)*F(4), i=1,4)]
      ! Integrate \(e^{-\bm{D}x}\bm{V}^{-1}\bm{F}\) to get \(\bm{C}\)
      C_r = [(this%integrator%solve_for(emDxVinvF_r(i), this%integrate_bound, &
              zero), i=1,4)]
      emDxVinvF_r(1) = C_r(4)%d_dx(1)
      C_i = [(this%integrator%solve_for(emDxVinvF_i(i), this%integrate_bound, &
              zero), i=1,4)]
      ! Get the values of \(e^{\bm{D}x}\bm{C}\) at the upper boundary
      do i=1,4
        tmp = C_r(i)%get_boundary(-this%integrate_bound, 1)
        rtmp = tmp%raw()
        tmp = C_i(i)%get_boundary(-this%integrate_bound, 1)
        ctmp = tmp%raw()
        C_bounds(i) = cmplx(rtmp(1), ctmp(1), r8)
      end do
      ! Calculate \(e^{\bm{D}x}\bm{C}\)
      eDxC_r = [(this%eDx_r(i)*C_r(i) - this%eDx_i(i)*C_i(i), i=1,4)]
      eDxC_i = [(this%eDx_r(i)*C_i(i) + this%eDx_i(i)*C_r(i), i=1,4)]
    end associate
    
    ! Compute RHS for the linear system satisfying the boundary conditions
    if (this%vel_bound_loc == this%integrate_bound) then
      rhs(1:2) = bound_vals(1:2)
    else
      rhs(1:2) = bound_vals(1:2) - matmul(this%bound_matrix(1:2,:), &
                                          C_bounds)
    end if
    if (this%dvel_bound_loc == this%integrate_bound) then
      rhs(3:4) = bound_vals(3:4)
    else
      rhs(3:4) = bound_vals(3:4) - matmul(this%bound_matrix(3:4,:), &
                                          C_bounds)
    end if

    ! Compute coefficients for inhomogeneous components of solution so
    ! that boundary conditions are satisfied
    call la_gesvx(this%bound_matrix_scaled, rhs, B, this%factored_matrix, &
                  this%pivots, 'F', trans, this%equed, this%r_scales,     &
                  this%c_scales, info=i)
!print*,this%D_r
!print*,this%D_i
!print*,B
    ! Calculate \(\bm{E} = e^{\bm{D}x}\bm{B} + e^{\bm{D}x}\bm{C}\)
    E_r = [(this%eDx_r(i)*real(B(i)) - this%eDx_i(i)*aimag(B(i)) + &
            eDxC_r(i), i=1,4)]
    E_i = [(this%eDx_i(i)*real(B(i)) + this%eDx_r(i)*aimag(B(i)) + &
            eDxC_i(i), i=1,4)]

    ! Transform back to proper basis to get proper solution
    F = [(this%V_r(i,1)*E_r(1) - this%V_i(i,1)*E_i(1) + &
          this%V_r(i,2)*E_r(2) - this%V_i(i,2)*E_i(2) + &
          this%V_r(i,3)*E_r(3) - this%V_i(i,3)*E_i(3) + &
          this%V_r(i,4)*E_r(4) - this%V_i(i,4)*E_i(4), i=1,4)]
    velocity = F(1:2)
    velocity_dx = F(3:4)
    call velocity%clean_temp(); call velocity_dx%clean_temp()
  end subroutine solve_for