C ALGORITHM 774, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 23,NO. 3, September, 1997, P. 448--450. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Doc # Drivers # Src # This archive created: Mon Sep 1 10:21:21 1997 export PATH; PATH=/bin:$PATH if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'make.sun' then echo shar: will not over-write existing file "'make.sun'" else cat << \SHAR_EOF > 'make.sun' # Makefile for BCP suite on a SUN FFLAGS = -u .f.o: f77 -c $< SUBS = Drivers/Dp/driver.o Src/Dp/src.o bcpver: $(SUBS) f77 $(FFLAGS) $(SUBS) -o SUN_bcpver clean: /bin/rm -f $(SUBS) core SHAR_EOF fi # end of overwriting check if test -f 'make.ibm' then echo shar: will not over-write existing file "'make.ibm'" else cat << \SHAR_EOF > 'make.ibm' # Makefile for BCP suite on a RISC FFLAGS = -u SUBS = Drivers/Dp/driver.o Src/Dp/src.o .f.o: xlf -c $< bcpver: $(SUBS) xlf $(FFLAGS) $(SUBS) -o RISC_bcpver clean: rm -f $(SUBS) core SHAR_EOF fi # end of overwriting check if test -f 'make.hp' then echo shar: will not over-write existing file "'make.hp'" else cat << \SHAR_EOF > 'make.hp' # Makefile for BCP suite on a HP Apollo Serie 700 FFLAGS = -u .f.o: f77 -c $< SUBS = Drivers/Dp/driver.o Src/Dp/src.o bcpver: $(SUBS) f77 $(FFLAGS) $(SUBS) -o HP_bcpver clean: /bin/rm -f $(SUBS) core SHAR_EOF fi # end of overwriting check if test -f 'README.tex' then echo shar: will not over-write existing file "'README.tex'" else cat << \SHAR_EOF > 'README.tex' % % File: README.tex % Contents: "Using the FORTRAN subroutines for generating box constrained % optimization problems" % Comments: this the README file to the suite of FORTRAN routines % Authors: F. Facchinei, J. Judice and J. Soares % Dates: Oct/95, revised Oct/96, final version Mar/97 % \documentstyle{esub2acm} \def\bkR{{\rm I\kern-.17em R}} \def\bkN{{\rm I\kern-.17em N}} \begin{document} \title{Using the FORTRAN subroutines for generating box constrained optimization problems} \author{Francisco Facchinei\\ Dipartimento di Informatica e Sistemistica, Universit\`a di Roma ``La Sapienza'', Italy. \and Joaquim J\'udice \and Jo\~ao Soares\\ Departamento de Matem\'atica, Universidade de Coimbra, Portugal} \begin{abstract} We explain how to use a set of FORTRAN subroutines for generating box constrained optimization problems as introduced in \cite{soares95,soares95a}. The subroutines made available are to be linked to the user main program. One example of such a main program, which not only exemplifies the usage of the suite of routines but also tests the installation, is also supplied. \end{abstract} \begin{bottomstuff} \begin{authinfo} \address{ Departamento de Matem\'atica, Universidade de Coimbra, 3000 Coimbra, Portugal. Email: {\tt jsoares@mat.uc.pt} } \end{authinfo} \end{bottomstuff} \markboth{Facchinei, J\'udice and Soares} {Using the FORTRAN subroutines} \maketitle \section{\label{software}Implementation and design} The following subroutines have been implemented in ANSI Standard FORTRAN 77. There are no specific instances of machine dependence, except for the maximum size of a common block, which we refer to in Section \ref{lastsub}. Check Appendix~\ref{parameters} for a description of all parameters involved in the only routines that should be explicitly employed by the user, namely, {\tt BCP01}, {\tt BCP02}, {\tt BCP03}, {\tt BCP10}, {\tt BCP11} and {\tt BCP12}. \begin{enumerate} \item subroutine {\tt BCP01}: sets the bounds and the data required to evaluate the objective function $f$ and its derivatives. \item subroutine {\tt BCP02}: prints an error message on the screen. \item subroutine {\tt BCP03}: prints information concerning the problem generated. \item subroutine {\tt BCP04}: alters randomly the position of the elements in a vector of integer components. \item functions {\tt HI0, HI1, HI2}: compute the value, the first derivative and the second order derivative of the function $h_i$, respectively. \item subroutine {\tt BCP10}: computes the objective function $f$ at a given point $x$. \item subroutine {\tt BCP11}: computes the gradient vector of the objective function $f$ at a given point $x$. \item subroutine {\tt BCP12}: computes the product of the Hessian at a given point $x$ by a vector $d$. \item subroutine {\tt rand\/}: is a portable random number generator \cite{schrage79}. \end{enumerate} The box constrained problem is generated by a call to subroutine {\tt BCP01}, where the user determines the characteristics of the problem generated through the choice of various parameters. However, the exact definition of the problem also depends on some random choices. For example, the user can determine exactly the number of lower bounds and the number of those that are active at the solution, but their placement is determined randomly. On exit from {\tt BCP01}, if the variable {\tt INFRM} is not zero then an error has occurred and the box constrained problem has not been generated. The user may get an appropriate error message by calling subroutine {\tt BCP02}. If the problem has been successfully generated ({\tt INFRM = 0}), the user has access to the vectors {\tt L} and {\tt U} containing the lower and upper bounds (infinity was set to $10^{20}$). Additional information may be obtained through a call to {\tt BCP03}. Subroutines {\tt BCP10}, {\tt BCP11} and {\tt BCP12} provide the value of the generated problem objective function and its derivatives. Subroutine {\tt BCP10} provides in {\tt F} the value of the function $f$ calculated at {\tt X}, {\tt BCP11} stores in {\tt GF} the gradient of $f$ calculated at {\tt X}, and {\tt BCP12} computes in {\tt HFD} the product of the Hessian of $f$ at {\tt X} times the vector {\tt D}. These are the quantities used by most existing codes for the solution of large-scale box constrained nonlinear programs. The generation of the box constrained problem is based on an underlying unconstrained problem. The data concerning the unconstrained problem must be user-supplied through the definition of the following four FORTRAN subroutines: \begin{verbatim} SUBROUTINE FUNCTG(N,X,G) SUBROUTINE GRADG(N,X,GG) SUBROUTINE HESGD(N,X,D,HGD) SUBROUTINE SOLUG(N,XBAR) \end{verbatim} Subroutine {\tt FUNCTG\/} gives the value of the function $g$ at a given point $x$, {\tt GRADG\/} gives the gradient of $g$ at $x$ and {\tt HESGD\/} performs the product of the Hessian of $g$ at $x$ by a vector $d$. Finally, subroutine {\tt SOLUG\/} provides a solution of the unconstrained problem. A detailed description of the parameters follows. \begin{itemize} \item {\tt N} (Input, {\tt INTEGER}) The dimension of the problem. \item {\tt X} (Input, {\tt DOUBLE PRECISION}) A vector of dimension {\tt N} representing a point $x \in \bkR^n$. \item {\tt G} (Output, {\tt DOUBLE PRECISION}) A variable that gives the value of the unconstrained objective function evaluated at {\tt X}. \item {\tt GG} (Output, {\tt DOUBLE PRECISION}) A vector of dimension {\tt N} corresponding to the value of the gradient of the unconstrained objective function evaluated in {\tt X}. \item {\tt D} (Input, {\tt DOUBLE PRECISION}) A vector of dimension {\tt N}. \item {\tt HGD} (Output, {\tt DOUBLE PRECISION}) A variable that provides the product of the Hessian of the unconstrained objective function $g$ evaluated at {\tt X} times the vector {\tt D}. \item {\tt XBAR} (Output, {\tt DOUBLE PRECISION}) A vector of dimension {\tt N} corresponding to the unconstrained minimizer. \end{itemize} \section{\label{install}Installation notes} \subsection{Files} The FORTRAN subroutines that the user should link to his/her code are in the file {\em src.f}. The program {\tt BCPVER} in {\em driver.f} performs a test of the installation and can also give insight on the usage of the subroutines. Since the execution of the test of installation requires the user to provide the four subroutines listed in the previous section, we provide a demo set of such routines in the file {\em bcpex.f}. The code has already been compiled on a SUN SPARC 10, an IBM RS/6000 and an HP Model 712/60. For completeness, we provide each of the {\em makefiles\/} in the files {\em make.hp}, {\em make.ibm} and {\em make.sun}. \subsection{Limitations\label{lastsub}} The FORTRAN constant {\tt NN}, defining the size of the vectors in the common blocks, controls the size of the largest problem that can be generated by the suite. The value of {\tt NN} is currently fixed at 10000 in a {\tt PARAMETER} statement in the subroutines {\tt BCP01}, {\tt BCP03}, {\tt BCP10}, {\tt BCP11} and {\tt BCP12}, and in the program {\tt BCPVER}. The user must increase its value in all the {\tt PARAMETER} statements if he/she wishes to generate larger problems. %\bibliographystyle{esub2acm} %\bibliography{abbrv,titles1,titles2,titles3} \begin{thebibliography}{} \bibitem[\protect\citeauthoryear{Facchinei, J\'udice, and Soares}{Facchinei et~al.}{a}]{soares95} \bibsc{Facchinei, F., J\'udice, J., and Soares, J.} \bibyear{??a}. \newblock {FORTRAN} subroutines for generating box constrained optimization problems. \newblock ACM Trans. Math. Software. \bibitem[\protect\citeauthoryear{Facchinei, J\'udice, and Soares}{Facchinei et~al.}{b}]{soares95a} \bibsc{Facchinei, F., J\'udice, J., and Soares, J.} \bibyear{??b}. \newblock Generating box constrained optimization problems. \newblock ACM Trans. Math. Software. \bibitem[\protect\citeauthoryear{Schrage}{Schrage}{1979}]{schrage79} \bibsc{Schrage, L.} \bibyear{1979}. \newblock A more portable fortran random generator. \newblock \bibemphic{ACM Trans. Math. Software}~\bibemph{5},~2 (June), 132--138. \end{thebibliography} \newpage \appendix \label{parameters} \section{BCP01 parameter list\label{parameters}} \begin{verbatim} SUBROUTINE BCP01( N , L , U , BND , LOWBND, 1 ACTBND, PAIR , DEG , LOWM1 , UPPM1 , 2 LOWM2 , UPPM2 , WIDTH1, WIDTH2, HIF , 3 ISEED , INFRM ) c N (input) INTEGER c The number of variables in the problem. (Must be <= NN, c see common block below) c c L, U (output) DOUBLE PRECISION vectors of dimension N c The bounds of the new problem. c c BND (input) INTEGER c The number of bound constraints (or simply, bounds), expressed as c the percentage of the 2*n possible constraints, that are really c existing (Must be in [0,100]). c Example: if N is 500 and BND=40 then the problem has 400 bounds. c c LOWBND (input) INTEGER c The number of lower bounds, expressed as the percentage of c the number of constraints as determined by BND (Must be in c [0,100]). This quantity also determines the number of upper c bounds. c Example: if N is 500, BND is 40 and LOWBND=75 then there c are 300 (100) lower (upper) bounds. c c ACTBND (input) INTEGER c The number of active bounds at the unconstrained solution, c expressed as the percentage of the number of constraints as c determined by BND (Must be in [0,100]). These active c bounds are distributed proportionally to LOWBND. c Example: Continuing the previous example, if ACTBND is 50 then c 150 (50) of the lower (upper) bounds should be active at the c unconstrained solution. c c PAIR (input) INTEGER c If this parameter is set to 1 then the constraints as c determined by BND must come in pairs, i.e. each variable c has either no constraint or both lower and upper bounds (in c which case LOWBND must be set to 50). If PAIR is c set to 0 then the indices of variables having upper bounds are c chosen randomly and independently of the lower bounds c (Must be 0 or 1). c c DEG (input) INTEGER c The number of Lagrange multipliers (at the unconstrained solution) c uniformly distributed in the interval [LOWM1, UPPM1] (see below), c expressed as the percentage of the active constraints c (Must be in [0,100]). c The remaining active constraints have Lagrange multipliers c uniformly distributed in [LOWM2, UPPM2] (see below). c Example: Continuing the example, if DEG=10 then 15 (5) active c lower (upper) bounds (at the unconstrained solution) have c Lagrange multipliers uniformly distributed in [LOWM1, UPPM1], c while the remaining active bounds have Lagrange multipliers c uniformly distributed in [LOWM2, UPPM2]. c c LOWM1, UPPM1, LOWM2, UPPM2 (input) DOUBLE PRECISION c These parameters specify the ranges for the c Lagrange multipliers according to what explained above. c (Must be >= 0) c c WIDTH1, WIDTH2 (input) DOUBLE PRECISION c These two parameters control the value of the bounds that c are not active at the unconstrained solution. More precisely, c if the i-th component has a lower bound not active at the c unconstrained solution xbar then this bound is randomly chosen in c the interval [xbar_i - WIDTH2, xbar_i - WIDTH1]. Analogously, c if the i-th component has an upper bound not active at c the unconstrained solution xbar then this bound is randomly chosen c in the interval [xbar_i + WIDTH1, xbar_i + WIDTH2] (Must be >= 0). c c HIF (input) INTEGER c This parameter determines which function h_i is to be used by c the generator. If HIF is equal to c c 1: h_i(x_i)= beta_i * ( x_i - xbar_i ); c c 2: h_i(x_i)=( x_i - xbar_i )^3 + beta_i*( x_i - xbar_i ); c c 3: h_i(x_i)=( x_i - xbar_i )^(7/3) + beta_i*( x_i - xbar_i ). c c (Must be 1, 2 or 3) c c ISEED (input) INTEGER c A seed for the random number generator. By running the box c constrained problems generator with the same choice of c parameters but two different seeds, two different box c constrained problems are generated with the same overall c characteristics (Must be >= 0). c c INFRM (output) INTEGER c Determines if a successful call to GENBP has occurred c (INFRM=0) or not (INFRM<>0). \end{verbatim} \section{BCP02 parameter list} \begin{verbatim} SUBROUTINE BCP02(INFRM) c INFRM (input) INTEGER (see above) c Determines if a successful call to BCP01 has occurred c (INFRM=0) or not (INFRM<>0). c If INFRM is equal to c 1: N is less than 1 or greater than NN c 2, 3, 4, 5, 6, 11: Some integer parameter is defined c out of its range. c 7: PAIR=1 but LOWBND is not 50. c 8: Either LOWM1 or LOWM2 is less than zero. c 9: Either LOWM1 is greater than UPPM1 or LOWM2 is greater c than UPPM2. c 10: Either WIDTH1 is less than zero or WIDTH1 is greater c than WIDTH2. c 12: It is impossible to make a selection because LOWBND c is too large. c 13: It is impossible to make a selection because LOWBND c is too small. c 14: It is impossible to make a selection because ACTBND c is too large. \end{verbatim} \section{BCP03 parameter list} \begin{verbatim} SUBROUTINE BCP03( N, L, U ) c N (input) INTEGER c The number of variables in the problem. c c L, U (input) DOUBLE PRECISION vectors of dimension N c The bounds of the new problem. \end{verbatim} \section{BCP10 parameter list} \begin{verbatim} SUBROUTINE BCP10( N, X, F) c N (input) INTEGER c The number of variables in the problem. c c X (input) DOUBLE PRECISION vector of dimension N c The point where the objective function is to be evaluated. c c F (output) DOUBLE PRECISION c The objective function value. \end{verbatim} \section{BCP11 parameter list} \begin{verbatim} SUBROUTINE BCP11( N, X, GF) c N (input) INTEGER c The number of variables in the problem. c c X (input) DOUBLE PRECISION vector of dimension N c The point where the gradient of the objective function is to c be evaluated. c c GF (output) DOUBLE PRECISION vector of dimension N c The gradient of objective function values. \end{verbatim} \section{BCP12 parameter list} \begin{verbatim} SUBROUTINE BCP12( N, X, D, HFD ) c N (input) INTEGER c The number of variables in the problem. c c X (input) DOUBLE PRECISION vector of dimension N c The point where the Hessian matrix of the objective function c is to be evaluated. c c D (input) DOUBLE PRECISION vector of dimension N c The vector that will be multiplied by the Hessian matrix. c c HFD (output) DOUBLE PRECISION vector of dimension N c The resulting vector. \end{verbatim} \end{document} SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << \SHAR_EOF > 'driver.f' C====================================================================== C c C This file contains the FORTRAN program which demonstrates the c C use of the suite of routines for generating box constrained c C nonlinear optimization problems. c C c C Another purpose of this program is to test the integrity of the c C routines in the suite after installation. It is by no means a c C thorough test. c C c C The other routines that should be linked to this main program c C are in files bcpgen.f and bcpex.f c C c C The user is referred to the documentation c C c C Soares, J., Judice, J. and Facchinei, F., "Generating box c C constrained optimization problems", c C ACM TOMS, (??) c C and c C Soares, J., Judice, J. and Facchinei, F., "FORTRAN subroutines c C for generating box constrained optimization problems", c C ACM TOMS, (??). c C c C====================================================================== PROGRAM VERIFY C Authors: F. Facchinei, J. Soares and J. Judice C Date Written: September 10, 1995 C Date Modified: September 3, 1996 C C When testing BCP01 (see comments in subroutine BCP01 for details) C we make use of the following variables: C =================== C Input parameters: C ( N, *, *, BND, LOWBND, ACTBND, PAIR, DEG, LOWM1, UPPM1, C LOWM2, UPPM2, WIDTH1, WIDTH2, HIF, ISEED, * ) C Output parameters: C ( *, L, U, *, *, *, *, *, *, *, *, *, *, *, *, *, INFRM ) C When testing BCP10 (see comments in subroutine BCP10 for details) C we make use of the following variables: C ==================== C Input parameters: C ( N, X, *) C Output parameters: C ( *, *, F) C When testing BCP11 (see comments in subroutine BCP11 for details) C we make use of the following variables: C ==================== C Input parameters: C ( N, X, *) C Output parameters: C ( *, *, GF) C When testing BCP12 (see comments in subroutine BCP12 for details) C we make use of the following variables: C ==================== C Input parameters: C ( N, X, D, *) C Output parameters: C ( *, *, *, HFD) C .. See the comments in BCP01 for a description of the variables C in the common blocks. C .. Parameters .. INTEGER MAXN PARAMETER (MAXN=100) INTEGER NN PARAMETER (NN=10000) C .. C .. Scalars in Common .. INTEGER ACTL,ACTU,HI,INL,INU,NM1L,NM1U,NM2L,NM2U C .. C .. Arrays in Common .. REAL BETA(NN),XBAR(NN) INTEGER PTACT(NN) C .. C .. Local Scalars .. REAL ACUM,F,LOWM1,LOWM2,UPPM1,UPPM2,WIDTH1,WIDTH2 INTEGER ACTBND,BND,DEG,HIF,I,INFRM,ISEED,LOWBND,N,PAIR C .. C .. Local Arrays .. REAL BETA2(MAXN),D(MAXN),GF(MAXN),HFD(MAXN),L(MAXN),L2(MAXN), + U(MAXN),U2(MAXN),X(MAXN) INTEGER PTACT2(MAXN) C .. C .. External Subroutines .. EXTERNAL BCP01,BCP02,BCP10,BCP11,BCP12 C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. C .. Common blocks .. COMMON /BCPB1/INL,INU,ACTL,ACTU,NM1L,NM1U,NM2L,NM2U,HI COMMON /BCPB2/BETA,XBAR,PTACT C .. C N = 20 meaning that the problems have 20 variables C (N need to be less or equal to NN) N = 20 C Initialize X and D for the calls to BCP10, BCP11 and BCP12 DO 10 I = 1,N X(I) = 1.0 + 1.0/I D(I) = 1.0/I 10 CONTINUE C ******************* Generate a problem with 10 lower bounds and C VERIFICATION TEST 1 10 upper bounds, distributed randomly. C ******************* C At the unconstrained solution, 8 bounds should C be active. Among these active bounds, 4 should C have a null multiplier. The remaining 4 C multipliers should be uniformly distributed in C [1.0,2.0]. C C Moreover, at the unconstrained solution the C distance between each nonactive bound and the C correspondent component should be C uniformly distributed in [3.0,50.0]. C C Use the simplest function h_i, C h_i(x_i)= beta_i * ( x_i - xbar_i ), C as described in the companion paper, to C define the objective function. PRINT *,'** VERIFICATION TEST 1 **' PRINT * C BND = 50, meaning that 50% of the all possible bounds (40) C are finite. BND = 50 C LOWBND = 50, meaning that lower and upper bounds are C equally distributed among the finite bounds. LOWBND = 50 C PAIR = 0, meaning that finite bounds need not exist in pairs. PAIR = 0 C ACTBND = 40, meaning that 40% of the finite bounds are active C at the unconstrained solution. ACTBND = 40 C Set the two ranges for the multipliers LOWM1 = 0.0 UPPM1 = 0.0 LOWM2 = 1.0 UPPM2 = 2.0 C DEG = 50, meaning that at the unconstrained solution 50% of the C multipliers associated with the active constraints are in the C first range while the others are in the second range. DEG = 50 C Set the range for the nonactive but finite bounds. WIDTH1 = 3.0 WIDTH2 = 50.0 C Select the function h_i HIF = 1 C Set ISEED ISEED = 1994 C Generate the problem CALL BCP01(N,L,U,BND,LOWBND,ACTBND,PAIR,DEG,LOWM1,UPPM1,LOWM2, + UPPM2,WIDTH1,WIDTH2,HIF,ISEED,INFRM) C Check INFRM IF (INFRM.NE.0) THEN CALL BCP02(INFRM) PRINT *,'** VERIFICATION TEST 1 INTERRUPTED **' STOP END IF C The output from BCP01 should have been the following in C single precision. L2(1) = -0.1000000E+21 L2(2) = -23.95055 L2(3) = 1.000000 L2(4) = -21.96506 L2(5) = 1.000000 L2(6) = -0.1000000E+21 L2(7) = -0.1000000E+21 L2(8) = -34.71118 L2(9) = -0.1000000E+21 L2(10) = -0.1000000E+21 L2(11) = -21.71301 L2(12) = -0.1000000E+21 L2(13) = -0.1000000E+21 L2(14) = -14.06864 L2(15) = -15.64741 L2(16) = 1.000000 L2(17) = -0.1000000E+21 L2(18) = -0.1000000E+21 L2(19) = -0.1000000E+21 L2(20) = 1.000000 U2(1) = 1.000000 U2(2) = 1.000000 U2(3) = 0.1000000E+21 U2(4) = 0.1000000E+21 U2(5) = 11.12918 U2(6) = 1.000000 U2(7) = 0.1000000E+21 U2(8) = 4.876650 U2(9) = 0.1000000E+21 U2(10) = 26.85025 U2(11) = 0.1000000E+21 U2(12) = 0.1000000E+21 U2(13) = 39.96765 U2(14) = 1.000000 U2(15) = 49.16487 U2(16) = 0.1000000E+21 U2(17) = 45.27851 U2(18) = 0.1000000E+21 U2(19) = 0.1000000E+21 U2(20) = 0.1000000E+21 PTACT2(1) = 3 PTACT2(2) = 20 PTACT2(3) = 16 PTACT2(4) = 5 PTACT2(5) = 2 PTACT2(6) = 14 PTACT2(7) = 1 PTACT2(8) = 6 PTACT2(9) = 10 PTACT2(10) = 5 PTACT2(11) = 4 PTACT2(12) = 20 PTACT2(13) = 12 PTACT2(14) = 9 PTACT2(15) = 7 PTACT2(16) = 19 PTACT2(17) = 11 PTACT2(18) = 16 PTACT2(19) = 3 PTACT2(20) = 18 BETA2(1) = 0.0 BETA2(2) = 0.0 BETA2(3) = 1.409947 BETA2(4) = 1.979443 BETA2(5) = 0.0 BETA2(6) = 0.0 BETA2(7) = 1.459439 BETA2(8) = 1.787866 BETA2(9) = 0.0 BETA2(10) = 0.0 BETA2(11) = 0.0 BETA2(12) = 0.0 BETA2(13) = 0.0 BETA2(14) = 0.0 BETA2(15) = 0.0 BETA2(16) = 0.0 BETA2(17) = 0.0 BETA2(18) = 0.0 BETA2(19) = 0.0 BETA2(20) = 0.0 C Determine discrepancy ACUM = 0.0 DO 20 I = 1,N ACUM = ACUM + ABS(L(I)-L2(I)) ACUM = ACUM + ABS(U(I)-U2(I)) ACUM = ACUM + ABS(PTACT(I)-PTACT2(I)) ACUM = ACUM + ABS(BETA(I)-BETA2(I)) 20 CONTINUE IF (ACUM.GT.1.0) THEN PRINT *,' Discrepancy in generation = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 1 INTERRUPTED **' STOP END IF C Evaluate f(X) and compute discrepancy. CALL BCP10(N,X,F) ACUM = ABS(F-670.363) IF (ACUM.GT.1.0) THEN PRINT *,' Discrepancy in f(X) = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 1 INTERRUPTED **' STOP END IF C Evaluate grad f(X) and compute discrepancy. CALL BCP11(N,X,GF) ACUM = 0.0 DO 30 I = 1,N ACUM = ACUM + GF(I)*GF(I) 30 CONTINUE ACUM = ABS(ACUM-4389455.7521496) IF (ACUM.GT.1.0) THEN PRINT *,' Discrepancy in || grad f(X) ||^2 = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 1 INTERRUPTED **' STOP END IF C Evaluation of the product Hessian of f at X by the vector D. C and compute discrepancy CALL BCP12(N,X,D,HFD) ACUM = 0.0 DO 40 I = 1,N ACUM = ACUM + HFD(I)*HFD(I) 40 CONTINUE ACUM = ABS(ACUM-15192289.986351) IF (ACUM.GT.1.0) THEN PRINT *,' Discrepancy in || H(X)d ||^2 = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 1 INTERRUPTED **' STOP END IF PRINT *,'** VERIFICATION TEST 1 ENDED SUCCESSFULLY **' PRINT * C ******************* Generate a problem with 9 lower bounds and C VERIFICATION TEST 2 3 upper bounds, distributed ramdomly. C ******************* C At the unconstrained solution, 6 bounds should C be active. All of them should have a null C multiplier. C C Moreover, at the unconstrained solution the C distance between each nonactive bound and the C correspondent component should be C uniformly distributed in [100.0,200.0]. C C Use the following function h_i, C h_i(x_i)= ( x_i - xbar_i )^3 + beta_i ( x_i - xbar_i ) C as described in the companion paper, to C define the objective function. PRINT *,'** VERIFICATION TEST 2 **' PRINT * C BND = 30, meaning that 30% of the all possible bounds (40) C are finite. BND = 30 C LOWBND = 75, meaning that 75% of the finite bounds are lower C bounds LOWBND = 75 C PAIR = 0, meaning that finite bounds need not exist in pairs. PAIR = 0 C ACTBND = 50, meaning that 50% of the finite bounds are active C at the unconstrained solution. ACTBND = 50 C Set the two ranges for the multipliers LOWM1 = 0.0 UPPM1 = 0.0 LOWM2 = 0.0 UPPM2 = 0.0 C DEG = 100, meaning that at the unconstrained solution all the C multipliers associated with the active constraints are in the C first range. DEG = 100 C Set the range for the nonactive but finite bounds. WIDTH1 = 100.0 WIDTH2 = 200.0 C Select the function h_i HIF = 2 C Set ISEED ISEED = 1994 C Generate the problem CALL BCP01(N,L,U,BND,LOWBND,ACTBND,PAIR,DEG,LOWM1,UPPM1,LOWM2, + UPPM2,WIDTH1,WIDTH2,HIF,ISEED,INFRM) C Check INFRM IF (INFRM.NE.0) THEN CALL BCP02(INFRM) PRINT *,'** VERIFICATION TEST 2 INTERRUPTED **' STOP END IF C The output from BCP01 should have been the following in C single precision. L2(1) = -0.1000000E+21 L2(2) = -145.7033 L2(3) = 1.000000 L2(4) = -141.4788 L2(5) = 1.000000 L2(6) = -0.1000000E+21 L2(7) = -0.1000000E+21 L2(8) = -0.1000000E+21 L2(9) = -0.1000000E+21 L2(10) = -0.1000000E+21 L2(11) = -140.9426 L2(12) = -0.1000000E+21 L2(13) = -0.1000000E+21 L2(14) = -124.6780 L2(15) = -128.0370 L2(16) = 1.000000 L2(17) = -0.1000000E+21 L2(18) = -0.1000000E+21 L2(19) = -0.1000000E+21 L2(20) = 1.000000 U2(1) = 0.1000000E+21 U2(2) = 0.1000000E+21 U2(3) = 0.1000000E+21 U2(4) = 0.1000000E+21 U2(5) = 0.1000000E+21 U2(6) = 177.5269 U2(7) = 0.1000000E+21 U2(8) = 0.1000000E+21 U2(9) = 0.1000000E+21 U2(10) = 0.1000000E+21 U2(11) = 0.1000000E+21 U2(12) = 0.1000000E+21 U2(13) = 0.1000000E+21 U2(14) = 1.000000 U2(15) = 1.000000 U2(16) = 0.1000000E+21 U2(17) = 0.1000000E+21 U2(18) = 0.1000000E+21 U2(19) = 0.1000000E+21 U2(20) = 0.1000000E+21 PTACT2(1) = 5 PTACT2(2) = 16 PTACT2(3) = 3 PTACT2(4) = 20 PTACT2(5) = 14 PTACT2(6) = 15 PTACT2(7) = 9 PTACT2(8) = 4 PTACT2(9) = 16 PTACT2(10) = 2 PTACT2(11) = 10 PTACT2(12) = 3 PTACT2(13) = 7 PTACT2(14) = 18 PTACT2(15) = 13 PTACT2(16) = 11 PTACT2(17) = 5 PTACT2(18) = 19 PTACT2(19) = 8 PTACT2(20) = 20 BETA2(1) = 0.0 BETA2(2) = 0.0 BETA2(3) = 0.0 BETA2(4) = 0.0 BETA2(5) = 0.0 BETA2(6) = 0.0 BETA2(7) = 1.459439 BETA2(8) = 1.787866 BETA2(9) = 0.0 BETA2(10) = 0.0 BETA2(11) = 0.0 BETA2(12) = 0.0 BETA2(13) = 0.0 BETA2(14) = 0.0 BETA2(15) = 0.0 BETA2(16) = 0.0 BETA2(17) = 0.0 BETA2(18) = 0.0 BETA2(19) = 0.0 BETA2(20) = 0.0 C Determine discrepancy ACUM = 0.0 DO 50 I = 1,N ACUM = ACUM + ABS(L(I)-L2(I)) ACUM = ACUM + ABS(U(I)-U2(I)) ACUM = ACUM + ABS(PTACT(I)-PTACT2(I)) ACUM = ACUM + ABS(BETA(I)-BETA2(I)) 50 CONTINUE IF (ACUM.GT.1.0) THEN PRINT *,' Discrepancy in generation = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 2 INTERRUPTED **' STOP END IF C Evaluate f(X) and compute discrepancy. CALL BCP10(N,X,F) ACUM = ABS(F-669.924) IF (ACUM.GT.1.0) THEN PRINT *,' Discrepancy in f(X) = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 2 INTERRUPTED **' STOP END IF C Evaluate grad f(X) and compute discrepancy. CALL BCP11(N,X,GF) ACUM = 0.0 DO 60 I = 1,N ACUM = ACUM + GF(I)*GF(I) 60 CONTINUE ACUM = ABS(ACUM-4389190.5930229) IF (ACUM.GT.1.0) THEN PRINT *,' Discrepancy in || grad f(X) ||^2 = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 2 INTERRUPTED **' STOP END IF C Evaluation of the product Hessian of f at X by the vector D. C and compute discrepancy CALL BCP12(N,X,D,HFD) ACUM = 0.0 DO 70 I = 1,N ACUM = ACUM + HFD(I)*HFD(I) 70 CONTINUE ACUM = ABS(ACUM-15192921.294180) IF (ACUM.GT.1.0) THEN PRINT *,' Discrepancy in || H(X)d ||^2 = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 2 INTERRUPTED **' STOP END IF PRINT *,'** VERIFICATION TEST 2 ENDED SUCCESSFULLY **' PRINT * C ******************* Generate a problem with 20 lower bounds and C VERIFICATION TEST 3 20 upper bounds. C ******************* C At the unconstrained solution, 20 bounds should C be active. Among these active bounds, 10 should C have a multiplier uniformly distributed in C [0.001,0.1]. The other 10 are uniformly C distributed in [10.0,20.0]. C C Moreover, at the unconstrained solution, the C distance between each nonactive bound and the C correspondent component of such point should be C uniformly distributed in [1.0,2.0]. C C Use the following function h_i, C h_i(x_i)= (x_i-xbar_i)^(7/3)+ beta_i(x_i-xbar_i) C as described in the companion paper. PRINT *,'** VERIFICATION TEST 3 **' PRINT * C BND = 100, meaning that all possible bounds (40) are finite. BND = 100 C LOWBND = 50, meaning that 50% of the finite bounds are lower C bounds. (Note: Any other value would result in an error) LOWBND = 50 C PAIR = 0, meaning that finite bounds need not exist in pairs. C (Note: In this example its value is indifferent) PAIR = 0 C ACTBND = 50, meaning that 50% of the finite bounds are active C at the unconstrained solution. ACTBND = 50 C Set the two ranges for the multipliers LOWM1 = 0.001 UPPM1 = 0.1 LOWM2 = 10.0 UPPM2 = 20.0 C DEG = 50, meaning that at the unconstrained solution 50% of the C multipliers associated with the active constraints are in the C first range. DEG = 50 C Set the range for the nonactive but finite bounds. WIDTH1 = 1.0 WIDTH2 = 2.0 C Select the function h_i HIF = 3 C Set ISEED ISEED = 1994 C Generate the problem CALL BCP01(N,L,U,BND,LOWBND,ACTBND,PAIR,DEG,LOWM1,UPPM1,LOWM2, + UPPM2,WIDTH1,WIDTH2,HIF,ISEED,INFRM) C Check INFRM IF (INFRM.NE.0) THEN CALL BCP02(INFRM) PRINT *,'** VERIFICATION TEST 3 INTERRUPTED **' STOP END IF C The output from BCP01 should have been the following in C single precision. L2(1) = -0.1740453 L2(2) = 1.000000 L2(3) = 1.000000 L2(4) = 1.000000 L2(5) = 1.000000 L2(6) = -0.5027176 L2(7) = -0.4247884 L2(8) = 1.000000 L2(9) = -0.4194258 L2(10) = -0.3779922 L2(11) = 1.000000 L2(12) = -0.2903705 L2(13) = -0.9144108 L2(14) = 1.000000 L2(15) = 1.000000 L2(16) = 1.000000 L2(17) = -0.2567797 L2(18) = -0.4670331 L2(19) = -0.6959825 L2(20) = 1.000000 U2(1) = 1.000000 U2(2) = 2.765269 U2(3) = 2.115662 U2(4) = 2.374369 U2(5) = 2.164240 U2(6) = 1.000000 U2(7) = 1.000000 U2(8) = 2.015358 U2(9) = 1.000000 U2(10) = 1.000000 U2(11) = 2.597139 U2(12) = 1.000000 U2(13) = 1.000000 U2(14) = 2.960955 U2(15) = 2.390587 U2(16) = 2.926635 U2(17) = 1.000000 U2(18) = 1.000000 U2(19) = 1.000000 U2(20) = 2.122221 PTACT2(1) = 20 PTACT2(2) = 14 PTACT2(3) = 15 PTACT2(4) = 11 PTACT2(5) = 8 PTACT2(6) = 4 PTACT2(7) = 3 PTACT2(8) = 5 PTACT2(9) = 2 PTACT2(10) = 16 PTACT2(11) = 9 PTACT2(12) = 7 PTACT2(13) = 10 PTACT2(14) = 1 PTACT2(15) = 19 PTACT2(16) = 6 PTACT2(17) = 17 PTACT2(18) = 12 PTACT2(19) = 18 PTACT2(20) = 13 BETA2(1) = 0.4158475E-01 BETA2(2) = 0.9796487E-01 BETA2(3) = 0.5060987E-01 BETA2(4) = 0.1604595E-01 BETA2(5) = 0.3231400E-01 BETA2(6) = 11.04975 BETA2(7) = 13.14028 BETA2(8) = 18.67464 BETA2(9) = 14.59439 BETA2(10) = 17.87866 BETA2(11) = 0.9906817E-01 BETA2(12) = 0.8076739E-01 BETA2(13) = 0.9251605E-01 BETA2(14) = 0.4730846E-01 BETA2(15) = 0.6836787E-01 BETA2(16) = 18.87231 BETA2(17) = 16.89056 BETA2(18) = 19.72215 BETA2(19) = 10.18641 BETA2(20) = 12.96863 C Determine discrepancy ACUM = 0.0 DO 80 I = 1,N ACUM = ACUM + ABS(L(I)-L2(I)) ACUM = ACUM + ABS(U(I)-U2(I)) ACUM = ACUM + ABS(PTACT(I)-PTACT2(I)) ACUM = ACUM + ABS(BETA(I)-BETA2(I)) 80 CONTINUE IF (ACUM.GT.1.0) THEN PRINT *,' Discrepancy in generation = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 3 INTERRUPTED **' STOP END IF C Evaluate f(X) and compute discrepancy. CALL BCP10(N,X,F) ACUM = ABS(F-689.545) IF (ACUM.GT.1.0) THEN PRINT *,' Discrepancy in f(X) = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 3 INTERRUPTED **' STOP END IF C Evaluate grad f(X) and compute discrepancy. CALL BCP11(N,X,GF) ACUM = 0.0 DO 90 I = 1,N ACUM = ACUM + GF(I)*GF(I) 90 CONTINUE ACUM = ABS(ACUM-4384443.0353111) IF (ACUM.GT.1.0) THEN PRINT *,' Discrepancy in || grad f(X) ||^2 = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 3 INTERRUPTED **' STOP END IF C Evaluation of the product Hessian of f at X by the vector D. C and compute discrepancy CALL BCP12(N,X,D,HFD) ACUM = 0.0 DO 100 I = 1,N ACUM = ACUM + HFD(I)*HFD(I) 100 CONTINUE ACUM = ABS(ACUM-15191153.5) IF (ACUM.GT.1.0) THEN PRINT *,' Discrepancy in || H(X)d ||^2 = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 3 INTERRUPTED **' STOP END IF PRINT *,'** VERIFICATION TEST 3 ENDED SUCCESSFULLY **' PRINT * PRINT *,' END OF VERIFICATION ' END C Extended Rosenbrock function. See More, J., Garbow, B. and C Hillstrom, K. "Testing Unconstrained Optimization Software", C ACM TOMS, Vol. 7, N. 1, March 1981, 17-41. C****************************************************** C EXTENDED ROSENBROCK C****************************************************** C -- UNCONSTRAINED SOLUTION SUBROUTINE SOLUTG(N,XBAR) C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. REAL XBAR(N) C .. C .. Local Scalars .. INTEGER J C .. DO 10 J = 1,N XBAR(J) = 1.0E+0 10 CONTINUE RETURN END C -- OBJECTIVE FUNCTION SUBROUTINE FUNCTG(N,X,G) C .. Scalar Arguments .. REAL G INTEGER N C .. C .. Array Arguments .. REAL X(N) C .. C .. Local Scalars .. REAL T1,T2 INTEGER J C .. G = 0.0 DO 10 J = 1,N,2 T1 = 1.0E+0 - X(J) T2 = 10.0E+0* (X(J+1)-X(J)*X(J)) G = G + T1*T1 + T2*T2 10 CONTINUE RETURN END C -- GRADIENT SUBROUTINE GRADG(N,X,GG) C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. REAL GG(N),X(N) C .. C .. Local Scalars .. REAL T1 INTEGER J C .. DO 10 J = 1,N,2 T1 = X(J+1) - X(J)*X(J) GG(J) = -400.0E+0*X(J)*T1 - 2.0E+0* (1.0E+0-X(J)) GG(J+1) = 200.0E+0*T1 10 CONTINUE RETURN END C -- HESSIAN BY VECTOR PRODUCT SUBROUTINE HESGD(N,X,D,HGD) C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. REAL D(N),HGD(N),X(N) C .. C .. Local Scalars .. REAL T1,T2 INTEGER J C .. DO 10 J = 1,N,2 T1 = -400.0E+0* (X(J+1)-3.0E+0*X(J)*X(J)) + 2.0E+0 T2 = -400.0E+0*X(J) HGD(J) = T1*D(J) + T2*D(J+1) HGD(J+1) = T2*D(J) + 200.0E+0*D(J+1) 10 CONTINUE RETURN END SHAR_EOF fi # end of overwriting check if test -f 'RES' then echo shar: will not over-write existing file "'RES'" else cat << \SHAR_EOF > 'RES' ** VERIFICATION TEST 1 ** ** VERIFICATION TEST 1 ENDED SUCCESSFULLY ** ** VERIFICATION TEST 2 ** ** VERIFICATION TEST 2 ENDED SUCCESSFULLY ** ** VERIFICATION TEST 3 ** ** VERIFICATION TEST 3 ENDED SUCCESSFULLY ** END OF VERIFICATION SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << \SHAR_EOF > 'driver.f' C====================================================================== C c C This file contains the FORTRAN program which demonstrates the c C use of the suite of routines for generating box constrained c C nonlinear optimization problems. c C c C Another purpose of this program is to test the integrity of the c C routines in the suite after installation. It is by no means a c C thorough test. c C c C The other routines that should be linked to this main program c C are in file src.f c C c C The user is referred to the documentation c C c C Soares, J., Judice, J. and Facchinei, F., "Generating box c C constrained optimization problems", c C ACM TOMS, (??) c C and c C Soares, J., Judice, J. and Facchinei, F., "FORTRAN subroutines c C for generating box constrained optimization problems", c C ACM TOMS, (??). c C c C====================================================================== PROGRAM VERIFY C Authors: F. Facchinei, J. Soares and J. Judice C Date Written: September 10, 1995 C Date Modified: September 3, 1996 C C When testing BCP01 (see comments in subroutine BCP01 for details) C we make use of the following variables: C =================== C Input parameters: C ( N, *, *, BND, LOWBND, ACTBND, PAIR, DEG, LOWM1, UPPM1, C LOWM2, UPPM2, WIDTH1, WIDTH2, HIF, ISEED, * ) C Output parameters: C ( *, L, U, *, *, *, *, *, *, *, *, *, *, *, *, *, INFRM ) C When testing BCP10 (see comments in subroutine BCP10 for details) C we make use of the following variables: C ==================== C Input parameters: C ( N, X, *) C Output parameters: C ( *, *, F) C When testing BCP11 (see comments in subroutine BCP11 for details) C we make use of the following variables: C ==================== C Input parameters: C ( N, X, *) C Output parameters: C ( *, *, GF) C When testing BCP12 (see comments in subroutine BCP12 for details) C we make use of the following variables: C ==================== C Input parameters: C ( N, X, D, *) C Output parameters: C ( *, *, *, HFD) C .. See the comments in BCP01 for a description of the variables C in the common blocks. C .. Parameters .. INTEGER MAXN PARAMETER (MAXN=100) INTEGER NN PARAMETER (NN=10000) C .. C .. Scalars in Common .. INTEGER ACTL,ACTU,HI,INL,INU,NM1L,NM1U,NM2L,NM2U C .. C .. Arrays in Common .. DOUBLE PRECISION BETA(NN),XBAR(NN) INTEGER PTACT(NN) C .. C .. Local Scalars .. DOUBLE PRECISION ACUM,F,LOWM1,LOWM2,UPPM1,UPPM2,WIDTH1,WIDTH2 INTEGER ACTBND,BND,DEG,HIF,I,INFRM,ISEED,LOWBND,N,PAIR C .. C .. Local Arrays .. DOUBLE PRECISION BETA2(MAXN),D(MAXN),GF(MAXN),HFD(MAXN),L(MAXN), + L2(MAXN),U(MAXN),U2(MAXN),X(MAXN) INTEGER PTACT2(MAXN) C .. C .. External Subroutines .. EXTERNAL BCP01,BCP02,BCP10,BCP11,BCP12 C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. C .. Common blocks .. COMMON /BCPB1/INL,INU,ACTL,ACTU,NM1L,NM1U,NM2L,NM2U,HI COMMON /BCPB2/BETA,XBAR,PTACT C .. C N = 20 meaning that the problems have 20 variables C (N need to be less or equal to NN) N = 20 C Initialize X and D for the calls to BCP10, BCP11 and BCP12 DO 10 I = 1,N X(I) = 1.0D0 + 1.0D0/I D(I) = 1.0D0/I 10 CONTINUE C ******************* Generate a problem with 10 lower bounds and C VERIFICATION TEST 1 10 upper bounds, distributed randomly. C ******************* C At the unconstrained solution, 8 bounds should C be active. Among these active bounds, 4 should C have a null multiplier. The remaining 4 C multipliers should be uniformly distributed in C [1.0,2.0]. C C Moreover, at the unconstrained solution the C distance between each nonactive bound and the C correspondent component should be C uniformly distributed in [3.0,50.0]. C C Use the simplest function h_i, C h_i(x_i)= beta_i * ( x_i - xbar_i ), C as described in the companion paper, to C define the objective function. PRINT *,'** VERIFICATION TEST 1 **' PRINT * C BND = 50, meaning that 50% of the all possible bounds (40) C are finite. BND = 50 C LOWBND = 50, meaning that lower and upper bounds are C equally distributed among the finite bounds. LOWBND = 50 C PAIR = 0, meaning that finite bounds need not exist in pairs. PAIR = 0 C ACTBND = 40, meaning that 40% of the finite bounds are active C at the unconstrained solution. ACTBND = 40 C Set the two ranges for the multipliers LOWM1 = 0.0D0 UPPM1 = 0.0D0 LOWM2 = 1.0D0 UPPM2 = 2.0D0 C DEG = 50, meaning that at the unconstrained solution 50% of the C multipliers associated with the active constraints are in the C first range while the others are in the second range. DEG = 50 C Set the range for the nonactive but finite bounds. WIDTH1 = 3.0D0 WIDTH2 = 50.0D0 C Select the function h_i HIF = 1 C Set ISEED ISEED = 1994 C Generate the problem CALL BCP01(N,L,U,BND,LOWBND,ACTBND,PAIR,DEG,LOWM1,UPPM1,LOWM2, + UPPM2,WIDTH1,WIDTH2,HIF,ISEED,INFRM) C Check INFRM IF (INFRM.NE.0) THEN CALL BCP02(INFRM) PRINT *,'** VERIFICATION TEST 1 INTERRUPTED **' STOP END IF C The output from BCP01 should have been the following in C single precision. L2(1) = -0.1000000D+21 L2(2) = -23.95055D0 L2(3) = 1.000000D0 L2(4) = -21.96506D0 L2(5) = 1.000000D0 L2(6) = -0.1000000D+21 L2(7) = -0.1000000D+21 L2(8) = -34.71118D0 L2(9) = -0.1000000D+21 L2(10) = -0.1000000D+21 L2(11) = -21.71301D0 L2(12) = -0.1000000D+21 L2(13) = -0.1000000D+21 L2(14) = -14.06864D0 L2(15) = -15.64741D0 L2(16) = 1.000000D0 L2(17) = -0.1000000D+21 L2(18) = -0.1000000D+21 L2(19) = -0.1000000D+21 L2(20) = 1.000000D0 U2(1) = 1.000000D0 U2(2) = 1.000000D0 U2(3) = 0.1000000D+21 U2(4) = 0.1000000D+21 U2(5) = 11.12918D0 U2(6) = 1.000000D0 U2(7) = 0.1000000D+21 U2(8) = 4.876650D0 U2(9) = 0.1000000D+21 U2(10) = 26.85025D0 U2(11) = 0.1000000D+21 U2(12) = 0.1000000D+21 U2(13) = 39.96765D0 U2(14) = 1.000000D0 U2(15) = 49.16487D0 U2(16) = 0.1000000D+21 U2(17) = 45.27851D0 U2(18) = 0.1000000D+21 U2(19) = 0.1000000D+21 U2(20) = 0.1000000D+21 PTACT2(1) = 3 PTACT2(2) = 20 PTACT2(3) = 16 PTACT2(4) = 5 PTACT2(5) = 2 PTACT2(6) = 14 PTACT2(7) = 1 PTACT2(8) = 6 PTACT2(9) = 10 PTACT2(10) = 5 PTACT2(11) = 4 PTACT2(12) = 20 PTACT2(13) = 12 PTACT2(14) = 9 PTACT2(15) = 7 PTACT2(16) = 19 PTACT2(17) = 11 PTACT2(18) = 16 PTACT2(19) = 3 PTACT2(20) = 18 BETA2(1) = 0.0D0 BETA2(2) = 0.0D0 BETA2(3) = 1.409947D0 BETA2(4) = 1.979443D0 BETA2(5) = 0.0D0 BETA2(6) = 0.0D0 BETA2(7) = 1.459439D0 BETA2(8) = 1.787866D0 BETA2(9) = 0.0D0 BETA2(10) = 0.0D0 BETA2(11) = 0.0D0 BETA2(12) = 0.0D0 BETA2(13) = 0.0D0 BETA2(14) = 0.0D0 BETA2(15) = 0.0D0 BETA2(16) = 0.0D0 BETA2(17) = 0.0D0 BETA2(18) = 0.0D0 BETA2(19) = 0.0D0 BETA2(20) = 0.0D0 C Determine discrepancy ACUM = 0.0D0 DO 20 I = 1,N ACUM = ACUM + ABS(L(I)-L2(I)) ACUM = ACUM + ABS(U(I)-U2(I)) ACUM = ACUM + ABS(PTACT(I)-PTACT2(I)) ACUM = ACUM + ABS(BETA(I)-BETA2(I)) 20 CONTINUE IF (ACUM.GT.1.0D0) THEN PRINT *,' Discrepancy in generation = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 1 INTERRUPTED **' STOP END IF C Evaluate f(X) and compute discrepancy. CALL BCP10(N,X,F) ACUM = ABS(F-670.363D0) IF (ACUM.GT.1.0D0) THEN PRINT *,' Discrepancy in f(X) = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 1 INTERRUPTED **' STOP END IF C Evaluate grad f(X) and compute discrepancy. CALL BCP11(N,X,GF) ACUM = 0.0D0 DO 30 I = 1,N ACUM = ACUM + GF(I)*GF(I) 30 CONTINUE ACUM = ABS(ACUM-4389455.7521496D0) IF (ACUM.GT.1.0D0) THEN PRINT *,' Discrepancy in || grad f(X) ||^2 = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 1 INTERRUPTED **' STOP END IF C Evaluation of the product Hessian of f at X by the vector D. C and compute discrepancy CALL BCP12(N,X,D,HFD) ACUM = 0.0D0 DO 40 I = 1,N ACUM = ACUM + HFD(I)*HFD(I) 40 CONTINUE ACUM = ABS(ACUM-15192289.986351D0) IF (ACUM.GT.1.0D0) THEN PRINT *,' Discrepancy in || H(X)d ||^2 = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 1 INTERRUPTED **' STOP END IF PRINT *,'** VERIFICATION TEST 1 ENDED SUCCESSFULLY **' PRINT * C ******************* Generate a problem with 9 lower bounds and C VERIFICATION TEST 2 3 upper bounds, distributed ramdomly. C ******************* C At the unconstrained solution, 6 bounds should C be active. All of them should have a null C multiplier. C C Moreover, at the unconstrained solution the C distance between each nonactive bound and the C correspondent component should be C uniformly distributed in [100.0,200.0]. C C Use the following function h_i, C h_i(x_i)= ( x_i - xbar_i )^3 + beta_i ( x_i - xbar_i ) C as described in the companion paper, to C define the objective function. PRINT *,'** VERIFICATION TEST 2 **' PRINT * C BND = 30, meaning that 30% of the all possible bounds (40) C are finite. BND = 30 C LOWBND = 75, meaning that 75% of the finite bounds are lower C bounds LOWBND = 75 C PAIR = 0, meaning that finite bounds need not exist in pairs. PAIR = 0 C ACTBND = 50, meaning that 50% of the finite bounds are active C at the unconstrained solution. ACTBND = 50 C Set the two ranges for the multipliers LOWM1 = 0.0D0 UPPM1 = 0.0D0 LOWM2 = 0.0D0 UPPM2 = 0.0D0 C DEG = 100, meaning that at the unconstrained solution all the C multipliers associated with the active constraints are in the C first range. DEG = 100 C Set the range for the nonactive but finite bounds. WIDTH1 = 100.0D0 WIDTH2 = 200.0D0 C Select the function h_i HIF = 2 C Set ISEED ISEED = 1994 C Generate the problem CALL BCP01(N,L,U,BND,LOWBND,ACTBND,PAIR,DEG,LOWM1,UPPM1,LOWM2, + UPPM2,WIDTH1,WIDTH2,HIF,ISEED,INFRM) C Check INFRM IF (INFRM.NE.0) THEN CALL BCP02(INFRM) PRINT *,'** VERIFICATION TEST 2 INTERRUPTED **' STOP END IF C The output from BCP01 should have been the following in C single precision. L2(1) = -0.1000000D+21 L2(2) = -145.7033D0 L2(3) = 1.000000D0 L2(4) = -141.4788D0 L2(5) = 1.000000D0 L2(6) = -0.1000000D+21 L2(7) = -0.1000000D+21 L2(8) = -0.1000000D+21 L2(9) = -0.1000000D+21 L2(10) = -0.1000000D+21 L2(11) = -140.9426D0 L2(12) = -0.1000000D+21 L2(13) = -0.1000000D+21 L2(14) = -124.6780D0 L2(15) = -128.0370D0 L2(16) = 1.000000D0 L2(17) = -0.1000000D+21 L2(18) = -0.1000000D+21 L2(19) = -0.1000000D+21 L2(20) = 1.000000D0 U2(1) = 0.1000000D+21 U2(2) = 0.1000000D+21 U2(3) = 0.1000000D+21 U2(4) = 0.1000000D+21 U2(5) = 0.1000000D+21 U2(6) = 177.5269D0 U2(7) = 0.1000000D+21 U2(8) = 0.1000000D+21 U2(9) = 0.1000000D+21 U2(10) = 0.1000000D+21 U2(11) = 0.1000000D+21 U2(12) = 0.1000000D+21 U2(13) = 0.1000000D+21 U2(14) = 1.000000D0 U2(15) = 1.000000D0 U2(16) = 0.1000000D+21 U2(17) = 0.1000000D+21 U2(18) = 0.1000000D+21 U2(19) = 0.1000000D+21 U2(20) = 0.1000000D+21 PTACT2(1) = 5 PTACT2(2) = 16 PTACT2(3) = 3 PTACT2(4) = 20 PTACT2(5) = 14 PTACT2(6) = 15 PTACT2(7) = 9 PTACT2(8) = 4 PTACT2(9) = 16 PTACT2(10) = 2 PTACT2(11) = 10 PTACT2(12) = 3 PTACT2(13) = 7 PTACT2(14) = 18 PTACT2(15) = 13 PTACT2(16) = 11 PTACT2(17) = 5 PTACT2(18) = 19 PTACT2(19) = 8 PTACT2(20) = 20 BETA2(1) = 0.0D0 BETA2(2) = 0.0D0 BETA2(3) = 0.0D0 BETA2(4) = 0.0D0 BETA2(5) = 0.0D0 BETA2(6) = 0.0D0 BETA2(7) = 1.459439D0 BETA2(8) = 1.787866D0 BETA2(9) = 0.0D0 BETA2(10) = 0.0D0 BETA2(11) = 0.0D0 BETA2(12) = 0.0D0 BETA2(13) = 0.0D0 BETA2(14) = 0.0D0 BETA2(15) = 0.0D0 BETA2(16) = 0.0D0 BETA2(17) = 0.0D0 BETA2(18) = 0.0D0 BETA2(19) = 0.0D0 BETA2(20) = 0.0D0 C Determine discrepancy ACUM = 0.0D0 DO 50 I = 1,N ACUM = ACUM + ABS(L(I)-L2(I)) ACUM = ACUM + ABS(U(I)-U2(I)) ACUM = ACUM + ABS(PTACT(I)-PTACT2(I)) ACUM = ACUM + ABS(BETA(I)-BETA2(I)) 50 CONTINUE IF (ACUM.GT.1.0D0) THEN PRINT *,' Discrepancy in generation = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 2 INTERRUPTED **' STOP END IF C Evaluate f(X) and compute discrepancy. CALL BCP10(N,X,F) ACUM = ABS(F-669.924D0) IF (ACUM.GT.1.0D0) THEN PRINT *,' Discrepancy in f(X) = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 2 INTERRUPTED **' STOP END IF C Evaluate grad f(X) and compute discrepancy. CALL BCP11(N,X,GF) ACUM = 0.0D0 DO 60 I = 1,N ACUM = ACUM + GF(I)*GF(I) 60 CONTINUE ACUM = ABS(ACUM-4389190.5930229D0) IF (ACUM.GT.1.0D0) THEN PRINT *,' Discrepancy in || grad f(X) ||^2 = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 2 INTERRUPTED **' STOP END IF C Evaluation of the product Hessian of f at X by the vector D. C and compute discrepancy CALL BCP12(N,X,D,HFD) ACUM = 0.0D0 DO 70 I = 1,N ACUM = ACUM + HFD(I)*HFD(I) 70 CONTINUE ACUM = ABS(ACUM-15192921.294180D0) IF (ACUM.GT.1.0D0) THEN PRINT *,' Discrepancy in || H(X)d ||^2 = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 2 INTERRUPTED **' STOP END IF PRINT *,'** VERIFICATION TEST 2 ENDED SUCCESSFULLY **' PRINT * C ******************* Generate a problem with 20 lower bounds and C VERIFICATION TEST 3 20 upper bounds. C ******************* C At the unconstrained solution, 20 bounds should C be active. Among these active bounds, 10 should C have a multiplier uniformly distributed in C [0.001,0.1]. The other 10 are uniformly C distributed in [10.0,20.0]. C C Moreover, at the unconstrained solution, the C distance between each nonactive bound and the C correspondent component of such point should be C uniformly distributed in [1.0,2.0]. C C Use the following function h_i, C h_i(x_i)= (x_i-xbar_i)^(7/3)+ beta_i(x_i-xbar_i) C as described in the companion paper. PRINT *,'** VERIFICATION TEST 3 **' PRINT * C BND = 100, meaning that all possible bounds (40) are finite. BND = 100 C LOWBND = 50, meaning that 50% of the finite bounds are lower C bounds. (Note: Any other value would result in an error) LOWBND = 50 C PAIR = 0, meaning that finite bounds need not exist in pairs. C (Note: In this example its value is indifferent) PAIR = 0 C ACTBND = 50, meaning that 50% of the finite bounds are active C at the unconstrained solution. ACTBND = 50 C Set the two ranges for the multipliers LOWM1 = 0.001D0 UPPM1 = 0.1D0 LOWM2 = 10.0D0 UPPM2 = 20.0D0 C DEG = 50, meaning that at the unconstrained solution 50% of the C multipliers associated with the active constraints are in the C first range. DEG = 50 C Set the range for the nonactive but finite bounds. WIDTH1 = 1.0D0 WIDTH2 = 2.0D0 C Select the function h_i HIF = 3 C Set ISEED ISEED = 1994 C Generate the problem CALL BCP01(N,L,U,BND,LOWBND,ACTBND,PAIR,DEG,LOWM1,UPPM1,LOWM2, + UPPM2,WIDTH1,WIDTH2,HIF,ISEED,INFRM) C Check INFRM IF (INFRM.NE.0) THEN CALL BCP02(INFRM) PRINT *,'** VERIFICATION TEST 3 INTERRUPTED **' STOP END IF C The output from BCP01 should have been the following in C single precision. L2(1) = -0.1740453D0 L2(2) = 1.000000D0 L2(3) = 1.000000D0 L2(4) = 1.000000D0 L2(5) = 1.000000D0 L2(6) = -0.5027176D0 L2(7) = -0.4247884D0 L2(8) = 1.000000D0 L2(9) = -0.4194258D0 L2(10) = -0.3779922D0 L2(11) = 1.000000D0 L2(12) = -0.2903705D0 L2(13) = -0.9144108D0 L2(14) = 1.000000D0 L2(15) = 1.000000D0 L2(16) = 1.000000D0 L2(17) = -0.2567797D0 L2(18) = -0.4670331D0 L2(19) = -0.6959825D0 L2(20) = 1.000000D0 U2(1) = 1.000000D0 U2(2) = 2.765269D0 U2(3) = 2.115662D0 U2(4) = 2.374369D0 U2(5) = 2.164240D0 U2(6) = 1.000000D0 U2(7) = 1.000000D0 U2(8) = 2.015358D0 U2(9) = 1.000000D0 U2(10) = 1.000000D0 U2(11) = 2.597139D0 U2(12) = 1.000000D0 U2(13) = 1.000000D0 U2(14) = 2.960955D0 U2(15) = 2.390587D0 U2(16) = 2.926635D0 U2(17) = 1.000000D0 U2(18) = 1.000000D0 U2(19) = 1.000000D0 U2(20) = 2.122221D0 PTACT2(1) = 20 PTACT2(2) = 14 PTACT2(3) = 15 PTACT2(4) = 11 PTACT2(5) = 8 PTACT2(6) = 4 PTACT2(7) = 3 PTACT2(8) = 5 PTACT2(9) = 2 PTACT2(10) = 16 PTACT2(11) = 9 PTACT2(12) = 7 PTACT2(13) = 10 PTACT2(14) = 1 PTACT2(15) = 19 PTACT2(16) = 6 PTACT2(17) = 17 PTACT2(18) = 12 PTACT2(19) = 18 PTACT2(20) = 13 BETA2(1) = 0.4158475D-01 BETA2(2) = 0.9796487D-01 BETA2(3) = 0.5060987D-01 BETA2(4) = 0.1604595D-01 BETA2(5) = 0.3231400D-01 BETA2(6) = 11.04975D0 BETA2(7) = 13.14028D0 BETA2(8) = 18.67464D0 BETA2(9) = 14.59439D0 BETA2(10) = 17.87866D0 BETA2(11) = 0.9906817D-01 BETA2(12) = 0.8076739D-01 BETA2(13) = 0.9251605D-01 BETA2(14) = 0.4730846D-01 BETA2(15) = 0.6836787D-01 BETA2(16) = 18.87231D0 BETA2(17) = 16.89056D0 BETA2(18) = 19.72215D0 BETA2(19) = 10.18641D0 BETA2(20) = 12.96863D0 C Determine discrepancy ACUM = 0.0D0 DO 80 I = 1,N ACUM = ACUM + ABS(L(I)-L2(I)) ACUM = ACUM + ABS(U(I)-U2(I)) ACUM = ACUM + ABS(PTACT(I)-PTACT2(I)) ACUM = ACUM + ABS(BETA(I)-BETA2(I)) 80 CONTINUE IF (ACUM.GT.1.0D0) THEN PRINT *,' Discrepancy in generation = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 3 INTERRUPTED **' STOP END IF C Evaluate f(X) and compute discrepancy. CALL BCP10(N,X,F) ACUM = ABS(F-689.545D0) IF (ACUM.GT.1.0D0) THEN PRINT *,' Discrepancy in f(X) = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 3 INTERRUPTED **' STOP END IF C Evaluate grad f(X) and compute discrepancy. CALL BCP11(N,X,GF) ACUM = 0.0D0 DO 90 I = 1,N ACUM = ACUM + GF(I)*GF(I) 90 CONTINUE ACUM = ABS(ACUM-4384443.0353111D0) IF (ACUM.GT.1.0D0) THEN PRINT *,' Discrepancy in || grad f(X) ||^2 = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 3 INTERRUPTED **' STOP END IF C Evaluation of the product Hessian of f at X by the vector D. C and compute discrepancy CALL BCP12(N,X,D,HFD) ACUM = 0.0D0 DO 100 I = 1,N ACUM = ACUM + HFD(I)*HFD(I) 100 CONTINUE ACUM = ABS(ACUM-15191153.5D0) IF (ACUM.GT.1.0D0) THEN PRINT *,' Discrepancy in || H(X)d ||^2 = ',ACUM PRINT * PRINT *,'** VERIFICATION TEST 3 INTERRUPTED **' STOP END IF PRINT *,'** VERIFICATION TEST 3 ENDED SUCCESSFULLY **' PRINT * PRINT *,' END OF VERIFICATION ' END C Extended Rosenbrock function. See More, J., Garbow, B. and C Hillstrom, K. "Testing Unconstrained Optimization Software", C ACM TOMS, Vol. 7, N. 1, March 1981, 17-41. C****************************************************** C EXTENDED ROSENBROCK C****************************************************** C -- UNCONSTRAINED SOLUTION SUBROUTINE SOLUTG(N,XBAR) C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. DOUBLE PRECISION XBAR(N) C .. C .. Local Scalars .. INTEGER J C .. DO 10 J = 1,N XBAR(J) = 1.0D+0 10 CONTINUE RETURN END C -- OBJECTIVE FUNCTION SUBROUTINE FUNCTG(N,X,G) C .. Scalar Arguments .. DOUBLE PRECISION G INTEGER N C .. C .. Array Arguments .. DOUBLE PRECISION X(N) C .. C .. Local Scalars .. DOUBLE PRECISION T1,T2 INTEGER J C .. G = 0.0D0 DO 10 J = 1,N,2 T1 = 1.0D+0 - X(J) T2 = 10.0D+0* (X(J+1)-X(J)*X(J)) G = G + T1*T1 + T2*T2 10 CONTINUE RETURN END C -- GRADIENT SUBROUTINE GRADG(N,X,GG) C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. DOUBLE PRECISION GG(N),X(N) C .. C .. Local Scalars .. DOUBLE PRECISION T1 INTEGER J C .. DO 10 J = 1,N,2 T1 = X(J+1) - X(J)*X(J) GG(J) = -400.0D+0*X(J)*T1 - 2.0D+0* (1.0D+0-X(J)) GG(J+1) = 200.0D+0*T1 10 CONTINUE RETURN END C -- HESSIAN BY VECTOR PRODUCT SUBROUTINE HESGD(N,X,D,HGD) C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. DOUBLE PRECISION D(N),HGD(N),X(N) C .. C .. Local Scalars .. DOUBLE PRECISION T1,T2 INTEGER J C .. DO 10 J = 1,N,2 T1 = -400.0D+0* (X(J+1)-3.0D+0*X(J)*X(J)) + 2.0D+0 T2 = -400.0D+0*X(J) HGD(J) = T1*D(J) + T2*D(J+1) HGD(J+1) = T2*D(J) + 200.0D+0*D(J+1) 10 CONTINUE RETURN END SHAR_EOF fi # end of overwriting check if test -f 'RES' then echo shar: will not over-write existing file "'RES'" else cat << \SHAR_EOF > 'RES' ** VERIFICATION TEST 1 ** ** VERIFICATION TEST 1 ENDED SUCCESSFULLY ** ** VERIFICATION TEST 2 ** ** VERIFICATION TEST 2 ENDED SUCCESSFULLY ** ** VERIFICATION TEST 3 ** ** VERIFICATION TEST 3 ENDED SUCCESSFULLY ** END OF VERIFICATION SHAR_EOF fi # end of overwriting check cd .. cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << \SHAR_EOF > 'src.f' C====================================================================== C c C This file contains the FORTRAN routines that the user should c C link to his/her code in order to employ a box constrained c C generated problem. c C Other routines that should be linked, besides the main program, c C include subroutines to evaluate the unconstrained function c C as well as its derivatives. c C c C The user is referred to the documentation c C c C Soares, J., Judice, J. and Facchinei, F., "Generating box c C constrained optimization problems", c C ACM TOMS, (??) c C and c C Soares, J., Judice, J. and Facchinei, F., "FORTRAN subroutines c C for generating box constrained optimization problems", c C ACM TOMS, (??). c C c C To change the precision of this set of subroutines the user c C should do the following: c C c C 1. From double to real: c C 1.1. Replace all 'DOUBLE PRECISION' with 'REAL' c C 1.2. >> 'D+' >> 'E+' c C c C 2. From real to double: c C 1.1. Replace all 'REAL' with 'DOUBLE PRECISION' c C 1.2. >> 'E+' >> 'D+' c C c C IMPORTANT EXCEPTION: Do not replace lower case 'real'. In case c C of doubt make sure that RAND in its definition and calls is c C always single precision. c C c C====================================================================== C .. C .. Begin of subroutine BCP01 .. SUBROUTINE BCP01(N,L,U,BND,LOWBND,ACTBND,PAIR,DEG,LOWM1,UPPM1, + LOWM2,UPPM2,WIDTH1,WIDTH2,HIF,ISEED,INFRM) C _Authors: F. Facchinei, J. Soares and J. Judice C _Date Written: September 10, 1995 C _Date Modified: September 3, 1996 C _Purpose: C C BCP01 generates box constrained nonlinear optimization C problems of the form C C Minimize f(x) C s.t. l <= x <= u C C where x belongs to R^n and f is a continuously differentiable C function R^n -> R. C C More precisely, BCP01 sets up the bounds, and the data required C to evaluate the objective function f and its derivatives. C Unsuccessful generation is detected by having a nonzero value C in INFRM as output. C C _Arguments in parameter list: C C N (input) INTEGER C The number of variables in the problem. (Must be <= NN, C see common block below) C C L, U (output) DOUBLE PRECISION vectors of dimension N C The bounds of the new problem. C C BND (input) INTEGER C The number of bound constraints (or simply, bounds), C expressed as the percentage of the 2*n possible constraints, C constraints,that are really existing (Must be in [0,100]). C Example: if N is 500 and BND=40 then the problem has 400 C bounds. C C LOWBND (input) INTEGER C The number of lower bounds, expressed as the percentage of C the number of constraints as determined by BND (Must be in C [0,100]). This quantity also determines the number of upper C bounds. C Example: if N is 500, BND is 40 and LOWBND=75 then there C are 300 (100) lower (upper) bounds. C C ACTBND (input) INTEGER C The number of active bounds at the unconstrained solution, C expressed as the percentage of the number of constraints as C determined by BND (Must be in [0,100]). These active C bounds are distributed proportionally to LOWBND. C Example: Continuing the previous example, if ACTBND C is 50 then 150 (50) of the lower (upper) bounds should C be active at the unconstrained solution. C C PAIR (input) INTEGER C If this parameter is set to 1 then the constraints as C determined by BND must come in pairs, i.e. each variable C has either no constraint or both lower and upper bounds (in C which case LOWBND must be set to 50). If PAIR is C set to 0 then the indices of variables having upper bounds C are chosen randomly and independently of the lower bounds C (Must be 0 or 1). C C DEG (input) INTEGER C The number of Lagrange multipliers (at the unconstrained C solution) uniformly distributed in the interval C [LOWM1, UPPM1] (see below), expressed as the percentage of C the active constraints (Must be in [0,100]). C The remaining active constraints have Lagrange multipliers C uniformly distributed in [LOWM2, UPPM2] (see below). C Example: C Continuing the example, if DEG=10 then 15 (5) active C lower (upper) bounds (at the unconstrained solution) have C Lagrange multipliers uniformly distributed in [LOWM1, C UPPM1], while the remaining active bounds have Lagrange C multipliers uniformly distributed in [LOWM2, UPPM2]. C C LOWM1, UPPM1, LOWM2, UPPM2 (input) DOUBLE PRECISION C These parameters specify the ranges for the C Lagrange multipliers according to what explained above. C (Must be >= 0) C C WIDTH1, WIDTH2 (input) DOUBLE PRECISION C These two parameters control the value of the bounds that C are not active at the unconstrained solution. More C precisely, if the i-th component has a lower bound not C active at the unconstrained solution xbar then this bound C is randomly chosen in the interval C [xbar_i - WIDTH2, xbar_i - WIDTH1]. Analogously, C if the i-th component has an upper bound not active at C the unconstrained solution xbar then this bound is C randomly chosen in the interval C [xbar_i + WIDTH1, xbar_i + WIDTH2] (Must be >= 0). C C HIF (input) INTEGER C This parameter determines which function h_i is to be used C by the generator. If HIF is equal to C C 1: h_i(x_i)= beta_i * ( x_i - xbar_i ); C C 2: h_i(x_i)=( x_i - xbar_i )^3 + beta_i*( x_i - xbar_i ); C C 3: h_i(x_i)=(x_i - xbar_i)^(7/3) + beta_i*(x_i - xbar_i). C C (Must be 1, 2 or 3) C C ISEED (input) INTEGER C A seed for the random number generator. By running the box C constrained problems generator with the same choice of C parameters but two different seeds, two different box C constrained problems are generated with the same overall C characteristics (Must be >= 0). C C INFRM (output) INTEGER C Determines if a successful call to BCP01 has occurred C (INFRM=0) or not (INFRM<>0). C If INFRM is equal to C 1: N is less than 1 or greater than NN C 2, 3, 4, 5, 6, 11: Some integer parameter is defined C out of its range. C 7: PAIR=1 but LOWBND is not 50. C 8: Either LOWM1 or LOWM2 is less than zero. C 9: Either LOWM1 is greater than UPPM1 or LOWM2 is greater C than UPPM2. C 10: Either WIDTH1 is less than zero or WIDTH1 is greater C than WIDTH2. C 12: It is impossible to make a selection because LOWBND C is too large. C 13: It is impossible to make a selection because LOWBND C is too small. C 14: It is impossible to make a selection because ACTBND C is too large. C C _Explanation of the variables used in the common blocks: C C NN INTEGER C Maximum N allowed. C C INL INTEGER C Number of lower bounds. C C INU INTEGER C Number of upper bounds. C C ACTL INTEGER C Number of lower bounds active at the unconstrained solution C C ACTU INTEGER C Number of upper bounds active at the unconstrained solution C C NM1L INTEGER C Number of lower bounds active at the unconstrained solution C which have a Lagrange multiplier randomly chosen in C [LOWM1,UPPM1]. C C NM1U INTEGER C Number of upper bounds active at the unconstrained solution C which have a Lagrange multiplier randomly chosen in C [LOWM1,UPPM1]. C C NM2L INTEGER C Number of lower bounds active at the unconstrained solution C which have a Lagrange multiplier randomly chosen in C [LOWM2,UPPM2]. C C NM2U INTEGER C Number of upper bounds active at the unconstrained solution C which have a Lagrange multiplier randomly chosen in C [LOWM2,UPPM2]. C C HI INTEGER C Characterizes the function h_i used. HI=HIF. C C PTACT INTEGER vector of dimension NN C A vector of pointers that characterizes the bounds that are C active at the unconstrained solution, i.e., C PTACT(II) (for II=1,ACTL) pinpoints an active lower bound C PTACT(II) (for II=ACTL+1,ACTL+ACTU) pinpoints an active C upper bound. C C BETA DOUBLE PRECISION vector of dimension NN C At the unconstrained solution, BETA characterizes the C multipliers of the bounds that are active, i.e., C BETA(II) (for II=1,ACTL) is the multiplier associated C with the lower bound in the PTBETA(II) variable; C BETA(II) (for II=ACTL+1,ACTL+ACTU) is the multiplier C associated with the upper bound in the PTBETA(II) variable. C C XBAR DOUBLE PRECISION vector of dimension NN C The unconstrained solution. C .. Parameters .. INTEGER NN PARAMETER (NN=10000) REAL BIGINF PARAMETER (BIGINF=1.0E+20) C .. C .. Scalar Arguments .. REAL LOWM1,LOWM2,UPPM1,UPPM2,WIDTH1,WIDTH2 INTEGER ACTBND,BND,DEG,HIF,INFRM,ISEED,LOWBND,N,PAIR C .. C .. Array Arguments .. REAL L(N),U(N) C .. C .. Scalars in Common .. INTEGER ACTL,ACTU,HI,INL,INU,NM1L,NM1U,NM2L,NM2U C .. C .. Arrays in Common .. REAL BETA(NN),XBAR(NN) INTEGER PTACT(NN) C .. C .. Local Scalars .. REAL AUX INTEGER I,II,INBND,J,J1,J2,JJ C .. C .. External Functions .. REAL RAND EXTERNAL RAND C .. C .. External Subroutines .. EXTERNAL BCP04,SOLUTG C .. C .. Intrinsic Functions .. INTRINSIC MIN C .. C .. Common blocks .. COMMON /BCPB1/INL,INU,ACTL,ACTU,NM1L,NM1U,NM2L,NM2U,HI COMMON /BCPB2/BETA,XBAR,PTACT C .. C Test the input parameters. INFRM = 0 IF ((N.GT.NN) .OR. (N.LE.0)) INFRM = 1 IF ((BND.LT.0) .OR. (BND.GT.100)) INFRM = 2 IF ((LOWBND.LT.0) .OR. (LOWBND.GT.100)) INFRM = 3 IF ((ACTBND.LT.0) .OR. (ACTBND.GT.100)) INFRM = 4 IF ((DEG.LT.0) .OR. (DEG.GT.100)) INFRM = 5 IF ((PAIR.NE.0) .AND. (PAIR.NE.1)) INFRM = 6 IF ((PAIR.EQ.1) .AND. (LOWBND.NE.50)) INFRM = 7 IF ((LOWM1.LT.0) .OR. (LOWM2.LT.0)) INFRM = 8 IF ((LOWM1.GT.UPPM1) .OR. (LOWM2.GT.UPPM2)) INFRM = 9 IF ((WIDTH1.LE.0) .OR. (WIDTH1.GT.WIDTH2)) INFRM = 10 IF ((HIF.LE.0) .OR. (HIF.GE.4)) INFRM = 11 IF (INFRM.NE.0) RETURN C Determine the number of constraints INBND = (BND*2*N)/100 C Determine the number of lower bounds INL = (LOWBND*INBND)/100 C Safeguard the maximum number of lower bounds IF (INL.GT.N) THEN INFRM = 12 RETURN END IF C Determine the number of upper bounds INU = INBND - INL C Safeguard the maximum number of upper bounds IF (INU.GT.N) THEN INFRM = 13 RETURN END IF C Determine the number of active constraints at XBAR I = (ACTBND*INBND)/100 ACTL = (ACTBND*INL)/100 ACTU = I - ACTL C Safeguard some settings related to the parameter PAIR .. IF (PAIR.EQ.0) THEN C .. Since the bounds are supposed to be chosen randomly, C safeguard the total number of constraints C active at XBAR. IF (ACTL+ACTU.GT.N) THEN INFRM = 14 RETURN END IF ELSE C .. Since the bounds are supposed to come in pairs, C safeguard that the number of lower bounds must C equal the number of upper bounds .. J = MIN(INL,INU) INL = J INU = J C .. in which case the number of bounds and active bounds at XBAR C should be corrected .. INBND = INL + INU ACTL = (ACTBND*INL)/100 ACTU = ACTL C .. there should not exist one pair of bounds that are both C active at XBAR (Note: INL=INU). IF (ACTL+ACTU.GT.INL) THEN INFRM = 14 RETURN END IF END IF C Read the optimal unconstrained point CALL SOLUTG(N,XBAR) C Set default values to nonexistent bounds and make some random C settings. For the moment, PTACT is used as an auxiliary vector. C After a call to BCP04 it will contain the set of integer values C 1..N in a random order. DO 10 I = 1,N L(I) = -BIGINF U(I) = BIGINF PTACT(I) = I 10 CONTINUE CALL BCP04(N,PTACT(1),ISEED) C Defining the lower bounds .. C .. actives at XBAR (Indices are in PTACT(1)..PTACT(ACTL)) .. J = 1 DO 20 II = 1,ACTL I = PTACT(J) L(I) = XBAR(I) J = J + 1 20 CONTINUE C .. and non actives at XBAR C (Indices are in PTACT(ACTL+1)..PTACT(INL)) .. DO 30 II = ACTL + 1,INL I = PTACT(J) AUX = WIDTH1 + RAND(ISEED)* (WIDTH2-WIDTH1) L(I) = XBAR(I) - AUX J = J + 1 30 CONTINUE C .. (Note that PTACT(INL+1)..PTACT(N) contain indices of variables C that do not have lower bounds). C Defining the upper bounds .. IF (PAIR.EQ.0) THEN C .. since the bounds are supposed to be chosen randomly, C disorder PTACT(ACTL+1)..PTACT(N) .. J = N - ACTL CALL BCP04(J,PTACT(ACTL+1),ISEED) C .. and determine the upper bounds that are active at XBAR C (Indices are in PTACT(ACTL+1)..PTACT(ACTL+ACTU) ) .. J = 1 JJ = ACTL + ACTU DO 40 II = 1,ACTU I = PTACT(JJ) U(I) = XBAR(I) PTACT(JJ) = PTACT(J) PTACT(J) = I J = J + 1 JJ = JJ - 1 40 CONTINUE C .. (PTACT(1)..PTACT(ACTU) contain the indices of upper bounds C that are active at XBAR) C Now, disorder PTACT(ACTU+1)..PTACT(N) .. J = N - ACTU CALL BCP04(J,PTACT(ACTU+1),ISEED) C C .. to define the upper bounds that are not active at XBAR C (Indices are in PTACT(ACTU+1)..PTACT(INU)) .. J = ACTU + 1 DO 50 II = ACTU + 1,INU I = PTACT(J) AUX = WIDTH1 + RAND(ISEED)* (WIDTH2-WIDTH1) U(I) = XBAR(I) + AUX J = J + 1 50 CONTINUE ELSE C .. since the bounds are supposed to come in pairs, C disorder PTACT(ACTL+1)..PTACT(INL) .. J = INL - ACTL CALL BCP04(J,PTACT(ACTL+1),ISEED) C .. and determine the upper bounds that are actives at XBAR C (Indices are in PTACT(ACTL+1)..PTACT(ACTL+ACTU). Also note that C ACTL+ACTU.LE.INL=INU) .. J = 1 JJ = ACTL + ACTU DO 60 II = 1,ACTU I = PTACT(JJ) U(I) = XBAR(I) PTACT(JJ) = PTACT(J) PTACT(J) = I J = J + 1 JJ = JJ - 1 60 CONTINUE C .. (PTACT(1)..PTACT(ACTU) contain the indices of upper bounds C that are active at XBAR). C Now, disorder PTACT(ACTU+1)..PTACT(INU) .. J = INU - ACTU CALL BCP04(J,PTACT(ACTU+1),ISEED) C .. to define the bounds that are not active at XBAR C (Indices are in PTACT(ACTU+1)..PTACT(INU)). J = ACTU + 1 DO 70 II = ACTU + 1,INU I = PTACT(J) AUX = WIDTH1 + RAND(ISEED)* (WIDTH2-WIDTH1) U(I) = XBAR(I) + AUX J = J + 1 70 CONTINUE END IF C Below PTACT is finally defined for which it was intended for. C That is .. C .. PTACT(1)..PTACT(ACTL) contains the indices of lower bounds C that are active at XBAR C .. PTACT(ACTL+1)..PTACT(ACTL+ACTU) contains the indices of upper C bounds that are active at XBAR J1 = 1 J2 = ACTL + 1 DO 80 I = 1,N IF (L(I).EQ.XBAR(I)) THEN PTACT(J1) = I J1 = J1 + 1 ELSE IF (U(I).EQ.XBAR(I)) THEN PTACT(J2) = I J2 = J2 + 1 END IF 80 CONTINUE C Define the Lagrange Multipliers at XBAR .. C .. How many active constraints will have Lagrange C Multipliers in the interval [LOWM1,UPPM1] .. NM1L = ACTL*DEG/100 NM1U = ACTU*DEG/100 C .. How many active constraints will have Lagrange C Multipliers in the interval [LOWM2,UPPM2] .. NM2L = ACTL - NM1L NM2U = ACTU - NM1U C .. Now, disorder PTACT(1)..PTACT(ACTL) .. J = ACTL CALL BCP04(J,PTACT(1),ISEED) C .. to determine randomly the values of the multipliers at XBAR C for the active lower bounds .. DO 90 J = 1,NM1L BETA(J) = LOWM1 + RAND(ISEED)* (UPPM1-LOWM1) 90 CONTINUE DO 100 J = NM1L + 1,NM1L + NM2L BETA(J) = LOWM2 + RAND(ISEED)* (UPPM2-LOWM2) 100 CONTINUE C .. Now, disorder PTACT(ACTL+1)..PTACT(ACTL+ACTU) .. J = ACTU CALL BCP04(J,PTACT(ACTL+1),ISEED) C .. to determine randomly the values of the multipliers at XBAR C for the active upper bounds .. DO 110 J = ACTL + 1,ACTL + NM1U BETA(J) = LOWM1 + RAND(ISEED)* (UPPM1-LOWM1) 110 CONTINUE DO 120 J = ACTL + NM1U + 1,ACTL + NM1U + NM2U BETA(J) = LOWM2 + RAND(ISEED)* (UPPM2-LOWM2) 120 CONTINUE C Copy HIF into HI HI = HIF RETURN END C .. C .. End of subroutine BCP01 .. C .. C .. Begin of subroutine BCP02 .. SUBROUTINE BCP02(INFRM) C _Authors: F. Facchinei, J. Soares and J. Judice C _Date Written: September 10, 1995 C _Date Modified: September 3, 1996 C _Purpose: C C BCP02 prints an adequate error message on the screen according C to the integer value in INFRM. C C .. Scalar Arguments .. INTEGER INFRM C .. IF (INFRM.EQ.0) THEN PRINT *,' ** INFRM = 0 ** EVERYTHING IS FINE **' ELSE IF (INFRM.EQ.1) THEN PRINT *,' ** INFRM = 1 ** N<=0 OR N > NN **' ELSE IF (INFRM.EQ.2) THEN PRINT *,' ** INFRM = 2 ** BND<0 OR BND>100 **' ELSE IF (INFRM.EQ.3) THEN PRINT *,' ** INFRM = 3 ** LOWBND<0 OR LOWBND>100 **' ELSE IF (INFRM.EQ.4) THEN PRINT *,' ** INFRM = 4 ** ACTBND<0 OR ACTBND>100 **' ELSE IF (INFRM.EQ.5) THEN PRINT *,' ** INFRM = 5 ** DEG<0 OR DEG>100 **' ELSE IF (INFRM.EQ.6) THEN PRINT *,' ** INFRM = 6 ** PAIR<> 0 AND 1 **' ELSE IF (INFRM.EQ.7) THEN PRINT *,' ** INFRM = 7 ** PAIR=1 AND LOWBND<>50 **' ELSE IF (INFRM.EQ.8) THEN PRINT *,' ** INFRM = 8 ** LOWM1<0 OR LOWM2<0 **' ELSE IF (INFRM.EQ.9) THEN PRINT *,' ** INFRM = 9 ** LOWM1>UPPM1 OR LOWM2>UPPM2 **' ELSE IF (INFRM.EQ.10) THEN PRINT *,' ** INFRM =10 **WIDTH1.LE.0 OR WIDTH1>WIDTH2 **' ELSE IF (INFRM.EQ.11) THEN PRINT *,' ** INFRM =11 ** HIF<1 OR HIF>3 **' ELSE IF (INFRM.EQ.12) THEN PRINT *,' ** INFRM =12 ** LOWBND TOO LARGE **' ELSE IF (INFRM.EQ.13) THEN PRINT *,' ** INFRM =13 ** LOWBND TOO LITTLE **' ELSE IF (INFRM.EQ.14) THEN PRINT *,' ** INFRM =14 ** ACTBND TOO LARGE **' ELSE PRINT *,' ** INFRM OUT OF RANGE **' END IF END C .. C .. End of subroutine BCP02 .. C .. C .. Begin of subroutine BCP03 .. SUBROUTINE BCP03(N,L,U) C _Authors: F. Facchinei, J. Soares and J. Judice C _Date Written: September 10, 1995 C _Date Modified: September 3, 1996 C _Purpose: C C BCP03 prints on the screen all the data that is in the C common blocks. Such information was obtained by a C previous call to BCP01. C C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. REAL L(N),U(N) C .. C .. Parameters .. C .. See the comments in BCP01 for a description of the variables C in the common blocks. INTEGER NN PARAMETER (NN=10000) C .. C .. Scalars in Common .. INTEGER ACTL,ACTU,HI,INL,INU,NM1L,NM1U,NM2L,NM2U C .. C .. Arrays in Common .. REAL BETA(NN),XBAR(NN) INTEGER PTACT(NN) C .. C .. Local Scalars .. INTEGER I C .. C .. Common blocks .. COMMON /BCPB1/INL,INU,ACTL,ACTU,NM1L,NM1U,NM2L,NM2U,HI COMMON /BCPB2/BETA,XBAR,PTACT C .. WRITE (*,FMT=9000) DO 10 I = 1,N WRITE (*,FMT=9010) I,L(I),XBAR(I),U(I),PTACT(I),BETA(I) 10 CONTINUE PRINT * PRINT *,' #Lower Bounds = ',INL PRINT *,' #Upper Bounds = ',INU PRINT *,' #Active Lower Bounds = ',ACTL PRINT *,' #Active Upper Bounds = ',ACTU PRINT *,' #Active Lower Bounds with mult in the', + ' first interval = ',NM1L PRINT *,' #Active Upper Bounds with mult in the', + ' first interval = ',NM1U PRINT *,' #Active Lower Bounds with mult in the', + ' second interval = ',NM2L PRINT *,' #Active Upper Bounds with mult in the', + ' second interval = ',NM2U PRINT *,' Function h_i: ',HI RETURN C .. Formats .. 9000 FORMAT (9X,'L',14X,'XBAR',12X,'U',8X,'PTACT',7X,'BETA') 9010 FORMAT (I3,1X,G14.7,1X,G14.7,1X,G14.7,1X,I3,3X,G14.7) END C .. C .. End of subroutine BCP03 .. C .. C .. Begin of subroutine BCP04 .. SUBROUTINE BCP04(N,VEC,ISEED) C _Authors: F. Facchinei, J. Soares and J. Judice C _Date Written: September 10, 1995 C _Date Modified: September 3, 1996 C _Purpose: C C Given a vector of integer components, BCP04 alters the position C of its elements randomly. C C _Arguments in parameter list: C C N (input) INTEGER C The dimension of the integer vector. C C VEC (input/output) INTEGER vector of dimension N C The vector whose components positions are to be randomly C changed. C C ISEED (input/output) INTEGER C A seed for the random number generator. It is modified. C .. Scalar Arguments .. INTEGER ISEED,N C .. C .. Array Arguments .. INTEGER VEC(N) C .. C .. External Functions .. REAL RAND EXTERNAL RAND C .. C .. Local Scalars .. INTEGER COPY,I,II C .. DO 10 I = N,1,-1 C VEC(J), for J=I+1,N will remain unaltered II = 1 + I*RAND(ISEED) C Exchange VEC(II) with VEC(I) COPY = VEC(I) VEC(I) = VEC(II) VEC(II) = COPY C VEC(J), for J=I,N will remain unaltered 10 CONTINUE RETURN END C .. C .. End of subroutine BCP04 .. C .. C .. Begin of subroutines HI0, HI1 and HI2 altogether .. C _Authors: F. Facchinei, J. Soares and J. Judice C _Date Written: September 10, 1995 C _Date Modified: September 3, 1996 C _Purpose: C C The routine: C C HI0 evaluates the function h_i C HI1 evaluates the derivative h'_i C HI2 evaluates the second derivative h''_i C C .. C .. begin of subroutine HI0 .. REAL FUNCTION HI0(XI,XIBAR,BETAI,HI) C .. Scalar Arguments .. REAL BETAI,XI,XIBAR INTEGER HI C .. C .. Local Scalars .. REAL AUX,AUX2 C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. IF (HI.EQ.1) THEN HI0 = BETAI* (XI-XIBAR) ELSE IF (HI.EQ.2) THEN AUX = XI - XIBAR HI0 = AUX*AUX*AUX + BETAI*AUX ELSE IF (HI.EQ.3) THEN AUX = XI - XIBAR IF (AUX.NE.0.0) THEN AUX2 = AUX* (ABS(AUX)** (4.0E+0/3.0E+0)) ELSE AUX2 = 0.0E+0 END IF HI0 = AUX2 + BETAI*AUX END IF END C .. C .. end of subroutine HI0 .. C .. C .. begin of subroutine HI1 .. REAL FUNCTION HI1(XI,XIBAR,BETAI,HI) C .. Scalar Arguments .. REAL BETAI,XI,XIBAR INTEGER HI C .. C .. Local Scalars .. REAL AUX,AUX2 C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. IF (HI.EQ.1) THEN HI1 = BETAI ELSE IF (HI.EQ.2) THEN AUX = XI - XIBAR HI1 = 3.0E+0*AUX*AUX + BETAI ELSE IF (HI.EQ.3) THEN AUX = XI - XIBAR IF (AUX.NE.0.0) THEN AUX2 = (ABS(AUX))** (4.0E+0/3.0E+0) ELSE AUX2 = 0.0E+0 END IF HI1 = 7.0E+0*AUX2/3.0E+0 + BETAI END IF END C .. C .. end of subroutine HI1 .. C .. C .. begin of subroutine HI2 .. REAL FUNCTION HI2(XI,XIBAR,HI) C .. Scalar Arguments .. REAL XI,XIBAR INTEGER HI C .. C .. Local Scalars .. REAL AUX,AUX2 C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. IF (HI.EQ.1) THEN HI2 = 0.0E+0 ELSE IF (HI.EQ.2) THEN HI2 = 6E+0* (XI-XIBAR) ELSE IF (HI.EQ.3) THEN AUX = XI - XIBAR IF (AUX.NE.0.0) THEN AUX2 = (ABS(AUX))** (4.0E+0/3.0E+0)/AUX ELSE AUX2 = 0.0E+0 END IF HI2 = 28.0E+0*AUX2/9.0E+0 END IF END C .. C .. end of subroutine HI2 .. C .. C .. No more routines HI. C .. C .. Begin of subroutine BCP10 .. SUBROUTINE BCP10(N,X,F) C _Authors: F. Facchinei, J. Soares and J. Judice C _Date Written: September 10, 1995 C _Date Modified: September 3, 1996 C _Purpose: C C BCP10 evaluates the objective function of the generated problem C at a given point. C C _Arguments in parameter list: C C N (input) INTEGER C The number of variables in the problem. C C X (input) DOUBLE PRECISION vector of dimension N C The point where the objective function is to be evaluated. C C F (output) DOUBLE PRECISION C The objective function value. C .. Scalar Arguments .. REAL F INTEGER N C .. C .. Array Arguments .. REAL X(N) C .. C .. External Subroutines .. EXTERNAL FUNCTG C .. C .. External Functions .. REAL HI0 EXTERNAL HI0 C .. C .. Parameters .. C .. See the comments in BCP01 for a description of the variables C in the common blocks. INTEGER NN PARAMETER (NN=10000) C .. C .. Scalars in Common .. INTEGER ACTL,ACTU,HI,INL,INU,NM1L,NM1U,NM2L,NM2U C .. C .. Arrays in Common .. REAL BETA(NN),XBAR(NN) INTEGER PTACT(NN) C .. C .. Local Scalars .. INTEGER I,II C .. C .. Common blocks .. COMMON /BCPB1/INL,INU,ACTL,ACTU,NM1L,NM1U,NM2L,NM2U,HI COMMON /BCPB2/BETA,XBAR,PTACT C .. C Evaluate the unconstrained objective function CALL FUNCTG(N,X,F) C Use the function h_i according to the selection in HI DO 10 II = 1,ACTL I = PTACT(II) F = F + HI0(X(I),XBAR(I),BETA(II),HI) 10 CONTINUE DO 20 II = ACTL + 1,ACTU I = PTACT(II) F = F - HI0(X(I),XBAR(I),BETA(II),HI) 20 CONTINUE RETURN END C .. C .. End of subroutine BCP10. C .. C .. Begin of subroutine BCP11 .. SUBROUTINE BCP11(N,X,GF) C _Authors: F. Facchinei, J. Soares and J. Judice C _Date Written: September 10, 1995 C _Date Modified: September 3, 1996 C _Purpose: C C BCP11 evaluates the gradient of the objective function of the C generated problem at a given point. C C _Arguments in parameter list: C C N (input) INTEGER C The number of variables in the problem. C C X (input) DOUBLE PRECISION vector of dimension N C The point where the gradient of the objective function is to C be evaluated. C C GF (output) DOUBLE PRECISION vector of dimension N C The gradient of objective function values. C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. REAL GF(N),X(N) C .. C .. External Subroutines .. EXTERNAL GRADG C .. C .. External Functions .. REAL HI1 EXTERNAL HI1 C .. C .. Parameters .. C .. See the comments in BCP01 for a description of the variables C in the common blocks. INTEGER NN PARAMETER (NN=10000) C .. C .. Scalars in Common .. INTEGER ACTL,ACTU,HI,INL,INU,NM1L,NM1U,NM2L,NM2U C .. C .. Arrays in Common .. REAL BETA(NN),XBAR(NN) INTEGER PTACT(NN) C .. C .. Local Scalars .. INTEGER I,II C .. C .. Common blocks .. COMMON /BCPB1/INL,INU,ACTL,ACTU,NM1L,NM1U,NM2L,NM2U,HI COMMON /BCPB2/BETA,XBAR,PTACT C .. C Evaluate the unconstrained objective function gradient CALL GRADG(N,X,GF) DO 10 II = 1,ACTL I = PTACT(II) GF(I) = GF(I) + HI1(X(I),XBAR(I),BETA(II),HI) 10 CONTINUE DO 20 II = ACTL + 1,ACTU I = PTACT(II) GF(I) = GF(I) - HI1(X(I),XBAR(I),BETA(II),HI) 20 CONTINUE RETURN END C .. C .. End of subroutine BCP11. C .. C .. Begin of subroutine BCP12 .. SUBROUTINE BCP12(N,X,D,HFD) C _Authors: F. Facchinei, J. Soares and J. Judice C _Date Written: September 10, 1995 C _Date Modified: September 3, 1996 C _Purpose: C C BCP12 evaluates the product between the Hessian matrix of C the objective function at a given point and a vector. C C _Arguments in parameter list: C C N (input) INTEGER C The number of variables in the problem. C C X (input) DOUBLE PRECISION vector of dimension N C The point where the Hessian matrix of the objective function C is to be evaluated. C C D (input) DOUBLE PRECISION vector of dimension N C The vector that will be multiplied by the Hessian matrix. C C HFD (output) DOUBLE PRECISION vector of dimension N C The resulting vector. C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. REAL D(N),HFD(N),X(N) C .. C .. External Subroutines .. EXTERNAL HESGD C .. C .. External Functions .. REAL HI2 EXTERNAL HI2 C .. C .. Parameters .. C .. See the comments in BCP01 for a description of the variables in C the common blocks. INTEGER NN PARAMETER (NN=10000) C .. C .. Scalars in Common .. INTEGER ACTL,ACTU,HI,INL,INU,NM1L,NM1U,NM2L,NM2U C .. C .. Arrays in Common .. REAL BETA(NN),XBAR(NN) INTEGER PTACT(NN) C .. C .. Local Scalars .. INTEGER I,II C .. C .. Common blocks .. COMMON /BCPB1/INL,INU,ACTL,ACTU,NM1L,NM1U,NM2L,NM2U,HI COMMON /BCPB2/BETA,XBAR,PTACT C .. C Evaluate the unconstrained objective function Hessian X d CALL HESGD(N,X,D,HFD) DO 10 II = 1,ACTL I = PTACT(II) HFD(I) = HFD(I) + HI2(X(I),XBAR(I),HI)*D(I) 10 CONTINUE DO 20 II = ACTL + 1,ACTU I = PTACT(II) HFD(I) = HFD(I) - HI2(X(I),XBAR(I),HI)*D(I) 20 CONTINUE RETURN END C .. C .. End of subroutine BCP12 C .. C .. Begin of subroutine RAND REAL FUNCTION RAND(ISEED) C ********** C C function rand C C Rand is the portable random number generator of L. Schrage. C C The generator is full cycle, that is, every integer from C 1 to 2**31 - 2 is generated exactly once in the cycle. C It is completely described in TOMS 5(1979),132-138. C C The function statement is C C real function rand(iseed) C C where C C iseed is a positive integer variable which specifies C the seed to the random number generator. Given the C input seed, rand returns a random number in the C open interval (0,1). On output the seed is updated. C C Argonne National Laboratory. MINPACK Project. March 1981. C C ********** C .. Scalar Arguments .. INTEGER ISEED C .. C .. Local Scalars .. REAL C INTEGER A,B15,B16,FHI,K,LEFTLO,P,XALO,XHI C .. C .. Intrinsic Functions .. INTRINSIC FLOAT C .. C .. Data statements .. C C Set a = 7**5, b15 = 2**15, b16 = 2**16, p = 2**31-1, c = 1/p. C DATA A/16807/,B15/32768/,B16/65536/,P/2147483647/ DATA C/4.656612875E-10/ C .. C C There are 8 steps in rand. C C 1. Get 15 hi order bits of iseed. C 2. Get 16 lo bits of iseed and form lo product. C 3. Get 15 hi order bits of lo product. C 4. Form the 31 highest bits of full product. C 5. Get overflo past 31st bit of full product. C 6. Assemble all the parts and pre-substract p. C The parentheses are essential. C 7. Add p back if necessary. C 8. Multiply by 1/(2**31-1). C XHI = ISEED/B16 XALO = (ISEED-XHI*B16)*A LEFTLO = XALO/B16 FHI = XHI*A + LEFTLO K = FHI/B15 ISEED = (((XALO-LEFTLO*B16)-P)+ (FHI-K*B15)*B16) + K IF (ISEED.LT.0) ISEED = ISEED + P RAND = C*FLOAT(ISEED) RETURN C C Last card of function rand. C END SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << \SHAR_EOF > 'src.f' C====================================================================== C c C This file contains the FORTRAN routines that the user should c C link to his/her code in order to employ a box constrained c C generated problem. c C Other routines that should be linked, besides the main program, c C include subroutines to evaluate the unconstrained function c C as well as its derivatives. c C c C The user is referred to the documentation c C c C Soares, J., Judice, J. and Facchinei, F., "Generating box c C constrained optimization problems", c C ACM TOMS, (??) c C and c C Soares, J., Judice, J. and Facchinei, F., "FORTRAN subroutines c C for generating box constrained optimization problems", c C ACM TOMS, (??). c C c C To change the precision of this set of subroutines the user c C should do the following: c C c C 1. From double to real: c C 1.1. Replace all 'DOUBLE PRECISION' with 'REAL' c C 1.2. >> 'D+' >> 'E+' c C c C 2. From real to double: c C 1.1. Replace all 'REAL' with 'DOUBLE PRECISION' c C 1.2. >> 'E+' >> 'D+' c C c C IMPORTANT EXCEPTION: Do not replace lower case 'real'. In case c C of doubt make sure that RAND in its definition and calls is c C always single precision. c C c C====================================================================== C .. C .. Begin of subroutine BCP01 .. SUBROUTINE BCP01(N,L,U,BND,LOWBND,ACTBND,PAIR,DEG,LOWM1,UPPM1, + LOWM2,UPPM2,WIDTH1,WIDTH2,HIF,ISEED,INFRM) C _Authors: F. Facchinei, J. Soares and J. Judice C _Date Written: September 10, 1995 C _Date Modified: September 3, 1996 C _Purpose: C C BCP01 generates box constrained nonlinear optimization C problems of the form C C Minimize f(x) C s.t. l <= x <= u C C where x belongs to R^n and f is a continuously differentiable C function R^n -> R. C C More precisely, BCP01 sets up the bounds, and the data required C to evaluate the objective function f and its derivatives. C Unsuccessful generation is detected by having a nonzero value C in INFRM as output. C C _Arguments in parameter list: C C N (input) INTEGER C The number of variables in the problem. (Must be <= NN, C see common block below) C C L, U (output) DOUBLE PRECISION vectors of dimension N C The bounds of the new problem. C C BND (input) INTEGER C The number of bound constraints (or simply, bounds), C expressed as the percentage of the 2*n possible constraints, C constraints,that are really existing (Must be in [0,100]). C Example: if N is 500 and BND=40 then the problem has 400 C bounds. C C LOWBND (input) INTEGER C The number of lower bounds, expressed as the percentage of C the number of constraints as determined by BND (Must be in C [0,100]). This quantity also determines the number of upper C bounds. C Example: if N is 500, BND is 40 and LOWBND=75 then there C are 300 (100) lower (upper) bounds. C C ACTBND (input) INTEGER C The number of active bounds at the unconstrained solution, C expressed as the percentage of the number of constraints as C determined by BND (Must be in [0,100]). These active C bounds are distributed proportionally to LOWBND. C Example: Continuing the previous example, if ACTBND C is 50 then 150 (50) of the lower (upper) bounds should C be active at the unconstrained solution. C C PAIR (input) INTEGER C If this parameter is set to 1 then the constraints as C determined by BND must come in pairs, i.e. each variable C has either no constraint or both lower and upper bounds (in C which case LOWBND must be set to 50). If PAIR is C set to 0 then the indices of variables having upper bounds C are chosen randomly and independently of the lower bounds C (Must be 0 or 1). C C DEG (input) INTEGER C The number of Lagrange multipliers (at the unconstrained C solution) uniformly distributed in the interval C [LOWM1, UPPM1] (see below), expressed as the percentage of C the active constraints (Must be in [0,100]). C The remaining active constraints have Lagrange multipliers C uniformly distributed in [LOWM2, UPPM2] (see below). C Example: C Continuing the example, if DEG=10 then 15 (5) active C lower (upper) bounds (at the unconstrained solution) have C Lagrange multipliers uniformly distributed in [LOWM1, C UPPM1], while the remaining active bounds have Lagrange C multipliers uniformly distributed in [LOWM2, UPPM2]. C C LOWM1, UPPM1, LOWM2, UPPM2 (input) DOUBLE PRECISION C These parameters specify the ranges for the C Lagrange multipliers according to what explained above. C (Must be >= 0) C C WIDTH1, WIDTH2 (input) DOUBLE PRECISION C These two parameters control the value of the bounds that C are not active at the unconstrained solution. More C precisely, if the i-th component has a lower bound not C active at the unconstrained solution xbar then this bound C is randomly chosen in the interval C [xbar_i - WIDTH2, xbar_i - WIDTH1]. Analogously, C if the i-th component has an upper bound not active at C the unconstrained solution xbar then this bound is C randomly chosen in the interval C [xbar_i + WIDTH1, xbar_i + WIDTH2] (Must be >= 0). C C HIF (input) INTEGER C This parameter determines which function h_i is to be used C by the generator. If HIF is equal to C C 1: h_i(x_i)= beta_i * ( x_i - xbar_i ); C C 2: h_i(x_i)=( x_i - xbar_i )^3 + beta_i*( x_i - xbar_i ); C C 3: h_i(x_i)=(x_i - xbar_i)^(7/3) + beta_i*(x_i - xbar_i). C C (Must be 1, 2 or 3) C C ISEED (input) INTEGER C A seed for the random number generator. By running the box C constrained problems generator with the same choice of C parameters but two different seeds, two different box C constrained problems are generated with the same overall C characteristics (Must be >= 0). C C INFRM (output) INTEGER C Determines if a successful call to BCP01 has occurred C (INFRM=0) or not (INFRM<>0). C If INFRM is equal to C 1: N is less than 1 or greater than NN C 2, 3, 4, 5, 6, 11: Some integer parameter is defined C out of its range. C 7: PAIR=1 but LOWBND is not 50. C 8: Either LOWM1 or LOWM2 is less than zero. C 9: Either LOWM1 is greater than UPPM1 or LOWM2 is greater C than UPPM2. C 10: Either WIDTH1 is less than zero or WIDTH1 is greater C than WIDTH2. C 12: It is impossible to make a selection because LOWBND C is too large. C 13: It is impossible to make a selection because LOWBND C is too small. C 14: It is impossible to make a selection because ACTBND C is too large. C C _Explanation of the variables used in the common blocks: C C NN INTEGER C Maximum N allowed. C C INL INTEGER C Number of lower bounds. C C INU INTEGER C Number of upper bounds. C C ACTL INTEGER C Number of lower bounds active at the unconstrained solution C C ACTU INTEGER C Number of upper bounds active at the unconstrained solution C C NM1L INTEGER C Number of lower bounds active at the unconstrained solution C which have a Lagrange multiplier randomly chosen in C [LOWM1,UPPM1]. C C NM1U INTEGER C Number of upper bounds active at the unconstrained solution C which have a Lagrange multiplier randomly chosen in C [LOWM1,UPPM1]. C C NM2L INTEGER C Number of lower bounds active at the unconstrained solution C which have a Lagrange multiplier randomly chosen in C [LOWM2,UPPM2]. C C NM2U INTEGER C Number of upper bounds active at the unconstrained solution C which have a Lagrange multiplier randomly chosen in C [LOWM2,UPPM2]. C C HI INTEGER C Characterizes the function h_i used. HI=HIF. C C PTACT INTEGER vector of dimension NN C A vector of pointers that characterizes the bounds that are C active at the unconstrained solution, i.e., C PTACT(II) (for II=1,ACTL) pinpoints an active lower bound C PTACT(II) (for II=ACTL+1,ACTL+ACTU) pinpoints an active C upper bound. C C BETA DOUBLE PRECISION vector of dimension NN C At the unconstrained solution, BETA characterizes the C multipliers of the bounds that are active, i.e., C BETA(II) (for II=1,ACTL) is the multiplier associated C with the lower bound in the PTBETA(II) variable; C BETA(II) (for II=ACTL+1,ACTL+ACTU) is the multiplier C associated with the upper bound in the PTBETA(II) variable. C C XBAR DOUBLE PRECISION vector of dimension NN C The unconstrained solution. C .. Parameters .. INTEGER NN PARAMETER (NN=10000) DOUBLE PRECISION BIGINF PARAMETER (BIGINF=1.0D+20) C .. C .. Scalar Arguments .. DOUBLE PRECISION LOWM1,LOWM2,UPPM1,UPPM2,WIDTH1,WIDTH2 INTEGER ACTBND,BND,DEG,HIF,INFRM,ISEED,LOWBND,N,PAIR C .. C .. Array Arguments .. DOUBLE PRECISION L(N),U(N) C .. C .. Scalars in Common .. INTEGER ACTL,ACTU,HI,INL,INU,NM1L,NM1U,NM2L,NM2U C .. C .. Arrays in Common .. DOUBLE PRECISION BETA(NN),XBAR(NN) INTEGER PTACT(NN) C .. C .. Local Scalars .. DOUBLE PRECISION AUX INTEGER I,II,INBND,J,J1,J2,JJ C .. C .. External Functions .. REAL RAND EXTERNAL RAND C .. C .. External Subroutines .. EXTERNAL BCP04,SOLUTG C .. C .. Intrinsic Functions .. INTRINSIC MIN C .. C .. Common blocks .. COMMON /BCPB1/INL,INU,ACTL,ACTU,NM1L,NM1U,NM2L,NM2U,HI COMMON /BCPB2/BETA,XBAR,PTACT C .. C Test the input parameters. INFRM = 0 IF ((N.GT.NN) .OR. (N.LE.0)) INFRM = 1 IF ((BND.LT.0) .OR. (BND.GT.100)) INFRM = 2 IF ((LOWBND.LT.0) .OR. (LOWBND.GT.100)) INFRM = 3 IF ((ACTBND.LT.0) .OR. (ACTBND.GT.100)) INFRM = 4 IF ((DEG.LT.0) .OR. (DEG.GT.100)) INFRM = 5 IF ((PAIR.NE.0) .AND. (PAIR.NE.1)) INFRM = 6 IF ((PAIR.EQ.1) .AND. (LOWBND.NE.50)) INFRM = 7 IF ((LOWM1.LT.0D0) .OR. (LOWM2.LT.0D0)) INFRM = 8 IF ((LOWM1.GT.UPPM1) .OR. (LOWM2.GT.UPPM2)) INFRM = 9 IF ((WIDTH1.LE.0D0) .OR. (WIDTH1.GT.WIDTH2)) INFRM = 10 IF ((HIF.LE.0) .OR. (HIF.GE.4)) INFRM = 11 IF (INFRM.NE.0) RETURN C Determine the number of constraints INBND = (BND*2*N)/100 C Determine the number of lower bounds INL = (LOWBND*INBND)/100 C Safeguard the maximum number of lower bounds IF (INL.GT.N) THEN INFRM = 12 RETURN END IF C Determine the number of upper bounds INU = INBND - INL C Safeguard the maximum number of upper bounds IF (INU.GT.N) THEN INFRM = 13 RETURN END IF C Determine the number of active constraints at XBAR I = (ACTBND*INBND)/100 ACTL = (ACTBND*INL)/100 ACTU = I - ACTL C Safeguard some settings related to the parameter PAIR .. IF (PAIR.EQ.0) THEN C .. Since the bounds are supposed to be chosen randomly, C safeguard the total number of constraints C active at XBAR. IF (ACTL+ACTU.GT.N) THEN INFRM = 14 RETURN END IF ELSE C .. Since the bounds are supposed to come in pairs, C safeguard that the number of lower bounds must C equal the number of upper bounds .. J = MIN(INL,INU) INL = J INU = J C .. in which case the number of bounds and active bounds at XBAR C should be corrected .. INBND = INL + INU ACTL = (ACTBND*INL)/100 ACTU = ACTL C .. there should not exist one pair of bounds that are both C active at XBAR (Note: INL=INU). IF (ACTL+ACTU.GT.INL) THEN INFRM = 14 RETURN END IF END IF C Read the optimal unconstrained point CALL SOLUTG(N,XBAR) C Set default values to nonexistent bounds and make some random C settings. For the moment, PTACT is used as an auxiliary vector. C After a call to BCP04 it will contain the set of integer values C 1..N in a random order. DO 10 I = 1,N L(I) = -BIGINF U(I) = BIGINF PTACT(I) = I 10 CONTINUE CALL BCP04(N,PTACT(1),ISEED) C Defining the lower bounds .. C .. actives at XBAR (Indices are in PTACT(1)..PTACT(ACTL)) .. J = 1 DO 20 II = 1,ACTL I = PTACT(J) L(I) = XBAR(I) J = J + 1 20 CONTINUE C .. and non actives at XBAR C (Indices are in PTACT(ACTL+1)..PTACT(INL)) .. DO 30 II = ACTL + 1,INL I = PTACT(J) AUX = WIDTH1 + RAND(ISEED)* (WIDTH2-WIDTH1) L(I) = XBAR(I) - AUX J = J + 1 30 CONTINUE C .. (Note that PTACT(INL+1)..PTACT(N) contain indices of variables C that do not have lower bounds). C Defining the upper bounds .. IF (PAIR.EQ.0) THEN C .. since the bounds are supposed to be chosen randomly, C disorder PTACT(ACTL+1)..PTACT(N) .. J = N - ACTL CALL BCP04(J,PTACT(ACTL+1),ISEED) C .. and determine the upper bounds that are active at XBAR C (Indices are in PTACT(ACTL+1)..PTACT(ACTL+ACTU) ) .. J = 1 JJ = ACTL + ACTU DO 40 II = 1,ACTU I = PTACT(JJ) U(I) = XBAR(I) PTACT(JJ) = PTACT(J) PTACT(J) = I J = J + 1 JJ = JJ - 1 40 CONTINUE C .. (PTACT(1)..PTACT(ACTU) contain the indices of upper bounds C that are active at XBAR) C Now, disorder PTACT(ACTU+1)..PTACT(N) .. J = N - ACTU CALL BCP04(J,PTACT(ACTU+1),ISEED) C C .. to define the upper bounds that are not active at XBAR C (Indices are in PTACT(ACTU+1)..PTACT(INU)) .. J = ACTU + 1 DO 50 II = ACTU + 1,INU I = PTACT(J) AUX = WIDTH1 + RAND(ISEED)* (WIDTH2-WIDTH1) U(I) = XBAR(I) + AUX J = J + 1 50 CONTINUE ELSE C .. since the bounds are supposed to come in pairs, C disorder PTACT(ACTL+1)..PTACT(INL) .. J = INL - ACTL CALL BCP04(J,PTACT(ACTL+1),ISEED) C .. and determine the upper bounds that are actives at XBAR C (Indices are in PTACT(ACTL+1)..PTACT(ACTL+ACTU). Also note that C ACTL+ACTU.LE.INL=INU) .. J = 1 JJ = ACTL + ACTU DO 60 II = 1,ACTU I = PTACT(JJ) U(I) = XBAR(I) PTACT(JJ) = PTACT(J) PTACT(J) = I J = J + 1 JJ = JJ - 1 60 CONTINUE C .. (PTACT(1)..PTACT(ACTU) contain the indices of upper bounds C that are active at XBAR). C Now, disorder PTACT(ACTU+1)..PTACT(INU) .. J = INU - ACTU CALL BCP04(J,PTACT(ACTU+1),ISEED) C .. to define the bounds that are not active at XBAR C (Indices are in PTACT(ACTU+1)..PTACT(INU)). J = ACTU + 1 DO 70 II = ACTU + 1,INU I = PTACT(J) AUX = WIDTH1 + RAND(ISEED)* (WIDTH2-WIDTH1) U(I) = XBAR(I) + AUX J = J + 1 70 CONTINUE END IF C Below PTACT is finally defined for which it was intended for. C That is .. C .. PTACT(1)..PTACT(ACTL) contains the indices of lower bounds C that are active at XBAR C .. PTACT(ACTL+1)..PTACT(ACTL+ACTU) contains the indices of upper C bounds that are active at XBAR J1 = 1 J2 = ACTL + 1 DO 80 I = 1,N IF (L(I).EQ.XBAR(I)) THEN PTACT(J1) = I J1 = J1 + 1 ELSE IF (U(I).EQ.XBAR(I)) THEN PTACT(J2) = I J2 = J2 + 1 END IF 80 CONTINUE C Define the Lagrange Multipliers at XBAR .. C .. How many active constraints will have Lagrange C Multipliers in the interval [LOWM1,UPPM1] .. NM1L = ACTL*DEG/100 NM1U = ACTU*DEG/100 C .. How many active constraints will have Lagrange C Multipliers in the interval [LOWM2,UPPM2] .. NM2L = ACTL - NM1L NM2U = ACTU - NM1U C .. Now, disorder PTACT(1)..PTACT(ACTL) .. J = ACTL CALL BCP04(J,PTACT(1),ISEED) C .. to determine randomly the values of the multipliers at XBAR C for the active lower bounds .. DO 90 J = 1,NM1L BETA(J) = LOWM1 + RAND(ISEED)* (UPPM1-LOWM1) 90 CONTINUE DO 100 J = NM1L + 1,NM1L + NM2L BETA(J) = LOWM2 + RAND(ISEED)* (UPPM2-LOWM2) 100 CONTINUE C .. Now, disorder PTACT(ACTL+1)..PTACT(ACTL+ACTU) .. J = ACTU CALL BCP04(J,PTACT(ACTL+1),ISEED) C .. to determine randomly the values of the multipliers at XBAR C for the active upper bounds .. DO 110 J = ACTL + 1,ACTL + NM1U BETA(J) = LOWM1 + RAND(ISEED)* (UPPM1-LOWM1) 110 CONTINUE DO 120 J = ACTL + NM1U + 1,ACTL + NM1U + NM2U BETA(J) = LOWM2 + RAND(ISEED)* (UPPM2-LOWM2) 120 CONTINUE C Copy HIF into HI HI = HIF RETURN END C .. C .. End of subroutine BCP01 .. C .. C .. Begin of subroutine BCP02 .. SUBROUTINE BCP02(INFRM) C _Authors: F. Facchinei, J. Soares and J. Judice C _Date Written: September 10, 1995 C _Date Modified: September 3, 1996 C _Purpose: C C BCP02 prints an adequate error message on the screen according C to the integer value in INFRM. C C .. Scalar Arguments .. INTEGER INFRM C .. IF (INFRM.EQ.0) THEN PRINT *,' ** INFRM = 0 ** EVERYTHING IS FINE **' ELSE IF (INFRM.EQ.1) THEN PRINT *,' ** INFRM = 1 ** N<=0 OR N > NN **' ELSE IF (INFRM.EQ.2) THEN PRINT *,' ** INFRM = 2 ** BND<0 OR BND>100 **' ELSE IF (INFRM.EQ.3) THEN PRINT *,' ** INFRM = 3 ** LOWBND<0 OR LOWBND>100 **' ELSE IF (INFRM.EQ.4) THEN PRINT *,' ** INFRM = 4 ** ACTBND<0 OR ACTBND>100 **' ELSE IF (INFRM.EQ.5) THEN PRINT *,' ** INFRM = 5 ** DEG<0 OR DEG>100 **' ELSE IF (INFRM.EQ.6) THEN PRINT *,' ** INFRM = 6 ** PAIR<> 0 AND 1 **' ELSE IF (INFRM.EQ.7) THEN PRINT *,' ** INFRM = 7 ** PAIR=1 AND LOWBND<>50 **' ELSE IF (INFRM.EQ.8) THEN PRINT *,' ** INFRM = 8 ** LOWM1<0 OR LOWM2<0 **' ELSE IF (INFRM.EQ.9) THEN PRINT *,' ** INFRM = 9 ** LOWM1>UPPM1 OR LOWM2>UPPM2 **' ELSE IF (INFRM.EQ.10) THEN PRINT *,' ** INFRM =10 **WIDTH1.LE.0 OR WIDTH1>WIDTH2 **' ELSE IF (INFRM.EQ.11) THEN PRINT *,' ** INFRM =11 ** HIF<1 OR HIF>3 **' ELSE IF (INFRM.EQ.12) THEN PRINT *,' ** INFRM =12 ** LOWBND TOO LARGE **' ELSE IF (INFRM.EQ.13) THEN PRINT *,' ** INFRM =13 ** LOWBND TOO LITTLE **' ELSE IF (INFRM.EQ.14) THEN PRINT *,' ** INFRM =14 ** ACTBND TOO LARGE **' ELSE PRINT *,' ** INFRM OUT OF RANGE **' END IF END C .. C .. End of subroutine BCP02 .. C .. C .. Begin of subroutine BCP03 .. SUBROUTINE BCP03(N,L,U) C _Authors: F. Facchinei, J. Soares and J. Judice C _Date Written: September 10, 1995 C _Date Modified: September 3, 1996 C _Purpose: C C BCP03 prints on the screen all the data that is in the C common blocks. Such information was obtained by a C previous call to BCP01. C C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. DOUBLE PRECISION L(N),U(N) C .. C .. Parameters .. C .. See the comments in BCP01 for a description of the variables C in the common blocks. INTEGER NN PARAMETER (NN=10000) C .. C .. Scalars in Common .. INTEGER ACTL,ACTU,HI,INL,INU,NM1L,NM1U,NM2L,NM2U C .. C .. Arrays in Common .. DOUBLE PRECISION BETA(NN),XBAR(NN) INTEGER PTACT(NN) C .. C .. Local Scalars .. INTEGER I C .. C .. Common blocks .. COMMON /BCPB1/INL,INU,ACTL,ACTU,NM1L,NM1U,NM2L,NM2U,HI COMMON /BCPB2/BETA,XBAR,PTACT C .. WRITE (*,FMT=9000) DO 10 I = 1,N WRITE (*,FMT=9010) I,L(I),XBAR(I),U(I),PTACT(I),BETA(I) 10 CONTINUE PRINT * PRINT *,' #Lower Bounds = ',INL PRINT *,' #Upper Bounds = ',INU PRINT *,' #Active Lower Bounds = ',ACTL PRINT *,' #Active Upper Bounds = ',ACTU PRINT *,' #Active Lower Bounds with mult in the', + ' first interval = ',NM1L PRINT *,' #Active Upper Bounds with mult in the', + ' first interval = ',NM1U PRINT *,' #Active Lower Bounds with mult in the', + ' second interval = ',NM2L PRINT *,' #Active Upper Bounds with mult in the', + ' second interval = ',NM2U PRINT *,' Function h_i: ',HI RETURN C .. Formats .. 9000 FORMAT (9X,'L',14X,'XBAR',12X,'U',8X,'PTACT',7X,'BETA') 9010 FORMAT (I3,1X,G14.7,1X,G14.7,1X,G14.7,1X,I3,3X,G14.7) END C .. C .. End of subroutine BCP03 .. C .. C .. Begin of subroutine BCP04 .. SUBROUTINE BCP04(N,VEC,ISEED) C _Authors: F. Facchinei, J. Soares and J. Judice C _Date Written: September 10, 1995 C _Date Modified: September 3, 1996 C _Purpose: C C Given a vector of integer components, BCP04 alters the position C of its elements randomly. C C _Arguments in parameter list: C C N (input) INTEGER C The dimension of the integer vector. C C VEC (input/output) INTEGER vector of dimension N C The vector whose components positions are to be randomly C changed. C C ISEED (input/output) INTEGER C A seed for the random number generator. It is modified. C .. Scalar Arguments .. INTEGER ISEED,N C .. C .. Array Arguments .. INTEGER VEC(N) C .. C .. External Functions .. REAL RAND EXTERNAL RAND C .. C .. Local Scalars .. INTEGER COPY,I,II C .. DO 10 I = N,1,-1 C VEC(J), for J=I+1,N will remain unaltered II = 1 + I*RAND(ISEED) C Exchange VEC(II) with VEC(I) COPY = VEC(I) VEC(I) = VEC(II) VEC(II) = COPY C VEC(J), for J=I,N will remain unaltered 10 CONTINUE RETURN END C .. C .. End of subroutine BCP04 .. C .. C .. Begin of subroutines HI0, HI1 and HI2 altogether .. C _Authors: F. Facchinei, J. Soares and J. Judice C _Date Written: September 10, 1995 C _Date Modified: September 3, 1996 C _Purpose: C C The routine: C C HI0 evaluates the function h_i C HI1 evaluates the derivative h'_i C HI2 evaluates the second derivative h''_i C C .. C .. begin of subroutine HI0 .. DOUBLE PRECISION FUNCTION HI0(XI,XIBAR,BETAI,HI) C .. Scalar Arguments .. DOUBLE PRECISION BETAI,XI,XIBAR INTEGER HI C .. C .. Local Scalars .. DOUBLE PRECISION AUX,AUX2 C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. IF (HI.EQ.1) THEN HI0 = BETAI* (XI-XIBAR) ELSE IF (HI.EQ.2) THEN AUX = XI - XIBAR HI0 = AUX*AUX*AUX + BETAI*AUX ELSE IF (HI.EQ.3) THEN AUX = XI - XIBAR IF (AUX.NE.0.0D0) THEN AUX2 = AUX* (ABS(AUX)** (4.0D+0/3.0D+0)) ELSE AUX2 = 0.0D+0 END IF HI0 = AUX2 + BETAI*AUX END IF END C .. C .. end of subroutine HI0 .. C .. C .. begin of subroutine HI1 .. DOUBLE PRECISION FUNCTION HI1(XI,XIBAR,BETAI,HI) C .. Scalar Arguments .. DOUBLE PRECISION BETAI,XI,XIBAR INTEGER HI C .. C .. Local Scalars .. DOUBLE PRECISION AUX,AUX2 C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. IF (HI.EQ.1) THEN HI1 = BETAI ELSE IF (HI.EQ.2) THEN AUX = XI - XIBAR HI1 = 3.0D+0*AUX*AUX + BETAI ELSE IF (HI.EQ.3) THEN AUX = XI - XIBAR IF (AUX.NE.0.0D0) THEN AUX2 = (ABS(AUX))** (4.0D+0/3.0D+0) ELSE AUX2 = 0.0D+0 END IF HI1 = 7.0D+0*AUX2/3.0D+0 + BETAI END IF END C .. C .. end of subroutine HI1 .. C .. C .. begin of subroutine HI2 .. DOUBLE PRECISION FUNCTION HI2(XI,XIBAR,HI) C .. Scalar Arguments .. DOUBLE PRECISION XI,XIBAR INTEGER HI C .. C .. Local Scalars .. DOUBLE PRECISION AUX,AUX2 C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. IF (HI.EQ.1) THEN HI2 = 0.0D+0 ELSE IF (HI.EQ.2) THEN HI2 = 6D+0* (XI-XIBAR) ELSE IF (HI.EQ.3) THEN AUX = XI - XIBAR IF (AUX.NE.0.0D0) THEN AUX2 = (ABS(AUX))** (4.0D+0/3.0D+0)/AUX ELSE AUX2 = 0.0D+0 END IF HI2 = 28.0D+0*AUX2/9.0D+0 END IF END C .. C .. end of subroutine HI2 .. C .. C .. No more routines HI. C .. C .. Begin of subroutine BCP10 .. SUBROUTINE BCP10(N,X,F) C _Authors: F. Facchinei, J. Soares and J. Judice C _Date Written: September 10, 1995 C _Date Modified: September 3, 1996 C _Purpose: C C BCP10 evaluates the objective function of the generated problem C at a given point. C C _Arguments in parameter list: C C N (input) INTEGER C The number of variables in the problem. C C X (input) DOUBLE PRECISION vector of dimension N C The point where the objective function is to be evaluated. C C F (output) DOUBLE PRECISION C The objective function value. C .. Scalar Arguments .. DOUBLE PRECISION F INTEGER N C .. C .. Array Arguments .. DOUBLE PRECISION X(N) C .. C .. External Subroutines .. EXTERNAL FUNCTG C .. C .. External Functions .. DOUBLE PRECISION HI0 EXTERNAL HI0 C .. C .. Parameters .. C .. See the comments in BCP01 for a description of the variables C in the common blocks. INTEGER NN PARAMETER (NN=10000) C .. C .. Scalars in Common .. INTEGER ACTL,ACTU,HI,INL,INU,NM1L,NM1U,NM2L,NM2U C .. C .. Arrays in Common .. DOUBLE PRECISION BETA(NN),XBAR(NN) INTEGER PTACT(NN) C .. C .. Local Scalars .. INTEGER I,II C .. C .. Common blocks .. COMMON /BCPB1/INL,INU,ACTL,ACTU,NM1L,NM1U,NM2L,NM2U,HI COMMON /BCPB2/BETA,XBAR,PTACT C .. C Evaluate the unconstrained objective function CALL FUNCTG(N,X,F) C Use the function h_i according to the selection in HI DO 10 II = 1,ACTL I = PTACT(II) F = F + HI0(X(I),XBAR(I),BETA(II),HI) 10 CONTINUE DO 20 II = ACTL + 1,ACTU I = PTACT(II) F = F - HI0(X(I),XBAR(I),BETA(II),HI) 20 CONTINUE RETURN END C .. C .. End of subroutine BCP10. C .. C .. Begin of subroutine BCP11 .. SUBROUTINE BCP11(N,X,GF) C _Authors: F. Facchinei, J. Soares and J. Judice C _Date Written: September 10, 1995 C _Date Modified: September 3, 1996 C _Purpose: C C BCP11 evaluates the gradient of the objective function of the C generated problem at a given point. C C _Arguments in parameter list: C C N (input) INTEGER C The number of variables in the problem. C C X (input) DOUBLE PRECISION vector of dimension N C The point where the gradient of the objective function is to C be evaluated. C C GF (output) DOUBLE PRECISION vector of dimension N C The gradient of objective function values. C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. DOUBLE PRECISION GF(N),X(N) C .. C .. External Subroutines .. EXTERNAL GRADG C .. C .. External Functions .. DOUBLE PRECISION HI1 EXTERNAL HI1 C .. C .. Parameters .. C .. See the comments in BCP01 for a description of the variables C in the common blocks. INTEGER NN PARAMETER (NN=10000) C .. C .. Scalars in Common .. INTEGER ACTL,ACTU,HI,INL,INU,NM1L,NM1U,NM2L,NM2U C .. C .. Arrays in Common .. DOUBLE PRECISION BETA(NN),XBAR(NN) INTEGER PTACT(NN) C .. C .. Local Scalars .. INTEGER I,II C .. C .. Common blocks .. COMMON /BCPB1/INL,INU,ACTL,ACTU,NM1L,NM1U,NM2L,NM2U,HI COMMON /BCPB2/BETA,XBAR,PTACT C .. C Evaluate the unconstrained objective function gradient CALL GRADG(N,X,GF) DO 10 II = 1,ACTL I = PTACT(II) GF(I) = GF(I) + HI1(X(I),XBAR(I),BETA(II),HI) 10 CONTINUE DO 20 II = ACTL + 1,ACTU I = PTACT(II) GF(I) = GF(I) - HI1(X(I),XBAR(I),BETA(II),HI) 20 CONTINUE RETURN END C .. C .. End of subroutine BCP11. C .. C .. Begin of subroutine BCP12 .. SUBROUTINE BCP12(N,X,D,HFD) C _Authors: F. Facchinei, J. Soares and J. Judice C _Date Written: September 10, 1995 C _Date Modified: September 3, 1996 C _Purpose: C C BCP12 evaluates the product between the Hessian matrix of C the objective function at a given point and a vector. C C _Arguments in parameter list: C C N (input) INTEGER C The number of variables in the problem. C C X (input) DOUBLE PRECISION vector of dimension N C The point where the Hessian matrix of the objective function C is to be evaluated. C C D (input) DOUBLE PRECISION vector of dimension N C The vector that will be multiplied by the Hessian matrix. C C HFD (output) DOUBLE PRECISION vector of dimension N C The resulting vector. C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. DOUBLE PRECISION D(N),HFD(N),X(N) C .. C .. External Subroutines .. EXTERNAL HESGD C .. C .. External Functions .. DOUBLE PRECISION HI2 EXTERNAL HI2 C .. C .. Parameters .. C .. See the comments in BCP01 for a description of the variables in C the common blocks. INTEGER NN PARAMETER (NN=10000) C .. C .. Scalars in Common .. INTEGER ACTL,ACTU,HI,INL,INU,NM1L,NM1U,NM2L,NM2U C .. C .. Arrays in Common .. DOUBLE PRECISION BETA(NN),XBAR(NN) INTEGER PTACT(NN) C .. C .. Local Scalars .. INTEGER I,II C .. C .. Common blocks .. COMMON /BCPB1/INL,INU,ACTL,ACTU,NM1L,NM1U,NM2L,NM2U,HI COMMON /BCPB2/BETA,XBAR,PTACT C .. C Evaluate the unconstrained objective function Hessian X d CALL HESGD(N,X,D,HFD) DO 10 II = 1,ACTL I = PTACT(II) HFD(I) = HFD(I) + HI2(X(I),XBAR(I),HI)*D(I) 10 CONTINUE DO 20 II = ACTL + 1,ACTU I = PTACT(II) HFD(I) = HFD(I) - HI2(X(I),XBAR(I),HI)*D(I) 20 CONTINUE RETURN END C .. C .. End of subroutine BCP12 C .. C .. Begin of subroutine RAND REAL FUNCTION RAND(ISEED) C ********** C C function rand C C Rand is the portable random number generator of L. Schrage. C C The generator is full cycle, that is, every integer from C 1 to 2**31 - 2 is generated exactly once in the cycle. C It is completely described in TOMS 5(1979),132-138. C C The function statement is C C real function rand(iseed) C C where C C iseed is a positive integer variable which specifies C the seed to the random number generator. Given the C input seed, rand returns a random number in the C open interval (0,1). On output the seed is updated. C C Argonne National Laboratory. MINPACK Project. March 1981. C C ********** C .. Scalar Arguments .. INTEGER ISEED C .. C .. Local Scalars .. REAL C INTEGER A,B15,B16,FHI,K,LEFTLO,P,XALO,XHI C .. C .. Intrinsic Functions .. INTRINSIC FLOAT C .. C .. Data statements .. C C Set a = 7**5, b15 = 2**15, b16 = 2**16, p = 2**31-1, c = 1/p. C DATA A/16807/,B15/32768/,B16/65536/,P/2147483647/ DATA C/4.656612875E-10/ C .. C C There are 8 steps in rand. C C 1. Get 15 hi order bits of iseed. C 2. Get 16 lo bits of iseed and form lo product. C 3. Get 15 hi order bits of lo product. C 4. Form the 31 highest bits of full product. C 5. Get overflo past 31st bit of full product. C 6. Assemble all the parts and pre-substract p. C The parentheses are essential. C 7. Add p back if necessary. C 8. Multiply by 1/(2**31-1). C XHI = ISEED/B16 XALO = (ISEED-XHI*B16)*A LEFTLO = XALO/B16 FHI = XHI*A + LEFTLO K = FHI/B15 ISEED = (((XALO-LEFTLO*B16)-P)+ (FHI-K*B15)*B16) + K IF (ISEED.LT.0) ISEED = ISEED + P RAND = C*FLOAT(ISEED) RETURN C C Last card of function rand. C END SHAR_EOF fi # end of overwriting check cd .. cd .. # End of shell archive exit 0