C ================================================================= C File: problem03b.f C ================================================================= C ================================================================= C Module: Subroutines that define the problem C ================================================================= C Last update of any of the component of this module: C November 6, 2006. C ================================================================= C Problem 3b from E. G. Birgin, C. A. Floudas and J. M. Martinez, C Global minimization using an Augmented Lagrangian method with C variable lower-level constraints, submitted. C ****************************************************************** C ****************************************************************** subroutine inip(n,x,l,u,m,lambda,rho,equatn,linear,ml,meq,alin, +acol,aval,annz,b) implicit none C SCALAR ARGUMENTS integer annz,m,ml,meq,n 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 = 2 C Lower and upper bounds l(1) = 1.0d-05 l(2) = 1.0d-05 u(1) = 16.0d0 u(2) = 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 = 1 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) = .false. C For each constraint i, set linear(i) = .true. if it is a linear C constraint, otherwise set linear(i) = .false. linear(1) = .false. C Number of lower-level linear constraints ml = 0 meq = 0 annz = 0 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. C COMMON ARRAYS double precision k(4) C LOCAL SCALARS double precision a,b,c,d,t1,t2 C COMMON BLOCKS common /probdata/ k flag = 0 a = 1.0d0 + k(1) * x(1) b = 1.0d0 + k(2) * x(2) c = 1.0d0 + k(3) * x(1) d = 1.0d0 + k(4) * x(2) t1 = a * c * d t2 = a * b * d f = - ( k(1) * x(1) / t1 + k(2) * x(2) / t2 ) 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 COMMON ARRAYS double precision k(4) C LOCAL SCALARS double precision a,b,c,d,t1,t2,q1,q2,q3,q4,q13,q24,p11,p12,p21,p22 C COMMON BLOCKS common /probdata/ k flag = 0 a = 1.0d0 + k(1) * x(1) b = 1.0d0 + k(2) * x(2) c = 1.0d0 + k(3) * x(1) d = 1.0d0 + k(4) * x(2) t1 = a * c * d t2 = a * b * d q1 = k(1) / a q2 = k(2) / b q3 = k(3) / c q4 = k(4) / d q13 = q1 + q3 q24 = q2 + q4 p11 = ( 1.0d0 - x(1) * q13 ) / t1 p22 = ( 1.0d0 - x(2) * q24 ) / t2 p21 = - x(2) * q1 / t2 p12 = - x(1) * q4 / t1 g(1) = - ( k(1) * p11 + k(2) * p21 ) g(2) = - ( k(1) * p12 + k(2) * p22 ) 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. C COMMON ARRAYS double precision k(4) C LOCAL SCALARS double precision a,b,c,d,t1,t2,q1,q2,q3,q4,q13,q24,p111,p222,p211, + p122,p112,p212 C COMMON BLOCKS common /probdata/ k flag = 0 a = 1.0d0 + k(1) * x(1) b = 1.0d0 + k(2) * x(2) c = 1.0d0 + k(3) * x(1) d = 1.0d0 + k(4) * x(2) t1 = a * c * d t2 = a * b * d q1 = k(1) / a q2 = k(2) / b q3 = k(3) / c q4 = k(4) / d q13 = q1 + q3 q24 = q2 + q4 p211 = 2.0d0 * x(2) * q1 ** 2 / t2 p122 = 2.0d0 * x(1) * q4 ** 2 / t1 p112 = - k(4) * ( 1.0d0 - x(1) * q13 ) / ( d * t1 ) p212 = - k(1) * ( 1.0d0 - x(2) * q24 ) / ( a * t2 ) p111 = - 2.0d0 * ( q13 + x(1) * ( q1 * q3 - q13 ** 2 ) ) / t1 p222 = - 2.0d0 * ( q24 + x(2) * ( q2 * q4 - q24 ** 2 ) ) / t2 hnnz = 3 hlin(1) = 1 hcol(1) = 1 hval(1) = - ( k(1) * p111 + k(2) * p211 ) hlin(2) = 2 hcol(2) = 1 hval(2) = - ( k(1) * p112 + k(2) * p212 ) hlin(3) = 2 hcol(3) = 2 hval(3) = - ( k(1) * p122 + k(2) * p222 ) 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. flag = 0 c = sqrt( x(1) ) + sqrt( x(2) ) - 4.0d0 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. flag = 0 jcnnz = 2 indjac(1) = 1 valjac(1) = 0.5d0 * x(1) ** ( - 0.5d0 ) indjac(2) = 2 valjac(2) = 0.5d0 * x(2) ** ( - 0.5d0 ) 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. flag = 0 hcnnz = 2 hclin(1) = 1 hccol(1) = 1 hcval(1) = - 0.25d0 * x(1) ** ( - 1.5d0 ) hclin(2) = 2 hccol(2) = 2 hcval(2) = - 0.25d0 * x(2) ** ( - 1.5d0 ) 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 COMMON ARRAYS C double precision k(4) C LOCAL SCALARS C integer i,flag C double precision c,f C COMMON BLOCKS C common /probdata/ k C write(*,*) 'Converting to the original formulation: ' C x(5) = x(1) C x(6) = x(2) C x(1) = 1.0d0 / ( 1.0d0 + k(1) * x(5) ) C x(2) = x(1) / ( 1.0d0 + k(2) * x(6) ) C x(3) = ( 1.0d0 - x(1) ) / ( 1.0d0 + k(3) * x(5) ) C x(4) = ( x(3) - x(2) + x(1) ) / ( 1.0d0 + K(4) * x(6) ) C n = 6 C m = 5 C write(*,*) 'Solution: ' C do i = 1,n C write(*,100) i,x(i) C end do C call evalfo(n,x,f,flag) C write(*,200) f C write(*,*) 'Value of the constraints: ' C do i = 1,m C call evalco(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 ****************************************************************** C subroutine evalfo(n,x,f,flag) C implicit none C SCALAR ARGUMENTS C integer flag,n C double precision f C ARRAY ARGUMENTS C 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. C flag = 0 C f = - x(4) C end C ****************************************************************** C ****************************************************************** C subroutine evalco(n,x,ind,c,flag) C implicit none C SCALAR ARGUMENTS C integer ind,flag,n C double precision c C ARRAY ARGUMENTS C 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 C double precision k(4) C COMMON BLOCKS C common /probdata/ k C flag = 0 C if ( ind .eq. 1 ) then C c = x(1) + k(1) * x(1) * x(5) - 1.0d0 C else if ( ind .eq. 2 ) then C c = x(2) - x(1) + k(2) * x(2) * x(6) C else if ( ind .eq. 3 ) then C c = x(3) + x(1) + k(3) * x(3) * x(5) - 1.0d0 C else if ( ind .eq. 4 ) then C c = x(4) - x(3) + x(2) - x(1) + k(4) * x(4) * x(6) C else ! if ( ind .eq. 5 ) then C c = sqrt( x(5) ) + sqrt( x(6) ) - 4.0d0 C end if C end C ****************************************************************** C ****************************************************************** subroutine addc(n,l,u,f,ml,meq,alin,acol,aval,annz,b,mla,annza, +ind) C SCALAR ARGUMENTS integer annz,annza,meq,ml,mla,n double precision f C ARRAY ARGUMENTS integer alin(annz),acol(annz),ind(ml) double precision aval(annz),b(ml),l(n),u(n) mla = ml annza = annz 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 COMMON ARRAYS double precision k(4) C LOCAL ARRAYS double precision a(2),b(2),c(2),d(2),t1(2),t2(2),tmp(2) C COMMON BLOCKS common /probdata/ k C EXTERNAL FUNCTIONS integer P flag = 0 C a = 1.0d0 + k(1) * x(1) call sclmlt(k(1),y(P(1)),a) call scladd(1.0d0,a,a) C b = 1.0d0 + k(2) * x(2) call sclmlt(k(2),y(P(2)),b) call scladd(1.0d0,b,b) C c = 1.0d0 + k(3) * x(1) call sclmlt(k(3),y(P(1)),c) call scladd(1.0d0,c,c) C d = 1.0d0 + k(4) * x(2) call sclmlt(k(4),y(P(2)),d) call scladd(1.0d0,d,d) C t1 = a * c * d call mult(a,c,t1) call mult(d,t1,t1) C t2 = a * b * d call mult(a,b,t2) call mult(d,t2,t2) C f = - ( k(1) * x(1) / t1 + k(2) * x(2) / t2 ) call idiv(y(P(1)),t1,f) call sclmlt(k(1),f,f) call idiv(y(P(2)),t2,tmp) call sclmlt(k(2),tmp,tmp) call add(f,tmp,f) call sclmlt(-1.0d0,f,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. C COMMON ARRAYS double precision k(4) C LOCAL ARRAYS double precision a(2),b(2),c(2),d(2),t1(2),t2(2),q1(2),q2(2), + q3(2),q4(2),q13(2),q24(2),p111(2),p112(2),p122(2),p211(2), + p212(2),p222(2),tmp(2) C COMMON BLOCKS common /probdata/ k C EXTERNAL FUNCTIONS integer P flag = 0 C a = 1.0d0 + k(1) * x(1) call sclmlt(k(1),y(P(1)),a) call scladd(1.0d0,a,a) C b = 1.0d0 + k(2) * x(2) call sclmlt(k(2),y(P(2)),b) call scladd(1.0d0,b,b) C c = 1.0d0 + k(3) * x(1) call sclmlt(k(3),y(P(1)),c) call scladd(1.0d0,c,c) C d = 1.0d0 + k(4) * x(2) call sclmlt(k(4),y(P(2)),d) call scladd(1.0d0,d,d) C t1 = a * c * d call mult(a,c,t1) call mult(d,t1,t1) C t2 = a * b * d call mult(a,b,t2) call mult(d,t2,t2) C q1 = k(1) / a call ivl1(k(1),q1) call idiv(q1,a,q1) C q2 = k(2) / b call ivl1(k(2),q2) call idiv(q2,b,q2) C q3 = k(3) / c call ivl1(k(3),q3) call idiv(q3,c,q3) C q4 = k(4) / d call ivl1(k(4),q4) call idiv(q4,d,q4) C q13 = q1 + q3 call add(q1,q3,q13) C q24 = q2 + q4 call add(q2,q4,q24) C p211 = 2.0d0 * x(2) * q1 ** 2 / t2 call power(q1,2,p211) call mult(y(P(2)),p211,p211) call idiv(p211,t2,p211) call sclmlt(2.0d0,p211,p211) C p122 = 2.0d0 * x(1) * q4 ** 2 / t1 call power(q4,2,p122) call mult(y(P(1)),p122,p122) call idiv(p122,t1,p122) call sclmlt(2.0d0,p122,p122) C p112 = - k(4) * ( 1.0d0 - x(1) * q13 ) / ( d * t1 ) call ivl1(1.0d0,p112) call mult(y(P(1)),q13,tmp) call sub(p112,tmp,p112) call sclmlt(-k(4),p112,p112) call mult(d,t1,tmp) call idiv(p112,tmp,p112) C p212 = - k(1) * ( 1.0d0 - x(2) * q24 ) / ( a * t2 ) call ivl1(1.0d0,p212) call mult(y(P(2)),q24,tmp) call sub(p212,tmp,p212) call sclmlt(-k(1),p212,p212) call mult(a,t2,tmp) call idiv(p212,tmp,p212) C p111 = - 2.0d0 * ( q13 + x(1) * ( q1 * q3 - q13 ** 2 ) ) / t1 call mult(q1,q3,p111) call power(q13,2,tmp) call sub(p111,tmp,p111) call mult(y(P(1)),p111,p111) call add(q13,p111,p111) call idiv(p111,t1,p111) call sclmlt(-2.0d0,p111,p111) C p222 = - 2.0d0 * ( q24 + x(2) * ( q2 * q4 - q24 ** 2 ) ) / t2 call mult(q2,q4,p222) call power(q24,2,tmp) call sub(p222,tmp,p222) call mult(y(P(2)),p222,p222) call add(q24,p222,p222) call idiv(p222,t2,p222) call sclmlt(-2.0d0,p222,p222) hnnz = 3 hlin(1) = 1 hcol(1) = 1 C hval(1) = - ( k(1) * p111 + k(2) * p211 ) call sclmlt(-k(1),p111,hval(P(1))) call sclmlt(-k(2),p211,tmp) call add(hval(P(1)),tmp,hval(P(1))) hlin(2) = 2 hcol(2) = 1 C hval(2) = - ( k(1) * p112 + k(2) * p212 ) call sclmlt(-k(1),p112,hval(P(2))) call sclmlt(-k(2),p212,tmp) call add(hval(P(2)),tmp,hval(P(2))) hlin(3) = 2 hcol(3) = 2 C hval(3) = - ( k(1) * p122 + k(2) * p222 ) call sclmlt(-k(1),p122,hval(P(3))) call sclmlt(-k(2),p222,tmp) call add(hval(P(3)),tmp,hval(P(3))) 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 LOCAL ARRAYS double precision tmp(2) C EXTERNAL FUNCTION integer P flag = 0 C c = sqrt( x(1) ) + sqrt( x(2) ) - 4.0d0 call isqrt(y(P(1)),c) call isqrt(y(P(2)),tmp) call add(c,tmp,c) call scladd(-4.0d0,c,c) 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 LOCAL ARRAYS double precision tmppow(2) C EXTERNAL FUNCTION integer P flag = 0 jcnnz = 2 indjac(1) = 1 C valjac(1) = 0.5d0 * x(1) ** ( - 0.5d0 ) call ivl1(-0.5d0,tmppow) call iipowr(y(P(1)),tmppow,valjac(P(1))) call sclmlt(0.5d0,valjac(P(1)),valjac(P(1))) indjac(2) = 2 C valjac(2) = 0.5d0 * x(2) ** ( - 0.5d0 ) call ivl1(-0.5d0,tmppow) call iipowr(y(P(2)),tmppow,valjac(P(2))) call sclmlt(0.5d0,valjac(P(2)),valjac(P(2))) 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 LOCAL ARRAYS double precision tmppow(2) C EXTERNAL FUNCTION integer P flag = 0 hcnnz = 2 hclin(1) = 1 hccol(1) = 1 C hcval(1) = - 0.25d0 * x(1) ** ( - 1.5d0 ) call ivl1(-1.5d0,tmppow) call iipowr(y(P(1)),tmppow,hcval(P(1))) call sclmlt(-0.25d0,hcval(P(1)),hcval(P(1))) hclin(2) = 2 hccol(2) = 2 C hcval(2) = - 0.25d0 * x(2) ** ( - 1.5d0 ) call ivl1(-1.5d0,tmppow) call iipowr(y(P(2)),tmppow,hcval(P(2))) call sclmlt(-0.25d0,hcval(P(2)),hcval(P(2))) end