formk – lbfgsb

private subroutine formk(n, Nsub, Ind, Nenter, Ileave, Indx2, Iupdat, Updatd, Wn, Wn1, m, Ws, Wy, Sy, Theta, Col, Head, Info)

This subroutine forms the LEL^T factorization of the indefinite matrix:

  K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
      [L_a -R_z           theta*S'AA'S ]

where:

  E = [-I  0]
      [ 0  I]

The matrix K can be shown to be equal to the matrix M^[-1]N occurring in section 5.1 of [1], as well as to the matrix Mbar^[-1] Nbar in section 5.3.

References

R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, "A limited memory algorithm for bound constrained optimization", SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, "L-BFGS-B: a limited memory FORTRAN code for solving bound constrained optimization problems", Tech. Report, NAM-11, EECS Department, Northwestern University, 1994.

Credits

NEOS, November 1994. (Latest revision June 1996.) Optimization Technology Center. Argonne National Laboratory and Northwestern University. Written by Ciyou Zhu in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.

Arguments

Type	Intent		Name
integer,	intent(in)	::	n	the dimension of the problem.
integer,	intent(in)	::	Nsub	the number of subspace variables in free set.
integer,	intent(in)	::	Ind(n)	specifies the indices of subspace variables.
integer,	intent(in)	::	Nenter	the number of variables entering the free set.
integer,	intent(in)	::	Ileave	`indx2(ileave),...,indx2(n)` are the variables leaving the free set.
integer,	intent(in)	::	Indx2(n)	On entry `indx2(1),...,indx2(nenter)` are the variables entering the free set, while `indx2(ileave),...,indx2(n)` are the variables leaving the free set.
integer,	intent(in)	::	Iupdat	the total number of BFGS updates made so far.
logical,	intent(in)	::	Updatd	true if the L-BFGS matrix is updatd.
real(kind=wp)		::	Wn(2m,2m)	On exit the upper triangle of `wn` stores the `LEL^T` factorization of the `2col x 2col` indefinite matrix `[-D -Y'ZZ'Y/theta L_a'-R_z' ] [L_a -R_z theta*S'AA'S ]`
real(kind=wp)		::	Wn1(2m,2m)	On entry `wn1` stores the lower triangular part of `[Y' ZZ'Y L_a'+R_z'] [L_a+R_z S'AA'S ]` in the previous iteration. On exit `wn1` stores the corresponding updated matrices. The purpose of `wn1` is just to store these inner products so they can be easily updated and inserted into wn.
integer,	intent(in)	::	m	the maximum number of variable metric corrections used to define the limited memory matrix.
real(kind=wp),	intent(in)	::	Ws(n,m)	information defining the limited memory BFGS matrix: `ws(n,m)` stores `S`, a set of s-vectors
real(kind=wp),	intent(in)	::	Wy(n,m)	information defining the limited memory BFGS matrix: `wy(n,m)` stores `Y`, a set of y-vectors
real(kind=wp),	intent(in)	::	Sy(m,m)	information defining the limited memory BFGS matrix: `sy(m,m)` stores `S'Y`
real(kind=wp),	intent(in)	::	Theta	information defining the limited memory BFGS matrix: the scaling factor specifying `B_0 = theta I`
integer,	intent(in)	::	Col	information defining the limited memory BFGS matrix: the number of variable metric corrections stored
integer,	intent(in)	::	Head	information defining the limited memory BFGS matrix: the location of the 1st s- (or y-) vector in S (or Y)
integer,	intent(out)	::	Info	On exit: `info = 0` for normal return; `info = -1` when the 1st Cholesky factorization failed; `info = -2` when the 2st Cholesky factorization failed.

Calls

Help

Called by

Help

