mp_elec2.f90

Go to the documentation of this file.

!> @file mp_elec2.f90
!! @ingroup FORTOPENMP_EXAMPLES
!!
!! @copydoc elec2.f90
 
#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 electron
 
   Use proginfop
   Use conopt
   Use hesscalls
   Use omp_lib
   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
!
!  The model shown here in GAMS format
!
!
!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));
!
   Integer :: ne ! Number of Electrons
 
   real*8 time0, time1, time2
!
!  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, 'mp_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 with enough space for the Hessian of the Lagrangian -- first time single thread
!
   coi_error = max( coi_error, coidef_hessfac( cntvect, 3.01d0*(ne+1)/4.d0    ) ) ! enough to allow Hessian to be use
   hcalls = 0
   seed = 12345
   time0 = omp_get_wtime()
   coi_error = coi_solve( cntvect )
   time1 = omp_get_wtime() - time0
 
   write(*,*)
   write(*,*) 'End of Electron example one thread. Return code=',coi_error
 
   If ( coi_error /= 0 ) then
      call flog( "One Thread: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "One Thread: Status or Solution routine was not called", 1 )
   elseif ( sstat /= 1 .or. mstat /= 2 ) then
      call flog( "One Thread: Solver and Model Status was not as expected (1,2)", 1 )
   elseif ( hcalls == 0 ) then
      call flog( "One Thread: Solve did not use 2nd derivatives but should have enough memory", 1 )
   endif
!
!  Start CONOPT again -- this time using multiple threads:
!
   seed = 12345 ! Make sure we get the same model the second time
   coi_error = max( coi_error, coidef_threads( cntvect, 0 ) )  ! 0 means use as many threads as you can
   time0 = omp_get_wtime()
   coi_error = coi_solve( cntvect )
   time2 = omp_get_wtime() - time0
 
   write(*,*)
   write(*,*) 'End of Electron example multi threads. Return code=',coi_error
 
   If ( coi_error /= 0 ) then
      call flog( "Multi thread: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Multi thread: Status or Solution routine was not called", 1 )
   elseif ( sstat /= 1 .or. mstat /= 2 ) then
      call flog( "Multi thread: Solver and Model Status was not as expected (1,2)", 1 )
   elseif ( hcalls == 0 ) then
      call flog( "Multi thread: Solve did not use 2nd derivatives but should have enough memory", 1 )
   endif
 
   if ( coi_free( cntvect ) /= 0 ) call flog( "Error while freeing control vector", 1 )
 
   write(*,*)
   write(*,"('Time for single thread',f10.3)") time1
   write(*,"('Time for multi  thread',f10.3)") time2
   write(*,"('Speedup               ',f10.3)") time1/time2
   write(*,"('Efficiency            ',f10.3)") time1/time2/omp_get_max_threads()
 
   call flog( "Successful Solve", 0 )
 
End Program electron
 
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