solve_regularized_qr Subroutine

public subroutine solve_regularized_qr(factorization, rhs, solution)

Arguments

Type IntentOptional Attributes Name
type(regularized_qr_type), intent(in) :: factorization
real(kind=dp), intent(in) :: rhs(:)
real(kind=dp), intent(out) :: solution(:)

Called by

proc~~solve_regularized_qr~~CalledByGraph proc~solve_regularized_qr solve_regularized_qr proc~precompute_periodic_root_operator precompute_periodic_root_operator proc~precompute_periodic_root_operator->proc~solve_regularized_qr proc~core_build_panel_plan_impl core_build_panel_plan_impl proc~core_build_panel_plan_impl->proc~precompute_periodic_root_operator proc~core_build_plan_impl core_build_plan_impl proc~core_build_plan_impl->proc~precompute_periodic_root_operator

Source Code

  subroutine solve_regularized_qr(factorization, rhs, solution)
    type(regularized_qr_type), intent(in) :: factorization
    real(dp), intent(in) :: rhs(:)
    real(dp), intent(out) :: solution(:)
    real(dp), allocatable :: augmented_rhs(:), qtb(:), scaled_solution(:)
    integer(i32) :: col_idx

    if (.not. allocated(factorization%q) .or. .not. allocated(factorization%r) .or. &
        .not. allocated(factorization%col_scale)) then
      error stop 'solve_regularized_qr requires a prepared factorization.'
    end if
    if (size(rhs) /= factorization%mrow .or. size(solution) /= factorization%ncol) then
      error stop 'solve_regularized_qr dimension mismatch.'
    end if

    allocate (augmented_rhs(factorization%mrow + factorization%ncol))
    allocate (qtb(factorization%ncol), scaled_solution(factorization%ncol))
    augmented_rhs = 0.0_dp
    augmented_rhs(1:factorization%mrow) = rhs
    qtb = matmul(transpose(factorization%q), augmented_rhs)
    call solve_upper_triangular_system(factorization%r, qtb, scaled_solution)
    do col_idx = 1_i32, factorization%ncol
      solution(col_idx) = scaled_solution(col_idx)/factorization%col_scale(col_idx)
    end do
  end subroutine solve_regularized_qr