elec2.f90

Go to the documentation of this file.

!> @file elec2.f90
!! @ingroup FORT1THREAD_EXAMPLES
!!
!! 
!! Electron model from COPS test set, this time with a Hessian.
!! 
!! The hessian for this model is extremely dense and in the first solve
!! CONOPT will not use it.
!! During the second solve we redefine HessFac so the hessian is used.
!! 
!! Since the model is non-convex we will several times experience
!! negative curvature in the SQP sub-solver. With the correct Hessian
!! this seems to be more of a problems and the extra information is
!! actually counter-productive -- not because it takes a long time to
!! evaluate the Hessian but because the use takes time and the number
!! of iterations grow.
!! 
!!
!!  This is a CONOPT implementation of the GAMS model:
!!
!!
!! @verbatim
!! Set i       electrons /i1 * i%np%/
!!     ut(i,i) upper triangular part;
!! 
!! Alias (i,j);
!! ut(i,j)$(ord(j) > ord(i)) = yes;
!! 
!! Variables   x(i) x-coordinate of the electron
!!             y(i) y-coordinate of the electron
!!             z(i) z-coordinate of the electron
!!             potential Coulomb potential;
!! 
!! Equations   obj     objective
!!             ball(i) points on unit ball;
!! 
!! obj.. potential =e=
!!   sum{ut(i,j), 1.0/sqrt(sqr(x[i]-x[j]) + sqr(y[i]-y[j]) + sqr(z[i]-z[j]))};
!! 
!! ball(i).. sqr(x(i)) + sqr(y(i)) + sqr(z(i)) =e= 1;
!! 
!! 
!! * Set the starting point to a quasi-uniform distribution
!! * of electrons on a unit sphere
!! 
!! scalar pi a famous constant;
!! pi = 2*arctan(inf);
!! 
!! parameter theta(i), phi(i);
!! theta(i) = 2*pi*uniform(0,1);
!! phi(i)   = pi*uniform(0,1);
!! 
!! x.l(i) = cos(theta(i))*sin(phi(i));
!! y.l(i) = sin(theta(i))*sin(phi(i));
!! z.l(i) =               cos(phi(i));
!! @endverbatim
!!
!!
!! For more information about the individual callbacks, please have a look at the source code.
 
#if defined(_WIN32) && !defined(_WIN64)
#define dec_directives_win32
#endif
 
Module hesscalls
   Integer :: Hcalls
End Module hesscalls
 
Module random
   Integer :: seed
End Module random
 
!> Main program. A simple setup and call of CONOPT
!!
Program tutorial
 
   Use proginfo
   Use conopt
   Use hesscalls
   Use random
   implicit None
!
!  Declare the user callback routines as Integer, External:
!
   Integer, External :: elec_readmatrix ! Mandatory Matrix definition routine defined below
   Integer, External :: elec_fdeval     ! Function and Derivative evaluation routine
                                       ! needed a nonlinear model.
   Integer, External :: std_status     ! Standard callback for displaying solution status
   Integer, External :: std_solution   ! Standard callback for displaying solution values
   Integer, External :: std_message    ! Standard callback for managing messages
   Integer, External :: std_errmsg     ! Standard callback for managing error messages
   Integer, External :: elec_2dlagrsize ! Callback for second derivatives size
   Integer, External :: elec_2dlagrstr  ! Callback for second derivatives structure
   Integer, External :: elec_2dlagrval  ! Callback for second derivatives values
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Elec_ReadMatrix
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Elec_FDEval
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_Status
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_Solution
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_Message
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_ErrMsg
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Elec_2DLagrSize
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Elec_2DLagrStr
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Elec_2DLagrVal
#endif
!
!  Control vector
!
   INTEGER, Dimension(:), Pointer :: cntvect
   INTEGER :: coi_error
   Integer :: ne ! Number of Electrons
!
!  Create and initialize a Control Vector
!
   call startup
 
   coi_error = coi_create( cntvect )