Source Code

      subroutine formk(n,Nsub,Ind,Nenter,Ileave,Indx2,Iupdat,Updatd,Wn, &
                       Wn1,m,Ws,Wy,Sy,Theta,Col,Head,Info)
      implicit none

      integer,intent(in) :: n !! the dimension of the problem.
      integer,intent(in) :: Nsub !! the number of subspace variables in free set.
      integer,intent(in) :: m !! the maximum number of variable metric corrections
                              !! used to define the limited memory matrix.
      integer,intent(in) :: Nenter !! the number of variables entering the free set.
      integer,intent(in) :: Ileave !! `indx2(ileave),...,indx2(n)` are the variables leaving the free set.
      integer,intent(in) :: Iupdat !! the total number of BFGS updates made so far.
      integer,intent(out) :: Info !!  On exit:
                                  !!
                                  !!  * `info =  0` for normal return;
                                  !!  * `info = -1` when the 1st Cholesky factorization failed;
                                  !!  * `info = -2` when the 2st Cholesky factorization failed.
      integer,intent(in) :: Ind(n) !! specifies the indices of subspace variables.
      integer,intent(in) :: Indx2(n) !! On entry `indx2(1),...,indx2(nenter)` are the variables entering
                                     !! the free set, while `indx2(ileave),...,indx2(n)` are the
                                     !! variables leaving the free set.
      real(wp),intent(in) :: Ws(n,m) !! information defining the limited memory BFGS matrix: `ws(n,m)` stores `S`, a set of s-vectors
      real(wp),intent(in) :: Wy(n,m) !! information defining the limited memory BFGS matrix: `wy(n,m)` stores `Y`, a set of y-vectors
      real(wp),intent(in) :: Sy(m,m) !! information defining the limited memory BFGS matrix: `sy(m,m)` stores `S'Y`
      real(wp),intent(in) :: Theta !! information defining the limited memory BFGS matrix: the scaling factor specifying `B_0 = theta I`
      integer,intent(in) :: Col !! information defining the limited memory BFGS matrix: the number of variable metric corrections stored
      integer,intent(in) :: Head !! information defining the limited memory BFGS matrix: the location of the 1st s- (or y-) vector in S (or Y)
      real(wp) :: Wn(2*m,2*m) !! On exit the upper triangle of `wn` stores the `LEL^T` factorization
                              !! of the `2*col x 2*col` indefinite matrix
                              !!```
                              !!  [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
                              !!  [L_a -R_z           theta*S'AA'S ]
                              !!```
      real(wp) :: Wn1(2*m,2*m) !! On entry `wn1` stores the lower triangular part of
                               !!```
                               !!  [Y' ZZ'Y   L_a'+R_z']
                               !!  [L_a+R_z   S'AA'S   ]
                               !!```
                               !! in the previous iteration.
                               !!
                               !! On exit `wn1` stores the corresponding updated matrices.
                               !! The purpose of `wn1` is just to store these inner products
                               !! so they can be easily updated and inserted into wn.
      logical,intent(in) :: Updatd !! true if the L-BFGS matrix is updatd.

      integer :: m2 , ipntr , jpntr , iy , is , jy , js , is1 , js1 , k1 , &
                 i , k , col2 , pbegin , pend , dbegin , dend , upcl
      real(wp) :: temp1 , temp2 , temp3 , temp4

      ! Form the lower triangular part of
      !           WN1 = [Y' ZZ'Y   L_a'+R_z']
      !                 [L_a+R_z   S'AA'S   ]
      !    where L_a is the strictly lower triangular part of S'AA'Y
      !          R_z is the upper triangular part of S'ZZ'Y.

      if ( Updatd ) then
         if ( Iupdat>m ) then
            ! shift old part of WN1.
            do jy = 1 , m - 1
               js = m + jy
               call dcopy(m-jy,Wn1(jy+1,jy+1),1,Wn1(jy,jy),1)
               call dcopy(m-jy,Wn1(js+1,js+1),1,Wn1(js,js),1)
               call dcopy(m-1,Wn1(m+2,jy+1),1,Wn1(m+1,jy),1)
            enddo
         endif

         ! put new rows in blocks (1,1), (2,1) and (2,2).
         pbegin = 1
         pend = Nsub
         dbegin = Nsub + 1
         dend = n
         iy = Col
         is = m + Col
         ipntr = Head + Col - 1
         if ( ipntr>m ) ipntr = ipntr - m
         jpntr = Head
         do jy = 1 , Col
            js = m + jy
            temp1 = zero
            temp2 = zero
            temp3 = zero
            ! compute element jy of row 'col' of Y'ZZ'Y
            do k = pbegin , pend
               k1 = Ind(k)
               temp1 = temp1 + Wy(k1,ipntr)*Wy(k1,jpntr)
            enddo
            ! compute elements jy of row 'col' of L_a and S'AA'S
            do k = dbegin , dend
               k1 = Ind(k)
               temp2 = temp2 + Ws(k1,ipntr)*Ws(k1,jpntr)
               temp3 = temp3 + Ws(k1,ipntr)*Wy(k1,jpntr)
            enddo
            Wn1(iy,jy) = temp1
            Wn1(is,js) = temp2
            Wn1(is,jy) = temp3
            jpntr = mod(jpntr,m) + 1
         enddo

         ! put new column in block (2,1).
         jy = Col
         jpntr = Head + Col - 1
         if ( jpntr>m ) jpntr = jpntr - m
         ipntr = Head
         do i = 1 , Col
            is = m + i
            temp3 = zero
            ! compute element i of column 'col' of R_z
            do k = pbegin , pend
               k1 = Ind(k)
               temp3 = temp3 + Ws(k1,ipntr)*Wy(k1,jpntr)
            enddo
            ipntr = mod(ipntr,m) + 1
            Wn1(is,jy) = temp3
         enddo
         upcl = Col - 1
      else
         upcl = Col
      endif

      ! modify the old parts in blocks (1,1) and (2,2) due to changes
      ! in the set of free variables.
      ipntr = Head
      do iy = 1 , upcl
         is = m + iy
         jpntr = Head
         do jy = 1 , iy
            js = m + jy
            temp1 = zero
            temp2 = zero
            temp3 = zero
            temp4 = zero
            do k = 1 , Nenter
               k1 = Indx2(k)
               temp1 = temp1 + Wy(k1,ipntr)*Wy(k1,jpntr)
               temp2 = temp2 + Ws(k1,ipntr)*Ws(k1,jpntr)
            enddo
            do k = Ileave , n
               k1 = Indx2(k)
               temp3 = temp3 + Wy(k1,ipntr)*Wy(k1,jpntr)
               temp4 = temp4 + Ws(k1,ipntr)*Ws(k1,jpntr)
            enddo
            Wn1(iy,jy) = Wn1(iy,jy) + temp1 - temp3
            Wn1(is,js) = Wn1(is,js) - temp2 + temp4
            jpntr = mod(jpntr,m) + 1
         enddo
         ipntr = mod(ipntr,m) + 1
      enddo

      ! modify the old parts in block (2,1).
      ipntr = Head
      do is = m + 1 , m + upcl
         jpntr = Head
         do jy = 1 , upcl
            temp1 = zero
            temp3 = zero
            do k = 1 , Nenter
               k1 = Indx2(k)
               temp1 = temp1 + Ws(k1,ipntr)*Wy(k1,jpntr)
            enddo
            do k = Ileave , n
               k1 = Indx2(k)
               temp3 = temp3 + Ws(k1,ipntr)*Wy(k1,jpntr)
            enddo
            if ( is<=jy+m ) then
               Wn1(is,jy) = Wn1(is,jy) + temp1 - temp3
            else
               Wn1(is,jy) = Wn1(is,jy) - temp1 + temp3
            endif
            jpntr = mod(jpntr,m) + 1
         enddo
         ipntr = mod(ipntr,m) + 1
      enddo

      ! Form the upper triangle of WN = [D+Y' ZZ'Y/theta   -L_a'+R_z' ]
      !                                 [-L_a +R_z        S'AA'S*theta]

      m2 = 2*m
      do iy = 1 , Col
         is = Col + iy
         is1 = m + iy
         do jy = 1 , iy
            js = Col + jy
            js1 = m + jy
            Wn(jy,iy) = Wn1(iy,jy)/Theta
            Wn(js,is) = Wn1(is1,js1)*Theta
         enddo
         do jy = 1 , iy - 1
            Wn(jy,is) = -Wn1(is1,jy)
         enddo
         do jy = iy , Col
            Wn(jy,is) = Wn1(is1,jy)
         enddo
         Wn(iy,iy) = Wn(iy,iy) + Sy(iy,iy)
      enddo

      ! Form the upper triangle of WN= [  LL'            L^-1(-L_a'+R_z')]
      !                                [(-L_a +R_z)L'^-1   S'AA'S*theta  ]

      ! first Cholesky factor (1,1) block of wn to get LL'
      ! with L' stored in the upper triangle of wn.
      call dpofa(Wn,m2,Col,Info)
      if ( Info/=0 ) then
         Info = -1
         return
      endif
      ! then form L^-1(-L_a'+R_z') in the (1,2) block.
      col2 = 2*Col
      do js = Col + 1 , col2
         call dtrsl(Wn,m2,Col,Wn(1,js),11,Info)
      enddo

      ! Form S'AA'S*theta + (L^-1(-L_a'+R_z'))'L^-1(-L_a'+R_z') in the
      ! upper triangle of (2,2) block of wn.

      do is = Col + 1 , col2
         do js = is , col2
            Wn(is,js) = Wn(is,js) + ddot(Col,Wn(1,is),1,Wn(1,js),1)
         enddo
      enddo

      ! Cholesky factorization of (2,2) block of wn.

      call dpofa(Wn(Col+1,Col+1),m2,Col,Info)
      if ( Info/=0 ) then
         Info = -2
         return
      endif

      end subroutine formk

formk Subroutine

Contents

Source Code

private subroutine formk(n, Nsub, Ind, Nenter, Ileave, Indx2, Iupdat, Updatd, Wn, Wn1, m, Ws, Wy, Sy, Theta, Col, Head, Info)

References

Credits

Arguments

Calls

Called by

Source Code