This routine orthogonalizes the vector vnew
against the previous kmp
vectors in the
v
matrix. It uses a modified Gram-Schmidt orthogonalization procedure with conditional
reorthogonalization.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=rk), | intent(inout) | :: | vnew(n) |
On entry, vector containing a scaled product of the Jacobian and the vector |
||
real(kind=rk), | intent(in) | :: | v(n,ll) |
Matrix containing the previous |
||
real(kind=rk), | intent(inout) | :: | hes(ldhes,ll) |
On entry, an upper Hessenberg matrix of shape |
||
integer, | intent(in) | :: | n |
Order of the matrix |
||
integer, | intent(in) | :: | ll |
Current order of the matrix |
||
integer, | intent(in) | :: | ldhes |
Leading dimension of |
||
integer, | intent(in) | :: | kmp |
Number of previous vectors the new vector |
||
real(kind=rk), | intent(out) | :: | snormw |
L-2 norm of |
Type | Visibility | Attributes | Name | Initial | |||
---|---|---|---|---|---|---|---|
integer, | public | :: | i | ||||
integer, | public | :: | i0 | ||||
real(kind=rk), | public | :: | arg | ||||
real(kind=rk), | public | :: | sumdsq | ||||
real(kind=rk), | public | :: | tem | ||||
real(kind=rk), | public | :: | vnorm |
pure subroutine dorth(vnew, v, hes, n, ll, ldhes, kmp, snormw) !! This routine orthogonalizes the vector `vnew` against the previous `kmp` vectors in the !! `v` matrix. It uses a modified Gram-Schmidt orthogonalization procedure with conditional !! reorthogonalization. use daskr_kinds, only: rk, zero use blas_interfaces, only: daxpy, ddot, dnrm2 implicit none real(rk), intent(inout) :: vnew(n) !! On entry, vector containing a scaled product of the Jacobian and the vector `v(:,ll)`. !! On return, the new vector orthogonal to `v(:,i0)`, where `i0 = max(1, ll - kmp + 1)`. real(rk), intent(in) :: v(n, ll) !! Matrix containing the previous `ll` orthogonal vectors `v(:,1)` to `v(:,ll)`. real(rk), intent(inout) :: hes(ldhes, ll) !! On entry, an upper Hessenberg matrix of shape `(ll, ll)` containing in `hes(i,k)`, !! for `k < ll`, the scaled inner products of `a*v(:,k)` and `v(:,i)`. !! On return, an upper Hessenberg matrix with column `ll` filled in with the scaled !! inner products of `a*v(:,ll)` and `v(:,i)`. integer, intent(in) :: n !! Order of the matrix `a`. integer, intent(in) :: ll !! Current order of the matrix `hes`. integer, intent(in) :: ldhes !! Leading dimension of `hes`. integer, intent(in) :: kmp !! Number of previous vectors the new vector `vnew` must be made orthogonal to !! (`kmp <= maxl`). real(rk), intent(out) :: snormw !! L-2 norm of `vnew`. integer :: i, i0 real(rk) :: arg, sumdsq, tem, vnorm ! Get norm of unaltered VNEW for later use. vnorm = dnrm2(n, vnew, 1) ! Do Modified Gram-Schmidt on VNEW = A*V(LL). ! Scaled inner products give new column of HES. ! Projections of earlier vectors are subtracted from VNEW. i0 = max(1, ll - kmp + 1) do i = i0, ll hes(i, ll) = ddot(n, v(1, i), 1, vnew, 1) tem = -hes(i, ll) call daxpy(n, tem, v(1, i), 1, vnew, 1) end do ! Compute SNORMW = norm of VNEW. ! If VNEW is small compared to its input value (in norm), then ! Reorthogonalize VNEW to V(*,1) through V(*,LL). ! Correct if relative correction exceeds 1000*(unit roundoff). ! Finally, correct SNORMW using the dot products involved. snormw = dnrm2(n, vnew, 1) if (vnorm + snormw/1000 /= vnorm) return ! @todo: fix this comparison sumdsq = zero do i = i0, ll tem = -ddot(n, v(1, i), 1, vnew, 1) if (hes(i, ll) + tem/1000 == hes(i, ll)) cycle ! @todo: fix this comparison hes(i, ll) = hes(i, ll) - tem call daxpy(n, tem, v(1, i), 1, vnew, 1) sumdsq = sumdsq + tem**2 end do if (sumdsq == zero) return arg = max(zero, snormw**2 - sumdsq) snormw = sqrt(arg) end subroutine dorth