pseudospec_block_solve_vector Function

  function pseudospec_block_solve_vector(this, rhs, bound_loc, bound_val, &
                                         good_bound) result(solution)
    !* Author: Chris MacMackin
    !  Date: September 2017
    !
    ! Solves the linear(ised) system represented by this finite
    ! difference block, for a given right hand side state vector
    ! (represented by a vector field).
    !
    ! @Warning Currently this is only implemented for a 1-D field.
    !
    ! @Bug For some reason, calls to the `vector_dimensions()` method
    ! produce a segfault when `rhs` is
    ! `class(vector_field)`. Everything works fine if it is
    ! `class(cheb1d_vector_field)`, so this is used as a workaround.
    !
    class(pseudospec_block), intent(inout) :: this
    class(cheb1d_vector_field), intent(in) :: rhs
      !! The right hand side of the linear(ised) system.
    integer, intent(in)                    :: bound_loc
      !! Which boundary is being set. The boundary will be the one
      !! normal to dimension of number `abs(boundary)`. If the
      !! argument is negative, then the lower boundary is returned. If
      !! positive, then the upper boundary is returned.
    class(vector_field), intent(in)        :: bound_val
      !! The value of the result at the specified boundary.
    integer, intent(in), optional          :: good_bound
      !! If provided, indicates which boundary contains trusted
      !! information from which to calculate the power of the highest
      !! frequency mode. Defaults to the opposite of `bound_loc`.
    class(vector_field), pointer           :: solution

    real(r8), dimension(:), allocatable :: sol_vector, bound_vec
    logical                             :: valid_bound
    integer                             :: bloc, n, i
    class(scalar_field), pointer        :: bcomponent
 
    call rhs%guard_temp(); call bound_val%guard_temp()
    valid_bound = .true.
    n = size(this%xvals)
    select case(bound_loc)
    case(1)
      bloc = 1
    case(-1)
      bloc = n
    case default
      valid_bound = .false.
      !call logger%warning('pseudospec_block%solve_for', '1-D field does not '// &
      !                    'have boundary in dimension '//str(bound_loc))
    end select
    sol_vector = rhs%raw()
    call bound_val%allocate_scalar_field(bcomponent)
    call bcomponent%guard_temp()
    do i = 1, rhs%vector_dimensions()
      bcomponent = bound_val%component(i)
      if (valid_bound) then
        bound_vec = bcomponent%raw()
        call integrate_1d(sol_vector((i-1)*n+1:i*n), this%xvals, &
                          bloc, bound_vec(1), good_bound)
      else
        call integrate_1d(sol_vector((i-1)*n+1:i*n), this%xvals, &
                          good_bound=good_bound)
      end if
    end do
    call bcomponent%clean_temp()
#ifdef DEBUG
    call logger%debug('pseudospec_block%solve_for', &
                      'Successfully performed Chebyshev pseudospectral integration.')
#endif
    call rhs%allocate_vector_field(solution)
    call solution%unset_temp()
    call solution%assign_meta_data(rhs)
    call solution%set_from_raw(sol_vector)
    call rhs%clean_temp(); call bound_val%clean_temp()
    call solution%set_temp()
  end function pseudospec_block_solve_vector