C ================================================================= C File: alabb.f C ================================================================= C ================================================================= C Module: Auxiliary subroutines C ================================================================= C Last update of any of the component of this module: C C March 27, 2006. C ****************************************************************** C ****************************************************************** subroutine solver(n,x,l,u,m,lambda,rho,equatn,linear,linvar,ml, +meq,alin,acol,aval,annz,b,fknown,checkder,epsgap,epsfeas,epsopt, +sparse,maxoutit,iprint,ncomp,f,cnorm,snorm,nalpsupn,ldh,H,outiter, +totnnlp,totnnodes,time,inform,wi1,wi2,wd1,wd2,wd3,wd4,wd5,wd6,wd7, +wd8,wd9,wd10,wd11,wd12) implicit none C This subroutine: C C (1) Computes some machine-dependent constants. C C (2) Process the optionally input specification file alabb.dat. C C (3) Check analytic derivatives if desired by the user. C C (4) Open and close the output file alabb.out. C C (5) Computes the CPU elapsed time used by the solver. C C (6) Calls the solver. C C (7) Write an output file called solution.txt with the solution. C SCALAR ARGUMENTS logical checkder,sparse integer annz,inform,iprint,ldh,m,maxoutit,meq,ml,n,ncomp,outiter, + totnnlp,totnnodes double precision epsgap,epsfeas,epsopt,f,fknown,gap,nalpsupn, + cnorm,snorm real time C ARRAY ARGUMENTS integer alin(annz),acol(annz),wi1(n),wi2(ml) logical equatn(m),linear(m),linvar(n) double precision aval(annz),b(ml),H(ldh,n),l(n),lambda(m),rho(m), + u(n),wd1(m),wd2(m),wd3(n),wd4(m),wd5(m),wd6(n),wd7(n), + wd8(n),wd9(n),wd10(n),wd11(n),wd12(n),x(n) C PARAMETERS include "dim.f" C LOCAL SCALARS integer i double precision bignum,macheps C DECLARATIONS RELATED TO THE TIME MEASUREMENT USING DTIME c real dtime c external dtime real dum(2) data dum/0.0,0.0/ C EXTERNAL FUNCTIONS double precision d1mach C COMPUTE MACHINE-DEPENDENT CONSTANTS bignum = 1.0d+99 macheps = d1mach(4) C OPEN THE OUTPUT FILE open(10,file='alabb.out') C SET SOME SOLVER ARGUMENTS USING THE SPECIFICATION FILE call fparam(checkder,epsgap,epsfeas,epsopt,sparse,maxoutit,iprint, +ncomp) epsgap = epsopt / 2.0d0 C TEST DERIVATIVES if ( checkder ) then call checkd(m,n,l,u,wi1,wd3,wd6,wd7,wd8,wd9,macheps,inform) if ( inform .lt. 0 ) then if ( iprint .ge. 1 ) then write(*, 2000) inform write(10,2000) inform end if go to 500 end if end if C CALL ALABB IF THE PROBLEM HAS CONSTRAINTS. OTHERWISE, CALL ABB time = dtime(dum) c time = 0.0d0 if ( m .gt. 0 ) then C CALL THE AUGMENTED LAGRANGIAN SOLVER ALABB call easyalabb(n,x,l,u,m,lambda,rho,equatn,linear,linvar,ml, + meq,alin,acol,aval,annz,b,fknown,epsgap,epsfeas,epsopt,sparse, + maxoutit,iprint,ncomp,macheps,bignum,f,cnorm,snorm,nalpsupn, + gap,ldh,H,outiter,totnnlp,totnnodes,inform,wi2,wd1,wd2,wd3, + wd4,wd5,wd6,wd7,wd8,wd9,wd10,wd11,wd12) else C CALL THE BOX-CONSTRAINTS GLOBAL SOLVER ABB call abb(n,l,u,m,lambda,rho,equatn,linear,linvar,ml,meq,alin, + acol,aval,annz,b,fknown,epsgap,epsfeas,epsopt,sparse,f,x,wd3, + nalpsupn,gap,totnnodes,inform) outiter = 0 cnorm = 0.0d0 snorm = 0.0d0 end if time = dtime(dum) time = dum(1) c time = 0.0d0 C CLOSE THE OUTPUT FILE close(10) C SAVE THE SOLUTION open(10,file='solution.txt') C Solution point write(10,9000) do i = 1,n write(10,9010) i,x(i) end do C Lagrange multipliers and penalty parameters if ( m .gt. 0 ) then write(10,9020) do i = 1,m write(10,9030) i,lambda(i),rho(i) end do end if close(10) C WRITE STATISTICS open(10,file='alabb-tabline.out') write(10,9050) time,inform,n,m,ml,outiter,totnnodes,f,cnorm,gap close(10) C RETURN 500 continue C NON-EXECUTABLE STATEMENTS 2000 format(/,' Flag of ALABB = ',I3,' Fatal Error',/, + /,' The following codes means: ',/, + /,' -90 : error in evalf subroutine', + /,' -91 : error in evalc subroutine', + /,' -92 : error in evalg subroutine', + /,' -93 : error in evaljac subroutine', + /,' -94 : error in evalh subroutine', + /,' -95 : error in evalhc subroutine',/) 9000 format(/,'FINAL POINT:',//,2X,'INDEX',16X,'X(INDEX)') 9010 format( I7,1P,D24.16) 9020 format(/,'FINAL ESTIMATION OF THE LAGRANGE MULTIPLIERS AND ', + 'PENALTY PARAMETERS: ',//,2X,'INDEX',11X, + 'LAMBDA(INDEX)',15X,'RHO(INDEX)') 9030 format(I7,1P,D24.16,1X,1P,D24.16) 9050 format(F16.2,1X,I3,1X,I6,1X,I6,1X,I6,1X,I3,1X,I7,1X, +1P,D24.16,1X,1P,D7.1,1X,1P,D8.1) end C ****************************************************************** C ****************************************************************** subroutine fparam(checkder,epsgap,epsfeas,epsopt,sparse,maxoutit, +iprint,ncomp) implicit none C This subroutine set some alabb arguments related to stopping C criteria and output. The setting values are taken fron file C alabb.dat. Nothing is done if the file does not exist. File C alabb.dat is an alternative to run ALABB several times C with different arguments without having to compile it again. C SCALAR ARGUMENTS logical checkder,sparse integer iprint,maxoutit,ncomp double precision epsfeas,epsgap,epsopt C PARAMETERS integer nwords parameter ( nwords = 8 ) C DATA BLOCKS character * 1 lower(26),upper(26) data lower /'a','b','c','d','e','f','g','h','i','j','k','l','m', + 'n','o','p','q','r','s','t','u','v','w','x','y','z'/ data upper /'A','B','C','D','E','F','G','H','I','J','K','L','M', + 'N','O','P','Q','R','S','T','U','V','W','X','Y','Z'/ character * 10 dictionary(nwords) data dictionary(1) /'CHECK-DERI'/ data dictionary(2) /'GAP-TOLERA'/ data dictionary(3) /'FEASIBILIT'/ data dictionary(4) /'OPTIMALITY'/ data dictionary(5) /'MAX-OUTER-'/ data dictionary(6) /'OUTPUT-DET'/ data dictionary(7) /'NCOMP-ARRA'/ data dictionary(8) /'SPARSE-INN'/ character * 32 description(nwords) data description(1) /'CHECK-DERIVATIVES '/ data description(2) /'GAP-TOLERANCE '/ data description(3) /'FEASIBILITY-TOLERANCE '/ data description(4) /'OPTIMALITY-TOLERANCE '/ data description(5) /'MAX-OUTER-ITERATIONS '/ data description(6) /'OUTPUT-DETAIL '/ data description(7) /'NCOMP-ARRAY '/ data description(8) /'SPARSE-INNER-SOLVER '/ character * 1 addinfo(nwords) data addinfo(1) /' '/ data addinfo(2) /'D'/ data addinfo(3) /'D'/ data addinfo(4) /'D'/ data addinfo(5) /'I'/ data addinfo(6) /'I'/ data addinfo(7) /'I'/ data addinfo(8) /' '/ C LOCAL SCALARS integer i,ifirst,ikey,ilast,inum,j double precision dnum C LOCAL ARRAYS character * 80 line character * 10 keyword C OPENING THE SPECIFICATION FILE open(20,err=300,file='alabb.dat',status='old') C MAIN LOOP write(*, 9005) write(10,9005) 100 continue C READING LINES read(20,fmt=1000,err=400,end=200) line C PROCESS LINES C Find first character i = 1 110 if ( i .le. 80 .and. line(i:i) .eq. ' ' ) then i = i + 1 go to 110 end if C Skip blank lines if ( i .gt. 80 ) then go to 100 end if ifirst = i C Skip comments if ( line(ifirst:ifirst) .eq. '*' .or. + line(ifirst:ifirst) .eq. '#' ) then go to 100 end if C Find the end of the keyword i = ifirst + 1 120 if ( i .le. 80 .and. line(i:i) .ne. ' ' ) then i = i + 1 go to 120 end if ilast = i - 1 C Obtain the first 10 characters and convert to upper-case keyword = ' ' do i = 1,min( 10, ilast - ifirst + 1 ) keyword(i:i) = line(ifirst+i-1:ifirst+i-1) do j = 1,26 if ( keyword(i:i) .eq. lower(j) ) then keyword(i:i) = upper(j) end if end do end do C Look up the keyword in the dictionary i = 1 130 if ( i .le. nwords .and. keyword .ne. dictionary(i) ) then i = i + 1 go to 130 end if C Ignore unknown keywords if ( i .gt. nwords ) then write(*, 9020) line(ifirst:ilast) write(10,9020) line(ifirst:ilast) go to 100 end if ikey = i C Read additional information if needed if ( addinfo(ikey) .ne. ' ' ) then C Skip blanks i = ilast + 1 140 if ( i .le. 80 .and. line(i:i) .eq. ' ' ) then i = i + 1 go to 140 end if C Ignore keywords without the required information if ( i .gt. 80 ) then write(*, 9030) description(ikey) write(10,9030) description(ikey) go to 100 end if C Read additional information if ( addinfo(ikey) .eq. 'I' ) then read(unit=line(i:80),fmt=2000) inum else if ( addinfo(ikey) .eq. 'D' ) then read(unit=line(i:80),fmt=3000) dnum end if end if C Process keyword if ( addinfo(ikey) .eq. ' ' ) then write(*, 9040) description(ikey) write(10,9040) description(ikey) else if ( addinfo(ikey) .eq. 'I' ) then write(*, 9041) description(ikey),inum write(10,9041) description(ikey),inum else if ( addinfo(ikey) .eq. 'D' ) then write(*, 9042) description(ikey),dnum write(10,9042) description(ikey),dnum end if C Set the corresponding alabb argument if ( ikey .eq. 1 ) then checkder = .true. else if ( ikey .eq. 2 ) then epsgap = dnum else if ( ikey .eq. 3 ) then epsfeas = dnum else if ( ikey .eq. 4 ) then epsopt = dnum else if ( ikey .eq. 5 ) then maxoutit = inum else if ( ikey .eq. 6 ) then iprint = inum else if ( ikey .eq. 7 ) then ncomp = inum else if ( ikey .eq. 8 ) then sparse = .true. end if C IIERATE go to 100 C END OF LOOP C TERMINATIONS C CLOSING SPECIFICATION FILE 200 continue close(20) go to 500 C NO SPECIFICATION FILE 300 continue write(*, 9000) write(10,9000) go to 500 C ERROR READING THE SPECIFICATION FILE 400 continue write(*, 9010) write(10,9010) go to 500 C PRINTING PARAMETERS VALUES 500 continue write(*, 4000)checkder,epsgap,epsfeas,epsopt,sparse,maxoutit, + iprint,ncomp write(10,4000)checkder,epsgap,epsfeas,epsopt,sparse,maxoutit, + iprint,ncomp C NON-EXECUTABLE STATEMENTS 1000 format(A80) 2000 format(BN,I20) 3000 format(BN,F24.0) 4000 format(/,' ALABB PARAMETERS:', + /,' checkder = ', 19X,L1, + /,' epsgap = ',8X,1P,D12.4, + /,' epsfeas = ',8X,1P,D12.4, + /,' epsopt = ',8X,1P,D12.4, + /,' sparse = ', 19X,L1, + /,' maxoutit = ', I20, + /,' iprint = ', I20, + /,' ncomp = ', I20) 9000 format(/,' Specification file alabb.dat was not found in the', + ' current directory.',/,' Using default values.') 9005 format(/,' Specification file alabb.dat is being used.') 9010 format(/,' Error reading specification file alabb.dat.') 9020 format( ' Ignoring unknown keyword ',A32) 9030 format( ' Ignoring incomplete keyword ',A32) 9040 format(1X,A32) 9041 format(1X,A32,5X,I20) 9042 format(1X,A32,1X,1P,D24.8) end C ***************************************************************** C ***************************************************************** subroutine checkd(m,n,l,u,wi,wd1,wd2,wd3,wd4,x,macheps,inform) implicit none C This subrotutine checks the user supplied first and second C derivatives subroutines (evalg, evalh, evaljac and evalhc) for C computing the objective function gradient and Hessian and the C constraints gradients and Hessians, respectively. C SCALAR ARGUMENTS integer inform,m,n double precision macheps C ARRAY ARGUMENTS integer wi(n) double precision l(n),u(n),wd1(n),wd2(n),wd3(n),wd4(n),x(n) C Parameters of the subroutine: C ============================= C C On Entry: C ========= C C M integer: Number of constraints. C --------- C C N integer: Number of variables. C --------- C C L double precision l(n): Lower bounds. C ----------------------- C C U double precision u(n): Upper bounds. C ----------------------- C C WI double precision wi(n): N-dimensional integer working vector. C ----------------------- C C WD1 ... WD4, X double precision wd1(n) ... wd4(n), x(n) C ------------------------------------------------------- C C N-dimensional double precision working vectors. C C MACHEPS double precision C ------------------------ C C Smallest positive number such that 1 + macheps is not equal to 1. C C On Return: C ========== C C INFORM integer: Output flag. C -------------- C LOCAL SCALARS character answer integer i,j double precision drand,seed,smalll,smallu inform = 0 C SET A RANDOM POINT seed = 17325.0d0 do i = 1,n smalll = max( l(i), - 1.0d0 ) smallu = min( u(i), smalll + 1.0d0 ) x(i) = smalll + ( smallu - smalll ) * drand(seed) end do write(*,900) do i = 1,n write(*,910) i,x(i) end do C CHECK OBJECTIVE FUNCTION GRADIENT write(*,920) read(*,*) answer if ( answer .eq. 'A' .or. answer .eq. 'a' ) then go to 500 else if ( answer .eq. 'Y' .or. answer .eq. 'y' ) then call checkg(n,x,wd1,macheps,inform) if ( inform .lt. 0 ) then return end if end if C CHECK CONSTRAINTS GRADIENTS j = 1 110 if ( j .le. m ) then write(*,930) j read(*,*) answer if ( answer .eq. 'A' .or. answer .eq. 'a' ) then go to 500 else if ( answer .eq. 'S' .or. answer .eq. 's' ) then go to 200 else if ( answer .eq. 'Y' .or. answer .eq. 'y' ) then call checkjac(n,x,j,wi,wd1,wd2,macheps,inform) if ( inform .lt. 0 ) then return end if end if j = j + 1 go to 110 end if C CHECK HESSIAN OF THE OBJECTIVE FUNCTION 200 continue write(*,940) read(*,*) answer if ( answer .eq. 'A' .or. answer .eq. 'a' ) then go to 500 else if ( answer .eq. 'Y' .or. answer .eq. 'y' ) then call checkh(n,x,wd1,wd2,wd3,macheps,inform) if ( inform .lt. 0 ) then return end if end if C CHECK HESSIANS OF THE CONSTRAINTS j = 1 310 if ( j .le. m ) then write(*,950) j read(*,*) answer if ( answer .eq. 'A' .or. answer .eq. 'a' ) then go to 500 else if ( answer .eq. 'S' .or. answer .eq. 's' ) then go to 500 else if ( answer .eq. 'Y' .or. answer .eq. 'y' ) then call checkhc(n,x,j,wi,wd1,wd2,wd3,wd4,macheps,inform) if ( inform .lt. 0 ) then return end if end if j = j + 1 go to 310 end if C RETURN 500 continue write(*,*) 'Hit any letter to continue.' read(*,*) answer 900 format(/,1X,'Derivatives will be tested at the random point: ') 910 format( 1X,'x(',I6,') = ',1P,D15.8) 920 format(/,1X,'Check the gradient of the objective function?', + /,1X,'Type Y(es), N(o) or A(bort checking): ') 930 format(/,1X,'Check the gradient of constraint ',I5,'?', + /,1X,'Type Y(es), N(o), A(bort checking) or ', + 'S(kip constraints gradients): ') 940 format(/,1X,'Check the Hessian matrix of the objective function?', + /,1X,'Type Y(es), N(o) or A(bort checking): ') 950 format(/,1X,'Check the Hessian of constraint ',I5,'?', + /,1X,'Type Y(es), N(o), A(bort checking) or ', + 'S(kip constraints gradients): ') end C ***************************************************************** C ***************************************************************** subroutine checkg(n,x,g,macheps,inform) implicit none C This subrotutine checks the user supplied subroutine evalg for C computing the gradient of the objective function using central C finite differences with two different discretization steps. C SCALAR ARGUMENTS integer inform,n double precision macheps C ARRAY ARGUMENTS double precision g(n),x(n) C Parameters of the subroutine: C ============================= C C On Entry: C ========= C C N integer: Number of variables. C --------- C C X double precision x(n): Point at which the gradient will be tested. C ----------------------- C C G double precision g(n) C ----------------------- C C N-dimensional double precision working vector. C C MACHEPS double precision C ------------------------ C C Smallest positive number such that 1 + macheps is not equal to 1. C C On Return: C ========== C C INFORM integer: Output flag. C -------------- C LOCAL SCALARS integer flag,i double precision fminus,fplus,gdiff1,gdiff2,maxerr,step1,step2,tmp inform = 0 call evalg(n,x,g,flag) if ( flag .ne. 0 ) then inform = - 92 return end if write(*,100) maxerr = 0.0d0 do i = 1,n tmp = x(i) step1 = macheps ** (1.0d0/3.0d0) * max( abs( tmp ), 1.0d0 ) x(i) = tmp + step1 call evalf(n,x,fplus,flag) if ( flag .ne. 0 ) then inform = - 90 return end if x(i) = tmp - step1 call evalf(n,x,fminus,flag) if ( flag .ne. 0 ) then inform = - 90 return end if gdiff1 = ( fplus - fminus ) / ( 2.0d0 * step1 ) step2 = macheps ** (1.0d0/3.0d0) * max( abs( tmp ), 1.0d-03 ) x(i) = tmp + step2 call evalf(n,x,fplus,flag) if ( flag .ne. 0 ) then inform = - 90 return end if x(i) = tmp - step2 call evalf(n,x,fminus,flag) if ( flag .ne. 0 ) then inform = - 90 return end if x(i) = tmp gdiff2 = ( fplus - fminus ) / ( 2.0d0 * step2 ) tmp = min( abs( g(i) - gdiff1 ), abs( g(i) - gdiff2 ) ) write(*,110) i,g(i),gdiff1,gdiff2,tmp maxerr = max( maxerr, tmp ) end do write(*,120) maxerr return 100 format(/,1X,'Gradient vector of the objective function.', + /,1X,'Index',13X,'evalg',2X,'Central diff (two different ', + 'steps)',4X,'Absolute error') 110 format( 1X,I5,4(3X,1P,D15.8)) 120 format( 1X,'Maximum absolute error = ',1P,D15.8) end C ***************************************************************** C ***************************************************************** subroutine checkh(n,x,g,gplus1,gplus2,macheps,inform) implicit none C This subrotutine checks the user supplied subroutine evalh for C computing the Hessian of the objective function using central C finite differences with two different discretization steps. C SCALAR ARGUMENTS integer inform,n double precision macheps C ARRAY ARGUMENTS double precision g(n),gplus1(n),gplus2(n),x(n) C Parameters of the subroutine: C ============================= C C On Entry: C ========= C C N integer: Number of variables. C --------- C C X double precision x(n): Point at which the Hessian will be tested. C ----------------------- C C G, GPLUS1, GPLUS2 double precision g(n), gplus1(n), gplus2(n) C ------------------------------------------------------------- C C N-dimensional double precision working vectors. C C MACHEPS double precision C ------------------------ C C Smallest positive number such that 1 + macheps is not equal to 1. C C On Return: C ========== C C INFORM integer: Output flag. C -------------- C PARAMETERS include "dim.f" C LOCAL SCALARS integer flag,i,j,nnzh double precision elem,hdiff1,hdiff2,step1,step2,tmp C LOCAL ARRAYS integer hlin(nmax**2),hcol(nmax**2) double precision H(nmax,nmax),hval(nmax**2),maxerr(nmax) inform = 0 C Test viability if ( n .gt. nmax ) then write(*,*) 'Subroutine CheckH uses dense matrices up to ', + 'dimension ',nmax,' times ',nmax,'. The Hessian ', + 'checking will be skipped.' go to 500 end if C Compute the gradient of the objective function at x call evalg(n,x,g,flag) if ( flag .ne. 0 ) then inform = - 92 return end if C Compute the Hessian of the objective function at x and save in a C dense matrix call evalh(n,x,hlin,hcol,hval,nnzh,flag) if ( flag .ne. 0 ) then inform = - 94 return end if do j = 1,n do i = 1,n H(i,j) = 0.0d0 end do end do do i = 1,nnzh H(hlin(i),hcol(i)) = H(hlin(i),hcol(i)) + hval(i) end do C Test column by column write(*,100) do j = 1,n tmp = x(j) step1 = sqrt( macheps ) * max( abs( tmp ), 1.0d0 ) x(j) = tmp + step1 call evalg(n,x,gplus1,flag) if ( flag .ne. 0 ) then inform = - 92 return end if step2 = sqrt( macheps ) * max( abs( tmp ), 1.0d-03 ) x(j) = tmp + step2 call evalg(n,x,gplus2,flag) if ( flag .ne. 0 ) then inform = - 92 return end if x(j) = tmp write(*,105) j maxerr(j) = 0.0d0 do i = 1,n if ( i .ge. j ) then elem = H(i,j) else elem = H(j,i) end if hdiff1 = ( gplus1(i) - g(i) ) / step1 hdiff2 = ( gplus2(i) - g(i) ) / step2 tmp = min( abs( elem - hdiff1 ), abs( elem - hdiff2 ) ) if ( elem .ne. 0.0d0 .or. + hdiff1 .ne. 0.0d0 .or. hdiff2 .ne. 0.0d0 ) then write(*,110) i,elem,hdiff1,hdiff2,tmp end if maxerr(j) = max( maxerr(j), tmp ) end do write(*,120) maxerr(j) end do write(*,*) do j = 1,n write(*,130) j,maxerr(j) end do 500 continue return 100 format(/,1X,'Hessian matrix of the objective function column by ', + 'column.') 105 format(/,1X,'Column: ',I6,/, + /,1X,'Index',13X,'evalh',3X,'Incr. Quoc. (two different ', + 'steps)',4X,'Absolute error') 110 format( 1X,I5,4(3X,1P,D15.8)) 120 format( 1X,'Maximum absolute error = ',1P,D15.8) 130 format( 1X,'Column ',I6,' Maximum absolute error = ',1P,D15.8) end C ***************************************************************** C ***************************************************************** subroutine checkjac(n,x,ind,indjac,valjac,g,macheps,inform) implicit none C This subrotutine checks the user supplied subroutine evaljac for C computing the gradients of the constraints using central finite C differences with two different discretization steps. C SCALAR ARGUMENTS integer ind,inform,n double precision macheps C ARRAY ARGUMENTS integer indjac(n) double precision g(n),valjac(n),x(n) C Parameters of the subroutine: C ============================= C C On Entry: C ========= C C N integer: Number of variables. C --------- C C X double precision x(n) C ----------------------- C C Point at which the gradient of the ind-th constraint will be C tested. C C IND integer C ----------- C C Index of the constraint whose gradient will be tested. C C INDJAC integer indjac(n): N-dimensional integer working vector. C ------------------------ C C VALJAC, G double precision valjac(n), g(n) C ------------------------------------------ C C N-dimensional double precision working vectors. C C MACHEPS double precision C ------------------------ C C Smallest positive number such that 1 + macheps is not equal to 1. C C On Return: C ========== C C INFORM integer: Output flag. C -------------- C LOCAL SCALARS integer flag,i,nnzjac double precision cminus,cplus,jacdiff1,jacdiff2,maxerr,step1, + step2,tmp inform = 0 C COMPUTE THE GRADIENT OF THE CONSTRAINT AND SAVE IT INTO A DENSE C VECTOR call evaljac(n,x,ind,indjac,valjac,nnzjac,flag) if ( flag .ne. 0 ) then inform = - 93 return end if do i = 1,n g(i) = 0.0d0 end do do i = 1,nnzjac g(indjac(i)) = g(indjac(i)) + valjac(i) end do C COMPARE WITH CENTRAL FINITE DIFFERENCES write(*,100) ind maxerr = 0.0d0 do i = 1,n tmp = x(i) step1 = macheps ** (1.0d0/3.0d0) * max( abs( tmp ), 1.0d0 ) x(i) = tmp + step1 call evalc(n,x,ind,cplus,flag) if ( flag .ne. 0 ) then inform = - 91 return end if x(i) = tmp - step1 call evalc(n,x,ind,cminus,flag) if ( flag .ne. 0 ) then inform = - 91 return end if jacdiff1 = ( cplus - cminus ) / ( 2.0d0 * step1 ) step2 = macheps * (1.0d0/3.0d0) * max( abs( tmp ), 1.0d-03 ) x(i) = tmp + step2 call evalc(n,x,ind,cplus,flag) if ( flag .ne. 0 ) then inform = - 91 return end if x(i) = tmp - step2 call evalc(n,x,ind,cminus,flag) if ( flag .ne. 0 ) then inform = - 91 return end if x(i) = tmp jacdiff2 = ( cplus - cminus ) / ( 2.0d0 * step2 ) tmp = min( abs( g(i) - jacdiff1 ), abs( g(i) - jacdiff2 ) ) if ( g(i) .ne. 0.0d0 .or. + jacdiff1 .ne. 0.0d0 .or. jacdiff2 .ne. 0.0d0 ) then write(*,110) i,g(i),jacdiff1,jacdiff2,tmp end if maxerr = max( maxerr, tmp ) end do write(*,120) maxerr return 100 format(/,1X,'Gradient vector of constraints ',I5,'.', + /,1X,'Index',11X,'evaljac',2X,'Central diff (two ', + 'different steps)',4X,'Absolute error') 110 format( 1X,I5,4(3X,1P,D15.8)) 120 format( 1X,'Maximum absolute error = ',1P,D15.8) end C ***************************************************************** C ***************************************************************** subroutine checkhc(n,x,ind,indjac,g,gplus1,gplus2,valjac,macheps, +inform) implicit none C This subrotutine checks the user supplied subroutine evalhc for C computing the Hessians of the constraints using finite C differences. C SCALAR ARGUMENTS integer ind,inform,n double precision macheps C ARRAY ARGUMENTS integer indjac(n) double precision g(n),gplus1(n),gplus2(n),valjac(n),x(n) C Parameters of the subroutine: C ============================= C C On Entry: C ========= C C N integer: Number of variables. C --------- C C X double precision x(n) C ----------------------- C C Point at which the Hessian of the ind-th constraint will be C tested. C C IND integer C ----------- C C Index of the constraint whose gradient will be tested. C C INDJAC integer indjac(n): N-dimensional integer working vector. C ------------------------ C C G double precision g(n) C GPLUS1 double precision gplus1(n) C GPLUS2 double precision gplus2(n) C VALJAC double precision valjac(n) C --------------------------------- C C N-dimensional double precision working vectors. C C MACHEPS double precision C ------------------------ C C Smallest positive number such that 1 + macheps is not equal to 1. C C On Return: C ========== C C INFORM integer: Output flag. C -------------- C PARAMETERS include "dim.f" C LOCAL SCALARS integer flag,i,j,nnzh,nnzjac double precision elem,hdiff1,hdiff2,step1,step2,tmp C LOCAL ARRAYS integer hlin(nmax**2),hcol(nmax**2) double precision H(nmax,nmax),hval(nmax**2),maxerr(nmax) inform = 0 C Test viability if ( n .gt. nmax ) then write(*,*) 'Subroutine CheckHc uses dense matrices up to ', + 'dimension ',nmax,' times ',nmax,'. The Hessian ', + 'checking will be skipped.' go to 500 end if C Compute the constraint gradient at x and save in a dense vector call evaljac(n,x,ind,indjac,valjac,nnzjac,flag) if ( flag .ne. 0 ) then inform = - 93 return end if do i = 1,n g(i) = 0.0d0 end do do i = 1,nnzjac g(indjac(i)) = g(indjac(i)) + valjac(i) end do C Compute the Hessian of the constraints at x and save in C dense matrix call evalhc(n,x,ind,hlin,hcol,hval,nnzh,flag) if ( flag .ne. 0 ) then inform = - 95 return end if do j = 1,n do i = 1,n H(i,j) = 0.0d0 end do end do do i = 1,nnzh H(hlin(i),hcol(i)) = H(hlin(i),hcol(i)) + hval(i) end do write(*,100) ind do j = 1,n tmp = x(j) C Compute the constraint gradient at xplus1 and save in a C dense vector step1 = sqrt( macheps ) * max( abs( tmp ), 1.0d0 ) x(j) = tmp + step1 call evaljac(n,x,ind,indjac,valjac,nnzjac,flag) if ( flag .ne. 0 ) then inform = - 93 return end if do i = 1,n gplus1(i) = 0.0d0 end do do i = 1,nnzjac gplus1(indjac(i)) = valjac(i) end do C Compute the constraint gradient at xplus2 and save in a C dense vector step2 = sqrt( macheps ) * max( abs( tmp ), 1.0d-03 ) x(j) = tmp + step2 call evaljac(n,x,ind,indjac,valjac,nnzjac,flag) if ( flag .ne. 0 ) then inform = - 93 return end if do i = 1,n gplus2(i) = 0.0d0 end do do i = 1,nnzjac gplus2(indjac(i)) = valjac(i) end do x(j) = tmp write(*,105) j maxerr(j) = 0.0d0 do i = 1,n if ( i .ge. j ) then elem = H(i,j) else elem = H(j,i) end if hdiff1 = ( gplus1(i) - g(i) ) / step1 hdiff2 = ( gplus2(i) - g(i) ) / step2 tmp = min( abs( elem - hdiff1 ), abs( elem - hdiff2 ) ) if ( elem .ne. 0.0d0 .or. + hdiff1 .ne. 0.0d0 .or. hdiff2 .ne. 0.0d0 ) then write(*,110) i,elem,hdiff1,hdiff2,tmp end if maxerr(j) = max( maxerr(j), tmp ) end do write(*,120) maxerr(j) end do write(*,*) do j = 1,n write(*,130) j,maxerr(j) end do 500 continue return 100 format(/,1X,'Hessian matrix of constraint ',I5,' column by ', + 'column.') 105 format(/,1X,'Column: ',I6,/, + /,1X,'Index',12X,'evalhc',3X,'Incr. Quoc. (two different ', + 'steps)',4X,'Absolute error') 110 format( 1X,I5,4(3X,1P,D15.8)) 120 format( 1X,'Maximum absolute error = ',1P,D15.8) 130 format( 1X,'Column ',I6,' Maximum absolute error = ',1P,D15.8) end