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