!
!  Tell CONOPT about the size of the model by populating the Control Vector:
!
   ne = 150 ! Number of electrons
   coi_error = max( coi_error, coidef_numvar( cntvect, 3*ne ) ) ! # variables
   coi_error = max( coi_error, coidef_numcon( cntvect, ne+1 ) ) ! # constraints
   coi_error = max( coi_error, coidef_numnz( cntvect, 6*ne ) ) ! # nonzeros in the Jacobian
   coi_error = max( coi_error, coidef_numnlnz( cntvect, 6*ne ) ) ! # of which are nonlinear
   coi_error = max( coi_error, coidef_optdir( cntvect, -1 ) )   ! Minimize
   coi_error = max( coi_error, coidef_objcon( cntvect, ne+1 ) ) ! Objective is constraint NE+1
   coi_error = max( coi_error, coidef_debug2d( cntvect, 0 ) )    ! turn 2nd derivative debugging on or off
   coi_error = max( coi_error, coidef_optfile( cntvect, 'elec2.opt'  ) )
!
!  Tell CONOPT about the callback routines:
!
   coi_error = max( coi_error, coidef_readmatrix( cntvect, elec_readmatrix ) )
   coi_error = max( coi_error, coidef_fdeval( cntvect, elec_fdeval     ) )
   coi_error = max( coi_error, coidef_status( cntvect, std_status     ) )
   coi_error = max( coi_error, coidef_solution( cntvect, std_solution   ) )
   coi_error = max( coi_error, coidef_message( cntvect, std_message    ) )
   coi_error = max( coi_error, coidef_errmsg( cntvect, std_errmsg     ) )
   coi_error = max( coi_error, coidef_2dlagrsize( cntvect, elec_2dlagrsize ) )
   coi_error = max( coi_error, coidef_2dlagrstr( cntvect, elec_2dlagrstr ) )
   coi_error = max( coi_error, coidef_2dlagrval( cntvect, elec_2dlagrval ) )
 
#if defined(CONOPT_LICENSE_INT_1) && defined(CONOPT_LICENSE_INT_2) && defined(CONOPT_LICENSE_INT_3) && defined(CONOPT_LICENSE_TEXT)
   coi_error = max( coi_error, coidef_license( cntvect, conopt_license_int_1, conopt_license_int_2, conopt_license_int_3, conopt_license_text) )
#endif
 
   If ( coi_error .ne. 0 ) THEN
      write(*,*)
      write(*,*) '**** Fatal Error while loading CONOPT Callback routines.'
      write(*,*)
      call flog( "Skipping Solve due to setup errors", 1 )
   ENDIF
!
!  Start CONOPT:
!
   coi_error = max( coi_error, coidef_hessfac( cntvect, 1.0d0 ) ) ! very little space to prevent Hessian to be used
   hcalls = 0
   seed = 12345
   coi_error = coi_solve( cntvect )
   If ( coi_error /= 0 ) then
      call flog( "Solve 1: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Solve 1: Status or Solution routine was not called", 1 )
   elseif ( sstat /= 1 .or. mstat /= 2 ) then
      call flog( "Solve 1: Solver and Model Status was not as expected (1,2)", 1 )
   elseif ( hcalls /= 0 ) then
      call flog( "Solve 1: Solve used 2nd derivatives but should reject on memory", 1 )
   endif
!
!  Start CONOPT again with more space for the Hessian of the Lagrangian
!
   coi_error = max( coi_error, coidef_hessfac( cntvect, 0.0d0 ) ) ! enough to allow Hessian to be used
   hcalls = 0
   seed = 12345 ! Make sure we get the same model the second time
   coi_error = coi_solve( cntvect )
   If ( coi_error /= 0 ) then
      call flog( "Solve 2: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Solve 2: Status or Solution routine was not called", 1 )
   elseif ( sstat /= 1 .or. mstat /= 2 ) then
      call flog( "Solve 2: Solver and Model Status was not as expected (1,2)", 1 )
   elseif ( hcalls == 0 ) then
      call flog( "Solve 2: Solve did not use 2nd derivatives but should have enough memory", 1 )
   endif
 
 
   write(*,*)
   write(*,*) 'End of Electron example. Return code=',coi_error
 
   if ( coi_free(cntvect) /= 0 ) call flog( "Error while freeing control vector",1)
 
   call flog( "Successful Solve", 0 )
 
End Program tutorial
 
REAL FUNCTION rndx( )
!
!  Defines a pseudo random number between 0 and 1
!
   Use random
   IMPLICIT NONE
 
   seed = mod(seed*1027+25,1048576)
   rndx = float(seed)/float(1048576)
 
END FUNCTION rndx
!
! ============================================================================
!  Define information about the model:
!
 
