prepare_regularized_qr Subroutine

  subroutine prepare_regularized_qr(factorization, matrix, ridge, qr_tolerance)
    type(regularized_qr_type), intent(inout) :: factorization
    real(dp), intent(in) :: matrix(:, :), ridge, qr_tolerance
    real(dp), allocatable :: augmented(:, :)
    real(dp) :: ridge_sqrt
    integer(i32) :: mrow, ncol, col_idx

    mrow = int(size(matrix, 1), i32)
    ncol = int(size(matrix, 2), i32)
    if (mrow <= 0_i32 .or. ncol <= 0_i32 .or. ridge < 0.0_dp .or. qr_tolerance <= 0.0_dp) then
      error stop 'prepare_regularized_qr received invalid dimensions or tolerances.'
    end if

    if (allocated(factorization%q)) deallocate (factorization%q)
    if (allocated(factorization%r)) deallocate (factorization%r)
    if (allocated(factorization%col_scale)) deallocate (factorization%col_scale)
    allocate (factorization%q(mrow + ncol, ncol), factorization%r(ncol, ncol))
    allocate (factorization%col_scale(ncol), augmented(mrow + ncol, ncol))

    augmented = 0.0_dp
    ridge_sqrt = sqrt(ridge)
    do col_idx = 1_i32, ncol
      factorization%col_scale(col_idx) = sqrt(sum(matrix(:, col_idx)*matrix(:, col_idx)))
      if (factorization%col_scale(col_idx) <= tiny(1.0_dp)) factorization%col_scale(col_idx) = 1.0_dp
      augmented(1:mrow, col_idx) = matrix(:, col_idx)/factorization%col_scale(col_idx)
      augmented(mrow + col_idx, col_idx) = ridge_sqrt
    end do
    call factor_tall_matrix_qr(augmented, factorization%q, factorization%r, qr_tolerance)
    factorization%mrow = mrow
    factorization%ncol = ncol
    factorization%preparation_count = factorization%preparation_count + 1_i32
  end subroutine prepare_regularized_qr