factor_tall_matrix_qr Subroutine

public subroutine factor_tall_matrix_qr(matrix, q, r, qr_tolerance)

Re-orthogonalised (two-pass modified Gram-Schmidt) QR of a tall matrix. A column whose residual norm falls below qr_tolerance times its norm (at least 1) gets that floor as its pivot and a zero q column.

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in) :: matrix(:,:)
real(kind=dp), intent(out) :: q(:,:)
real(kind=dp), intent(out) :: r(:,:)
real(kind=dp), intent(in) :: qr_tolerance

Called by

proc~~factor_tall_matrix_qr~~CalledByGraph proc~factor_tall_matrix_qr factor_tall_matrix_qr proc~precompute_periodic_root_operator precompute_periodic_root_operator proc~precompute_periodic_root_operator->proc~factor_tall_matrix_qr proc~prepare_regularized_qr prepare_regularized_qr proc~precompute_periodic_root_operator->proc~prepare_regularized_qr proc~prepare_regularized_qr->proc~factor_tall_matrix_qr proc~core_complete_plan_impl core_complete_plan_impl proc~core_complete_plan_impl->proc~precompute_periodic_root_operator proc~core_build_panel_plan_impl core_build_panel_plan_impl proc~core_build_panel_plan_impl->proc~core_complete_plan_impl proc~core_build_plan_impl core_build_plan_impl proc~core_build_plan_impl->proc~core_complete_plan_impl

Source Code

  subroutine factor_tall_matrix_qr(matrix, q, r, qr_tolerance)
    real(dp), intent(in) :: matrix(:, :), qr_tolerance
    real(dp), intent(out) :: q(:, :), r(:, :)
    integer(i32) :: mrow, ncol, col_idx, basis_idx
    real(dp), allocatable :: v(:)
    real(dp) :: norm_v, corr, base_norm

    mrow = int(size(matrix, 1), i32)
    ncol = int(size(matrix, 2), i32)
    if (size(q, 1) /= mrow .or. size(q, 2) /= ncol) error stop 'factor_tall_matrix_qr q dimension mismatch.'
    if (size(r, 1) /= ncol .or. size(r, 2) /= ncol) error stop 'factor_tall_matrix_qr r dimension mismatch.'

    q = 0.0_dp
    r = 0.0_dp
    allocate (v(mrow))
    do col_idx = 1_i32, ncol
      v = matrix(:, col_idx)
      base_norm = max(sqrt(sum(v*v)), 1.0_dp)
      do basis_idx = 1_i32, col_idx - 1_i32
        r(basis_idx, col_idx) = dot_product(q(:, basis_idx), v)
        v = v - r(basis_idx, col_idx)*q(:, basis_idx)
      end do
      do basis_idx = 1_i32, col_idx - 1_i32
        corr = dot_product(q(:, basis_idx), v)
        r(basis_idx, col_idx) = r(basis_idx, col_idx) + corr
        v = v - corr*q(:, basis_idx)
      end do
      norm_v = sqrt(sum(v*v))
      if (norm_v <= qr_tolerance*base_norm) then
        r(col_idx, col_idx) = qr_tolerance*base_norm
      else
        r(col_idx, col_idx) = norm_v
        q(:, col_idx) = v/norm_v
      end if
    end do
  end subroutine factor_tall_matrix_qr