!>  Define information about the model
!!
!! @include{doc} readMatrix_params.dox
Integer Function elec_readmatrix( lower, curr, upper, vsta, type, rhs, esta,      &
                                 colsta, rowno, value, nlflag, n, m, nz,  &
                                 usrmem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Elec_ReadMatrix
#endif
   implicit none
   integer, intent (in)                      :: n      ! number of variables
   integer, intent (in)                      :: m      ! number of constraints
   integer, intent (in)                      :: nz     ! number of nonzeros
   real*8,  intent (in out), dimension(n)    :: lower  ! vector of lower bounds
   real*8,  intent (in out), dimension(n)    :: curr   ! vector of initial values
   real*8,  intent (in out), dimension(n)    :: upper  ! vector of upper bounds
   integer, intent (in out), dimension(n)    :: vsta   ! vector of initial variable status
                                                       ! (not defined here)
   integer, intent (out),    dimension(m)    :: type   ! vector of equation types
   integer, intent (in out), dimension(m)    :: esta   ! vector of initial equation status
                                                       ! (not defined here)
   real*8,  intent (in out), dimension(m)    :: rhs    ! vector of right hand sides
   integer, intent (in out), dimension(n+1)  :: colsta ! vector with start of column indices
   integer, intent (out),    dimension(nz)   :: rowno  ! vector of row numbers
   integer, intent (in out), dimension(nz)   :: nlflag ! vector of nonlinearity flags
   real*8,  intent (in out), dimension(nz)   :: value  ! vector of matrix values
   real*8   usrmem(*)                                  ! optional user memory
 
   Integer :: ne
   Integer :: i, k
   real*8, parameter :: pi = 3.141592
   real*8 :: theta, phi
   Real, External :: rndx
 
   ne = n / 3
!
!  Information about Variables:
!  Default: Lower = -Inf, Curr = 0, and Upper = +inf.
!  Default: the status information in Vsta is not used.
!
   DO i = 1, ne
      theta = rndx()
      phi   = rndx()
      curr(i     ) = cos(theta)*sin(phi)
      curr(i+  ne) = sin(theta)*sin(phi)
      curr(i+2*ne) =            cos(phi)
   enddo
!
!  Information about Constraints:
!  Default: Rhs = 0
!  Default: the status information in Esta and the function
!           value in FV are not used.
!  Default: Type: There is no default.
!           0 = Equality,
!           1 = Greater than or equal,
!           2 = Less than or equal,
!           3 = Non binding.
!
!  Constraint 1 (Objective)
!  Rhs = 0.0 and type Non binding
!
   DO i = 1, ne
      Type(i) = 0
      rhs(i)  = 1.d0
   Enddo
   type(ne+1) = 3
!
!  Information about the Jacobian. We have to define Rowno, Value,
!  Nlflag and Colsta.
!
!  Colsta = Start of column indices (No Defaults):
!  Rowno  = Row indices
!  Value  = Value of derivative (by default only linear
!                                derivatives are used)
!  Nlflag = 0 for linear and 1 for nonlinear derivative
!           (not needed for completely linear models)
!
!  Indices
!       x(1)  x(2)  .. y(1)   y(2) ..  z(1)   z(2) ..
!  1:    1            2*Ne+1          3*Ne+1
!  2:          3             2*Ne+3          3*Ne+3
!  ..
!  Obj:  2     4      2*Ne+2 2*Ne+4   3*Ne+2 3*Ne+4
!
!  Nonlinearity Structure: L = 0 are linear and NL = 1 are nonlinear
!  All nonzeros are nonlinear
!
   DO i = 1, n+1
      colsta(i) = 2*i-1
   enddo
   k = 0
   DO i = 1, ne         ! x variablex
      k = k + 1
      rowno(k)  = i    ! Ball constraint
      nlflag(k) = 1
      k = k + 1
      rowno(k)  = ne+1 ! Objective
      nlflag(k) = 1
   enddo
   DO i = 1, ne         ! y variablex
      k = k + 1
      rowno(k)  = i    ! Ball constraint
      nlflag(k) = 1
      k = k + 1
      rowno(k)  = ne+1 ! Objective
      nlflag(k) = 1
   enddo
   DO i = 1, ne         ! z variablex
      k = k + 1
      rowno(k)  = i    ! Ball constraint
      nlflag(k) = 1
      k = k + 1
      rowno(k)  = ne+1 ! Objective
      nlflag(k) = 1
   enddo
 
   elec_readmatrix = 0 ! Return value means OK
 
end Function elec_readmatrix
!
!==========================================================================
!  Compute nonlinear terms and non-constant Jacobian elements
!
 
!>  Compute nonlinear terms and non-constant Jacobian elements
!!
!! @include{doc} fdeval_params.dox
Integer Function elec_fdeval( x, g, jac, rowno, jcnm, mode, ignerr, errcnt, &
                             n, nz, thread, usrmem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Elec_FDEval
#endif
   implicit none
   integer, intent (in)                      :: n      ! number of variables
   integer, intent (in)                      :: rowno  ! number of the row to be evaluated
   integer, intent (in)                      :: nz     ! number of nonzeros in this row
   real*8,  intent (in),     dimension(n)    :: x      ! vector of current solution values
   real*8,  intent (in out)                  :: g      ! constraint value
   real*8,  intent (in out), dimension(n)    :: jac    ! vector of derivatives for current constraint
   integer, intent (in),     dimension(nz)   :: jcnm   ! list of variables that appear nonlinearly
                                                       ! in this row. Ffor information only.
   integer, intent (in)                      :: mode   ! evaluation mode: 1 = function value
                                                       ! 2 = derivatives, 3 = both
   integer, intent (in)                      :: ignerr ! if 1 then errors can be ignored as long
                                                       ! as errcnt is incremented
   integer, intent (in out)                  :: errcnt ! error counter to be incremented in case
                                                       ! of function evaluation errors.
   integer, intent (in)                      :: thread
   real*8   usrmem(*)                                  ! optional user memory
 
   Integer :: ne
   Integer :: i, j
   real*8 :: dist, dist32, tx, ty, tz
 
   ne = n / 3
!
!  Row NE+1: the objective function is nonlinear
!
   if ( rowno == ne+1 ) then
!
!  Mode = 1 or 3. Function value:
!
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         g = 0.d0
         do i = 1, ne-1
            do j = i+1,ne
               dist = (x(i)-x(j))**2 + (x(i+ne)-x(j+ne))**2 + (x(i+2*ne)-x(j+2*ne))**2
               g    = g + 1.d0/sqrt(dist)
            enddo
         enddo
      endif
!
!  Mode = 2 or 3: Derivative values: w.r.t. x: 2*(x(i)-x(j))*(-1/2)*dist(-3/2)
!
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         do i = 1, 3*ne
            jac(i) = 0.d0
         enddo
         do i = 1, ne-1
            do j = i+1,ne
               dist = (x(i)-x(j))**2 + (x(i+ne)-x(j+ne))**2 + (x(i+2*ne)-x(j+2*ne))**2
               dist32 = (1.d0/sqrt(dist))**3
               tx = -(x(i)-x(j))*dist32
               ty = -(x(i+ne)-x(j+ne))*dist32
               tz = -(x(i+2*ne)-x(j+2*ne))*dist32
               jac(i) = jac(i) + tx
               jac(j) = jac(j) - tx
               jac(i+ne) = jac(i+ne) + ty
               jac(j+ne) = jac(j+ne) - ty
               jac(i+2*ne) = jac(i+2*ne) + tz
               jac(j+2*ne) = jac(j+2*ne) - tz
            enddo
         enddo
      endif
!
   Else  ! this is ball constraint rowno
!
!  Mode = 1 or 3. Function value:
!
      i = rowno
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         g = x(i)**2 + x(i+ne)**2 + x(i+2*ne)**2
      endif
!
!  Mode = 2 or 3: Derivative values:
!
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         jac(i)      = 2.d0*x(i)
         jac(i+ne)   = 2.d0*x(i+ne)
         jac(i+2*ne) = 2.d0*x(i+2*ne)
      endif
 
   endif
   elec_fdeval = 0
 
end Function elec_fdeval
 
Integer Function elec_2dlagrsize( NODRV, NumVar, NumCon, NumHess, MaxHess, UsrMem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Elec_2DLagrSize
#endif
   Use hesscalls
   Implicit None
 
   Integer, Intent (IN)               :: numvar, numcon, maxhess
   Integer, Intent (IN OUT)           :: nodrv, numhess
   real*8,  Intent(IN OUT)            :: usrmem(*)
!
!  Sizing call -- the hessian is completely dense
!
   numhess = numvar * (numvar+1)/2
   elec_2dlagrsize = 0
 
End Function elec_2dlagrsize
 
 
!>  Specify the structure of the Lagrangian of the Hessian
!!
!! @include{doc} 2DLagrStr_params.dox
Integer Function elec_2dlagrstr( HSRW, HSCL, NODRV, NumVar, NumCon, NumHess, UsrMem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Elec_2DLagrStr
#endif
   Use hesscalls
   Implicit None
 
   Integer, Intent (IN)               :: numvar, numcon, numhess
   Integer, Intent (IN)               :: nodrv
   Integer, Dimension(*)              :: hsrw, hscl
   real*8,  Intent(IN OUT)            :: usrmem(*)
 
   integer           :: i, j, k
!
!  Define structure of Hessian: A dense lower-triangular matrix
!
   k = 0
   do i = 1, numvar
      do j = i, numvar
         k = k + 1
         hsrw(k) = j
         hscl(k) = i
      enddo
   enddo
   elec_2dlagrstr = 0
 
End Function elec_2dlagrstr
 
 
!>  Compute the Lagrangian of the Hessian
!!
!! @include{doc} 2DLagrVal_params.dox
Integer Function elec_2dlagrval( X, U, HSRW, HSCL, HSVL, NODRV, NumVar, NumCon, NumHess, UsrMem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Elec_2DLagrVal
#endif
   Use hesscalls
   Implicit None
 
   Integer, Intent (IN)                      :: numvar, numcon, numhess
   Integer, Intent (IN OUT)                  :: nodrv
   real*8,  Dimension(NumVar), Intent (IN)   :: x
   real*8,  Dimension(NumCon), Intent (IN)   :: u
   Integer, Dimension(NumHess), Intent (IN)  :: hsrw, hscl
   real*8,  Dimension(NumHess), Intent (Out) :: hsvl
   real*8,  Intent(IN OUT)                   :: usrmem(*)
 
   Integer :: ne ! Number of electrons
   integer :: k, l, ix, iy, iz, jx, jy, jz
   real*8 :: tx, ty, tz, dist, distsq, dist3, dist5
!
!  Evaluation of 2nd derivatives. Hsvl is zero on entry:
!
   ne = numvar / 3
   hcalls = hcalls + 1
   k = 0
!   do i = 1, N
!      do j = i, N
!         k = k + 1
!         Hsvl(k) = 0.d0
!      enddo
!   enddo
   do l = 1, ne ! constraint no l with 3 elements on the diagonal, each with value 2.0
      k       = pos(l,l)           ! 1+(l-1)*(2*N+2-l)/2 ! variable x(l)
      hsvl(k) = hsvl(k) + 2.d0*u(l)
      k       = pos(l+ne,l+ne)     ! 1+(l+NE-1)*(2*N+2-(l+NE))/2 ! variable y(l)
      hsvl(k) = hsvl(k) + 2.d0*u(l)
      k       = pos(l+2*ne,l+2*ne) ! 1+(l+2*NE-1)*(2*N+2-(l+2*NE))/2 ! variable z(l)
      hsvl(k) = hsvl(k) + 2.d0*u(l)
   enddo
   if ( u(ne+1) /= 0.0d0 ) Then ! Nonzero dual on the objective
      do ix = 1, ne-1
         iy = ix + ne
         iz = iy + ne
         do jx = ix+1, ne
            jy = jx + ne
            jz = jy + ne
!  Objective computation:
!           dist = (x(i)-x(j))**2 + (x(i+NE)-x(j+NE))**2 + (x(i+2*NE)-x(j+2*NE))**2
!           g    = g + 1.d0/sqrt(dist)
            tx = x(ix)-x(jx) ! x1 - x2
            ty = x(iy)-x(jy) ! y1 - y2
            tz = x(iz)-x(jz) ! z1 - z2
! The first derivatives of tx, ty, and tz w.r.t. the x variables are +1 or -1
! The second derivatives of tx, ty, and tz are zero
            dist   = tx**2 + ty**2 + tz**2
! The first derivatives of dist are
!           d(dist)/dx = d(dist)/d(tx)*d(tx)/dx = 2*tx*(+/-)   (+1 for x1 and -1 for x2)
! and similarly for d(dist)/dy = 2*ty*(+/-1) and d(dist)/dz = 2*tz*(+/-1)
! The second derivatives of dist are
!           d2(dist)/dxidxj = 2 if xi=xj and 0 otherwise
!           d2(dist)/dyidyj = 2 if yi=yj and 0 otherwise
!           d2(dist)/dzidzj = 2 if zi=zj and 0 otherwise
! All cross 2nd derivatives dxdy, dxdz, and dydz are zero
!
            distsq  = 1/sqrt(dist) ! = dist**(-0.5)
            dist3  = distsq**3
            dist5  = distsq**5
! The first derivatives of distsq are
!           d(distsq)/dx = -0.5*dist**(-1.5)*d(dist)/dx = -dist**(-1.5)*tx*(+/-1) = -distsq**3*tx*(+/-1)
!           d(distsq)/dy = -0.5*dist**(-1.5)*d(dist)/dy = -dist**(-1.5)*ty*(+/-1) = -distsq**3*ty*(+/-1)
!           d(distsq)/dz = -0.5*dist**(-1.5)*d(dist)/dz = -dist**(-1.5)*tz*(+/-1) = -distsq**3*tz*(+/-1)
! The second derivatives of distsq are
!           d2(distsq)/dxidxj = -1.5*dist**(-2.5)*d(dist)/dxj*tx*(+/-1) - dist**(-1.5)*d(tx)/dxj*(+/-1)
!                             = -1.5*distsq**5*2*tx*tx*(+/-1) - distsq**3*(+/-1)
!           d2(distsq)/dxidyj = -1.5*dist**(-2.5)*d(dist)/dyj*tx*(+/-1) - dist**(-1.5)*d(tx)/dyj*(+/-1)
!                             = -1.5*distsq**5*2*ty*tx*(+/-1)
            k = pos(ix,ix)
            hsvl(k) = hsvl(k) + 3.0d0*dist5*tx**2 - dist3
            k = pos(jx,jx)
            hsvl(k) = hsvl(k) + 3.0d0*dist5*tx**2 - dist3
            k = pos(ix,jx)
            hsvl(k) = hsvl(k) - 3.0d0*dist5*tx**2 + dist3
            k = pos(iy,iy)
            hsvl(k) = hsvl(k) + 3.0d0*dist5*ty**2 - dist3
            k = pos(jy,jy)
            hsvl(k) = hsvl(k) + 3.0d0*dist5*ty**2 - dist3
            k = pos(iy,jy)
            hsvl(k) = hsvl(k) - 3.0d0*dist5*ty**2 + dist3
            k = pos(iz,iz)
            hsvl(k) = hsvl(k) + 3.0d0*dist5*tz**2 - dist3
            k = pos(jz,jz)
            hsvl(k) = hsvl(k) + 3.0d0*dist5*tz**2 - dist3
            k = pos(iz,jz)
            hsvl(k) = hsvl(k) - 3.0d0*dist5*tz**2 + dist3
!
!  Off-diagonal blocks
!
            k = pos(ix,iy)
            hsvl(k) = hsvl(k) + 3.0d0*dist5*tx*ty
            k = pos(ix,jy)
            hsvl(k) = hsvl(k) - 3.0d0*dist5*tx*ty
            k = pos(jx,iy)
            hsvl(k) = hsvl(k) - 3.0d0*dist5*tx*ty
            k = pos(jx,jy)
            hsvl(k) = hsvl(k) + 3.0d0*dist5*tx*ty
            k = pos(ix,iz)
            hsvl(k) = hsvl(k) + 3.0d0*dist5*tx*tz
            k = pos(ix,jz)
            hsvl(k) = hsvl(k) - 3.0d0*dist5*tx*tz
            k = pos(jx,iz)
            hsvl(k) = hsvl(k) - 3.0d0*dist5*tx*tz
            k = pos(jx,jz)
            hsvl(k) = hsvl(k) + 3.0d0*dist5*tx*tz
            k = pos(iy,iz)
            hsvl(k) = hsvl(k) + 3.0d0*dist5*ty*tz
            k = pos(iy,jz)
            hsvl(k) = hsvl(k) - 3.0d0*dist5*ty*tz
            k = pos(jy,iz)
            hsvl(k) = hsvl(k) - 3.0d0*dist5*ty*tz
            k = pos(jy,jz)
            hsvl(k) = hsvl(k) + 3.0d0*dist5*ty*tz
         enddo
      enddo
   endif
 
   elec_2dlagrval = 0
 
   Contains
 
   Integer function pos(i,j)
   Integer, Intent(in) :: i, j
   pos = 1+(i-1)*(2*numvar+2-i)/2 + (j-i)  ! Hessian position for variable (j,i), j >= i
   End function pos
 
End Function elec_2dlagrval