mxdpgf Subroutine

public pure subroutine mxdpgf(n, a, inf, alf, tau)

factorization a+e=l*d*trans(l) of a dense symmetric positive definite matrix a+e where d and e are diagonal positive definite matrices and l is a lower triangular matrix. if a is sufficiently positive definite then e=0.

Method

  • p.e.gill, w.murray : newton type methods for unconstrained and linearly constrained optimization, math. programming 28 (1974) pp. 311-350.

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) on input a given dense symmetric (usually positive definite) matrix a stored in the packed form. on output the computed factorization a+e=l*d*trans(l).
integer, intent(out) :: inf

an information obtained in the factorization process. if:

  • inf=0 then a is sufficiently positive definite and e=0. if
  • inf<0 then a is not sufficiently positive definite and e>0. if
  • inf>0 then a is indefinite and inf is an index of the most negative diagonal element used in the factorization process.
real(kind=wp), intent(inout) :: alf

on input a desired tolerance for positive definiteness. on output the most negative diagonal element used in the factorization process (if inf>0).
real(kind=wp), intent(out) :: tau

maximum diagonal element of the matrix e.

Called by

proc~~mxdpgf~~CalledByGraph proc~mxdpgf mxdpgf proc~dual_range_space_quad_prog psqp_class%dual_range_space_quad_prog proc~dual_range_space_quad_prog->proc~mxdpgf proc~psqp psqp_class%psqp proc~psqp->proc~dual_range_space_quad_prog proc~psqpn psqp_class%psqpn proc~psqpn->proc~psqp

Source Code

      pure subroutine mxdpgf(n,a,inf,alf,tau)

      real(wp),intent(inout) :: alf  !! on input a desired tolerance for positive definiteness. on
                                     !! output the most negative diagonal element used in the
                                     !! factorization process (if inf>0).
      real(wp),intent(out) :: tau  !! maximum diagonal element of the matrix e.
      integer,intent(out) :: inf  !! an information obtained in the factorization process. if:
                                  !!
                                  !!  * `inf=0` then a is sufficiently positive definite and e=0. if
                                  !!  * `inf<0` then a is not sufficiently positive definite and e>0. if
                                  !!  * `inf>0` then a is indefinite and inf is an index of the
                                  !!    most negative diagonal element used in the factorization
                                  !!    process.
      integer,intent(in) :: n  !! order of the matrix a.
      real(wp),intent(inout) :: a(*)  !! `a(n*(n+1)/2)`  on input a given dense symmetric (usually positive
                                      !! definite) matrix a stored in the packed form. on output the
                                      !! computed factorization `a+e=l*d*trans(l)`.

      real(wp) :: bet , del , gam , rho , sig , tol
      integer :: i , ij , ik , j , k , kj , kk , l

      l = 0
      inf = 0
      tol = alf

      ! estimation of the matrix norm

      alf = 0.0_wp
      bet = 0.0_wp
      gam = 0.0_wp
      tau = 0.0_wp
      kk = 0
      do k = 1 , n
         kk = kk + k
         bet = max(bet,abs(a(kk)))
         kj = kk
         do j = k + 1 , n
            kj = kj + j - 1
            gam = max(gam,abs(a(kj)))
         end do
      end do
      bet = max(tol,bet,gam/n)
      ! del = tol*bet
      del = tol*max(bet,1.0_wp)
      kk = 0
      do k = 1 , n
         kk = kk + k
         ! determination of a diagonal correction
         sig = a(kk)
         if ( alf>sig ) then
            alf = sig
            l = k
         endif
         gam = 0.0_wp
         kj = kk
         do j = k + 1 , n
            kj = kj + j - 1
            gam = max(gam,abs(a(kj)))
         end do
         gam = gam*gam
         rho = max(abs(sig),gam/bet,del)
         if ( tau<rho-sig ) then
            tau = rho - sig
            inf = -1
         endif
         ! gaussian elimination
         a(kk) = rho
         kj = kk
         do j = k + 1 , n
            kj = kj + j - 1
            gam = a(kj)
            a(kj) = gam/rho
            ik = kk
            ij = kj
            do i = k + 1 , j
               ik = ik + i - 1
               ij = ij + 1
               a(ij) = a(ij) - a(ik)*gam
            end do
         end do
      end do
      if ( l>0 .and. abs(alf)>del ) inf = l
      end subroutine mxdpgf