mxdpgu Subroutine

public pure subroutine mxdpgu(n, a, alf, x, y)

correction of a dense symmetric positive definite matrix a+e in the factored form a+e=l*d*trans(l) obtained by the subroutine mxdpgf. the correction is defined as a+e:=a+e+alf*x*trans(x) where alf is a given scaling factor and x is a given vector.

Method

  • p.e.gill, w.murray, m.saunders: methods for computing and modifying the ldv factors of a matrix, math. of comp. 29 (1974) pp. 1051-1077.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: n

order of the matrix a.
real(kind=wp), intent(inout) :: a(*)

a(n*(n+1)/2) factorization a+e=l*d*trans(l) obtained by the subroutine mxdpgf.
real(kind=wp), intent(in) :: alf

scaling factor in the correction term.
real(kind=wp), intent(in) :: x(*)

vector in the correction term.
real(kind=wp), intent(out) :: y(*)

auxiliary vector.

Calls

proc~~mxdpgu~~CallsGraph proc~mxdpgu mxdpgu proc~mxvscl mxvscl proc~mxdpgu->proc~mxvscl

Called by

proc~~mxdpgu~~CalledByGraph proc~mxdpgu mxdpgu proc~bfgs_variable_metric_update bfgs_variable_metric_update proc~bfgs_variable_metric_update->proc~mxdpgu proc~psqp psqp_class%psqp proc~psqp->proc~bfgs_variable_metric_update proc~psqpn psqp_class%psqpn proc~psqpn->proc~psqp

Source Code

      pure subroutine mxdpgu(n,a,alf,x,y)

      integer,intent(in) :: n  !! order of the matrix a.
      real(wp),intent(in) :: alf  !! scaling factor in the correction term.
      real(wp),intent(inout) :: a(*)  !! `a(n*(n+1)/2)` factorization `a+e=l*d*trans(l)`
                                      !! obtained by the subroutine [[mxdpgf]].
      real(wp),intent(in) :: x(*)  !! vector in the correction term.
      real(wp),intent(out) :: y(*)  !! auxiliary vector.

      real(wp) :: alfr
      real(wp) :: b , d , p , r , t , to
      integer :: i , ii , ij , j

      real(wp), parameter :: zero = 0.0_wp
      real(wp), parameter :: one = 1.0_wp
      real(wp), parameter :: four = 4.0_wp
      real(wp), parameter :: con = 1.0e-8_wp

      if ( alf>=zero ) then
         ! forward correction in case when the scaling factor is nonnegative
         alfr = sqrt(alf)
         call mxvscl(n,alfr,x,y)
         to = one
         ii = 0
         do i = 1 , n
            ii = ii + i
            d = a(ii)
            p = y(i)
            t = to + p*p/d
            r = to/t
            a(ii) = d/r
            b = p/(d*t)
            if ( a(ii)<=four*d ) then
               ! an easy formula for limited diagonal element
               ij = ii
               do j = i + 1 , n
                  ij = ij + j - 1
                  d = a(ij)
                  y(j) = y(j) - p*d
                  a(ij) = d + b*y(j)
               end do
            else
               ! a more complicate but numerically stable formula for unlimited
               ! diagonal element
               ij = ii
               do j = i + 1 , n
                  ij = ij + j - 1
                  d = a(ij)
                  a(ij) = r*d + b*y(j)
                  y(j) = y(j) - p*d
               end do
            endif
            to = t
         end do
      else
         ! backward correction in case when the scaling factor is negative
         alfr = sqrt(-alf)
         call mxvscl(n,alfr,x,y)
         to = one
         ij = 0
         do i = 1 , n
            d = y(i)
            do j = 1 , i - 1
               ij = ij + 1
               d = d - a(ij)*y(j)
            end do
            y(i) = d
            ij = ij + 1
            to = to - d*d/a(ij)
         end do
         if ( to<=zero ) to = con
         ii = n*(n+1)/2
         do i = n , 1 , -1
            d = a(ii)
            p = y(i)
            t = to + p*p/d
            a(ii) = d*to/t
            b = -p/(d*to)
            to = t
            ij = ii
            do j = i + 1 , n
               ij = ij + j - 1
               d = a(ij)
               a(ij) = d + b*y(j)
               y(j) = y(j) + p*d
            end do
            ii = ii - i
         end do
      endif
      end subroutine mxdpgu