dorth Subroutine

pure subroutine dorth(vnew, v, hes, n, ll, ldhes, kmp, snormw)

Uses

  • proc~~dorth~~UsesGraph proc~dorth dorth module~blas_interfaces blas_interfaces proc~dorth->module~blas_interfaces module~daskr_kinds daskr_kinds proc~dorth->module~daskr_kinds module~blas_interfaces->module~daskr_kinds iso_fortran_env iso_fortran_env module~daskr_kinds->iso_fortran_env

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.

Arguments

Type IntentOptional Attributes Name
real(kind=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(kind=rk), intent(in) :: v(n,ll)

Matrix containing the previous ll orthogonal vectors v(:,1) to v(:,ll).

real(kind=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(kind=rk), intent(out) :: snormw

L-2 norm of vnew.


Calls

proc~~dorth~~CallsGraph proc~dorth dorth interface~daxpy daxpy proc~dorth->interface~daxpy interface~ddot ddot proc~dorth->interface~ddot interface~dnrm2 dnrm2 proc~dorth->interface~dnrm2

Variables

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

Source Code

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