C ================================================================= C File: problem03a.f C ================================================================= C ****************************************************************** C ****************************************************************** subroutine inip(n,x,l,u,m,lambda,rho,equatn,linear,ml,meq,alin, +acol,aval,annz,b,fknown) implicit none C SCALAR ARGUMENTS integer annz,m,ml,meq,n double precision fknown C ARRAY ARGUMENTS logical equatn(*),linear(*) integer alin(*),acol(*) double precision aval(*),b(*),l(*),lambda(*),rho(*),u(*),x(*) C This subroutine must set some problem data. For achieving this C objective YOU MUST MODIFY it according to your problem. See below C where your modifications must be inserted. C C Parameters of the subroutine: C C On Entry: C C This subroutine has no input parameters. C C On Return C C n integer, C number of variables, C C x double precision x(n), C initial point, C C l double precision l(n), C lower bounds on x, C C u double precision u(n), C upper bounds on x, C C m integer, C number of constraints (excluding the bounds), C C lambda double precision lambda(m), C initial estimation of the Lagrange multipliers, C C rho double precision rho(m), C initial penalty parameters. C C equatn logical equatn(m) C for each constraint j, set equatn(j) = .true. if it is an C equality constraint of the form c_j(x) = 0, and set C equatn(j) = .false. if it is an inequality constraint of C the form c_j(x) <= 0, C C linear logical linear(m) C for each constraint j, set linear(j) = .true. if it is a C linear constraint, and set linear(j) = .false. if it is a C nonlinear constraint. C COMMON ARRAYS double precision k(4) C LOCAL SCALARS integer i double precision seed C COMMON BLOCKS common /probdata/ k C EXTERNAL FUNCTIONS double precision drand C Set problem data C write(*,*) 'Seed for the random initial point: ' C read(*,*) seed seed = 123456.0d0 k(1) = 9.755988d-02 k(2) = 0.99d0 * k(1) k(3) = 3.919080d-02 k(4) = 0.90d0 * k(3) C Number of variables n = 6 C Lower and upper bounds l(1) = 0.0d0 l(2) = 0.0d0 l(3) = 0.0d0 l(4) = 0.0d0 l(5) = 1.0d-05 l(6) = 1.0d-05 u(1) = 1.0d0 u(2) = 1.0d0 u(3) = 1.0d0 u(4) = 1.0d0 u(5) = 16.0d0 u(6) = 16.0d0 C Initial point do i = 1,n x(i) = l(i) + ( u(i) - l(i) ) * drand(seed) end do C Number of constraints (equalities plus inequalities) m = 5 C Lagrange multipliers approximation. do i = 1,m lambda(i) = 0.0d0 end do C Initial penalty parameters do i = 1,m rho(i) = 10.0d0 end do C For each constraint i, set equatn(i) = .true. if it is an equality C constraint of the form c_i(x) = 0, and set equatn(i) = .false. if C it is an inequality constraint of the form c_i(x) <= 0. equatn(1) = .true. equatn(2) = .true. equatn(3) = .true. equatn(4) = .true. equatn(5) = .false. C For each constraint i, set linear(i) = .true. if it is a linear C constraint, otherwise set linear(i) = .false. do i = 1,m linear(i) = .false. end do C Number of lower-level linear constraints ml = 0 meq = 0 annz = 0 fknown = -3.880823580776311d-01 end C ****************************************************************** C ****************************************************************** subroutine evalf(n,x,f,flag) implicit none C SCALAR ARGUMENTS integer flag,n double precision f C ARRAY ARGUMENTS double precision x(n) C This subroutine must compute the objective function. For achieving C this objective YOU MUST MODIFY it according to your problem. See C below where your modifications must be inserted. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C On Return C C f double precision, C objective function value at x, C C flag integer, C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the objective C function. (For example, trying to compute the square root C of a negative number, dividing by zero or a very small C number, etc.) If everything was o.k. you must set it C equal to zero. flag = 0 f = - x(4) end C ****************************************************************** C ****************************************************************** subroutine evalg(n,x,g,flag) implicit none C SCALAR ARGUMENTS integer flag,n C ARRAY ARGUMENTS double precision g(n),x(n) C This subroutine must compute the gradient vector of the objective C function. For achieving these objective YOU MUST MODIFY it in the C way specified below. However, if you decide to use numerical C derivatives (we dont encourage this option at all!) you dont need C to modify evalg. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C On Return C C g double precision g(n), C gradient vector of the objective function evaluated at x, C C flag integer, C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of any component C of the gradient vector. (For example, trying to compute C the square root of a negative number, dividing by zero or C a very small number, etc.) If everything was o.k. you C must set it equal to zero. C LOCAL SCALAR integer i flag = 0 do i = 1,n g(i) = 0.0d0 end do g(4) = - 1.0d0 end C ****************************************************************** C ****************************************************************** subroutine evalh(n,x,hlin,hcol,hval,hnnz,flag) implicit none C SCALAR ARGUMENTS integer flag,n,hnnz C ARRAY ARGUMENTS integer hcol(*),hlin(*) double precision hval(*),x(n) C This subroutine might compute the Hessian matrix of the objective C function. For achieving this objective YOU MAY MODIFY it according C to your problem. To modify this subroutine IS NOT MANDATORY. See C below where your modifications must be inserted. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C On Return C C hnnz integer, C number of perhaps-non-null elements of the computed C Hessian, C C hlin integer hlin(hnnz), C see below, C C hcol integer hcol(hnnz), C see below, C C hval double precision hval(hnnz), C the non-null value of the (hlin(k),hcol(k)) position C of the Hessian matrix of the objective function must C be saved at hval(k). Just the lower triangular part of C Hessian matrix must be computed, C C flag integer, C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the Hessian C matrix of the objective funtion. (For example, trying C to compute the square root of a negative number, C dividing by zero or a very small number, etc.) If C everything was o.k. you must set it equal to zero. flag = 0 hnnz = 0 end C ****************************************************************** C ****************************************************************** subroutine evalc(n,x,ind,c,flag) implicit none C SCALAR ARGUMENTS integer ind,flag,n double precision c C ARRAY ARGUMENTS double precision x(n) C This subroutine must compute the ind-th constraint of your C problem. For achieving this objective YOU MUST MOFIFY it C according to your problem. See below the places where your C modifications must be inserted. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C ind integer, C index of the constraint to be computed, C C On Return C C c double precision, C ind-th constraint evaluated at x, C C flag integer C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the C constraint. (For example, trying to compute the square C root of a negative number, dividing by zero or a very C small number, etc.) If everything was o.k. you must set C it equal to zero. C COMMON ARRAYS double precision k(4) C COMMON BLOCKS common /probdata/ k flag = 0 if ( ind .eq. 1 ) then c = x(1) + k(1) * x(1) * x(5) - 1.0d0 else if ( ind .eq. 2 ) then c = x(2) - x(1) + k(2) * x(2) * x(6) else if ( ind .eq. 3 ) then c = x(3) + x(1) + k(3) * x(3) * x(5) - 1.0d0 else if ( ind .eq. 4 ) then c = x(4) - x(3) + x(2) - x(1) + k(4) * x(4) * x(6) else ! if ( ind .eq. 5 ) then c = sqrt( x(5) ) + sqrt( x(6) ) - 4.0d0 end if end C ****************************************************************** C ****************************************************************** subroutine evaljac(n,x,ind,indjac,valjac,jcnnz,flag) implicit none C SCALAR ARGUMENTS integer flag,ind,n,jcnnz C ARRAY ARGUMENTS integer indjac(n) double precision x(n),valjac(n) C This subroutine must compute the gradient of the constraint ind. C For achieving these objective YOU MUST MODIFY it in the way C specified below. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C ind integer, C index of the constraint whose gradient will be computed, C C On Return C C jcnnz integer, C number of perhaps-non-null elements of the computed C gradient, C C indjac integer indjac(jcnnz), C see below, C C valjac double precision valjac(jcnnz), C the non-null value of the partial derivative of the C ind-th constraint with respect to the indjac(k)-th C variable must be saved at valjac(k). C C flag integer C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the C constraint. (For example, trying to compute the square C root of a negative number, dividing by zero or a very C small number, etc.) If everything was o.k. you must set C it equal to zero. C COMMON ARRAYS double precision k(4) C COMMON BLOCKS common /probdata/ k flag = 0 if ( ind .eq. 1 ) then jcnnz = 2 indjac(1) = 1 valjac(1) = 1.0d0 + k(1) * x(5) indjac(2) = 5 valjac(2) = k(1) * x(1) else if ( ind .eq. 2 ) then jcnnz = 3 indjac(1) = 1 valjac(1) = - 1.0d0 indjac(2) = 2 valjac(2) = 1.0d0 + k(2) * x(6) indjac(3) = 6 valjac(3) = k(2) * x(2) else if ( ind .eq. 3 ) then jcnnz = 3 indjac(1) = 1 valjac(1) = 1.0d0 indjac(2) = 3 valjac(2) = 1.0d0 + k(3) * x(5) indjac(3) = 5 valjac(3) = k(3) * x(3) else if ( ind .eq. 4 ) then jcnnz = 5 indjac(1) = 1 valjac(1) = - 1.0d0 indjac(2) = 2 valjac(2) = 1.0d0 indjac(3) = 3 valjac(3) = - 1.0d0 indjac(4) = 4 valjac(4) = 1.0d0 + k(4) * x(6) indjac(5) = 6 valjac(5) = k(4) * x(4) else ! if ( ind .eq. 5 ) then jcnnz = 2 indjac(1) = 5 valjac(1) = 0.5d0 * x(5) ** ( - 0.5d0 ) indjac(2) = 6 valjac(2) = 0.5d0 * x(6) ** ( - 0.5d0 ) end if end C ****************************************************************** C ****************************************************************** subroutine evalhc(n,x,ind,hclin,hccol,hcval,hcnnz,flag) implicit none C SCALAR ARGUMENTS integer flag,ind,n,hcnnz C ARRAY ARGUMENTS integer hccol(*),hclin(*) double precision hcval(*),x(n) C This subroutine might compute the Hessian matrix of the ind-th C constraint. For achieving this objective YOU MAY MODIFY it C according to your problem. To modify this subroutine IS NOT C MANDATORY. See below where your modifications must be inserted. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C ind integer, C index of the constraint whose Hessian will be computed, C C On Return C C hcnnz integer, C number of perhaps-non-null elements of the computed C Hessian, C C hclin integer hclin(hcnnz), C see below, C C hccol integer hccol(hcnnz), C see below, C C hcval double precision hcval(hcnnz), C the non-null value of the (hclin(k),hccol(k)) position C of the Hessian matrix of the ind-th constraint must C be saved at hcval(k). Just the lower triangular part of C Hessian matrix must be computed, C C flag integer, C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the Hessian C matrix of the ind-th constraint. (For example, trying C to compute the square root of a negative number, C dividing by zero or a very small number, etc.) If C everything was o.k. you must set it equal to zero. C COMMON ARRAYS double precision k(4) C COMMON BLOCKS common /probdata/ k flag = 0 if ( ind .eq. 1 ) then hcnnz = 1 hclin(1) = 5 hccol(1) = 1 hcval(1) = k(1) else if ( ind .eq. 2 ) then hcnnz = 1 hclin(1) = 6 hccol(1) = 2 hcval(1) = k(2) else if ( ind .eq. 3 ) then hcnnz = 1 hclin(1) = 5 hccol(1) = 3 hcval(1) = k(3) else if ( ind .eq. 4 ) then hcnnz = 1 hclin(1) = 6 hccol(1) = 4 hcval(1) = k(4) else ! if ( ind .eq. 5 ) then hcnnz = 2 hclin(1) = 5 hccol(1) = 5 hcval(1) = - 0.25d0 * x(5) ** ( - 1.5d0 ) hclin(2) = 6 hccol(2) = 6 hcval(2) = - 0.25d0 * x(6) ** ( - 1.5d0 ) end if end C ****************************************************************** C ****************************************************************** subroutine evalhlp(n,x,m,lambda,p,hp,goth,flag) implicit none C SCALAR ARGUMENTS logical goth integer flag,m,n C ARRAY ARGUMENTS double precision hp(n),lambda(m),p(n),x(n) C This subroutine might compute the product of the Hessian of the C Lagrangian times vector p (just the Hessian of the objective C function in the unconstrained or bound-constrained case). C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C m integer, C number of constraints, C C lambda double precision lambda(m), C vector of Lagrange multipliers, C C p double precision p(n), C vector of the matrix-vector product, C C goth logical, C can be used to indicate if the Hessian matrices were C computed at the current point. It is set to .false. C by the optimization method every time the current C point is modified. Sugestion: if its value is .false. C then compute the Hessians, save them in a common C structure and set goth to .true.. Otherwise, just use C the Hessians saved in the common block structure, C C On Return C C hp double precision hp(n), C Hessian-vector product, C C goth logical, C see above, C C flag integer, C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the C Hessian-vector product. (For example, trying to compute C the square root of a negative number, dividing by zero C or a very small number, etc.) If everything was o.k. you C must set it equal to zero. flag = - 1 end C ****************************************************************** C ****************************************************************** subroutine endp(n,x,l,u,m,lambda,rho,equatn,linear) implicit none C SCALAR ARGUMENTS integer m,n C ARRAY ARGUMENTS logical equatn(m),linear(m) double precision l(n),lambda(m),rho(m),u(n),x(n) C This subroutine can be used to do some extra job after the solver C has found the solution,like some extra statistics, or to save the C solution in some special format or to draw some graphical C representation of the solution. If the information given by the C solver is enough for you then leave the body of this subroutine C empty. C C Parameters of the subroutine: C C The paraemters of this subroutine are the same parameters of C subroutine inip. But in this subroutine there are not output C parameter. All the parameters are input parameters. C LOCAL SCALARS C integer i,flag C double precision c,f C write(*,*) 'Solution: ' C do i = 1,n C write(*,100) i,x(i) C end do C call evalf(n,x,f,flag) C write(*,200) f C write(*,*) 'Value of the constraints: ' C do i = 1,m C call evalc(n,x,i,c,flag) C write(*,300) i,c C end do C write(*,*) 'Data from H.S. Ryoo and N.V. Sahinidis (2005): ' C x(1) = 0.771462d0 C x(2) = 0.516997d0 C x(3) = 0.204234d0 C x(4) = 0.388812d0 C x(5) = 3.036504d0 C x(6) = 5.096052d0 C write(*,*) 'Reported solution: ' C do i = 1,n C write(*,100) i,x(i) C end do C call evalf(n,x,f,flag) C write(*,200) f C write(*,*) 'Value of the constraints: ' C do i = 1,m C call evalc(n,x,i,c,flag) C write(*,300) i,c C end do C NON EXECUTABLE STATEMENTS C100 format(1X,'x(',I1,') = ',1P,D24.16) C200 format(1X,'Objective function value = ',1P,D24.16) C300 format(1X,'c(',I1,') = ',1P,D24.16) end C ****************************************************************** C ****************************************************************** subroutine addc(n,l,u,f,ml,meq,alin,acol,aval,annz,b,fknown,mla, +annza,ind) C SCALAR ARGUMENTS integer annz,annza,meq,ml,mla,n double precision f,fknown C ARRAY ARGUMENTS integer alin(annz),acol(annz),ind(ml) double precision aval(annz),b(ml),l(n),u(n) C COMMON ARRAYS double precision k(4) C LOCAL SCALARS double precision tmp5,tmp6 C COMMON BLOCKS common /probdata/ k C Under and over linear estimations of equality constraints 1 C to 4 plus linear understimation of inequality constraint 5. mla = ml annza = annz C ( 1.0d0 + k(1) * l(5) ) * x(1) + k(1) * l(1) * x(5) - C ( 1.0d0 + k(1) * l(1) * l(5) ) <= 0 ind(mla+1) = 1 alin(annza+1) = mla + 1 acol(annza+1) = 1 aval(annza+1) = - ( 1.0d0 + k(1) * l(5) ) alin(annza+2) = mla + 1 acol(annza+2) = 5 aval(annza+2) = - k(1) * l(1) b(mla+1) = - ( 1.0d0 + k(1) * l(1) * l(5) ) mla = mla + 1 annza = annza + 2 C ( 1.0d0 + k(1) * u(5) ) * x(1) + k(1) * u(1) * x(5) - C ( 1.0d0 + k(1) * u(1) * u(5) ) <= 0 ind(mla+1) = 1 alin(annza+1) = mla + 1 acol(annza+1) = 1 aval(annza+1) = - ( 1.0d0 + k(1) * u(5) ) alin(annza+2) = mla + 1 acol(annza+2) = 5 aval(annza+2) = - k(1) * u(1) b(mla+1) = - ( 1.0d0 + k(1) * u(1) * u(5) ) mla = mla + 1 annza = annza + 2 C - ( 1.0d0 + k(1) * l(5) ) * x(1) - k(1) * u(1) * x(5) + C ( 1.0d0 + k(1) * u(1) * l(5) ) <= 0 ind(mla+1) = - 1 alin(annza+1) = mla + 1 acol(annza+1) = 1 aval(annza+1) = 1.0d0 + k(1) * l(5) alin(annza+2) = mla + 1 acol(annza+2) = 5 aval(annza+2) = k(1) * u(1) b(mla+1) = 1.0d0 + k(1) * u(1) * l(5) mla = mla + 1 annza = annza + 2 C - ( 1.0d0 + k(1) * u(5) ) * x(1) - k(1) * l(1) * x(5) + C ( 1.0d0 + k(1) * l(1) * u(5) ) <= 0 ind(mla+1) = - 1 alin(annza+1) = mla + 1 acol(annza+1) = 1 aval(annza+1) = 1.0d0 + k(1) * u(5) alin(annza+2) = mla + 1 acol(annza+2) = 5 aval(annza+2) = k(1) * l(1) b(mla+1) = 1.0d0 + k(1) * l(1) * u(5) mla = mla + 1 annza = annza + 2 C - x(1) + ( 1.0d0 + k(2) * l(6) ) * x(2) + ( k(2) * l(2) ) * x(6) - C k(2) * l(2) * l(6) <= 0 ind(mla+1) = 2 alin(annza+1) = mla + 1 acol(annza+1) = 1 aval(annza+1) = 1.0d0 alin(annza+2) = mla + 1 acol(annza+2) = 2 aval(annza+2) = - ( 1.0d0 + k(2) * l(6) ) alin(annza+3) = mla + 1 acol(annza+3) = 6 aval(annza+3) = - k(2) * l(2) b(mla+1) = - k(2) * l(2) * l(6) mla = mla + 1 annza = annza + 3 C - x(1) + ( 1.0d0 + k(2) * u(6) ) * x(2) + ( k(2) * u(2) ) * x(6) - C k(2) * u(2) * u(6) <= 0 ind(mla+1) = 2 alin(annza+1) = mla + 1 acol(annza+1) = 1 aval(annza+1) = 1.0d0 alin(annza+2) = mla + 1 acol(annza+2) = 2 aval(annza+2) = - ( 1.0d0 + k(2) * u(6) ) alin(annza+3) = mla + 1 acol(annza+3) = 6 aval(annza+3) = - k(2) * u(2) b(mla+1) = - k(2) * u(2) * u(6) mla = mla + 1 annza = annza + 3 C x(1) - ( 1.0d0 + k(2) * l(6) ) * x(2) - ( k(2) * u(2) ) * x(6) + C k(2) * u(2) * l(6) <= 0 ind(mla+1) = - 2 alin(annza+1) = mla + 1 acol(annza+1) = 1 aval(annza+1) = - 1.0d0 alin(annza+2) = mla + 1 acol(annza+2) = 2 aval(annza+2) = 1.0d0 + k(2) * l(6) alin(annza+3) = mla + 1 acol(annza+3) = 6 aval(annza+3) = k(2) * u(2) b(mla+1) = k(2) * u(2) * l(6) mla = mla + 1 annza = annza + 3 C x(1) - ( 1.0d0 + k(2) * u(6) ) * x(2) - ( k(2) * l(2) ) * x(6) + C k(2) * l(2) * u(6) <= 0 ind(mla+1) = - 2 alin(annza+1) = mla + 1 acol(annza+1) = 1 aval(annza+1) = - 1.0d0 alin(annza+2) = mla + 1 acol(annza+2) = 2 aval(annza+2) = 1.0d0 + k(2) * u(6) alin(annza+3) = mla + 1 acol(annza+3) = 6 aval(annza+3) = k(2) * l(2) b(mla+1) = k(2) * l(2) * u(6) mla = mla + 1 annza = annza + 3 C x(1) + ( 1.0d0 + k(3) * l(5) ) * x(3) + ( k(3) * l(3) ) * x(5) - C ( 1.0d0 + k(3) * l(3) * l(5) ) <= 0 ind(mla+1) = 3 alin(annza+1) = mla + 1 acol(annza+1) = 1 aval(annza+1) = - 1.0d0 alin(annza+2) = mla + 1 acol(annza+2) = 3 aval(annza+2) = - ( 1.0d0 + k(3) * l(5) ) alin(annza+3) = mla + 1 acol(annza+3) = 5 aval(annza+3) = - k(3) * l(3) b(mla+1) = - ( 1.0d0 + k(3) * l(3) * l(5) ) mla = mla + 1 annza = annza + 3 C x(1) + ( 1.0d0 + k(3) * u(5) ) * x(3) + ( k(3) * u(3) ) * x(5) - C ( 1.0d0 + k(3) * u(3) * u(5) ) <= 0 ind(mla+1) = 3 alin(annza+1) = mla + 1 acol(annza+1) = 1 aval(annza+1) = - 1.0d0 alin(annza+2) = mla + 1 acol(annza+2) = 3 aval(annza+2) = - ( 1.0d0 + k(3) * u(5) ) alin(annza+3) = mla + 1 acol(annza+3) = 5 aval(annza+3) = - k(3) * u(3) b(mla+1) = - ( 1.0d0 + k(3) * u(3) * u(5) ) mla = mla + 1 annza = annza + 3 C - x(1) - ( 1.0d0 + k(3) * l(5) ) * x(3) - ( k(3) * u(3) ) * x(5) + C 1.0d0 + k(3) * u(3) * l(5) <= 0 ind(mla+1) = - 3 alin(annza+1) = mla + 1 acol(annza+1) = 1 aval(annza+1) = 1.0d0 alin(annza+2) = mla + 1 acol(annza+2) = 3 aval(annza+2) = 1.0d0 + k(3) * l(5) alin(annza+3) = mla + 1 acol(annza+3) = 5 aval(annza+3) = k(3) * u(3) b(mla+1) = 1.0d0 + k(3) * u(3) * l(5) mla = mla + 1 annza = annza + 3 C - x(1) - ( 1.0d0 + k(3) * u(5) ) * x(3) - ( k(3) * l(3) ) * x(5) + C 1.0d0 + k(3) * l(3) * u(5) <= 0 ind(mla+1) = - 3 alin(annza+1) = mla + 1 acol(annza+1) = 1 aval(annza+1) = 1.0d0 alin(annza+2) = mla + 1 acol(annza+2) = 3 aval(annza+2) = 1.0d0 + k(3) * u(5) alin(annza+3) = mla + 1 acol(annza+3) = 5 aval(annza+3) = k(3) * l(3) b(mla+1) = 1.0d0 + k(3) * l(3) * u(5) mla = mla + 1 annza = annza + 3 C - x(1) + x(2) - x(3) + x(4) + C k(4) * ( l(4) * x(6) + l(6) * x(4) - l(4) * l(6) ) <= 0 ind(mla+1) = 4 alin(annza+1) = mla + 1 acol(annza+1) = 1 aval(annza+1) = 1.0d0 alin(annza+2) = mla + 1 acol(annza+2) = 2 aval(annza+2) = - 1.0d0 alin(annza+3) = mla + 1 acol(annza+3) = 3 aval(annza+3) = 1.0d0 alin(annza+4) = mla + 1 acol(annza+4) = 4 aval(annza+4) = - ( 1.0d0 + k(4) * l(6) ) alin(annza+5) = mla + 1 acol(annza+5) = 6 aval(annza+5) = - k(4) * l(4) b(mla+1) = - k(4) * l(4) * l(6) mla = mla + 1 annza = annza + 5 C - x(1) + x(2) - x(3) + x(4) + C k(4) * ( u(4) * x(6) + u(6) * x(4) - u(4) * u(6) ) <= 0 ind(mla+1) = 4 alin(annza+1) = mla + 1 acol(annza+1) = 1 aval(annza+1) = 1.0d0 alin(annza+2) = mla + 1 acol(annza+2) = 2 aval(annza+2) = - 1.0d0 alin(annza+3) = mla + 1 acol(annza+3) = 3 aval(annza+3) = 1.0d0 alin(annza+4) = mla + 1 acol(annza+4) = 4 aval(annza+4) = - ( 1.0d0 + k(4) * u(6) ) alin(annza+5) = mla + 1 acol(annza+5) = 6 aval(annza+5) = - k(4) * u(4) b(mla+1) = - k(4) * u(4) * u(6) mla = mla + 1 annza = annza + 5 C x(1) - x(2) + x(3) - x(4) - C k(4) * ( u(4) * x(6) + l(6) * x(4) - u(4) * l(6) ) <= 0 ind(mla+1) = - 4 alin(annza+1) = mla + 1 acol(annza+1) = 1 aval(annza+1) = - 1.0d0 alin(annza+2) = mla + 1 acol(annza+2) = 2 aval(annza+2) = 1.0d0 alin(annza+3) = mla + 1 acol(annza+3) = 3 aval(annza+3) = - 1.0d0 alin(annza+4) = mla + 1 acol(annza+4) = 4 aval(annza+4) = 1.0d0 + k(4) * l(6) alin(annza+5) = mla + 1 acol(annza+5) = 6 aval(annza+5) = k(4) * u(4) b(mla+1) = k(4) * u(4) * l(6) mla = mla + 1 annza = annza + 5 C x(1) - x(2) + x(3) - x(4) - C k(4) * ( l(4) * x(6) + u(6) * x(4) - l(4) * u(6) ) <= 0 ind(mla+1) = - 4 alin(annza+1) = mla + 1 acol(annza+1) = 1 aval(annza+1) = - 1.0d0 alin(annza+2) = mla + 1 acol(annza+2) = 2 aval(annza+2) = 1.0d0 alin(annza+3) = mla + 1 acol(annza+3) = 3 aval(annza+3) = - 1.0d0 alin(annza+4) = mla + 1 acol(annza+4) = 4 aval(annza+4) = 1.0d0 + k(4) * u(6) alin(annza+5) = mla + 1 acol(annza+5) = 6 aval(annza+5) = k(4) * l(4) b(mla+1) = k(4) * l(4) * u(6) mla = mla + 1 annza = annza + 5 C Understimation of the concave functions sqrt( x(5) ) and sqrt( x(6) ) tmp5 = ( sqrt( u(5) ) - sqrt( l(5) ) ) / ( u(5) - l(5) ) tmp6 = ( sqrt( u(6) ) - sqrt( l(6) ) ) / ( u(6) - l(6) ) ind(mla+1) = 5 alin(annza+1) = mla + 1 acol(annza+1) = 5 aval(annza+1) = - tmp5 alin(annza+2) = mla + 1 acol(annza+2) = 6 aval(annza+2) = - tmp6 b(mla+1) = sqrt( l(5) ) - l(5) * tmp5 + + sqrt( l(6) ) - l(6) * tmp6 - 4.0d0 mla = mla + 1 annza = annza + 2 C - x(4) - f <= 0 ind(mla+1) = 0 alin(annza+1) = mla + 1 acol(annza+1) = 4 aval(annza+1) = 1.0d0 b(mla+1) = - min( f, fknown ) mla = mla + 1 annza = annza + 1 end C ****************************************************************** C ****************************************************************** subroutine setlv(n,linvar) implicit none C SCALAR ARGUMENTS integer n C ARRAY ARGUMENTS logical linvar(n) C LOCAL SCALARS integer i do i = 1,n linvar(i) = .false. end do end C ****************************************************************** C ****************************************************************** subroutine evalfi(n,y,f,flag) implicit none C SCALAR ARGUMENTS integer flag,n double precision f(2) C ARRAY ARGUMENTS double precision y(2*n) C This subroutine computes the intervalar objective function. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C y double precision y(2*n), C bounds on x variables, y(2*i-1) <= x(i) <= y(2*i), C C On Return C C f double precision, C objective function value at x, C C flag integer, C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the objective C function. (For example, trying to compute the square root C of a negative number, dividing by zero or a very small C number, etc.) If everything was o.k. you must set it C equal to zero. C EXTERNAL FUNCTIONS integer P flag = 0 C f = - x(4) call sclmlt(-1.0d0,y(P(4)),f) end C ****************************************************************** C ****************************************************************** subroutine evalhi(n,y,hlin,hcol,hval,hnnz,flag) implicit none C SCALAR ARGUMENTS integer flag,n,hnnz C ARRAY ARGUMENTS integer hcol(*),hlin(*) double precision hval(*),y(2*n) C This subroutine computes the intervalar Hessian matrix of the C objective function. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C y double precision y(2*n), C bounds on x variables, y(2*i-1) <= x(i) <= y(2*i), C C On Return C C hnnz integer, C number of perhaps-non-null elements of the computed C Hessian, C C hlin integer hlin(hnnz), C see below, C C hcol integer hcol(hnnz), C see below, C C hval double precision hval(hnnz), C the interval of the non-null value of the (hlin(k),hcol(k)) C position of the Hessian matrix of the objective function must C be saved at hval(P(k)). Just the lower triangular part of c Hessian matrix must be computed, C C flag integer, C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the Hessian C matrix of the objective funtion. (For example, trying C to compute the square root of a negative number, C dividing by zero or a very small number, etc.) If C everything was o.k. you must set it equal to zero. flag = 0 hnnz = 0 end C ****************************************************************** C ****************************************************************** subroutine evalci(n,y,ind,c,flag) implicit none C SCALAR ARGUMENTS integer ind,flag,n double precision c(2) C ARRAY ARGUMENTS double precision y(2*n) C This subroutine computes the interval of the ind-th constraint of C your problem. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C y double precision y(2*n), C bounds on x variables, y(2*i-1) <= x(i) <= y(2*i), C C ind integer, C index of the constraint to be computed, C C On Return C C c double precision, C interval of the ind-th constraint "evaluated at y", C C flag integer C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the C constraint. (For example, trying to compute the square C root of a negative number, dividing by zero or a very C small number, etc.) If everything was o.k. you must set C it equal to zero. C COMMON ARRAYS double precision k(4) C LOCAL ARRAYS double precision tmp(2) C COMMON BLOCKS common /probdata/ k C EXTERNAL FUNCTION integer P flag = 0 if ( ind .eq. 1 ) then C c = x(1) + k(1) * x(1) * x(5) - 1.0d0 call ivli(y(P(1)),c) call mult(y(P(1)),y(P(5)),tmp) call sclmlt(k(1),tmp,tmp) call add(c,tmp,c) call scladd(-1.0d0,c,c) else if ( ind .eq. 2 ) then C c = x(2) - x(1) + k(2) * x(2) * x(6) call ivli(y(P(2)),c) call sub(c,y(P(1)),c) call mult(y(P(2)),y(P(6)),tmp) call sclmlt(k(2),tmp,tmp) call add(c,tmp,c) else if ( ind .eq. 3 ) then C c = x(3) + x(1) + k(3) * x(3) * x(5) - 1.0d0 call ivli(y(P(3)),c) call add(c,y(P(1)),c) call mult(y(P(3)),y(P(5)),tmp) call sclmlt(k(3),tmp,tmp) call add(c,tmp,c) call scladd(-1.0d0,c,c) else if ( ind .eq. 4 ) then C c = x(4) - x(3) + x(2) - x(1) + k(4) * x(4) * x(6) call ivli(y(P(4)),c) call sub(c,y(P(3)),c) call add(c,y(P(2)),c) call sub(c,y(P(1)),c) call mult(y(P(4)),y(P(6)),tmp) call sclmlt(k(4),tmp,tmp) call add(c,tmp,c) else ! if ( ind .eq. 5 ) then C c = sqrt( x(5) ) + sqrt( x(6) ) - 4.0d0 call isqrt(y(P(5)),c) call isqrt(y(P(6)),tmp) call add(c,tmp,c) call scladd(-4.0d0,c,c) end if end C ****************************************************************** C ****************************************************************** subroutine evaljaci(n,y,ind,indjac,valjac,jcnnz,flag) implicit none C SCALAR ARGUMENTS integer flag,ind,n,jcnnz C ARRAY ARGUMENTS integer indjac(n) double precision y(2*n),valjac(2*n) C This subroutine must compute the intervalar gradient of the C constraint ind. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C y double precision y(2*n), C bounds on x variables, y(2*i-1) <= x(i) <= y(2*i), C C ind integer, C index of the constraint whose intervalar gradient will be C computed, C C On Return C C jcnnz integer, C number of perhaps-non-null elements of the computed C gradient, C C indjac integer indjac(jcnnz), C see below, C C valjac double precision valjac(jcnnz), C the interval of the non-null value of the partial C derivative of the ind-th constraint with respect to the C indjac(k)-th variable must be saved at valjac(P(k)). C C flag integer C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the C constraint. (For example, trying to compute the square C root of a negative number, dividing by zero or a very C small number, etc.) If everything was o.k. you must set C it equal to zero. C COMMON ARRAYS double precision k(4) C LOCAL ARRAYS double precision tmppow(2) C COMMON BLOCKS common /probdata/ k C EXTERNAL FUNCTION integer P flag = 0 if ( ind .eq. 1 ) then jcnnz = 2 indjac(1) = 1 C valjac(1) = 1.0d0 + k(1) * x(5) call sclmlt(k(1),y(P(5)),valjac(P(1))) call scladd(1.0d0,valjac(P(1)),valjac(P(1))) indjac(2) = 5 C valjac(2) = k(1) * x(1) call sclmlt(k(1),y(P(1)),valjac(P(2))) else if ( ind .eq. 2 ) then jcnnz = 3 indjac(1) = 1 C valjac(1) = - 1.0d0 call ivl1(-1.0d0,valjac(P(1))) indjac(2) = 2 C valjac(2) = 1.0d0 + k(2) * x(6) call sclmlt(k(2),y(P(6)),valjac(P(2))) call scladd(1.0d0,valjac(P(2)),valjac(P(2))) indjac(3) = 6 C valjac(3) = k(2) * x(2) call sclmlt(k(2),y(P(2)),valjac(P(3))) else if ( ind .eq. 3 ) then jcnnz = 3 indjac(1) = 1 C valjac(1) = 1.0d0 call ivl1(1.0d0,valjac(P(1))) indjac(2) = 3 C valjac(2) = 1.0d0 + k(3) * x(5) call sclmlt(k(3),y(P(5)),valjac(P(2))) call scladd(1.0d0,valjac(P(2)),valjac(P(2))) indjac(3) = 5 C valjac(3) = k(3) * x(3) call sclmlt(k(3),y(P(3)),valjac(P(3))) else if ( ind .eq. 4 ) then jcnnz = 5 indjac(1) = 1 C valjac(1) = - 1.0d0 call ivl1(-1.0d0,valjac(P(1))) indjac(2) = 2 C valjac(2) = 1.0d0 call ivl1(1.0d0,valjac(P(2))) indjac(3) = 3 C valjac(3) = - 1.0d0 call ivl1(-1.0d0,valjac(P(3))) indjac(4) = 4 C valjac(4) = 1.0d0 + k(4) * x(6) call sclmlt(k(4),y(P(6)),valjac(P(4))) call scladd(1.0d0,valjac(P(4)),valjac(P(4))) indjac(5) = 6 C valjac(5) = k(4) * x(4) call sclmlt(k(4),y(P(4)),valjac(P(5))) else ! if ( ind .eq. 5 ) then jcnnz = 2 indjac(1) = 5 C valjac(1) = 0.5d0 * x(5) ** ( - 0.5d0 ) call ivl1(-0.5d0,tmppow) call iipowr(y(P(5)),tmppow,valjac(P(1))) call sclmlt(0.5d0,valjac(P(1)),valjac(P(1))) indjac(2) = 6 C valjac(2) = 0.5d0 * x(6) ** ( - 0.5d0 ) call ivl1(-0.5d0,tmppow) call iipowr(y(P(6)),tmppow,valjac(P(2))) call sclmlt(0.5d0,valjac(P(2)),valjac(P(2))) end if end C ****************************************************************** C ****************************************************************** subroutine evalhci(n,y,ind,hclin,hccol,hcval,hcnnz,flag) implicit none C SCALAR ARGUMENTS integer flag,ind,n,hcnnz C ARRAY ARGUMENTS integer hccol(*),hclin(*) double precision hcval(*),y(2*n) C This subroutine might compute the intervalar Hessian matrix of C the ind-th constraint. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C y double precision y(2*n), C bounds on x variables, y(2*i-1) <= x(i) <= y(2*i), C C ind integer, C index of the constraint whose Hessian will be computed, C C On Return C C hcnnz integer, C number of perhaps-non-null elements of the computed C Hessian, C C hclin integer hclin(hcnnz), C see below, C C hccol integer hccol(hcnnz), C see below, C C hcval double precision hcval(hcnnz), C the interval of the non-null value of the (hclin(k),hccol(k)) C position of the Hessian matrix of the ind-th constraint must C be saved at hcval(P(k)). Just the lower triangular part of C Hessian matrix must be computed, C C flag integer, C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the Hessian C matrix of the ind-th constraint. (For example, trying C to compute the square root of a negative number, C dividing by zero or a very small number, etc.) If C everything was o.k. you must set it equal to zero. C COMMON ARRAYS double precision k(4) C LOCAL ARRAYS double precision tmppow(2) C COMMON BLOCKS common /probdata/ k C EXTERNAL FUNCTION integer P flag = 0 if ( ind .eq. 1 ) then hcnnz = 1 hclin(1) = 5 hccol(1) = 1 C hcval(1) = k(1) call ivl1(k(1),hcval(P(1))) else if ( ind .eq. 2 ) then hcnnz = 1 hclin(1) = 6 hccol(1) = 2 C hcval(1) = k(2) call ivl1(k(2),hcval(P(1))) else if ( ind .eq. 3 ) then hcnnz = 1 hclin(1) = 5 hccol(1) = 3 C hcval(1) = k(3) call ivl1(k(3),hcval(P(1))) else if ( ind .eq. 4 ) then hcnnz = 1 hclin(1) = 6 hccol(1) = 4 C hcval(1) = k(4) call ivl1(k(4),hcval(P(1))) else ! if ( ind .eq. 5 ) then hcnnz = 2 hclin(1) = 5 hccol(1) = 5 C hcval(1) = - 0.25d0 * x(5) ** ( - 1.5d0 ) call ivl1(-1.5d0,tmppow) call iipowr(y(P(5)),tmppow,hcval(P(1))) call sclmlt(-0.25d0,hcval(P(1)),hcval(P(1))) hclin(2) = 6 hccol(2) = 6 C hcval(2) = - 0.25d0 * x(6) ** ( - 1.5d0 ) call ivl1(-1.5d0,tmppow) call iipowr(y(P(6)),tmppow,hcval(P(2))) call sclmlt(-0.25d0,hcval(P(2)),hcval(P(2))) end if end