C SAFE CALL TO THE USER PROVIDED SUBROUTINES C ****************************************************************** C ****************************************************************** subroutine vinip(n,x,l,u,m,lambda,equatn,linear,coded,checkder, +inform) implicit none C SCALAR ARGUMENTS logical checkder integer inform,m,n C ARRAY ARGUMENTS logical coded(6),equatn(m),linear(m) double precision l(n),lambda(m),u(n),x(n) #include "dim.par" #include "algparam.com" C LOCAL SCALARS integer i C Avoid huge bounds do i = 1,n l(i) = max( l(i), - 1.0d+20 ) u(i) = min( u(i), 1.0d+20 ) end do C Project initial guess do i = 1,n x(i) = max( l(i), min( x(i), u(i) ) ) end do C Check derivatives if ( checkder ) then call checkd(n,l,u,m,inform) if ( inform .lt. 0 ) return end if end C ****************************************************************** C ****************************************************************** subroutine vendp(n,x,l,u,m,lambda,equatn,linear,inform) implicit none C SCALAR ARGUMENTS integer inform,m,n C ARRAY ARGUMENTS logical equatn(m),linear(m) double precision l(n),lambda(m),u(n),x(n) C LOCAL SCALARS integer i C Save solution open(10,file='solution.txt') C Point write(10,100) do i = 1,n write(10,200) i,x(i) end do C Lagrange multipliers if ( m .gt. 0 ) then write(10,300) do i = 1,m write(10,200) i,lambda(i) end do end if close(10) C NON-EXECUTABLE STATEMENTS 100 format(/,'FINAL POINT:',//,2X,'INDEX',16X,'X(INDEX)') 200 format(I7,1P,D24.16) 300 format(/,'FINAL ESTIMATION OF THE LAGRANGE MULTIPLIERS: ', + //,2X,'INDEX',11X,'LAMBDA(INDEX)') end C ****************************************************************** C ****************************************************************** subroutine vevalf(n,x,f,inform) implicit none C SCALAR ARGUMENTS integer inform,n double precision f C ARRAY ARGUMENTS double precision x(n) #include "dim.par" #include "counters.com" C LOCAL SCALARS integer flag call evalf(n,x,f,flag) fcnt = fcnt + 1 if ( flag .ne. 0 ) then inform = - 90 call reperr(inform) return end if end C ****************************************************************** C ****************************************************************** subroutine vevalg(n,x,g,inform) implicit none C SCALAR ARGUMENTS integer inform,n C ARRAY ARGUMENTS double precision g(n),x(n) #include "dim.par" #include "algparam.com" #include "counters.com" C LOCAL SCALARS integer flag if ( gcoded ) then call evalg(n,x,g,flag) gcnt = gcnt + 1 if ( flag .ne. 0 ) then inform = - 92 call reperr(inform) return end if else call ievalg(n,x,g,inform) if ( inform .lt. 0 ) return end if end C ****************************************************************** C ****************************************************************** subroutine ievalg(n,x,g,inform) implicit none C SCALAR ARGUMENTS integer inform,n C ARRAY ARGUMENTS double precision g(n),x(n) #include "machconst.com" C LOCAL SCALARS integer j double precision fminus,fplus,step,tmp do j = 1,n tmp = x(j) step = macheps13 * max( 1.0d0, abs( tmp ) ) x(j) = tmp + step call setp(n,x) call vevalf(n,x,fplus,inform) if ( inform .lt. 0 ) return x(j) = tmp - step call setp(n,x) call vevalf(n,x,fminus,inform) if ( inform .lt. 0 ) return g(j) = ( fplus - fminus ) / ( 2.0d0 * step ) x(j) = tmp end do end C ****************************************************************** C ****************************************************************** subroutine vevalh(n,x,hlin,hcol,hval,hnnz,inform) implicit none C SCALAR ARGUMENTS integer inform,n,hnnz C ARRAY ARGUMENTS integer hcol(*),hlin(*) double precision hval(*),x(n) #include "dim.par" #include "counters.com" C LOCAL SCALARS integer flag,i call evalh(n,x,hlin,hcol,hval,hnnz,flag) hcnt = hcnt + 1 if ( flag .ne. 0 ) then inform = - 94 call reperr(inform) return end if do i = 1,hnnz if ( hlin(i) .lt. 1 .or. hlin(i) .gt. n .or. + hcol(i) .lt. 1 .or. hcol(i) .gt. n .or. + hcol(i) .gt. hlin(i) ) then write(* ,100) i,hlin(i),hcol(i),hval(i) write(10,100) i,hlin(i),hcol(i),hval(i) hlin(i) = 1 hcol(i) = 1 hval(i) = 0.0d0 end if end do C NON-EXECUTABLE STATEMENTS 100 format(/,1X,'WARNING: There is an element out of range, or in ', + 'the upper triangle, of the Hessian of the objetive ', + 'function computed by the user-supplied subroutine ', + 'EVALH. It will be ignored.', + /,1X,'Position: ',I16, + /,1X,'Row : ',I16, + /,1X,'Column : ',I16, + /,1X,'Value : ',1P,D24.16) end C ****************************************************************** C ****************************************************************** subroutine vevalc(n,x,ind,c,inform) implicit none C SCALAR ARGUMENTS integer ind,inform,n double precision c C ARRAY ARGUMENTS double precision x(n) #include "dim.par" #include "counters.com" C LOCAL SCALARS integer flag call evalc(n,x,ind,c,flag) ccnt(ind) = ccnt(ind) + 1 if ( flag .ne. 0 ) then inform = - 91 call reperr(inform) return end if end C ****************************************************************** C ****************************************************************** subroutine vevaljac(n,x,ind,jcvar,jcval,jcnnz,inform) implicit none C SCALAR ARGUMENTS integer inform,ind,n,jcnnz C ARRAY ARGUMENTS integer jcvar(n) double precision x(n),jcval(n) #include "dim.par" #include "algparam.com" #include "counters.com" C LOCAL SCALARS integer flag,i if ( jaccoded ) then call evaljac(n,x,ind,jcvar,jcval,jcnnz,flag) jccnt(ind) = jccnt(ind) + 1 if ( flag .ne. 0 ) then inform = - 93 call reperr(inform) return end if do i = 1,jcnnz if ( jcvar(i) .lt. 1 .or. jcvar(i) .gt. n ) then write(* ,100) ind,i,jcvar(i),jcval(i) write(10,100) ind,i,jcvar(i),jcval(i) jcvar(i) = 1 jcval(i) = 0.0d0 end if end do else call ievaljac(n,x,ind,jcvar,jcval,jcnnz,inform) if ( inform .lt. 0 ) return end if C NON-EXECUTABLE STATEMENTS 100 format(/,1X,'WARNING: There is an element out of range in ', + 'gradient of constraint ',I16,' computed by the ', + 'user-supplied subroutine EVALJAC. It will be ', + 'ignored.', + /,1X,'Position: ',I16, + /,1X,'Variable: ',I16, + /,1X,'Value : ',1P,D24.16) end C ****************************************************************** C ****************************************************************** subroutine ievaljac(n,x,ind,jcvar,jcval,jcnnz,inform) implicit none C SCALAR ARGUMENTS integer ind,inform,n,jcnnz C ARRAY ARGUMENTS integer jcvar(n) double precision jcval(n),x(n) #include "dim.par" #include "machconst.com" C LOCAL SCALARS integer j double precision cminus,cplus,step,tmp jcnnz = 0 do j = 1,n tmp = x(j) step = macheps13 * max( 1.0d0, abs( tmp ) ) x(j) = tmp + step call setp(n,x) call vevalc(n,x,ind,cplus,inform) if ( inform .lt. 0 ) return x(j) = tmp - step call setp(n,x) call vevalc(n,x,ind,cminus,inform) if ( inform .lt. 0 ) return jcvar(jcnnz + 1) = j jcval(jcnnz + 1) = ( cplus - cminus ) / ( 2.0d0 * step ) if ( abs( jcval(jcnnz + 1) ) .gt. 0.0d0 ) then jcnnz = jcnnz + 1 end if x(j) = tmp end do end C ****************************************************************** C ****************************************************************** subroutine vevalhc(n,x,ind,hlin,hcol,hval,hnnz,inform) implicit none C SCALAR ARGUMENTS integer inform,ind,n,hnnz C ARRAY ARGUMENTS integer hcol(*),hlin(*) double precision hval(*),x(n) #include "dim.par" #include "counters.com" C LOCAL SCALARS integer flag,i call evalhc(n,x,ind,hlin,hcol,hval,hnnz,flag) hccnt(ind) = hccnt(ind) + 1 if ( flag .ne. 0 ) then inform = - 95 call reperr(inform) return end if do i = 1,hnnz if ( hlin(i) .lt. 1 .or. hlin(i) .gt. n .or. + hcol(i) .lt. 1 .or. hcol(i) .gt. n .or. + hcol(i) .gt. hlin(i) ) then write(* ,100) ind,i,hlin(i),hcol(i),hval(i) write(10,100) ind,i,hlin(i),hcol(i),hval(i) hlin(i) = 1 hcol(i) = 1 hval(i) = 0.0d0 end if end do C NON-EXECUTABLE STATEMENTS 100 format(/,1X,'WARNING: There is an element out of range, or in ', + 'the upper triangle, of the Hessian of constraint ', + I16,' computed by the user-supplied subroutine ', + 'EVALHC. It will be ignored.', + /,1X,'Position: ',I16, + /,1X,'Row : ',I16, + /,1X,'Column : ',I16, + /,1X,'Value : ',1P,D24.16) end C ****************************************************************** C ****************************************************************** subroutine vevalhl(n,x,m,lambda,sf,sc,hlin,hcol,hval,hnnz,inform) implicit none C SCALAR ARGUMENTS integer hnnz,inform,m,n double precision sf C ARRAY ARGUMENTS integer hlin(*),hcol(*) double precision hval(*),lambda(m),sc(m),x(n) #include "dim.par" #include "algparam.com" #include "counters.com" C LOCAL SCALARS integer flag,i if ( hlcoded ) then call evalhl(n,x,m,lambda,sf,sc,hlin,hcol,hval,hnnz,flag) hlcnt = hlcnt + 1 if ( flag .ne. 0 ) then inform = - 96 call reperr(inform) return end if do i = 1,hnnz if ( hlin(i) .lt. 1 .or. hlin(i) .gt. n .or. + hcol(i) .lt. 1 .or. hcol(i) .gt. n .or. + hcol(i) .gt. hlin(i) ) then write(* ,100) i,hlin(i),hcol(i),hval(i) write(10,100) i,hlin(i),hcol(i),hval(i) hlin(i) = 1 hcol(i) = 1 hval(i) = 0.0d0 end if end do else if ( hcoded .and. ( hccoded .or. m .eq. 0 ) ) then call ievalhl(n,x,m,lambda,sf,sc,hlin,hcol,hval,hnnz,inform) if ( inform .lt. 0 ) return end if C NON-EXECUTABLE STATEMENTS 100 format(/,1X,'WARNING: There is an element out of range, or in ', + 'the upper triangle, of the Hessian of the ', + 'Lagrangian computed by the user-supplied ', + 'subroutine EVALHL. It will be ignored.', + /,1X,'Position: ',I16, + /,1X,'Row : ',I16, + /,1X,'Column : ',I16, + /,1X,'Value : ',1P,D24.16) end C ****************************************************************** C ****************************************************************** subroutine ievalhl(n,x,m,lambda,sf,sc,hlin,hcol,hval,hnnz,inform) implicit none C SCALAR ARGUMENTS integer hnnz,inform,m,n double precision sf C ARRAY ARGUMENTS integer hlin(*),hcol(*) double precision hval(*),lambda(m),sc(m),x(n) #include "dim.par" C LOCAL SCALARS integer col,con,hnnztmp,i,ind,itmp,j,lin,nextj,rnnz double precision val C LOCAL ARRAYS integer hcon(hnnzmax),pos(nmax),rind(nmax),stlin(nmax) double precision rval(nmax) C ================================================================== C COMPUTE HESSIANS C ================================================================== C COMPUTE HESSIAN OF THE OBJECTIVE FUNCTION call vevalh(n,x,hlin,hcol,hval,hnnz,inform) if ( inform .lt. 0 ) return C For each element of the Hessian of the objective function, C set constraint index as zero do i = 1,hnnz hval(i) = hval(i) / sf hcon(i) = 0 end do if ( m .eq. 0 ) return C COMPUTE HESSIANS OF THE CONSTRAINTS ind = 0 do j = 1,m C Compute Hessian of constraint j call vevalhc(n,x,j,hlin(hnnz+ind+1),hcol(hnnz+ind+1), + hval(hnnz+ind+1),hnnztmp,inform) if ( inform .lt. 0 ) return C For each element of the Hessian, set constraint as j do i = hnnz + ind + 1,hnnz + ind + hnnztmp hval(i) = hval(i) / sc(j) hcon(i) = j end do ind = ind + hnnztmp end do if ( ind .eq. 0 ) return hnnz = hnnz + ind C ================================================================== C SET ROW LINKED LISTS C ================================================================== C Initialize pointers to the first element of each row do i = 1,n stlin(i) = 0 end do C Set row linked lists do i = 1,hnnz lin = hlin(i) itmp = stlin(lin) stlin(lin) = i hlin(i) = itmp end do C ================================================================== C BUILD HESSIAN OF THE LAGRANGIAN ROW BY ROW C ================================================================== C Initialize array pos do i = 1,n pos(i) = 0 end do do i = 1,n C Initialize the i-th row of the Hessian of the Lagrangian rnnz = 0 C Process the i-th row of all the Hessians j = stlin(i) 10 if ( j .ne. 0 ) then C Process element (i,hcol(j)) of the Hessian of constraint C hcon(j) (Hessian of the objective function if hcon(j)=0) col = hcol(j) con = hcon(j) if ( con .eq. 0 ) then val = hval(j) else val = hval(j) * lambda(con) end if if ( pos(col) .ne. 0 ) then rval(pos(col)) = rval(pos(col)) + val else rnnz = rnnz + 1 pos(col) = rnnz rind(pos(col)) = col rval(pos(col)) = val end if C Get next element in the i-th row linked list j = hlin(j) go to 10 end if C Clean array pos do j = 1,rnnz pos(rind(j)) = 0 end do C Set i-th row of hl (over the i-th rows of the Hessians) C and mark remaining elements to be deleted j = stlin(i) 20 if ( j .ne. 0 ) then nextj = hlin(j) if ( rnnz .ne. 0 ) then hlin(j) = i hcol(j) = rind(rnnz) hval(j) = rval(rnnz) rnnz = rnnz - 1 else hlin(j) = 0 end if j = nextj go to 20 end if end do C Eliminate remaining elements (marked with hlin(j)=0) j = 1 30 if ( j .le. hnnz ) then if ( hlin(j) .eq. 0 ) then if ( j .ne. hnnz ) then hlin(j) = hlin(hnnz) hcol(j) = hcol(hnnz) hval(j) = hval(hnnz) end if hnnz = hnnz - 1 else j = j + 1 end if go to 30 end if end C ****************************************************************** C ****************************************************************** subroutine vevalhlp(n,x,m,lambda,sf,sc,p,hp,gothl,inform) implicit none C SCALAR ARGUMENTS logical gothl integer inform,m,n double precision sf C ARRAY ARGUMENTS double precision hp(n),lambda(m),p(n),sc(m),x(n) #include "dim.par" #include "algparam.com" #include "counters.com" C LOCAL SCALARS integer flag if ( hlpcoded ) then call evalhlp(n,x,m,lambda,sf,sc,p,hp,gothl,flag) hlpcnt = hlpcnt + 1 if ( flag .ne. 0 ) then inform = - 97 call reperr(inform) return end if else if ( truehl ) then call ievalhlp(n,x,m,lambda,sf,sc,p,hp,gothl,inform) if ( inform .lt. 0 ) return end if end C ****************************************************************** C ****************************************************************** subroutine ievalhlp(n,x,m,lambda,sf,sc,p,hp,gothl,inform) implicit none C SCALAR ARGUMENTS logical gothl integer inform,m,n double precision sf C ARRAY ARGUMENTS double precision hp(n),lambda(m),p(n),sc(m),x(n) #include "dim.par" #include "hessdat.com" C LOCAL SCALARS integer col,i,lin double precision val if ( .not. gothl ) then gothl = .true. call vevalhl(n,x,m,lambda,sf,sc,hlin,hcol,hval,hnnz,inform) if ( inform .lt. 0 ) return end if do i = 1,n hp(i) = 0.0d0 end do do i = 1,hnnz lin = hlin(i) col = hcol(i) val = hval(i) hp(lin) = hp(lin) + p(col) * val if ( lin .ne. col ) then hp(col) = hp(col) + p(lin) * val end if end do end C ****************************************************************** C ****************************************************************** subroutine vevalfc(n,x,f,m,c,inform) implicit none C SCALAR ARGUMENTS integer inform,m,n double precision f C ARRAY ARGUMENTS double precision c(m),x(n) #include "dim.par" #include "counters.com" C LOCAL SCALARS integer flag call evalfc(n,x,f,m,c,flag) fccnt = fccnt + 1 if ( flag .ne. 0 ) then inform = - 98 call reperr(inform) return end if end C ****************************************************************** C ****************************************************************** subroutine vevalgjac(n,x,g,m,jcfun,jcvar,jcval,jcnnz,inform) implicit none C SCALAR ARGUMENTS integer inform,jcnnz,m,n C ARRAY ARGUMENTS integer jcfun(*),jcvar(*) double precision g(n),jcval(*),x(n) #include "dim.par" #include "algparam.com" #include "counters.com" C LOCAL SCALARS integer flag,i if ( gjaccoded ) then call evalgjac(n,x,g,m,jcfun,jcvar,jcval,jcnnz,flag) gjccnt = gjccnt + 1 if ( flag .ne. 0 ) then inform = - 99 call reperr(inform) return end if do i = 1,jcnnz if ( jcfun(i) .lt. 1 .or. jcfun(i) .gt. m .or. + jcvar(i) .lt. 1 .or. jcvar(i) .gt. n ) then write(* ,100) i,jcfun(i),jcvar(i),jcval(i) write(10,100) i,jcfun(i),jcvar(i),jcval(i) jcfun(i) = 1 jcvar(i) = 1 jcval(i) = 0.0d0 end if end do else call ievalgjac(n,x,g,m,jcfun,jcvar,jcval,jcnnz,inform) if ( inform .lt. 0 ) return end if C NON-EXECUTABLE STATEMENTS 100 format(/,1X,'WARNING: There is an element out of range in the ', + 'Jacobian of the constraints computed by the ', + 'user-supplied subroutine EVALGJAC. It will be ', + 'ignored.', + /,1X,'Position : ',I16, + /,1X,'Constraint: ',I16, + /,1X,'Variable : ',I16, + /,1X,'Value : ',1P,D24.16) end C ****************************************************************** C ****************************************************************** subroutine ievalgjac(n,x,g,m,jcfun,jcvar,jcval,jcnnz,inform) implicit none C SCALAR ARGUMENTS integer inform,jcnnz,m,n C ARRAY ARGUMENTS integer jcfun(*),jcvar(*) double precision g(n),jcval(*),x(n) #include "dim.par" #include "machconst.com" C LOCAL SCALARS integer i,j double precision fminus,fplus,step,tmp C LOCAL ARRAYS double precision cminus(mmax),cplus(mmax) jcnnz = 0 do i = 1,n tmp = x(i) step = macheps13 * max( 1.0d0, abs( tmp ) ) x(i) = tmp + step call setp(n,x) call vevalfc(n,x,fplus,m,cplus,inform) if ( inform .lt. 0 ) return x(i) = tmp - step call setp(n,x) call vevalfc(n,x,fminus,m,cminus,inform) if ( inform .lt. 0 ) return do j = 1,m jcfun(jcnnz + 1) = j jcvar(jcnnz + 1) = i jcval(jcnnz + 1) = ( cplus(j) - cminus(j) ) / + ( 2.0d0 * step ) if ( abs( jcval(jcnnz + 1) ) .gt. 0.0d0 ) then jcnnz = jcnnz + 1 end if end do g(i) = ( fplus - fminus ) / ( 2.0d0 * step ) x(i) = tmp end do end C ****************************************************************** C ****************************************************************** subroutine reperr(inform) C SCALAR ARGUMENTS integer inform if ( inform .eq. -90 ) then write(* ,100) 'EVALF' write(10,100) 'EVALF' else if ( inform .eq. -91 ) then write(* ,100) 'EVALC' write(10,100) 'EVALC' else if ( inform .eq. -92 ) then write(* ,100) 'EVALG' write(10,100) 'EVALG' else if ( inform .eq. -93 ) then write(* ,100) 'EVALJAC' write(10,100) 'EVALJAC' else if ( inform .eq. -94 ) then write(* ,100) 'EVALH' write(10,100) 'EVALH' else if ( inform .eq. -95 ) then write(* ,100) 'EVALHC' write(10,100) 'EVALHC' else if ( inform .eq. -96 ) then write(* ,100) 'EVALHL' write(10,100) 'EVALHL' else if ( inform .eq. -97 ) then write(* ,100) 'EVALHLP' write(10,100) 'EVALHLP' else if ( inform .eq. -98 ) then write(* ,100) 'EVALFC' write(10,100) 'EVALFC' else if ( inform .eq. -99 ) then write(* ,100) 'EVALGJAC' write(10,100) 'EVALGJAC' end if C NON-EXECUTABLE STATEMENTS 100 format(/,1X,'*** There was an error in the user supplied ', + 'subroutine ',A10,' ***',/) end C ****************************************************************** C ****************************************************************** subroutine vsetp(n,x) implicit none C SCALAR ARGUMENTS integer n C ARRAY ARGUMENTS double precision x(n) #include "dim.par" #include "graddat.com" gotc = .false. call setp(n,x) end C ****************************************************************** C ****************************************************************** subroutine vunsetp() implicit none #include "dim.par" #include "graddat.com" gotc = .false. call unsetp() end