sed > Makefile <<'#####' 's/^-//' -# Makefile for blp suite - -FFLAGS = -ansi -v -C -g -u - -BLAS = dot.o copy.o scal.o axpy.o nrm2.o - -SUBS = blpgen.o blpmkm.o blpmx.o blpdot.o blpax.o\ - blpfun.o blpchk.o blpflp.o dload.o rand.o - -change: change.o - f77 change.o -o change - -usage: usage.o $(SUBS) $(BLAS) - f77 $(FFLAGS) usage.o $(SUBS) $(BLAS) -o SUN_usage - -verify: verify.o $(SUBS) $(BLAS) - f77 $(FFLAGS) verify.o $(SUBS) $(BLAS) -o SUN_verify - -check: $(SUBS) $(BLAS) - f77 -c $(FFLAGS) $(SUBS) $(BLAS) - -archive: $(SUBS) $(BLAS) - ar rcv blp_archive $(SUBS) $(BLAS) - -clean: - /bin/rm -f usage.o verify.o $(SUBS) $(BLAS) core - - ##### sed > README <<'#####' 's/^-//' -Included with this distribution are the suite programs: - -blpax.f blpflp.f blpmkm.f dload.f -blpchk.f blpfun.f blpmx.f rand.f -blpdot.f blpgen.f - -and utility programs: - -usage.f verify.f change.f - -and the makefile: - -Makefile - -In addition, the following modified level-1 BLAS routines are included -for completeness: - -axpy.f dot.f scal.f -copy.f nrm2.f - -This contribution provides software for generating test problems -for disjointly constrained bilinear programming. -The algorithm constructs problems with a number of favorable properties -that can be selected and controlled by the user. -The intention is to provide a set of FORTRAN 77 routines that can be -used for testing, verifying and comparing solution techniques for -these problems. - -In {Vicente, L.N., Calamai, P.H. and Judice, J.J., "Generation of Disjointly -Constrained Bilinear Programming Test Problems", Computational Optimization -and Applications, (1993)} we describe a technique for generating random -disjointly constrained bilinear programming problems. -Our hope is that these codes will be used for testing, verifying and -comparing new (and established) solution methods for these problems. - - -Implementation and Design. - -This code was implemented entirely in portable ANSI standard FORTRAN 77 with -one exception: both uppercase and lowercase were used to improve readability. - -There are several routines in the suite: - -subroutine blpgen - -This is the central routine of the suite and the one that -controls the properties of the problems that are generated. -See Usage Notes. - -subroutine blpfun - -This routine evaluates the functions and constraints generated -by subroutine blpgen. -See Usage Notes. - -subroutine blpchk - -This service routine checks a candidate solution against all minimizers -of the generated problem. -See Usage Notes. - -subroutine blpax - -This routine is used to evaluate the left-hand-side of the constraints. - -subroutine blpmx - -This routine applies the transformation constructed by -subroutine blpmkm. - -subroutine blpmkm - -This routine constructs the transformation. - -subroutine blpdot - -This routine splits the inner-product of two (potentially) sparse -vectors between two scalars. - -subroutine blpflp - -This routine computes a vector minus a constant times a vector. - -subroutine dload - -This routine copies a scalar into a vector. - -function rand - -This portable random number generator of L. Schrage is completely -described in ACM TOMS 5(1979), 132-138. - -modified level-1 BLAS routines which have been renamed to -axpy, copy, dot, scal and nrm2. - - -Usage Notes. - -Each test problem that is generated (ie. each call to subroutine -blpgen) has properties that depend on the choice of values given to the parameters passed to subroutine blpgen: - -the number of x and y variables is controlled by scalars k1 and k2; - -the solution characteristics are controlled by scalars k11, k12, -k13 and k14 and vectors delta1 and delta3; and - -the characteristics of the transformation are controlled by scalars lcondm -and non0. - - -The influence of these parameters is described in the code comments -and in the cited manuscript. -Example calls to subroutine blpgen appear in program usage -that accompanies this code. -An annotated description of these examples can be obtained by -stripping out the comments in this program. - -Two additional subroutines that may be called by the user -are subroutine blpfun and subroutine blpchk. -Such calls are demonstrated in program usage. - -A call to subroutine blpfun - -subroutine blpfun(v,bilfun,res) - -evaluates the bilinear objective (scalar parameter bilfun) -and the constraint residuals (vector parameter -res) at the point passed (vector parameter v). - -A call to subroutine blpchk - -subroutine blpchk(v,mr,errl2,type) - -compares the candidate solution (passed in vector parameter v) -against all minimizers. -It is important to note that this check is made in the -untransformed space (to do so in the transformed space would be a -combinatorial process) using the point wv (the point in the -untransformed space corresponding to v). -On termination, blpchk returns the vector mr} (the point in the -transformed space corresponding to r, where r is the minimizer -closest to wv), the Euclidean distance between wv and -r (in scalar parameter errl2) and the classification of the -minimizer mr (in scalar parameter type). - - -Installation Notes. - -\subsection{Preliminaries.} - -The precision of the programs in this suite can be changed -in one of two ways. -The easiest approach is to obtain the FORTRAN 77 program -change.f (see {Grcar, J.F., "The change tool for changing programs and -scripts." Sandia Report SAND92-8225 UC-405, (Sept. 1992).) and use it to -manipulate the precision change blocks in all routines in the suite. -This procedure is described in the Sandia report. -Alternatively, a text editor can be used to make a single precision -version of the suite by appending a FORTRAN comment character to the -start of all lines between the precision double change blocks -(i.e., those lines between sequential {precision > double} -and {end precision > double} comments) and deleting -the FORTRAN comment character at the start of all lines between -the precision single change blocks (i.e., those lines between sequential -{precision > single} and {end precision > single} comments). -A double precision version of this suite is obtained, in this fashion, -by interchanging the append and delete operations in the above instructions. - - -Linking. - -This suite makes extensive use of some level-1 BLAS routines. -A modified version of these routines, that incorporate the -precision change blocks described above, have been included with the suite. -Users may wish to strip these modified routines from the suite, -modify the calls to the modified BLAS routines to conform with -the appropriate (i.e., correct precision) BLAS calls, and link -to the appropriate BLAS library. - - -Verification. - -A rudimentary test of the installation can be performed by -compiling program verify that is included with the suite, -linking to the suite and running the resulting module. - - -Limitations. - -The parameter maxk2 controls the size of the largest -problem that can be generated by the suite. -The value of maxk2 is currently fixed at 2000 in -a parameter statement in subroutines blpgen, blpfun,blpax, -blpchk, blpmx and blpmkm and in program verify. -A user wishing to generate larger problems must change -these parameter statements accordingly. ##### sed > axpy.f <<'#####' 's/^-//' - subroutine axpy(n,a,x,incx,y,incy) - -c constant times a vector plus a vector. -c uses unrolled loops for increments equal to one. -c jack dongarra, linpack, 3/11/78. -c modified to include precision changes, paul calamai, 20/10/92. - -c To change the precision of this program, run the change -c program on axpy.f -c Alternatively, to make a single precision version append a -c comment character to the start of all lines between sequential -c precision > double -c and -c end precision > double -c comments and delete the comment character at the start of all -c lines between sequential -c precision > single -c and -c end precision > single -c comments. To make a double precision version interchange -c the append and delete operations in these instructions. - -c*****precision > double - double precision x(*),y(*),a - double precision zero - parameter (zero=0.d0) -c*****end precision > double -c*****precision > single -c real x(*),y(*),a -c real zero -c parameter (zero=0.0) -c*****end precision > single - integer i,incx,incy,ix,iy,m,mp1,n - - if(n.le.0)return - if (a .eq. zero) return - if(incx.eq.1.and.incy.eq.1)go to 20 - -c code for unequal increments or equal increments -c not equal to 1 - - ix = 1 - iy = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 - if(incy.lt.0)iy = (-n+1)*incy + 1 - do 10 i = 1,n - y(iy) = y(iy) + a*x(ix) - ix = ix + incx - iy = iy + incy - 10 continue - return - -c code for both increments equal to 1 - - -c clean-up loop - - 20 m = mod(n,4) - if( m .eq. 0 ) go to 40 - do 30 i = 1,m - y(i) = y(i) + a*x(i) - 30 continue - if( n .lt. 4 ) return - 40 mp1 = m + 1 - do 50 i = mp1,n,4 - y(i) = y(i) + a*x(i) - y(i + 1) = y(i + 1) + a*x(i + 1) - y(i + 2) = y(i + 2) + a*x(i + 2) - y(i + 3) = y(i + 3) + a*x(i + 3) - 50 continue - return - end ##### sed > blpax.f <<'#####' 's/^-//' - subroutine blpax(x,y) - integer maxk2 - parameter (maxk2=2000) - -c ********** -c -c See comments in subroutine blpgen -c -c To change the precision of this program, run the change -c program on blpax.f -c Alternatively, to make a single precision version append a -c comment character to the start of all lines between sequential -c precision > double -c and -c end precision > double -c comments and delete the comment character at the start of all -c lines between sequential -c precision > single -c and -c end precision > single -c comments. To make a double precision version interchange -c the append and delete operations in these instructions. -c -c Authors: Luis N. Vicente and Paul H. Calamai -c Date Written: December, 1992 -c -c ********** - -c declarations for passed parameters - -c*****precision > double - double precision x(*),y(*) -c*****end precision > double -c*****precision > single -c real x(*),y(*) -c*****end precision > single - -c declarations for local parameters - - integer i,istart,iend,ic1,ic2,ic3,ic4,ic5,j4,j3,j2,j1,k,na - -c parameter declarations - -c*****precision > double - double precision minus3,minus2,two,three,threhf,mfihf - parameter (minus3=-3.d0,minus2=-2.d0) - parameter (two=2.d0,three=3.d0,threhf=1.5d0,mfihf=-2.5d0) -c*****end precision > double -c*****precision > single -c real minus3,minus2,two,three,threhf,mfihf -c parameter (minus3=-3.0,minus2=-2.0) -c parameter (two=2.0,three=3.0,threhf=1.5,mfihf=-2.5) -c*****end precision > single - -c declarations for named common /dim/ - - integer ik1,ik2,nx,ny,ik11,ik12,ik13,ik14,n - common/dim/ik1,ik2,nx,ny,ik11,ik12,ik13,ik14,n - -c declarations for named common /param/ - -c*****precision > double - double precision delt1(maxk2),delt3(maxk2) -c*****end precision > double -c*****precision > single -c real delt1(maxk2),delt3(maxk2) -c*****end precision > single - common/param/delt1,delt3 - - istart = 1 - iend = ik11 - na = 3*ik2 - do 10 i=istart,iend - ic2=2*i - ic1=ic2-1 - ic4=ic2+nx - ic3=ic4-1 - j3=3*i - j2=j3-1 - j1=j3-2 - y(j1) = x(ic2) - y(j2) = minus2*x(ic1) - x(ic2) - y(j3) = two*x(ic1) - x(ic2) - y(j1+na) = -delt1(i)*x(ic3) + x(ic4) - y(j2+na) = delt1(i)*x(ic3) + x(ic4) - y(j3+na) = minus2*x(ic4) - 10 continue - - istart = istart + ik11 - iend = iend + ik12 - do 20 i=istart,iend - ic2=2*i - ic1=ic2-1 - ic4=ic2+nx - ic3=ic4-1 - j3=3*i - j2=j3-1 - j1=j3-2 - y(j1) = x(ic2) - y(j2) = minus2*x(ic1) - x(ic2) - y(j3) = two*x(ic1) - x(ic2) - y(j1+na) = minus3*x(ic3) + x(ic4) - y(j2+na) = three*x(ic3) + x(ic4) - y(j3+na) = + minus2*x(ic4) - 20 continue - - istart = istart + ik12 - iend = iend + ik13 - k=1 - do 30 i=istart,iend - ic2=2*i - ic1=ic2-1 - ic4=ic2+nx - ic3=ic4-1 - j3=3*i - j2=j3-1 - j1=j3-2 - y(j1) = x(ic2) - y(j2) = minus2*x(ic1) - x(ic2) - y(j3) = two*x(ic1) - x(ic2) - y(j1+na) = -delt3(k)*x(ic3) + x(ic4) - y(j2+na) = delt3(k)*x(ic3) + x(ic4) - y(j3+na) = minus2*x(ic4) - k=k+1 - 30 continue - - istart = istart + ik13 - iend = ik1 - do 40 i=istart,iend - ic2=2*i - ic1=ic2-1 - ic4=ic2+nx - ic3=ic4-1 - j3=3*i - j2=j3-1 - j1=j3-2 - y(j1) = x(ic2) - y(j2) = minus2*x(ic1) - x(ic2) - y(j3) = two*x(ic1) - x(ic2) - y(j1+na) = mfihf*x(ic3) + x(ic4) - y(j2+na) = x(ic3) + x(ic4) - y(j3+na) = threhf*x(ic3) + minus2*x(ic4) - 40 continue - - istart = 1 - iend = ik2-ik1 - ic5 = 2*(ik1+ik2) - j4 = 3*(ik1+ik2) - do 50 i=istart,iend - ic4=2*(i+ik1) - ic3=ic4-1 - j3=3*(i+ik1) - j2=j3-1 - j1=j3-2 - y(j1) = x(ic4) - y(j2) = minus2*x(ic3) - x(ic4) - y(j3) = two*x(ic3) - x(ic4) - y(j4+i) = x(ic5+i) - 50 continue - - return - -c last card of subroutine blpax - - end ##### sed > blpchk.f <<'#####' 's/^-//' - subroutine blpchk(v,mr,errl2,type) - integer maxk2 - parameter (maxk2=2000) - -c ********** -c -c This subroutine checks Wv (the point corresponding to v -c in the untransformed space) and finds the closest (local -c or global) solution r (in the untransformed space) to -c the problem generated by subroutine blpgen and returns: -c -c errl2 - The l2 distance from Wv to r (and an upper-bound -c on the distance from v to Mr) -c mr - The point corresponding to r in the transformed -c space (ie. Mr) -c type - A classifier for the point mr -c -c value classification -c -c 1 mr is the unique global minimizer of the -c untransformed problem -c 2 mr is a nonunique global minimizer of the -c untransformed problem -c 3 mr is a local minimizer of the untransformed -c problem -c -c To change the precision of this program, run the change -c program on qbpchk.f -c Alternatively, to make a single precision version append a -c comment character to the start of all lines between sequential -c precision > double -c and -c end precision > double -c comments and delete the comment character at the start of all -c lines between sequential -c precision > single -c and -c end precision > single -c comments. To make a double precision version interchange -c the append and delete operations in these instructions. -c -c ******************** -c Subprograms called -c Fortran supplied: abs, sqrt -c Local suite: blpmx -c ******************** -c -c Authors: Luis N. Vicente and Paul H. Calamai -c Date Written: December, 1992 -c -c ******************** - -c declarations for passed parameters - - integer type -c*****precision > double - double precision v(*),mr(*) - double precision errl2 -c*****end precision > double -c*****precision > single -c real v(*),mr(*) -c real errl2 -c*****end precision > single - -c declarations for local parameters - - integer i,istart,iend,ibest,ic1,ic2,ic3,ic4,ic5 -c*****precision > double - double precision errto1,errto2,errto3,errto4,errbes - double precision x1to,x2to,y1to,y2to - double precision wv(4*maxk2) -c*****end precision > double -c*****precision > single -c real errto1,errto2,errto3,errto4,errbes -c real x1to,x2to,y1to,y2to -c real wv(4*maxk2) -c*****end precision > single - -c function declarations - - intrinsic abs,sqrt - -c parameter declarations - -c*****precision > double - double precision minus3,minus2,minus1,zero,one,two,three, - 1 fivehf,mfihf - parameter (minus3=-3.d0,minus2=-2.d0,minus1=-1.d0,zero=0.d0) - parameter (one=1.d0,two=2.d0,three=3.d0, - 1 fivehf=2.5d0,mfihf=-2.5d0) -c*****end precision > double -c*****precision > single -c real minus3,minus2,minus1,zero,one,two,three, -c 1 fivehf,mfihf -c parameter (minus3=-3.0,minus2=-2.0,minus1=-1.0,zero=0.0) -c parameter (one=1.0,two=2.0,three=3.0, -c 1 fivehf=2.5,mfihf=-2.5) -c*****end precision > single - -c declarations for named common /obj/ - -c*****precision > double - double precision cd(4*maxk2),ab(6*maxk2) -c*****end precision > double -c*****precision > single -c real cd(4*maxk2),ab(6*maxk2) -c*****end precision > single - common/obj/cd,ab - -c declarations for named common /param/ - -c*****precision > double - double precision delt1(maxk2),delt3(maxk2) -c*****end precision > double -c*****precision > single -c real delt1(maxk2),delt3(maxk2) -c*****end precision > single - common/param/delt1,delt3 - -c declarations for named common /dim/ - - integer ik1,ik2,nx,ny,ik11,ik12,ik13,ik14,n - common/dim/ik1,ik2,nx,ny,ik11,ik12,ik13,ik14,n - -c obtain Wv - - call blpmx(v,wv,.true.) - -c find the closest solution (in the untransformed space) -c to Wv (ie. find r) - - type=1 - - if (ik11.gt.0 .or. ik12.gt.0) then - type = 2 - endif - - errl2=zero - -c check the components 1..ik11+ik12+ik13 - - istart = 1 - iend = ik11+ik12+ik13 - if (iend.gt.0) then - do 10 i=istart,iend - - ic2=2*i - ic1=ic2-1 - ic4=ic2+nx - ic3=ic4-1 - - x1to = abs(wv(ic1)) - x2to = abs(wv(ic2)+minus2) - y1to = abs(wv(ic3)+minus2) - y2to = abs(wv(ic4)) - errto1 = x1to**2+y1to**2+x2to**2+y2to**2 - ibest = 1 - errbes = errto1 - - x1to = abs(wv(ic1)+minus2) - x2to = abs(wv(ic2)+minus2) - y1to = abs(wv(ic3)) - y2to = abs(wv(ic4)) - errto2 = x1to**2+y1to**2+x2to**2+y2to**2 - if (errto2 .lt. errbes) then - errbes = errto2 - ibest = 2 - endif - - x1to = abs(wv(ic1)+minus1) - x2to = abs(wv(ic2)) - y1to = abs(wv(ic3)+minus1) - if (i.le.ik11) then - y2to = abs(wv(ic4)-delt1(i)) - endif - if (ik11.lt.i .and. i.le.(ik11+ik12)) then - y2to = abs(wv(ic4)+minus3) - endif - if ((ik11+ik12).lt.i .and. i.le.(ik11+ik12+ik13)) then - y2to = abs(wv(ic4)-delt3(i-ik11-ik12)) - endif - errto3 = x1to**2+y1to**2+x2to**2+y2to**2 - if (errto3 .lt. errbes) then - errbes = errto3 - ibest = 3 - endif - - x1to = abs(wv(ic1)+minus1) - x2to = abs(wv(ic2)+minus2) - y1to = abs(wv(ic3)+minus1) - y2to = abs(wv(ic4)) - errto4 = x1to**2+y1to**2+x2to**2+y2to**2 - if (errto4 .lt. errbes) then - errbes = errto4 - ibest = 4 - endif - - if (ibest.eq.1) then - if ((ik11+ik12).lt.i .and. i.le.(ik11+ik12+ik13)) then - type = 3 - endif - errl2 = errl2 + errto1 - mr(ic1) = zero - mr(ic2) = two - mr(ic3) = two - mr(ic4) = zero - endif - if (ibest.eq.2) then - if ((ik11+ik12).lt.i .and. i.le.(ik11+ik12+ik13)) then - type = 3 - endif - errl2 = errl2 + errto2 - mr(ic1) = two - mr(ic2) = two - mr(ic3) = zero - mr(ic4) = zero - endif - if (ibest.eq.3) then - errl2 = errl2 + errto3 - mr(ic1) = one - mr(ic2) = zero - mr(ic3) = one - if (i.le.ik11) then - type = 3 - mr(ic4) = delt1(i) - endif - if (ik11.lt.i .and. i.le.(ik11+ik12)) then - mr(ic4) = three - endif - if ((ik11+ik12).lt.i .and. i.le.(ik11+ik12+ik13)) then - mr(ic4) = delt3(i-ik11-ik12) - endif - endif - if (ibest.eq.4) then - type = 3 - errl2 = errl2 + errto4 - mr(ic1) = one - mr(ic2) = two - mr(ic3) = one - mr(ic4) = zero - endif - - 10 continue - endif - -c check the components ik11+ik12+ik13+1...ik1 - - if (ik14.gt.0) then - istart = iend+1 - iend = ik1 - do 20 i=istart,iend - - ic2=2*i - ic1=ic2-1 - ic4=ic2+nx - ic3=ic4-1 - - x1to = abs(wv(ic1)+minus2) - x2to = abs(wv(ic2)+minus2) - y1to = abs(wv(ic3)) - y2to = abs(wv(ic4)) - errto1 = x1to**2+y1to**2+x2to**2+y2to**2 - ibest = 1 - errbes = errto1 - - x1to = abs(wv(ic1)+minus1) - x2to = abs(wv(ic2)) - y1to = abs(wv(ic3)+minus1) - y2to = abs(wv(ic4)+mfihf) - errto2 = x1to**2+y1to**2+x2to**2+y2to**2 - if (errto2 .lt. errbes) then - errbes = errto2 - ibest = 2 - endif - - if (ibest.eq.1) then - errl2 = errl2 + errto1 - mr(ic1) = two - mr(ic2) = two - mr(ic3) = zero - mr(ic4) = zero - endif - if (ibest.eq.2) then - type = 3 - errl2 = errl2 + errto2 - mr(ic1) = one - mr(ic2) = zero - mr(ic3) = one - mr(ic4) = fivehf - endif - - 20 continue - endif - -c check the components ik1+1...ik2 - - istart = ik1+1 - iend = ik2 - if (istart.le.iend) then - do 30 i=istart,iend - - ic2=2*i - ic1=ic2-1 - ic5=ik1+nx+i - - x1to = abs(wv(ic1)+minus1) - x2to = abs(wv(ic2)) - y1to = abs(wv(ic5)+minus2) - errto1 = x1to**2+y1to**2+x2to**2 - errl2 = errl2 + errto1 - mr(ic1) = one - mr(ic2) = zero - mr(ic5) = two - - 30 continue - endif - - errl2=sqrt(errl2) - -c transform mr to get point Mr in transformed space - - call blpmx(mr,mr,.false.) - -c end of subroutine blpchk - - end ##### sed > blpdot.f <<'#####' 's/^-//' - subroutine blpdot(non0,x,y,ind,produ,prodl) -c -c ******************** -c -c splits the sparse dot product of x and y between -c produ and prodl -c -c To change the precision of this program, run the change -c program on blpdot.f -c Alternatively, to make a single precision version append a -c comment character to the start of all lines between sequential -c precision > double -c and -c end precision > double -c comments and delete the comment character at the start of all -c lines between sequential -c precision > single -c and -c end precision > single -c comments. To make a double precision version interchange -c the append and delete operations in these instructions. -c -c ******************** -c -c Authors: Luis N. Vicente and Paul H. Calamai -c Date Written: December, 1992 -c -c ******************** - -c*****precision > double - double precision x(*),y(*),produ,prodl,zero -c*****end precision > double -c*****precision > single -c real x(*),y(*),produ,prodl,zero -c*****end precision > single - integer i,ind(*),ixy,non0 - - integer ik1,ik2,nx,ny,ik11,ik12,ik13,ik14,n - common/dim/ik1,ik2,nx,ny,ik11,ik12,ik13,ik14,n - -c*****precision > double - parameter (zero=0.d0) -c*****end precision > double -c*****precision > single -c parameter (zero=0.0) -c*****end precision > single - - produ = zero - prodl = zero - if(n.le.0)return - do 10 i = 1,non0 - ixy = ind(i) - if (ixy.le.nx) then - produ=produ + x(ixy)*y(ixy) - else - prodl=prodl + x(ixy)*y(ixy) - endif - 10 continue - return - -c Last card of subroutine blpdot - - end ##### sed > blpflp.f <<'#####' 's/^-//' - subroutine blpflp(non0,da,db,dx,dy,inc) - -c ********* -c -c vector minus constant times a vector. -c -c To change the precision of this program, run the change -c program on blpflp.f -c Alternatively, to make a single precision version append a -c comment character to the start of all lines between sequential -c precision > double -c and -c end precision > double -c comments and delete the comment character at the start of all -c lines between sequential -c precision > single -c and -c end precision > single -c comments. To make a double precision version interchange -c the append and delete operations in these instructions. -c -c ********** -c -c Authors: Luis N. Vicente and Paul H. Calamai -c Date Written: December, 1992 -c -c ********** - -c*****precision > double - double precision dx(*),dy(*),da,db -c*****end precision > double -c*****precision > single -c real dx(*),dy(*),da,db -c*****end precision > single - integer i,inc(*),ixy,non0 - - integer ik1,ik2,nx,ny,ik11,ik12,ik13,ik14,n - common/dim/ik1,ik2,nx,ny,ik11,ik12,ik13,ik14,n - -c*****precision > double - double precision zero - parameter (zero=0.d0) -c*****end precision > double -c*****precision > single -c real zero -c parameter (zero=0.0) -c*****end precision > single - - if(non0.le.0)return - if (da .eq. zero .and. db .eq. zero) return - - do 10 i = 1,non0 - ixy=inc(i) - if (ixy.le.nx) then - dy(ixy) = dy(ixy) - da*dx(ixy) - else - dy(ixy) = dy(ixy) - db*dx(ixy) - endif - 10 continue - return - -c last card of subroutine blpflp - - end ##### sed > blpfun.f <<'#####' 's/^-//' - subroutine blpfun(v,bilfun,res) - integer maxk2 - parameter (maxk2=2000) - -c ********** -c -c This subroutine computes the value of the bilinear -c objective function and the residuals of the constraints -c of the bilinear disjointly constrained -c programming problem generated using subroutine -c blpgen -c -c Details of this problem generation can be found -c in the comments of subroutine blpgen -c -c To change the precision of this program, run the change -c program on blpfun.f -c Alternatively, to make a single precision version append a -c comment character to the start of all lines between sequential -c precision > double -c and -c end precision > double -c comments and delete the comment character at the start of all -c lines between sequential -c precision > single -c and -c end precision > single -c comments. To make a double precision version interchange -c the append and delete operations in these instructions. -c -c ********** -c Subprograms called -c Local suite: blpax,blpmx -c ********** -c -c Authors: Luis N. Vicente and Paul H. Calamai -c Date Written: December, 1992 -c -c ********** - -c declarations for passed parameters - -c*****precision > double - double precision v(*),res(*),bilfun -c*****end precision > double -c*****precision > single -c real v(*),res(*),bilfun -c*****end precision > single - -c declarations for local parameters - -c*****precision > double - double precision mv(4*maxk2) -c*****end precision > double -c*****precision > single -c real mv(4*maxk2) -c*****end precision > single - integer i,ipnx,ipik1,j - -c declarations for named common /m/ - -c*****precision > double - double precision d(4*maxk2),z(4*maxk2) -c*****end precision > double -c*****precision > single -c real d(4*maxk2),z(4*maxk2) -c*****end precision > single - integer num0 - integer pz(4*maxk2) - common/m/d,z,pz,num0 - -c declarations for named common /obj/ - -c*****precision > double - double precision cd(4*maxk2),ab(6*maxk2) -c*****end precision > double -c*****precision > single -c real cd(4*maxk2),ab(6*maxk2) -c*****end precision > single - common/obj/cd,ab - -c declarations for named common /dim/ - - integer ik1,ik2,nx,ny,ik11,ik12,ik13,ik14,n - common/dim/ik1,ik2,nx,ny,ik11,ik12,ik13,ik14,n - -c store Mv in mv - - call blpmx(v,mv,.false.) - -c compute the objective function - -c all the components of the first ik1 kernel problems - - do 10 i=1,2*ik1 - ipnx=i+nx - bilfun = bilfun + mv(i)*(mv(ipnx) + cd(i)) + - 1 mv(ipnx)*cd(ipnx) - 10 continue - -c all the components of the remmaining ik2-ik1 kernel problems - - do 20 i=1,ik2-ik1 - ipik1=i+ik1 - j=2*ipik1 - bilfun = bilfun + mv(2*ny+i)*(mv(j-1) + mv(j) + cd(2*ny+i)) + - 1 mv(j-1)*cd(j-1) + mv(j)*cd(j) - - 20 continue - -c compute the residuals - - call blpax(mv,res) - do 30 i=1,4*ik2+2*ik1 - res(i) = ab(i) - res(i) - 30 continue - return - -c last card of subroutine blpfun - - end ##### sed > blpgen.f <<'#####' 's/^-//' - subroutine blpgen(delta1,delta3,k1,k2,k11,k12,k13,k14,lcondm,non0, - 1 info,iseed) - integer maxk2 - parameter (maxk2=2000) - -c ********** -c -c To change the precision of this program, run the change -c program on blpgen.f -c Alternatively, to make a single precision version append a -c comment character to the start of all lines between sequential -c precision > double -c and -c end precision > double -c comments and delete the comment character at the start of all -c lines between sequential -c precision > single -c and -c end precision > single -c comments. To make a double precision version interchange -c the append and delete operations in these instructions. -c -c This routine generates random disjointly constrained bilinear -c programs of the form: -c -c Minimize f(x,y) = (c'Mx)x + x'(Mx Q My)y + (d'My)y -c -c subject to (A Mx)x <= a , (B My)y <= b -c -c -c where c = (c' ... c' ) , d = (d' ... d' ) , -c 1 k2 1 k2 -c -c Q = diag (Q ... Q ) , -c 1 k2 -c -c A = diag (A ... A ) , B = diag (B ... B ) , -c 1 k2 1 k2 -c -c a = (a' ... a' ) and b = (b' ... b' ) -c 1 k2 1 k2 -c -c and, for k = 1 , ... , k1, -c -c ( 1 0 ) -c c' = d' = (-1, -1) , Q = ( ) , -c k k k ( 0 1 ) -c -c -c ( 0 1 ) ( s1(k) 1 ) -c A = ( -2 -1 ) , B = ( s2(k) 1 ) , -c k ( 2 -1 ) k ( s3(k) -2 ) -c -c -c a' = (2, -2, 2) , b' = (0, s4(k), 0) -c k k -c -c where -c -c k1 = abs(k11) + k12 + abs(k13) + k14, -c -c s1' = (-delta1', -3*1(k12)', -delta3', -(5/2)*1(k14)') -c -c s2' = (delta1', 3*1(k12)', delta3', 1(k14)') -c -c s3' = (0(abs(k11)+k12+abs(k13)), (3/2)*1(k14)) -c -c s4' = (2*delta1, 6*1(k12)', 2*delta3', (7/2)*1(k14)) -c -c and, for k = k1 + 1 , ... , k2, -c -c ( 1 ) -c c' = (-1, -1) , d' = -2 , Q = ( ) , -c k k k ( 1 ) -c -c -c ( 0 1 ) -c A = ( -2 -1 ) , B = ( 1 ) , -c k ( 2 -1 ) k -c -c -c a' = (2, -2, 2) , b' = 2. -c k k -c -c -c Here -c -c k1 and k2 are nonnegative integers satisfying -c 0 <= k1 <= k2 and k2 > 0 -c nx = 2*k2 is the number of x variables -c ny = k1 + k2 is the number of y variables -c n is the total number of variables (i.e. n = nx + ny) -c -c 0(m) is the zero vector of order m -c 1(m) is the ones vector of order m -c -c Mx are My are (sparse) positive definite matrices -c defined by Mx = Dx Hx and My = Dy Hy -c -c Hx and Hy are order nx and ny Householder matrices generated -c by the (Householder) unit 2-norm vectors zx and zy -c respectively, z'=(zx' zy') has non0 nonzero components in (-1,1) -c in rows pz(1) through pz(non0), and D=diag(d(1) ,..., d(n)) -c is a positive definite diagonal matrix with 2-norm condition -c 10**lcondm, where D = diag (Dx, Dy) -c -c H = diag (Hx, Hy) , M = diag (Mx, My) -c W is the inverse of M -c -c k11 is as described below -c k12 is as described below -c k13 is as described below -c k14 is as described below -c delta1 is as described below -c delta3 is as described below -c -c This problem has -c -c 4^{abs(k11)+k12+abs(k13)} * 2^{k14} -c -c local solutions including -c -c 2^{abs(k11)} * 3^{k12} -c -c global solutions. -c -c The user is referred to: -c Vicente, L.N., Calamai, P.H. and Judice, J.J., -c "Generation of Disjointly Constrained Bilinear Programming -c Test Problems", -c Computational Optimization and Applications, (1993) -c for a description of how this is accomplished. -c -c ********** -c -c INTEGER INPUT PARAMETERS (unchanged on exit) -c -c -c k1 and k2 - 2*k2 is the number of x variables and k1 + k2 is the -c number of y variables -c -c k1 and k2 should verify -c 0 <= k1 <= k2 and k2 > 0 -c -c k11 - in the untransformed space the components -c (x(2*i-1),x(2*i),y(2*i-1),y(2*i)) -c have their global minimizers at -c (0,2,2,0), (2,2,0,0) -c and their local minimizers at -c (1,0,1,delta1(i)), (1,2,1,0) -c for i=1,...,ik11 where ik11 = abs(k11) -c -c k12 - in the untransformed space the components -c (x(2*i-1),x(2*i),y(2*i-1),y(2*i)) -c have their global minimizers at -c (0,2,2,0), (2,2,0,0), (1,0,1,3) -c and their local minimizer at -c (1,2,1,0) -c for i=ik11+1,...,ik11+k12 -c -c k13 - in the untransformed space the components -c (x(2*i-1),x(2*i),y(2*i-1),y(2*i)) -c have their global minimizer at -c (1,0,1,delta3(i-ik11-k12)) -c and their local minimizers at -c (0,2,2,0), (2,2,0,0), (1,2,1,0) -c for i=ik11+k12+1,...,ik11+k12+ik13 where ik13 = abs(k13) -c -c k14 - in the untransformed space the components -c (x(2*i-1),x(2*i),y(2*i-1),y(2*i)) -c have their global minimizer at -c (2,2,0,0) -c and their local minimizer at -c (1,0,1,5/2) -c for i=ik11+k12+ik13+1,...,k1 -c -c ***************************************************************** -c NOTE - ik11 + k12 + ik13 + k14 must equal k1 in this version -c ***************************************************************** -c -c lcondm - log (base 10) of the condition number of M (lcondm >= 0) -c -c non0 - (maximum) number of nonzeros in z (0 < non0 <= n) -c -c INTEGER INPUT PARAMETERS (changed on exit) -c -c iseed - see comments in subprogram rand -c -c REAL INPUT PARAMETERS -c -c delta1 - real vector of length ik11 that controls the minima -c corresponding to components 1 through 2*ik11 of x -c and y as described above under k11. -c The use of delta1 depends on the sign of k11 as follows: -c -c value use -c -c k11 < 0 input values are ignored and the ik11 components of -c delta1 are generated randomly to fall in the open -c interval (1,3) -c -c k11 > 0 if -3 < delta1(1) < -1 then the ik11 components of -c delta1 are all set to -delta1(1); otherwise, the -c values that are input are used and these must all -c lie in the open interval (1,3) -c -c delta3 - real vector of length ik13 that controls the minima -c corresponding to components 2*(ik11+k12)+1 through 2*(ik11+ -c k12+ik13) of x and y as described above under k13. -c The use of delta3 depends on the sign of k13 as follows: -c -c value use -c -c k13 < 0 input values are ignored and the ik13 components of -c delta3 are generated randomly to fall in the open -c interval (3,delta3(1)). (NOTE: delta3(1) must be -c greater than 3 when k13 < 0) -c -c k13 > 0 if delta3(1) < -3 then the ik13 components of -c delta3 are all set to -delta3(1); otherwise, the -c values that are input are used and these must all -c be greater than 3 -c -c DIAGNOSTIC PARAMETERS -c -c info - on exit contains a value indicating the state of program -c blpgen when it terminated -c -c value definition -c -c 0 program terminated normally -c -c >0 the state is given by that subset of the following -c conditions whose value sum to the value of info: -c -c value condition -c -c 1 k1 and/or k2 outside permissible range -c 2 lcondm outside permissible range -c 4 non0 outside permissible range -c 8 k11 and/or k12 and/or k13 and/or k14 outside -c permissible range -c 16 components of delta1 outside permissible range -c 32 components of delta3 outside permissible range -c -c PARAMETERS OF COMMON /M/ -c -c d - the diagonals of D -c -c z - the Householder vectors that generate Hx and Hy -c -c pz - row indices of (possible) nonzeros in z -c -c num0 - same as non0 -c -c PARAMETERS OF COMMON /OBJ/ -c -c cd - the coefficients of the linear term of f(x,y) -c -c ab - the upper-bounds on the linear constraints -c -c PARAMETERS OF COMMON /PARAM/ -c -c delt1 - same as delta1 -c -c delt3 - same as delta3 -c -c PARAMETERS OF COMMON /DIM/ -c -c ik1 - same as k1 -c -c ik2 - same as k2 -c -c dimx - same as nx -c -c dimy - same as ny -c -c ik11 - abs(k11) -c -c ik12 - same as k12 -c -c ik13 - abs(k13) -c -c ik14 - same as k14 -c -c n - nx+ny -c -c ******************** -c Subprograms called -c Modified Linpack: copy,scal -c TOMS (MINPACK project) supplied: rand -c Fortran supplied: abs -c Local suite: blpmkm,dload -c ******************** -c -c Authors: Luis N. Vicente and Paul H. Calamai -c Date Written: December, 1992 -c -c ******************** - -c declarations for passed parameters - - integer k1,k2,k11,k12,k13,k14,lcondm,non0,info,iseed -c*****precision > double - double precision delta1(*),delta3(*) -c*****end precision > double -c*****precision > single -c real delta1(*),delta3(*) -c*****end precision > single - -c declarations for local parameters - - logical fault - integer i,k1sum -c*****precision > double - double precision delta -c*****end precision > double -c*****precision > single -c real delta -c*****end precision > single - -c function declarations - - real rand - -c parameter declarations - -c*****precision > double - double precision minus3,minus2,minus1,zero,one,two,three, - 1 six,sevehf - parameter (minus3=-3.d0,minus2=-2.d0,minus1=-1.d0) - parameter (zero=0.d0,one=1.d0,two=2.d0,three=3.d0, - 1 six=6.d0,sevehf=3.5d0) -c*****end precision > double -c*****precision > single -c real minus3,minus2,minus1,zero,one,two,three, -c 1 six,sevehf -c parameter (minus3=-3.0,minus2=-2.0,minus1=-1.0) -c parameter (zero=0.0,one=1.0,two=2.0,three=3.0, -c 1 six=6.0,sevehf=3.5) -c*****end precision > single - -c declarations for named common /m/ - - integer num0 - integer pz(4*maxk2) -c*****precision > double - double precision d(4*maxk2),z(4*maxk2) -c*****end precision > double -c*****precision > single -c real d(4*maxk2),z(4*maxk2) -c*****end precision > single - common/m/d,z,pz,num0 - -c declarations for named common /obj/ - -c*****precision > double - double precision cd(4*maxk2),ab(6*maxk2) -c*****end precision > double -c*****precision > single -c real cd(4*maxk2),ab(6*maxk2) -c*****end precision > single - common/obj/cd,ab - -c declarations for named common /param/ - -c*****precision > double - double precision delt1(maxk2),delt3(maxk2) -c*****end precision > double -c*****precision > single -c real delt1(maxk2),delt3(maxk2) -c*****end precision > single - common/param/delt1,delt3 - -c declarations for named common /dim/ - - integer ik1,ik2,dimx,dimy,ik11,ik12,ik13,ik14,n - common/dim/ik1,ik2,dimx,dimy,ik11,ik12,ik13,ik14,n - -c store information for commons - - ik1 = k1 - ik2 = k2 - dimx = 2*k2 - dimy = k1+k2 - ik11 = abs(k11) - ik12 = k12 - ik13 = abs(k13) - ik14 = k14 - n = dimx+dimy - num0 = non0 - - do 10 i=1,ik11 - delt1(i) = delta1(i) -10 continue - - do 20 i=1,ik13 - delt3(i) = delta3(i) -20 continue - -c check validity of integer input parameters and set info - - k1sum = ik11+k12+ik13+k14 - info = 0 - if (k1.lt.0 .or. k2.lt.1 .or. k2.gt.maxk2 .or. k1.gt.k2) then - info=1 - endif - if (lcondm.lt.0) then - info=info+2 - endif - if (non0.lt.1 .or. non0.gt.(dimx+dimy)) then - info=info+4 - endif - if (k12.lt.0 .or. k14.lt.0 .or. k1sum.ne.k1) then - info=info+8 - endif - if (info.ne.0) return - -c check and adjust delt1 when appropriate - - if (k11.ne.0) then - if (k11.gt.0) then - if (minus3.lt.delt1(1) .and. delt1(1).lt.minus1) then - delta=-delt1(1) - call dload(ik11,delta,delt1,1) - else - fault=.false. - do 30 i=1,ik11 - if (delt1(i).lt.one .or. delt1(i).gt.three) fault=.true. - 30 continue - if (fault) info=info+16 - endif - else - do 40 i=1,ik11 - delt1(i)=two*rand(iseed)+one - 40 continue - endif - endif - -c check and adjust delt3 when appropriate - - if (k13.ne.0) then - if (k13.gt.0) then - if (delt3(1).lt.minus3) then - delta=-delt3(1) - call dload(ik13,delta,delt3,1) - else - fault=.false. - do 50 i=1,ik13 - if (delt3(i).lt.three) fault=.true. - 50 continue - if (fault) info=info+32 - endif - else - delta=delt3(1) - if (delta.le.three) then - info=info+32 - else - do 60 i=1,ik13 - delt3(i)=(delta-three)*rand(iseed)+three - 60 continue - endif - endif - endif - if (info.ne.0) return - -c store data for constructing M - - call blpmkm(lcondm,iseed) - -c construct cd - - call dload(2*ik1+2*ik2,minus1,cd,1) - call dload(ik2-ik1,minus2,cd(2*ik1+2*ik2+1),1) - -c construct ab - - call dload(k2,two,ab(1),3) - call dload(k2,minus2,ab(2),3) - call dload(k2,two,ab(3),3) - - call dload(ik11,zero,ab(3*k2+1),3) - call copy(ik11,delt1,1,ab(3*k2+2),3) - call scal(ik11,two,ab(3*k2+2),3) - call dload(ik11,zero,ab(3*k2+3),3) - - call dload(k12,zero,ab(3*k2+3*ik11+1),3) - call dload(k12,six,ab(3*k2+3*ik11+2),3) - call dload(k12,zero,ab(3*k2+3*ik11+3),3) - - call dload(ik13,zero,ab(3*k2+3*ik11+3*k12+1),3) - call copy(ik13,delt3,1,ab(3*k2+3*ik11+3*k12+2),3) - call scal(ik13,two,ab(3*k2+3*ik11+3*k12+2),3) - call dload(ik13,zero,ab(3*k2+3*ik11+3*k12+3),3) - - call dload(k14,zero,ab(3*k2+3*ik11+3*k12+3*ik13+1),3) - call dload(k14,sevehf,ab(3*k2+3*ik11+3*k12+3*ik13+2),3) - call dload(k14,zero,ab(3*k2+3*ik11+3*k12+3*ik13+3),3) - - call dload(k2-k1,two,ab(3*k2+3*k1+1),1) - -c end of subroutine blpgen - - end ##### sed > blpmkm.f <<'#####' 's/^-//' - subroutine blpmkm(condm,iseed) - integer maxk2 - parameter (maxk2=2000) - -c ********** -c -c This subroutine defines the n by n matrices M = HDH, -c where H is a block-diagonal matrix defined by H=diag(Hx,Hy), -c Hx and Hy are order nx and ny Householder matrices generated -c by the (Householder) unit 2-norm vectors zx and zy -c respectively, z'=[zx' zy'] has non0 nonzero components from (-1,1) -c in rows pz(1) through pz(non0), and D=diag(d(1) ,..., d(n)) -c is a positive definite diagonal matrix with 2-norm condition -c 10**lcondm -c -c To change the precision of this program, run the change -c program on blpmkm.f -c Alternatively, to make a single precision version append a -c comment character to the start of all lines between sequential -c precision > double -c and -c end precision > double -c comments and delete the comment character at the start of all -c lines between sequential -c precision > single -c and -c end precision > single -c comments. To make a double precision version interchange -c the append and delete operations in these instructions. -c -c ********** -c Subprograms called -c Modified Linpack: nrm2,scal -c TOMS (MINPACK project) supplied: rand -c Fortran supplied: int,dble -c Local suite: dload -c ********** -c -c Authors: Luis N. Vicente and Paul H. Calamai -c Date Written: December, 1992 -c -c ********** - -c declarations for passed parameters - - integer condm,iseed - -c declarations for local parameters - - integer indx,k -c*****precision > double - double precision znorm,tempz,factor -c*****end precision > double -c*****precision > single -c real znorm,tempz,factor -c*****end precision > single - -c function declarations - -c*****precision > double - double precision nrm2 -c*****end precision > double -c*****precision > single -c real nrm2 -c*****end precision > single - real rand - intrinsic int - -c parameter declarations - -c*****precision > double - double precision zero,one - parameter (zero=0.d0,one=1.d0) -c*****end precision > double -c*****precision > single -c real zero,one -c parameter (zero=0.0,one=1.0) -c*****end precision > single - -c declarations for common /m/ - -c*****precision > double - double precision d(4*maxk2),z(4*maxk2) -c*****end precision > double -c*****precision > single -c real d(4*maxk2),z(4*maxk2) -c*****end precision > single - integer pz(4*maxk2) - integer num0 - common/m/d,z,pz,num0 - -c declarations for common /dim/ - - integer ik1,ik2,nx,ny,ik11,ik12,ik13,ik14,n - common/dim/ik1,ik2,nx,ny,ik11,ik12,ik13,ik14,n - -c put num0 random elements in (-1,1) into random locations -c of z and store locations in pz(1) through pz(num0) - - call dload(n,zero,z,1) - do 20 k = 1, num0 - tempz = 2*rand(iseed) - one - 10 continue - indx=int(n*rand(iseed)+1) - if (indx.gt.n) go to 10 - if (z(indx).ne.zero) go to 10 - z(indx)=tempz - pz(k)=indx - 20 continue - -c make zx and zy unit 2-norm vectors - - znorm = nrm2(nx,z,1) - if (znorm.ne.zero) then - call scal(nx,one/znorm,z,1) - endif - znorm = nrm2(ny,z(nx+1),1) - if (znorm.ne.zero) then - call scal(ny,one/znorm,z(nx+1),1) - endif - -c store diagonal of D - - do 30 k = 1,n -c*****precision > double - factor = (dble(k-1)/dble(n-1)) -c*****end precision > double -c*****precision > single -c factor = (real(k-1)/real(n-1)) -c*****end precision > single - d(k) = 10**(-condm*factor) - 30 continue - return - -c Last card of subroutine blpmkm - - end ##### sed > blpmx.f <<'#####' 's/^-//' - subroutine blpmx(x,y,invert) - integer maxk2 - parameter (maxk2=2000) - -c ********** -c -c This subroutine forms the product Mx (or M**(-1)x when invert -c is true), where M is the transformation constructed in -c subroutine blpmkm (see comments in subroutine blpmkm and in -c subroutine blpgen). -c -c To change the precision of this program, run the change -c program on blpmx.f -c Alternatively, to make a single precision version append a -c comment character to the start of all lines between sequential -c precision > double -c and -c end precision > double -c comments and delete the comment character at the start of all -c lines between sequential -c precision > single -c and -c end precision > single -c comments. To make a double precision version interchange -c the append and delete operations in these instructions. -c -c ********** -c Subprograms called -c Modified Linpack: copy -c Local suite: blpdot, blpflp -c ********** -c -c Authors: Luis N. Vicente and Paul H. Calamai -c Date Written: December, 1992 -c -c ********** - -c declarations for passed parameters - - logical invert -c*****precision > double - double precision x(*),y(*) -c*****end precision > double -c*****precision > single -c real x(*),y(*) -c*****end precision > single - -c declarations for local parameters - - integer indx,k -c*****precision > double - double precision produ,prodl -c*****end precision > double -c*****precision > single -c real produ,prodl -c*****end precision > single - -c parameter declarations - -c*****precision > double - double precision two - parameter (two=2.d0) -c*****end precision > double -c*****precision > single -c real two -c parameter (two=2.0) -c*****end precision > single - -c declarations for named common /m/ - -c*****precision > double - double precision d(4*maxk2),z(4*maxk2) -c*****end precision > double -c*****precision > single -c real d(4*maxk2),z(4*maxk2) -c*****end precision > single - integer pz(4*maxk2) - integer num0 - common/m/d,z,pz,num0 - -c declarations for named common /dim/ - - integer ik1,ik2,nx,ny,ik11,ik12,ik13,ik14,n - common/dim/ik1,ik2,nx,ny,ik11,ik12,ik13,ik14,n - -c store Householder inner-products in produ and prodl - - call blpdot(num0,z,x,pz,produ,prodl) - produ=two*produ - prodl=two*prodl - -c store Householder reflections in y - - call copy(n,x,1,y,1) - call blpflp(num0,produ,prodl,z,y,indx) - -c store Dy (or D**(-1)y when invert) in y - - if (invert) then - do 30 k=1,n - y(k)=y(k)/d(k) - 30 continue - else - do 40 k=1,n - y(k)=y(k)*d(k) - 40 continue - endif - -c store Householder inner-products in produ and prodl - - call blpdot(num0,z,y,pz,produ,prodl) - produ=two*produ - prodl=two*prodl - -c store Householder reflections in y - - call blpflp(num0,produ,prodl,z,y,indx) - return - -c last card of subroutine blpmx - - end ##### sed > change.f <<'#####' 's/^-//' -C/////////////////////////////////////////////////////////////////////// -C -C C H A N G E -C -C VERSION 2.02 OF AUGUST 1992 -C -C COMMENT-OUT OR UN-COMMENT-IN LINES OF A FORTRAN PROGRAM OR A UNIX -C SHELL SCRIPT -C -C WRITTEN BY: DR. JOSEPH F. GRCAR -C SANDIA NATIONAL LABORATORIES -C DEPARTMENT 8745 -C LIVERMORE, CA 94551-0969 USA -C -C (510) 294-2662 -C -C na.grcar@na-net.ornl.gov -C sepp@snll-arpagw.llnl.gov -C -C/////////////////////////////////////////////////////////////////////// -C -C DOCUMENTATION: -C -C J. F. Grcar, "The Change Tool for Changing Programs and Scripts," -C Sandia National Laboratories Report SAND92-8225, Livermore, -C California, September 1992. -C -C/////////////////////////////////////////////////////////////////////// -C -C CHANGES FROM THE PREVIOUS VERSION: -C -C 1) ALTER ERROR MESSAGES. -C -C 2) REMOVE EXTERNAL STATEMENT GIVING TROUBLE TO SOME COMPILERS. -C -C/////////////////////////////////////////////////////////////////////// - - PROGRAM CHANGE - - IMPLICIT COMPLEX (A - P, R - Z), INTEGER (Q) - CHARACTER - + FILE1*80, FILE2*80, ITEM*8, LIST*80, MARK*8, REMARK*6, SMARK*8 -C THE MICROWAY NDP FORTRAN COMPILER VERSION 4.0.2 DOES PASSES -C CHARACTER STRINGS INCORRECTLY WITH THE FOLLOWING STATEMENT. -C EXTERNAL EXTENT, SQUEEZ - INTEGER - + CMND, CMNTS, CODE, COUNT, INDEX, ITEMP, J, K, LEN1, LEN2, LEN3, - + LEN4, LEN5, LENGTH, LINES, LONG, MAXL, MAXN, MAXS, NAMES, - + NUMBER, REPORT, ROUTE, SAVNUM, STATUS, TOPNUM, TYPE, UNIT1, - + UNIT2 - LOGICAL - + ACTIVE, ERROR, EXIST, FLAG, FOUND, INSIDE, MANUAL, MAY, NEW, - + OPEN, SPECL, TOOBIG - - PARAMETER (REPORT = 1000, MAXN = 100) - PARAMETER (QADD = 1, QCOPY = 2, QSTRIP = 3) - PARAMETER (QFORT = 1, QSHELL = 2) - - PARAMETER (MAXS = 200) -C THE FOLLOWING CHARACTER STRINGS SHOULD HAVE LENGTH EQUAL TO MAXS - CHARACTER - + LINE*200, NAME*200, SAVLIN*200, TEMP*200, TEMP2*200, TEMP3*200, - + TOP*200 - - DIMENSION - + ACTIVE(MAXN, 2), CMND(0 : MAXN), CODE(MAXN), LONG(MAXN), - + NAME(MAXN), REMARK(MAXN, 2), SMARK(MAXN), TYPE(0 : MAXN) - -C FORCE STATIC ALLOCATION OF VARIABLES. THIS IS NOT NEEDED BY THE -C CHANGE PROGRAM, BUT IT SEEMS TO HELP WITH TROUBLESOME COMPILERS. - SAVE - -C/////////////////////////////////////////////////////////////////////// -C -C PROLOGUE. -C -C/////////////////////////////////////////////////////////////////////// - -C/// CHOOSE THE UNIT NUMBERS. - - UNIT1 = 20 - UNIT2 = 21 - -C/// DECIDE WHETHER THE COPY MAY OVERWRITE THE ORIGINAL. - - MAY = .FALSE. -C*****COPY MAY OVERWRITE ORIGINAL (RECOMMENDED) - MAY = .TRUE. -C*****END COPY MAY OVERWRITE ORIGINAL (RECOMMENDED) - -C/// DECIDE WHETHER TO LIMIT LINES TO 72 CHARACTERS. - - MAXL = MAXS -C*****LINE LENGTH RESTRICTED TO 72 CHARACTERS - MAXL = 72 -C*****END LINE LENGTH RESTRICTED TO 72 CHARACTERS - -C/// WRITE THE HEADER. - - WRITE (*, 10001) - -C/// PRINT ERROR MESSAGES FOR THE USER MANUAL. - - MANUAL = .FALSE. -C MANUAL = .TRUE. - - IF (MANUAL) THEN - LINE = 'C*****ABRAHAM LINCOLN' - NAME(1) = 'GEORGE WASHINGTON' - NUMBER = 8701 - SAVLIN = ' + (/13X, ''ERROR. IN THE ACTIVE CHANGE BLOCK ' - + // 'WHICH BEGINS ON LINE '', A, '',''' - SAVNUM = 4528 - STATUS = 43 - TOP = 'C*****GEORGE WASHINGTON' - TOPNUM = 4521 - INDEX = 1 - NAMES = 1 - GO TO 9101 - END IF - -C/////////////////////////////////////////////////////////////////////// -C -C (1) OPEN THE ORIGINAL FILE. -C -C/////////////////////////////////////////////////////////////////////// - - ASSIGN 1010 TO ROUTE -1010 CONTINUE - -C/// PROMPT FOR THE FILE NAME. - - WRITE (*, 10002) - READ - + (END = 99999, - + ERR = 9101, - + FMT = '(A)', - + IOSTAT = STATUS, - + UNIT = *) - + FILE1 - - CALL SQUEEZ (LENGTH, FILE1) - IF (FILE1 .EQ. ' ') GO TO 99999 - -C/// CHECK THE FILE. - - ERROR = .TRUE. - INQUIRE - + (ERR = 9102, - + EXIST = EXIST, - + FILE = FILE1, - + IOSTAT = STATUS, - + OPENED = OPEN) - ERROR = .FALSE. - - ERROR = .NOT. EXIST - IF (ERROR) GO TO 9103 - - ERROR = EXIST .AND. OPEN - IF (ERROR) GO TO 9104 - -C/// OPEN THE FILE. - - ERROR = .TRUE. - OPEN - + (ACCESS = 'SEQUENTIAL', - + ERR = 9105, - + FILE = FILE1, - + FORM = 'FORMATTED', - + IOSTAT = STATUS, - + STATUS = 'OLD', - + UNIT = UNIT1) - ERROR = .FALSE. - -C/////////////////////////////////////////////////////////////////////// -C -C (2) LOOP THROUGH THE ORIGINAL FILE. -C -C/////////////////////////////////////////////////////////////////////// - - INSIDE = .FALSE. - MARK = ' ' - NUMBER = 0 - NAMES = 0 - -C/// TOP OF THE LOOP. READ THE NEXT LINE. - -2010 CONTINUE - READ - + (END = 2050, - + ERR = 9201, - + FMT = '(A)', - + IOSTAT = STATUS, - + UNIT = UNIT1) - + LINE - NUMBER = NUMBER + 1 - -C/// GATHER INFORMATION ABOUT THE LINE. - - LENGTH = 0 - DO 2020 J = MAXS, 1, - 1 - IF (LINE (J : J) .NE. ' ') THEN - LENGTH = J - GO TO 2030 - END IF -2020 CONTINUE -2030 CONTINUE - ERROR = MAXL .LT. LENGTH - IF (ERROR) GO TO 9202 - - SPECL = LINE (1 : 6) .EQ. 'C*****' - + .OR. LINE (1 : 6) .EQ. 'c*****' - + .OR. LINE (1 : 6) .EQ. '******' - + .OR. LINE (1 : 6) .EQ. '#*****' - -C/// TREAT SPECIAL LINES FOUND FROM OUTSIDE. - - IF (SPECL) THEN - IF (LINE (1 : 1) .EQ. '#') THEN - TYPE(0) = QSHELL - ELSE - TYPE(0) = QFORT - END IF - - IF (.NOT. INSIDE) THEN - FOUND = .FALSE. - DO 2040 J = 1, NAMES - IF (TYPE(0) .EQ. TYPE(J) .AND. - + LINE (7 : ) .EQ. NAME(J)) THEN - FOUND = .TRUE. - INDEX = J - END IF -2040 CONTINUE - - IF (.NOT. FOUND) THEN - NAMES = NAMES + 1 - ERROR = NAMES .GT. MAXN - IF (ERROR) GO TO 9203 - - NEW = .TRUE. - INDEX = NAMES - CODE(INDEX) = INDEX - NAME(INDEX) = LINE (7 : ) - TYPE(INDEX) = TYPE(0) - END IF - - CMNTS = 0 - INSIDE = .TRUE. - LINES = 0 - MARK = LINE (1 : 1) - TOOBIG = .FALSE. - TOP = LINE - TOPNUM = NUMBER - -C/// TREAT SPECIAL LINES FOUND FROM INSIDE. - - ELSE - ERROR = .NOT. ( - + (LINE (7 : 7) .EQ. 'E' .OR. LINE (7 : 7) .EQ. 'e') .AND. - + (LINE (8 : 8) .EQ. 'N' .OR. LINE (8 : 8) .EQ. 'n') .AND. - + (LINE (9 : 9) .EQ. 'D' .OR. LINE (9 : 9) .EQ. 'd') .AND. - + LINE (10 : 10) .EQ. ' ') - IF (ERROR) GO TO 9204 - - ERROR = .NOT. (NAME(INDEX) .EQ. LINE (11 : )) - IF (ERROR) GO TO 9205 - - ERROR = .NOT. (TOP (1 : 1) .EQ. LINE (1 : 1)) - IF (ERROR) GO TO 9206 - - ERROR = LINES .EQ. 0 - IF (ERROR) GO TO 9207 - - FLAG = CMNTS .LT. LINES - IF (NEW) THEN - ACTIVE(INDEX, 1) = FLAG - NEW = .FALSE. - ELSE - ERROR = ACTIVE(INDEX, 1) .NEQV. FLAG - IF (ERROR) GO TO 9208 - END IF - - ERROR = ACTIVE(INDEX, 1) .AND. TOOBIG - IF (ERROR) GO TO 9209 - - INDEX = 0 - INSIDE = .FALSE. - END IF - -C/// TREAT SIMPLE LINES FROM OUTSIDE. - - ELSE - IF (.NOT. INSIDE) THEN - -C/// TREAT SIMPLE LINES FROM INSIDE. - - ELSE - LINES = LINES + 1 - - IF (TYPE (INDEX) .EQ. QFORT) THEN - IF (LINE (1 : 1) .EQ. 'C' .OR. LINE (1 : 1) .EQ. 'c' - + .OR. LINE (1 : 1) .EQ. '*') CMNTS = CMNTS + 1 - ELSE IF (TYPE (INDEX) .EQ. QSHELL) THEN - IF (LINE (1 : 1) .EQ. '#') CMNTS = CMNTS + 1 - END IF - - IF (.NOT. TOOBIG .AND. MAXL .LE. LENGTH) THEN - SAVLIN = LINE - SAVNUM = NUMBER - TOOBIG = .TRUE. - END IF - END IF - END IF - -C/// BOTTOM OF THE LOOP THROUGH THE FILE. - - IF (MOD (NUMBER, REPORT) .EQ. 0) THEN - IF (NUMBER .EQ. REPORT) WRITE (*, 10004) - WRITE (TEMP, 10005) NUMBER - CALL SQUEEZ (LENGTH, TEMP) - WRITE (*, 10006) TEMP (1 : LENGTH) - END IF - - GO TO 2010 -2050 CONTINUE - - IF (MOD (NUMBER, REPORT) .NE. 0) THEN - IF (NUMBER .LT. REPORT) WRITE (*, 10004) - WRITE (TEMP, 10005) NUMBER - CALL SQUEEZ (LENGTH, TEMP) - WRITE (*, 10006) TEMP (1 : LENGTH) - END IF - - ERROR = INSIDE - IF (ERROR) GO TO 9210 - -C/////////////////////////////////////////////////////////////////////// -C -C (3) DETERMINE WHICH BLOCKS SHOULD BE ACTIVE. -C -C/////////////////////////////////////////////////////////////////////// - -C/// DECIDE WHETHER TO WRITE A COPY. - - IF (NAMES .EQ. 0) THEN - WRITE (*, 10010) - GO TO 99999 - END IF - -C/// SORT THE CHANGE BLOCKS. - - DO 3010 INDEX = 1, NAMES - IF (TYPE(INDEX) .EQ. QFORT) THEN - SMARK(INDEX) = ')' - ELSE - SMARK(INDEX) = '#' - END IF -3010 CONTINUE - - DO 3020 K = 1, NAMES - DO 3020 J = K + 1, NAMES - IF (NAME(J) .LT. NAME(K)) THEN - ITEMP = CODE(J) - CODE(J) = CODE(K) - CODE(K) = ITEMP - - TEMP = SMARK(J) - SMARK(J) = SMARK(K) - SMARK(K) = TEMP - - TEMP = NAME(J) - NAME(J) = NAME(K) - NAME(K) = TEMP - - FLAG = ACTIVE(J, 1) - ACTIVE(J, 1) = ACTIVE(K, 1) - ACTIVE(K, 1) = FLAG - END IF -3020 CONTINUE - -C/// INITIALIZE THE CHOICES AND LENGTHS. - - DO 3030 J = 1, NAMES - ACTIVE(J, 2) = ACTIVE(J, 1) - CALL EXTENT (LONG(J), NAME(J)) -3030 CONTINUE - -C/// TOP OF THE QUERY LOOP. - - LIST = ' ' - WRITE (*, 10007) -3040 CONTINUE - -C/// DISPLAY THE NAMES OF THE CHANGE BLOCKS. ASK WHICH TO CHANGE. - - IF (LIST .EQ. ' ') THEN - DO 3050 K = 1, 2 - DO 3050 J = 1, NAMES - IF (ACTIVE(J, K)) THEN - REMARK(J, K) = 'ACTIVE' - ELSE - REMARK(J, K) = ' -- ' - END IF -3050 CONTINUE - - WRITE (*, 10008) (J, SMARK(J), REMARK(J, 1), REMARK(J, 2), - + NAME(J) (1 : LONG(J)), J = 1, NAMES) - - WRITE (*, 10009) - READ - + (END = 3090, - + ERR = 9101, - + FMT = '(A)', - + IOSTAT = STATUS, - + UNIT = *) - + LIST - IF (LIST .EQ. ' ') GO TO 3090 - END IF - -C/// GET THE NUMBER OF THE NEXT BLOCK TO CHANGE. - - CALL SQUEEZ (LENGTH, LIST) - COUNT = 1 - DO 3060 J = 1, LEN (LIST) - IF (LIST (J : J) .EQ. ' ') GO TO 3070 - COUNT = J -3060 CONTINUE -3070 CONTINUE - - IF (COUNT .LE. 8) THEN - ITEM = ' ' - ITEM (9 - COUNT : 8) = LIST (1 : COUNT) - LIST (1 : COUNT) = ' ' - READ - + (ERR = 3040, - + FMT = '(I8)', - + IOSTAT = STATUS, - + UNIT = ITEM) - + J - ELSE - LIST (1 : COUNT) = ' ' - GO TO 3040 - END IF - -C/// CHANGE THE BLOCK. - - IF (1 .LE. J .AND. J .LE. NAMES) ACTIVE(J, 2) = .NOT. ACTIVE(J, 2) - -C/// BOTTOM OF THE QUERY LOOP. - - GO TO 3040 -3090 CONTINUE - -C/// DETERMINE THE ACTION TO BE TAKEN ON EACH BLOCK. - - DO 3100 J = 1, NAMES - IF (ACTIVE(J, 1) .EQV. ACTIVE(J, 2)) THEN - CMND(J) = QCOPY - ELSE IF (.NOT. ACTIVE(J, 1) .AND. ACTIVE(J, 2)) THEN - CMND(J) = QSTRIP - ELSE IF (ACTIVE(J, 1) .AND. .NOT. ACTIVE(J, 2)) THEN - CMND(J) = QADD - END IF -3100 CONTINUE - -C/// DECIDE WHETHER TO WRITE A COPY. - - FLAG = .FALSE. - DO 3110 J = 1, NAMES - FLAG = FLAG .OR. (CMND(J) .NE. QCOPY) -3110 CONTINUE - - IF (.NOT. FLAG) THEN - WRITE (*, 10011) - GO TO 99999 - END IF - -C/////////////////////////////////////////////////////////////////////// -C -C (4) OPEN THE COPY FILE. -C -C/////////////////////////////////////////////////////////////////////// - -C/// TOP OF THE BLOCK. - - ASSIGN 4010 TO ROUTE -4010 CONTINUE - -C/// PROMPT FOR THE FILE NAME. - - CALL SQUEEZ (LENGTH, FILE1) - WRITE (*, 10003) FILE1 (1 : LENGTH) - READ - + (END = 99999, - + ERR = 9101, - + FMT = '(A)', - + IOSTAT = STATUS, - + UNIT = *) - + FILE2 - - CALL SQUEEZ (LENGTH, FILE2) - IF (FILE2 .EQ. ' ') GO TO 99999 - -C/// USE THE ORIGINAL FILE. - - IF (FILE1 .EQ. FILE2) THEN - IF (.NOT. MAY) GO TO 9401 - -C/// OPEN THE SCRATCH FILE. - - ERROR = .TRUE. - OPEN - + (ACCESS = 'SEQUENTIAL', - + ERR = 9402, - + FORM = 'FORMATTED', - + IOSTAT = STATUS, - + STATUS = 'SCRATCH', - + UNIT = UNIT2) - ERROR = .FALSE. - -C/// COPY FROM THE ORIGINAL FILE TO THE SCRATCH FILE. - - NUMBER = 0 - REWIND UNIT1 -4020 CONTINUE - READ - + (END = 4050, - + ERR = 9201, - + FMT = '(A)', - + IOSTAT = STATUS, - + UNIT = UNIT1) - + LINE - NUMBER = NUMBER + 1 - - LENGTH = 1 - DO 4030 J = MAXS, 1, - 1 - IF (LINE (J : J) .NE. ' ') THEN - LENGTH = J - GO TO 4040 - END IF -4030 CONTINUE -4040 CONTINUE - ERROR = MAXL .LT. LENGTH - IF (ERROR) GO TO 9202 - - WRITE (UNIT2, '(A)') LINE (1 : LENGTH) - GO TO 4020 -4050 CONTINUE - -C/// SWITCH THE UNIT NUMBERS. - - ITEMP = UNIT1 - UNIT1 = UNIT2 - UNIT2 = ITEMP - -C/// REWIND THE FILES. - - REWIND UNIT1 - REWIND UNIT2 - -C/// USE A NEW FILE. - - ELSE - -C/// CHECK THE FILE. - - ERROR = .TRUE. - INQUIRE - + (ERR = 9102, - + EXIST = EXIST, - + FILE = FILE2, - + IOSTAT = STATUS, - + OPENED = OPEN) - ERROR = .FALSE. - - ERROR = EXIST .AND. OPEN - IF (ERROR) GO TO 9104 - -C/// OPEN THE FILE. - - IF (EXIST) THEN - ERROR = .TRUE. - OPEN - + (ACCESS = 'SEQUENTIAL', - + ERR = 9105, - + FILE = FILE2, - + FORM = 'FORMATTED', - + IOSTAT = STATUS, - + STATUS = 'OLD', - + UNIT = UNIT2) - ERROR = .FALSE. - -C/// CREATE THE FILE. - - ELSE - ERROR = .TRUE. - OPEN - + (ACCESS = 'SEQUENTIAL', - + ERR = 9107, - + FILE = FILE2, - + FORM = 'FORMATTED', - + IOSTAT = STATUS, - + STATUS = 'NEW', - + UNIT = UNIT2) - ERROR = .FALSE. - END IF - -C/// REWIND THE FILES. - - REWIND UNIT1 - END IF - -C/////////////////////////////////////////////////////////////////////// -C -C (5) MAKE THE COPY. -C -C/////////////////////////////////////////////////////////////////////// - - NUMBER = 0 - -C/// TOP OF THE LOOP. READ THE NEXT LINE. - -5010 CONTINUE - READ - + (END = 5050, - + ERR = 9201, - + FMT = '(A)', - + IOSTAT = STATUS, - + UNIT = UNIT1) - + LINE - NUMBER = NUMBER + 1 - -C/// GATHER INFORMATION ABOUT THE LINE. - - LENGTH = 0 - DO 5020 J = MAXS, 1, - 1 - IF (LINE (J : J) .NE. ' ') THEN - LENGTH = J - GO TO 5030 - END IF -5020 CONTINUE -5030 CONTINUE - ERROR = MAXL .LT. LENGTH - IF (ERROR) GO TO 9202 - - SPECL = LINE (1 : 6) .EQ. 'C*****' - + .OR. LINE (1 : 6) .EQ. 'c*****' - + .OR. LINE (1 : 6) .EQ. '******' - + .OR. LINE (1 : 6) .EQ. '#*****' - -C/// TREAT SPECIAL LINES FOUND FROM OUTSIDE. - - IF (SPECL) THEN - IF (LINE (1 : 1) .EQ. '#') THEN - TYPE(0) = QSHELL - ELSE - TYPE(0) = QFORT - END IF - - IF (.NOT. INSIDE) THEN - FOUND = .FALSE. - DO 5040 J = 1, NAMES - IF (TYPE(0) .EQ. TYPE(J) .AND. - + LINE (7 : ) .EQ. NAME(J)) THEN - FOUND = .TRUE. - INDEX = J - END IF -5040 CONTINUE - ERROR = .NOT. FOUND - IF (ERROR) GO TO 9501 - - INSIDE = .TRUE. - MARK = LINE (1 : 1) - QCMND = QCOPY - -C/// TREAT SPECIAL LINES FOUND FROM INSIDE. - - ELSE - INDEX = 0 - INSIDE = .FALSE. - QCMND = QCOPY - END IF - -C/// TREAT SIMPLE LINES FROM OUTSIDE. - - ELSE - IF (.NOT. INSIDE) THEN - QCMND = QCOPY - -C/// TREAT SIMPLE LINES FROM INSIDE. - - ELSE - LINES = LINES + 1 - QCMND = CMND(INDEX) - END IF - END IF - -C/// WRITE THE LINE. - - IF (QCMND .EQ. QADD) THEN - IF (0 .LT. LENGTH) THEN - TEMP = MARK (1 : 1) // LINE (1 : LENGTH) - ELSE - TEMP = MARK (1 : 1) - END IF - LINE = TEMP - LENGTH = LENGTH + 1 - ELSE IF (QCMND .EQ. QSTRIP) THEN - LENGTH = LENGTH - 1 - TEMP = LINE (2 : ) - LINE = TEMP - END IF - - IF (1 .LE. LENGTH) THEN - WRITE (UNIT2, '(A)') LINE (1 : LENGTH) - ELSE - WRITE (UNIT2, '()') - END IF - -C/// BOTTOM OF THE LOOP THROUGH THE FILE. - - IF (MOD (NUMBER, REPORT) .EQ. 0) THEN - IF (NUMBER .EQ. REPORT) WRITE (*, 10004) - WRITE (TEMP, 10005) NUMBER - CALL SQUEEZ (LENGTH, TEMP) - WRITE (*, 10006) TEMP (1 : LENGTH) - END IF - - GO TO 5010 -5050 CONTINUE - - IF (MOD (NUMBER, REPORT) .NE. 0) THEN - IF (NUMBER .LT. REPORT) WRITE (*, 10004) - WRITE (TEMP, 10005) NUMBER - CALL SQUEEZ (LENGTH, TEMP) - WRITE (*, 10006) TEMP (1 : LENGTH) - END IF - -C/////////////////////////////////////////////////////////////////////// -C -C INFORMATIVE MESSAGES. -C -C/////////////////////////////////////////////////////////////////////// - - GO TO 99999 - -10001 FORMAT - + (/ ' CHANGE: VERSION 2.02 OF AUGUST 1992', - + ' BY JOSEPH GRCAR.') - -C*****$ EDIT DESCRIPTOR > YES (RECOMMENDED) -10002 FORMAT - + (/13X, 'ENTER THE NAME OF THE ORIGINAL FILE.' - + //1X, ' FILE? ', $) - -10003 FORMAT - + (/13X, 'ENTER THE NAME FOR THE COPY,' - + /13X, 'OR ENTER BLANK TO QUIT.' - + //13X, A, ' (NAME OF THE ORIGINAL FILE)' - + //1X, ' FILE? ', $) -C*****END $ EDIT DESCRIPTOR > YES (RECOMMENDED) - -C*****$ EDIT DESCRIPTOR > NO (FORTRAN 77 STANDARD) -C10002 FORMAT -C + (/13X, 'ENTER THE NAME OF THE ORIGINAL FILE.' -C + /) -C -C10003 FORMAT -C + (/13X, 'ENTER THE NAME FOR THE COPY,' -C + /13X, 'OR ENTER BLANK TO QUIT.' -C + //1X, A, ' (NAME OF THE ORIGINAL FILE)' -C + /) -C*****END $ EDIT DESCRIPTOR > NO (FORTRAN 77 STANDARD) - -10004 FORMAT - + () - -10005 FORMAT - + (I10, ' LINES') - -10006 FORMAT - + (13X, A) - -10007 FORMAT - + (/13X, 'STATUS OF THE CHANGE BLOCKS.') - -10008 FORMAT - + (/12X, 'ORIGINAL COPY NAME OF THE CHANGE BLOCK' - + /12X, ' ------ ------ --------------------------' - + /(1X, I9, A1, ' ', A6, 2X, A6, 2X, A)) - -C*****$ EDIT DESCRIPTOR > YES (RECOMMENDED) -10009 FORMAT - + (/13X, 'ENTER THE NUMBERS OF BLOCKS TO CHANGE,' - + /13X, 'OR ENTER BLANK TO ACCEPT THINGS AS THEY ARE.' - + //1X, ' NUMBER? ', $) -C*****END $ EDIT DESCRIPTOR > YES (RECOMMENDED) - -C*****$ EDIT DESCRIPTOR > NO (FORTRAN 77 STANDARD) -C10009 FORMAT -C + (/13X, 'ENTER THE NUMBERS OF BLOCKS TO CHANGE,' -C + /13X, 'OR ENTER BLANK TO ACCEPT THINGS AS THEY ARE.' -C + /) -C*****END $ EDIT DESCRIPTOR > NO (FORTRAN 77 STANDARD) - -10010 FORMAT - + (/13X, 'NO CHANGE BLOCKS.') - -10011 FORMAT - + (/13X, 'NO CHANGES.') - -C*****PAUSE AT FINISH > NO (RECOMMENDED) -10012 FORMAT - + (/13X, 'FINISH.' - + /) -C*****END PAUSE AT FINISH > NO (RECOMMENDED) - -C*****PAUSE AT FINISH > YES (SOME WINDOW ENVIRONMENTS) -C10012 FORMAT -C + (/13X, 'PRESS RETURN TO FINISH.') -C*****END PAUSE AT FINISH > YES (SOME WINDOW ENVIRONMENTS) - -C/////////////////////////////////////////////////////////////////////// -C -C ERROR MESSAGES. -C -C/////////////////////////////////////////////////////////////////////// - -9101 WRITE (TEMP, '(I20)') STATUS - CALL SQUEEZ (LENGTH, TEMP) - WRITE (*, 99101) TEMP (1 : LENGTH) - IF (.NOT. MANUAL) GO TO 99999 - -9102 WRITE (TEMP, '(I20)') STATUS - CALL SQUEEZ (LENGTH, TEMP) - WRITE (*, 99102) TEMP (1 : LENGTH) - IF (.NOT. MANUAL) GO TO 99999 - -9103 WRITE (*, 99103) - IF (.NOT. MANUAL) GO TO ROUTE - -9104 WRITE (*, 99104) - IF (.NOT. MANUAL) GO TO ROUTE - -9105 WRITE (TEMP, '(I20)') STATUS - CALL SQUEEZ (LENGTH, TEMP) - WRITE (*, 99105) TEMP (1 : LENGTH) - IF (.NOT. MANUAL) GO TO 99999 - -9107 WRITE (TEMP, '(I20)') STATUS - CALL SQUEEZ (LENGTH, TEMP) - WRITE (*, 99107) TEMP (1 : LENGTH) - IF (.NOT. MANUAL) GO TO 99999 - -9201 IF (NUMBER .EQ. 1) THEN - WRITE (TEMP, '(I20, A)') NUMBER, ' LINE' - ELSE - WRITE (TEMP, '(I20, A)') NUMBER, ' LINES' - END IF - CALL SQUEEZ (LEN1, TEMP) - WRITE (TEMP2, '(I20)') STATUS - CALL SQUEEZ (LEN2, TEMP2) - WRITE (*, 99201) TEMP (1 : LEN1), TEMP2 (1 : LEN2) - IF (.NOT. MANUAL) GO TO 99999 - -9202 WRITE (TEMP, '(I20)') NUMBER - WRITE (TEMP2, '(I20)') MAXL - CALL SQUEEZ (LEN1, TEMP) - CALL SQUEEZ (LEN2, TEMP2) - WRITE (*, 99202) TEMP (1 : LEN1), TEMP2 (1 : LEN2) - IF (.NOT. MANUAL) GO TO 99999 - -9203 WRITE (TEMP, '(I20)') MAXN - CALL SQUEEZ (LENGTH, TEMP) - WRITE (*, 99203) TEMP (1 : LENGTH) - IF (.NOT. MANUAL) GO TO 99999 - -9204 CALL EXTENT (LEN1, TOP) - CALL EXTENT (LEN2, LINE) - WRITE (*, 99204) - + TOPNUM, TOP (1 : LEN1), NUMBER, LINE (1 : LEN2) - IF (.NOT. MANUAL) GO TO 99999 - -9205 CALL EXTENT (LEN1, TOP) - CALL EXTENT (LEN2, LINE) - WRITE (*, 99205) - + TOPNUM, TOP (1 : LEN1), NUMBER, LINE (1 : LEN2) - IF (.NOT. MANUAL) GO TO 99999 - -9206 CALL EXTENT (LEN1, TOP) - CALL EXTENT (LEN2, LINE) - WRITE (*, 99206) - + TOPNUM, TOP (1 : LEN1), NUMBER, LINE (1 : LEN2) - IF (.NOT. MANUAL) GO TO 99999 - -9207 CALL EXTENT (LEN1, TOP) - CALL EXTENT (LEN2, LINE) - WRITE (*, 99207) - + TOPNUM, TOP (1 : LEN1), NUMBER, LINE (1 : LEN2) - IF (.NOT. MANUAL) GO TO 99999 - -9208 CALL SQUEEZ (LEN1, TOP) - IF (FLAG) THEN - WRITE (*, 99208) 'ACTIVE', 'INACTIVE', TOPNUM, TOP (1 : LEN1) - ELSE - WRITE (*, 99208) 'INACTIVE', 'ACTIVE', TOPNUM, TOP (1 : LEN1) - END IF - IF (.NOT. MANUAL) GO TO 99999 - -9209 WRITE (TEMP, '(I20)') TOPNUM - CALL SQUEEZ (LEN1, TEMP) - WRITE (TEMP2, '(I20)') SAVNUM - CALL SQUEEZ (LEN2, TEMP2) - WRITE (TEMP3, '(I20)') MAXL - 1 - CALL SQUEEZ (LEN3, TEMP3) - CALL EXTENT (LEN4, TOP) - SAVLIN (58 : MAXS) = ' ' - CALL EXTENT (LEN5, SAVLIN) - WRITE (*, 99209) - + TEMP (1 : LEN1), TEMP2 (1 : LEN2), TEMP3 (1 : LEN3), - + TOPNUM, TOP (1 : LEN4), SAVNUM, SAVLIN (1 : LEN5) // '...' - IF (.NOT. MANUAL) GO TO 99999 - -9210 WRITE (*, 99210) - IF (.NOT. MANUAL) GO TO 99999 - -9401 WRITE (*, 99401) - IF (.NOT. MANUAL) GO TO 4010 - -9402 WRITE (TEMP, '(I20)') STATUS - CALL SQUEEZ (LENGTH, TEMP) - WRITE (*, 99402) TEMP (1 : LENGTH) - IF (.NOT. MANUAL) GO TO 99999 - -9501 WRITE (*, 99501) - IF (.NOT. MANUAL) GO TO 5010 - -99101 FORMAT - + (/13X, 'ERROR. READ FAILS WITH I/O STATUS ', A, '.') - -99102 FORMAT - + (/13X, 'ERROR. INQUIRE FAILS WITH I/O STATUS ', A, '.') - -99103 FORMAT - + (/13X, 'ERROR. THE FILE DOES NOT EXIST.') - -99104 FORMAT - + (/13X, 'ERROR. THE FILE IS OPEN TO SOMEONE ELSE.') - -99105 FORMAT - + (/13X, 'ERROR. OPEN FAILS WITH I/O STATUS ', A, '.') - -99107 FORMAT - + (/13X, 'ERROR. FILE CREATE FAILS WITH I/O STATUS ', A, '.') - -99201 FORMAT - + (/13X, 'ERROR. AFTER ', A, ' SUCCESSFULLY READ,' - + /21X, 'A FILE READ FAILS WITH I/O STATUS ', A, '.') - -99202 FORMAT - + (/13X, 'ERROR. LINE NUMBER ', A, ' HAS OVER ', A, - + ' CHARACTERS.') - -99203 FORMAT - + (/13X, 'ERROR. THERE ARE OVER ', A, ' CHANGE BLOCK NAMES.') - -99204 FORMAT - + (/13X, 'ERROR. THE LINE AT THE BOTTOM OF A CHANGE BLOCK' - + /21X, 'DOES NOT SAY "END".' - + //1X, I9, ': ', A - + /1X, I9, ': ', A) - -99205 FORMAT - + (/13X, 'ERROR. THE LINES AT THE TOP AND BOTTOM OF A CHANGE' - + /21X, 'BLOCK HAVE DIFFERENT NAMES.' - + //1X, I9, ': ', A - + /1X, I9, ': ', A) - -99206 FORMAT - + (/13X, 'ERROR. THE LINES AT THE TOP AND BOTTOM OF A CHANGE' - + /21X, 'BLOCK HAVE DIFFERENT COMMENT MARKS.' - + //1X, I9, ': ', A - + /1X, I9, ': ', A) - -99207 FORMAT - + (/13X, 'ERROR. A CHANGE BLOCK HAS NO LINES BETWEEN THE TOP AND' - + /21X, 'BOTTOM LINE.' - + //1X, I9, ': ', A - + /1X, I9, ': ', A) - -99208 FORMAT - + (/13X, 'ERROR. THE CHANGE BLOCKS ARE INCONSISTENT. AMONG' - + /21X, 'BLOCKS WITH THE NAME BELOW, THE FIRST IS ', A, ',' - + /21X, 'BUT THE BLOCK BEGINNING ON THE LINE BELOW IS ', A, '.' - + //1X, I9, ': ', A) - -99209 FORMAT - + (/13X, 'ERROR. IN THE ACTIVE CHANGE BLOCK WHICH BEGINS ON ', - + 'LINE ', A, ',' - + /21X, 'LINE NUMBER ', A, ' IS LONGER THAN ', A, ' CHARACTERS.' - + //1X, I9, ': ', A - + /1X, I9, ': ', A) - -99210 FORMAT - + (/13X, 'ERROR. THE FILE ENDS IN THE MIDDLE OF A CHANGE', - + ' BLOCK.') - -99401 FORMAT - + (/13X, 'ERROR. THE COPY MAY NOT OVERWRITE THE ORIGINAL FILE.') - -99402 FORMAT - + (/13X, 'ERROR. SCRATCH OPEN FAILS WITH I/O STATUS ', A, '.') - -99501 FORMAT - + (/13X, 'ERROR. A CHANGE BLOCK NAME IS NOT RECOGNIZED.') - -C/// EXIT. - -99999 CONTINUE - WRITE (*, 10012) -C*****PAUSE AT FINISH > YES (SOME WINDOW ENVIRONMENTS) -C PAUSE ' ' -C*****END PAUSE AT FINISH > YES (SOME WINDOW ENVIRONMENTS) - END - - SUBROUTINE EXTENT (LENGTH, STRING) - -C/////////////////////////////////////////////////////////////////////// -C -C T O O L S -C -C EXTENT -C -C RETURN THE LENGTH OF A STRING. THAT IS, RETURN THE POSITION OF -C THE LAST NONBLANK CHARACTER, OR 1 IF ALL ARE BLANK. -C -C/////////////////////////////////////////////////////////////////////// - - IMPLICIT COMPLEX (A - Z) - CHARACTER STRING*(*) - INTEGER J, LENGTH - -C FORCE STATIC ALLOCATION OF VARIABLES. THIS IS NOT NEEDED BY THE -C CHANGE PROGRAM, BUT IT SEEMS TO HELP WITH TROUBLESOME COMPILERS. - SAVE - - LENGTH = 1 - DO 0100 J = 1, LEN (STRING) - IF (STRING(J : J) .NE. ' ') LENGTH = J -0100 CONTINUE - - RETURN - END - - SUBROUTINE SQUEEZ (LENGTH, STRING) - -C/////////////////////////////////////////////////////////////////////// -C -C T O O L S -C -C SQUEEZ -C -C SQUEEZE LEADING BLANKS AND MULTIPLE BLANKS FROM A CHARACTER -C STRING. RETURN THE LENGTH OF THE STRING (OR 1 IF ALL BLANK). -C -C/////////////////////////////////////////////////////////////////////// - - IMPLICIT COMPLEX (A - P, R - Z), INTEGER (Q) - CHARACTER CHAR*1, STRING*(*) - INTEGER END, J, LENGTH - LOGICAL BLANK - -C FORCE STATIC ALLOCATION OF VARIABLES. THIS IS NOT NEEDED BY THE -C CHANGE PROGRAM, BUT IT SEEMS TO HELP WITH TROUBLESOME COMPILERS. - SAVE - -C/// SQUEEZE THE STRING. - - LENGTH = 0 - BLANK = .TRUE. - DO 0100 J = 1, LEN (STRING) - CHAR = STRING (J : J) - IF (.NOT. BLANK .OR. CHAR .NE. ' ') THEN - BLANK = CHAR .EQ. ' ' - LENGTH = LENGTH + 1 - STRING (LENGTH : LENGTH) = CHAR - END IF -0100 CONTINUE - -C/// ADJUST THE LENGTH AND PAD THE STRING. - - IF (0 .LT. LENGTH) THEN - IF (STRING (LENGTH : LENGTH) .EQ. ' ') LENGTH = LENGTH - 1 -C THE ABSOFT FORTRAN 77 COMPILER VERSION 3.1.2 FOR THE MACINTOSH -C MIS-COMPILES THE FOLLOWING STATEMENT. -C IF (LENGTH .LT. LEN (STRING)) STRING (LENGTH + 1 : ) = ' ' - END = LEN (STRING) - IF (LENGTH .LT. LEN (STRING)) STRING (LENGTH + 1 : END) = ' ' - ELSE - LENGTH = 1 - END IF - -C/// EXIT. - - RETURN - END ##### sed > copy.f <<'#####' 's/^-//' - subroutine copy(n,x,incx,y,incy) - -c copies a vector, x, to a vector, y. -c uses unrolled loops for increments equal to 1. -c jack dongarra, linpack, 3/11/78. -c modified to include precision changes, paul calamai, 20/10/92. - -c To change the precision of this program, run the change -c program on copy.f -c Alternatively, to make a single precision version append a -c comment character to the start of all lines between sequential -c precision > double -c and -c end precision > double -c comments and delete the comment character at the start of all -c lines between sequential -c precision > single -c and -c end precision > single -c comments. To make a double precision version interchange -c the append and delete operations in these instructions. - -c*****precision > double - double precision x(*),y(*) -c*****end precision > double -c*****precision > single -c real x(*),y(*) -c*****end precision > single - integer i,incx,incy,ix,iy,m,mp1,n - - if(n.le.0)return - if(incx.eq.1.and.incy.eq.1)go to 20 - -c code for unequal increments or equal increments -c not equal to 1 - - ix = 1 - iy = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 - if(incy.lt.0)iy = (-n+1)*incy + 1 - do 10 i = 1,n - y(iy) = x(ix) - ix = ix + incx - iy = iy + incy - 10 continue - return - -c code for both increments equal to 1 - - -c clean-up loop - - 20 m = mod(n,7) - if( m .eq. 0 ) go to 40 - do 30 i = 1,m - y(i) = x(i) - 30 continue - if( n .lt. 7 ) return - 40 mp1 = m + 1 - do 50 i = mp1,n,7 - y(i) = x(i) - y(i + 1) = x(i + 1) - y(i + 2) = x(i + 2) - y(i + 3) = x(i + 3) - y(i + 4) = x(i + 4) - y(i + 5) = x(i + 5) - y(i + 6) = x(i + 6) - 50 continue - return - end ##### sed > dload.f <<'#####' 's/^-//' - subroutine dload(n,da,dx,incx) - -c copies a scalar, a, to a vector, x. -c uses unrolled loops when incx equals one. - -c To change the precision of this program, run the change -c program on dload.f -c Alternatively, to make a single precision version append a -c comment character to the start of all lines between sequential -c precision > double -c and -c end precision > double -c comments and delete the comment character at the start of all -c lines between sequential -c precision > single -c and -c end precision > single -c comments. To make a double precision version interchange -c the append and delete operations in these instructions. - -c*****precision > double - double precision da,dx(*) -c*****end precision > double -c*****precision > single -c real da,dx(*) -c*****end precision > single - integer i,incx,ix,m,mp1,n - - if(n.le.0)return - if(incx.eq.1)go to 20 - -c code for incx not equal to 1 - - ix = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 - do 10 i = 1,n - dx(ix) = da - ix = ix + incx - 10 continue - return - -c code for incx equal to 1 and clean-up loop - - 20 m = mod(n,7) - if( m .eq. 0 ) go to 40 - do 30 i = 1,m - dx(i) = da - 30 continue - if( n .lt. 7 ) return - 40 mp1 = m + 1 - do 50 i = mp1,n,7 - dx(i) = da - dx(i + 1) = da - dx(i + 2) = da - dx(i + 3) = da - dx(i + 4) = da - dx(i + 5) = da - dx(i + 6) = da - 50 continue - return - end ##### sed > dot.f <<'#####' 's/^-//' -c*****precision > double - double precision function dot(n,x,incx,y,incy) -c*****end precision > double -c*****precision > single -c real function dot(n,x,incx,y,incy) -c*****end precision > single - -c forms the dot product of two vectors. -c uses unrolled loops for increments equal to one. -c jack dongarra, linpack, 3/11/78. -c modified to include precision changes, paul calamai, 20/10/92. - -c To change the precision of this program, run the change -c program on dot.f -c Alternatively, to make a single precision version append a -c comment character to the start of all lines between sequential -c precision > double -c and -c end precision > double -c comments and delete the comment character at the start of all -c lines between sequential -c precision > single -c and -c end precision > single -c comments. To make a double precision version interchange -c the append and delete operations in these instructions. - -c*****precision > double - double precision x(*),y(*),temp - double precision zero - parameter (zero=0.d0) -c*****end precision > double -c*****precision > single -c real x(*),y(*),temp -c real zero -c parameter (zero=0.0) -c*****end precision > single - integer i,incx,incy,ix,iy,m,mp1,n - - temp = zero - dot = zero - if(n.le.0)return - if(incx.eq.1.and.incy.eq.1)go to 20 - -c code for unequal increments or equal increments -c not equal to 1 - - ix = 1 - iy = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 - if(incy.lt.0)iy = (-n+1)*incy + 1 - do 10 i = 1,n - temp = temp + x(ix)*y(iy) - ix = ix + incx - iy = iy + incy - 10 continue - dot = temp - return - -c code for both increments equal to 1 - - -c clean-up loop - - 20 m = mod(n,5) - if( m .eq. 0 ) go to 40 - do 30 i = 1,m - temp = temp + x(i)*y(i) - 30 continue - if( n .lt. 5 ) go to 60 - 40 mp1 = m + 1 - do 50 i = mp1,n,5 - temp = temp + x(i)*y(i) + x(i + 1)*y(i + 1) + - * x(i + 2)*y(i + 2) + x(i + 3)*y(i + 3) + x(i + 4)*y(i + 4) - 50 continue - 60 dot = temp - return - end ##### sed > nrm2.f <<'#####' 's/^-//' -c*****precision > double - double precision function nrm2 ( n, x, incx) -c*****end precision > double -c*****precision > single -c real function nrm2 ( n, x, incx) -c*****end precision > single - -c To change the precision of this program, run the change -c program on nrm2.f -c Alternatively, to make a single precision version append a -c comment character to the start of all lines between sequential -c precision > double -c and -c end precision > double -c comments and delete the comment character at the start of all -c lines between sequential -c precision > single -c and -c end precision > single -c comments. To make a double precision version interchange -c the append and delete operations in these instructions. - - integer next,j,nn,n,incx,i -c*****precision > double - double precision x(*),cutlo,cuthi,hitest,sum,xmax,zero,one - parameter (zero=0.d0,one=1.d0) -c*****end precision > double -c*****precision > single -c real x(*),cutlo,cuthi,hitest,sum,xmax,zero,one -c parameter (zero=0.0e0,one=1.0e0) -c*****end precision > single - intrinsic abs,sqrt - -c euclidean norm of the n-vector stored in x() with storage -c increment incx . -c if n .le. 0 return with result = 0. -c if n .ge. 1 then incx must be .ge. 1 -c -c c.l.lawson, 1978 jan 08 -c modified, paul calamai, 20/10/92. -c -c four phase method using two built-in constants that are -c hopefully applicable to all machines. -c cutlo = maximum of sqrt(u/eps) over all known machines. -c cuthi = minimum of sqrt(v) over all known machines. -c where -c eps = smallest no. such that eps + 1. .gt. 1. -c u = smallest positive no. (underflow limit) -c v = largest no. (overflow limit) -c -c brief outline of algorithm.. -c -c phase 1 scans zero components. -c move to phase 2 when a component is nonzero and .le. cutlo -c move to phase 3 when a component is .gt. cutlo -c move to phase 4 when a component is .ge. cuthi/m -c where m = n for x() real and m = 2*n for complex. -c -c values for cutlo and cuthi.. -c from the environmental parameters listed in the imsl converter -c document the limiting values are as follows.. -c cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are -c univac and dec at 2**(-103) -c thus cutlo = 2**(-51) = 4.44089e-16 -c cuthi, s.p. v = 2**127 for univac, honeywell, and dec. -c thus cuthi = 2**(63.5) = 1.30438e19 -c cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. -c thus cutlo = 2**(-33.5) = 8.23181d-11 -c cuthi, d.p. same as s.p. cuthi = 1.30438d19 -c*****precision > double - parameter (cutlo=8.232d-11,cuthi=1.304d19) -c*****end precision > double -c*****precision > single -c parameter (cutlo=4.441e-16,cuthi=1.304e19) -c*****end precision > single - - if(n .gt. 0) go to 10 - nrm2 = zero - go to 300 - - 10 assign 30 to next - sum = zero - nn = n * incx -c begin main loop - i = 1 - 20 go to next,(30, 50, 70, 110) - 30 if( abs(x(i)) .gt. cutlo) go to 85 - assign 50 to next - xmax = zero - -c phase 1. sum is zero - - 50 if( x(i) .eq. zero) go to 200 - if( abs(x(i)) .gt. cutlo) go to 85 - -c prepare for phase 2. - assign 70 to next - go to 105 - -c prepare for phase 4. - - 100 i = j - assign 110 to next - sum = (sum / x(i)) / x(i) - 105 xmax = abs(x(i)) - go to 115 - -c phase 2. sum is small. -c scale to avoid destructive underflow. - - 70 if( abs(x(i)) .gt. cutlo ) go to 75 - -c common code for phases 2 and 4. -c in phase 4 sum is large. scale to avoid overflow. - - 110 if( abs(x(i)) .le. xmax ) go to 115 - sum = one + sum * (xmax / x(i))**2 - xmax = abs(x(i)) - go to 200 - - 115 sum = sum + (x(i)/xmax)**2 - go to 200 - - -c prepare for phase 3. - - 75 sum = (sum * xmax) * xmax - - -c for real or d.p. set hitest = cuthi/n -c for complex set hitest = cuthi/(2*n) - - 85 hitest = cuthi/float( n ) - -c phase 3. sum is mid-range. no scaling. - - do 95 j =i,nn,incx - if(abs(x(j)) .ge. hitest) go to 100 - 95 sum = sum + x(j)**2 - nrm2 = sqrt( sum ) - go to 300 - - 200 continue - i = i + incx - if ( i .le. nn ) go to 20 - -c end of main loop. - -c compute square root and adjust for scaling. - - nrm2 = xmax * sqrt(sum) - 300 continue - return - end ##### sed > rand.f <<'#####' 's/^-//' - real function rand(iseed) - integer 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 ********** - integer a,b15,b16,fhi,k,leftlo,p,xhi,xalo - real c -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 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 ##### sed > scal.f <<'#####' 's/^-//' - subroutine scal(n,a,x,incx) - -c scales a vector by a constant. -c uses unrolled loops for increment equal to 1. -c jack dongarra, linpack, 3/11/78. -c modified to include precision changes, paul calamai, 20/10/92. - -c To change the precision of this program, run the change -c program on scal.f -c Alternatively, to make a single precision version append a -c comment character to the start of all lines between sequential -c precision > double -c and -c end precision > double -c comments and delete the comment character at the start of all -c lines between sequential -c precision > single -c and -c end precision > single -c comments. To make a double precision version interchange -c the append and delete operations in these instructions. - -c*****precision > double - double precision a,x(*) -c*****end precision > double -c*****precision > single -c real a,x(*) -c*****end precision > single - integer i,incx,m,mp1,n,nincx - - if(n.le.0)return - if(incx.eq.1)go to 20 - -c code for increment not equal to 1 - - nincx = n*incx - do 10 i = 1,nincx,incx - x(i) = a*x(i) - 10 continue - return - -c code for increment equal to 1 - - -c clean-up loop - - 20 m = mod(n,5) - if( m .eq. 0 ) go to 40 - do 30 i = 1,m - x(i) = a*x(i) - 30 continue - if( n .lt. 5 ) return - 40 mp1 = m + 1 - do 50 i = mp1,n,5 - x(i) = a*x(i) - x(i + 1) = a*x(i + 1) - x(i + 2) = a*x(i + 2) - x(i + 3) = a*x(i + 3) - x(i + 4) = a*x(i + 4) - 50 continue - return - end ##### sed > usage.f <<'#####' 's/^-//' - program usage - -c The purpose of this program is to demonstrate the use -c of a suite of routines for generating disjointly constrained -c bilinear programming problems. -c -c The comments in subroutine blpgen describe the form -c of the problems generated. -c -c The user is referred to "Generation of disjointly constrained -c bilinear programming test problems", Computational Optimization -c and Applications, 1993 for further details. -c -c To change the precision of this program, run the change -c program on usage.f. -c Alternatively, to make a single precision version append a -c comment character to the start of all lines between sequential -c precision > double -c and -c end precision > double -c comments and delete the comment character at the start of all -c lines between sequential -c precision > single -c and -c end precision > single -c comments. To make a double precision version interchange -c the append and delete operations in these instructions. -c -c ******************** -c Subprograms called -c Local suite: blpchk, blpgen, blpfun -c ******************** -c -c Authors: Luis N. Vicente and Paul H. Calamai -c Date Written: December, 1992 -c -c ******************** - -c declarations for parameters to be passed to subroutine blpgen -c (ie. delta1,delta3,k1,k2,k11,k12,k13,k14,lcondm,non0,info and iseed) -c see comments in subroutine blpgen for details - -c Note: the dimension of delta1 and delta3 need be no greater -c than k1 - - integer k1,k2,k11,k12,k13,k14,lcondm,non0,info,iseed -c*****precision > double - double precision delta1(2),delta3(2) -c*****end precision > double -c*****precision > single -c real delta1(2),delta3(2) -c*****end precision > single - -c declarations for parameters to be passed to subroutines -c blpchk and blpfun - -c Note: this would typically be the current solution produced -c by whatever solver was being tested by the user - -c*****precision > double - double precision soln(11) -c*****end precision > double -c*****precision > single -c real soln(11) -c*****end precision > single - -c declarations for parameters returned by subroutine blpchk -c (ie. mr,errl2,type) - -c Note: the dimension of mr need be no greater than -c k1+3*k2 - - integer type -c*****precision > double - double precision mr(11),errl2 -c*****end precision > double -c*****precision > single -c real mr(11),errl2 -c*****end precision > single - -c declarations for parameters returned by subroutine blpfun -c (ie. bilfun,res) - -c Note: the dimension of res need be no greater than -c 2*k1+4*k2 - -c*****precision > double - double precision bilfun,res(16) -c*****end precision > double -c*****precision > single -c real bilfun,res(16) -c*****end precision > single - -c parameter declarations - -c*****precision > double - double precision minus4,two,four - parameter (minus4=-4.d0,two=2.0d0,four=4.d0) -c*****end precision > double -c*****precision > single -c real minus4,two,four -c parameter (minus4=-4.0,two=2.0,four=4.0) -c*****end precision > single - -c To generate disjointly constrained bilinear -c programming problems with 6 x variables and 5 -c y variables we set k1=2 and k2=3. - - k1=2 - k2=3 - -c We will demonstrate the usage of the suite of routines -c using three different versions of this 11 dimensional -c problem. Subroutines blpchk and blpfun will only be called -c in the first version of this problem. - -c ********* -c VERSION 1 -c ********* - -c For this version we want to generate a problem with 2 global -c solutions. To accomplish this we must set abs(k11)=1 and k12=0. -c We also want an additional 14 local solutions. This can be -c accomplished, among other ways, by setting k14=0 and -c abs(k11)+k12+abs(k13)=2 --> abs(k13)=1. -c If we also want to provide values for delta1 and delta3 -c we must set k11>0 and k13>0. - - k11=1 - k12=0 - k13=1 - k14=0 - -c Since k11=1 we must provide a value for delta1(1) in -c the open interval (1,3). Similarily, since k13=1 -c we must provide a value for delta3(1) in the open -c interval (3,+infinity). -c We choose delta1(1)=2.0 (to coincide with the diagrams -c corresponding to Class 1 in the COA article) -c and delta3(1)=4.0. - - delta1(1)=two - delta3(1)=four - -c Finally, to obtain a sparse (but ideally conditioned) -c transformation matrix (see comments in subroutines -c blpgen and blpmkm) we set non0=1 and lcondm=0. - -c NOTE: when lcondm=0 the transformation matrix should -c (with exact arithmetic) be the identity matrix -c (We do not, however, use this fact explicitly -c in the code). - - non0=1 - lcondm=0 - -c We now need simply to initialize the seed for the -c random number generator (see subroutine rand) and -c then call subroutine blpgen. - - iseed=1 - call blpgen(delta1,delta3,k1,k2,k11,k12,k13,k14, - 1 lcondm,non0,info,iseed) - -c The value of info returned by subroutine blpgen -c should be checked to make sure that termination -c was normal (ie. info=0). - - if (info.ne.0) then - print *,'** Something is wrong (info=',info,')' - print *,'** VERSION 1 run abandoned **' - else - -c Much of the data defining the problem generated by -c subroutine blpgen is stored in named common blocks -c (see comments in subroutine blpgen). However, the -c user may only need to call subroutine blpchk and -c subroutine blpfun in order to access this information. - -c To demonstrate this concept suppose that the user -c wishes to evaluate the functions and constraints -c of the generated bilinear problem. This can easily -c be accomplished by calling subroutine blpfun. - -c Suppose, for example, that the user wishes to evaluate -c the bilinear problem at the point stored in vector soln. -c Here we set soln to a feasible point of the problem we -c have just generated. - - soln(1) = 0.25d0 - soln(2) = 2.00d0 - soln(3) = 1.75d0 - soln(4) = 2.00d0 - soln(5) = 1.00d0 - soln(6) = 0.00d0 - soln(7) = 1.75d0 - soln(8) = 0.00d0 - soln(9) = 0.25d0 - soln(10)= 0.25d0 - soln(11)= 2.00d0 - -c The following call to subroutine blpfun - - call blpfun(soln,bilfun,res) - print *, 'bilfun = ',bilfun - print *, 'res = ' - print 10, res - -c would return the value of the objective -c evaluated at soln and the residuals of the -c constraints (ie. their slacks) in the vector res. -c In this particular case the values returned should -c (approximately) be: - -c bilfun = -9.8750000000000 -c res = -c 0.00000000000000 0.50000000000000 -c 3.50000000000000 0.00000000000000 -c 3.50000000000000 0.50000000000000 -c 2.00000000000000 0.00000000000000 -c 0.00000000000000 3.50000000000000 -c 0.50000000000000 0.00000000000000 -c 0.75000000000000 6.75000000000000 -c 0.50000000000000 0.00000000000000 - -c Suppose, in addition, that the user wishes to see -c if the point (in soln) is a minimizer of the generated -c problem. To an extent this question can be answered -c using the following call to subroutine blpchk. - - call blpchk(soln,mr,errl2,type) - print *,'errl2 = ',errl2 - print *,'type = ',type - print *,'mr = ' - print 10, mr - -c As described in the comments of subroutine blpchk -c this call actually results in a comparison against -c all minimizers in the untransformed space. If the value -c returned in errl2 is "small" then the point stored in -c soln is "close" to a minimizer of the type specified by -c the value returned in type (see comments in subroutine -c blpchk). The point returned in vector mr corresponds to -c the closest minimizer (in the untransformed space). -c In this particular case the values returned should -c (approximately) be: - -c errl2 = 0.55901699437495 -c type = 3 -c mr = -c 0.00000000000000 2.00000000000000 -c 2.00000000000000 2.00000000000000 -c 1.00000000000000 0.00000000000000 -c 2.00000000000000 0.00000000000000 -c 0.00000000000000 0.00000000000000 -c 2.00000000000000 - - print *,'** VERSION 1 **' - print *,' Check output above against values in' - print *,' program usage.' - endif - -c ********* -c VERSION 2 -c ********* - -c For this version we want to generate a problem with a -c unique global solution. To accomplish this we set -c k11=k12=0. - - k11=0 - k12=0 - -c In order to generate an additional 15 local solutions (see -c the comments in subroutine blpgen or the COA article) we -c must have 2*abs(k13)+k14=4. This can be accomplished, among -c other ways, by setting abs(k13)=2 and k14=0. -c To get subroutine blpgen to set both components -c of delta3 to 4.0 we set k13=2 and delta3(1)=-4.0 -c as per the instructions in subroutine blpgen. - - k13=2 - k14=0 - delta3(1)=minus4 - -c Finally, to obtain a (relatively) dense transformation matrix -c with two-norm condition 10 we set non0=6 and lcondm=1. - - non0=6 - lcondm=1 - -c We now need simply to call subroutine blpgen (iseed has -c already been initialized). - - call blpgen(delta1,delta3,k1,k2,k11,k12,k13,k14, - 1 lcondm,non0,info,iseed) - if (info.ne.0) then - print *,'** Something is wrong (info=',info,')' - print *,'** VERSION 2 **' - endif - -c ********* -c VERSION 3 -c ********* - -c For this version we want to generate a problem with a -c unique global solution. To accomplish this we set -c k11=k12=0. - - k11=0 - k12=0 - -c In order to generate an additional 15 local solutions (see -c the comments in subroutine blpgen or the COA article) we -c must have 2*abs(k13)+k14=4. This can be accomplished, among -c other ways, by setting abs(k13)=2 and k14=0. -c To get subroutine blpgen to select both components -c of delta3 in (three,delta3(1)) we set k13=-2 and delta3(1) -c to some value greater than three ( for instance -c delta3(1)=4.0) as per the instructions in subroutine blpgen. - - k13=-2 - k14=0 - delta3(1)=four - -c Finally, to obtain a (relatively) sparse transformation matrix -c with two-norm condition 100 we set non0=2 and lcondm=2. - - non0=2 - lcondm=2 - -c We now need simply to call subroutine blpgen (iseed has -c already been initialized). - - call blpgen(delta1,delta3,k1,k2,k11,k12,k13,k14, - 1 lcondm,non0,info,iseed) - if (info.ne.0) then - print *,'** Something is wrong (info=',info,')' - print *,'** VERSION 3 **' - endif - -c format specifications - - 10 format(9x,2f19.14) - stop -c last card of usage - end ##### sed > verify.f <<'#####' 's/^-//' - program verify - integer maxk2 - parameter (maxk2=2000) - -c The purpose of this program is to test (albeit minimally) the -c integrity of the routines in this suite after they are -c installed. It is by no means a thorough test. For usage -c examples see program usage. -c -c The user is referred to the article "Generation of disjointly -c constrained bilinear programming test problems", Computational -c Optimization and Applications, 1993 for -c details about the suite, and to the comments in -c subroutines blpgen, blpchk and blpfun for details about -c the parameters. -c -c To change the precision of this program, run the change -c program on verify.f. -c Alternatively, to make a single precision version append a -c comment character to the start of all lines between sequential -c precision > double -c and -c end precision > double -c comments and delete the comment character at the start of all -c lines between sequential -c precision > single -c and -c end precision > single -c comments. To make a double precision version interchange -c the append and delete operations in these instructions. -c -c ******************** -c Subprograms called -c Local suite: blpax,blpchk,blpgen,blpmkm,blpmx,blpfun,dload -c Modified Linpack: axpy,nrm2 -c ******************** -c -c Authors: Paul H. Calamai and Luis N. Vicente -c Date Written: December, 1992 -c -c ******************** - -c declarations for parameters to be passed to subroutine blpgen -c (ie. delta1,delta3,lcondm,non0,info and iseed) -c see comments in subroutine blqgen for details - -c Note: the dimension of delta1 and delta3 need be no greater -c than maxk2 - - integer lcondm,non0,info,iseed -c*****precision > double - double precision delta1(maxk2),delta3(maxk2) -c*****end precision > double -c*****precision > single -c real delta1(maxk2),delta3(maxk2) -c*****end precision > single - -c declarations for parameters returned by subroutine blpchk -c (ie. mr,errl2,type) - -c Note: the dimension of mr need be no greater than -c 4*maxk2 - - integer type -c*****precision > double - double precision mr(4*maxk2),errl2 -c*****end precision > double -c*****precision > single -c real mr(4*maxk2),errl2 -c*****end precision > single - -c declarations for parameters returned by subroutine blpfun -c (ie. bilfun,res) - -c Note: the dimension of res need be no greater than -c 6*maxk2 - -c*****precision > double - double precision bilfun,res(6*maxk2) -c*****end precision > double -c*****precision > single -c real bilfun,res(6*maxk2) -c*****end precision > single - -c declarations for local parameters - - integer i -c*****precision > double - double precision error,x(4*maxk2),y(4*maxk2) -c*****end precision > double -c*****precision > single -c real error,x(4*maxk2),y(4*maxk2) -c*****end precision > single - -c function declarations - -c*****precision > double - double precision nrm2 -c*****end precision > double -c*****precision > single -c real nrm2 -c*****end precision > single - -c parameter declarations - -c*****precision > double - double precision minus1,one,two,four - parameter (minus1=-1.d0,one=1.d0,two=2.0d0,four=4.d0) -c*****end precision > double -c*****precision > single -c real minus1,one,two,four -c parameter (minus1=-1.0,one=1.0,two=2.0,four=4.0) -c*****end precision > single - -c declarations for named common /m/ - - integer num0 - integer pz(4*maxk2) -c*****precision > double - double precision d(4*maxk2),z(4*maxk2) -c*****end precision > double -c*****precision > single -c real d(4*maxk2),z(4*maxk2) -c*****end precision > single - common/m/d,z,pz,num0 - -c declarations for named common /dim/ - - integer ik1,ik2,nx,ny,ik11,ik12,ik13,ik14,n - common/dim/ik1,ik2,nx,ny,ik11,ik12,ik13,ik14,n - -c ******************* -c VERIFICATION TEST 1 -c ******************* - -c Subroutines tested: -c blpchk, blpgen, blpfun -c ******************* - -c The first test will be conducted using VERSION 1 of the -c problem described in program usage. See the comments -c in program usage for details. - - delta1(1)=two - delta3(1)=four - ik1=2 - ik2=3 - ik11=1 - ik12=0 - ik13=1 - ik14=0 - lcondm=0 - non0=1 - iseed=1 - -c To conduct a test of subroutines blpfun and blpchk we begin -c by calling the problem generator (ie. subroutine blpgen) -c using the above data. The fact that lcondm=0 is essential -c since this yields an identity transformation. - - call blpgen(delta1,delta3,ik1,ik2,ik11,ik12,ik13,ik14, - 1 lcondm,non0,info,iseed) - -c The value of info returned by subroutine blpgen should -c be checked to make sure that termination was normal -c (ie. info=0). - - if (info.ne.0) then - print *,'** Something is wrong (info=',info,')' - print *,'** VERIFICATION TEST 1 abandoned **' - else - -c Since the transformed and untransformed problem are identical -c we have complete knowledge of all minimizers of the generated -c problem. To test these subroutines we will pass them one of -c the known global solutions using vector parameter x. - - x(1) = 0.00d0 - x(2) = 2.00d0 - x(3) = 1.00d0 - x(4) = 0.00d0 - x(5) = 1.00d0 - x(6) = 0.00d0 - x(7) = 2.00d0 - x(8) = 0.00d0 - x(9) = 1.00d0 - x(10) = 4.00d0 - x(11) = 2.00d0 - -c The following call to subroutine blpfun - - call blpfun(x,bilfun,res) - print *, 'bilfun = ',bilfun - print *, 'res = ' - print 10, (res(i),i=1,2*ik1+4*ik2) - -c should return the value of the objective and the -c residuals of the constraints (ie. their slacks) -c in the vector res. -c In this particular case the values returned should -c (approximately) be: - -c bilfun = -12.00000000000000 -c res = -c 0.00000000000000 0.00000000000000 -c 4.00000000000000 2.00000000000000 -c 0.00000000000000 0.00000000000000 -c 2.00000000000000 0.00000000000000 -c 0.00000000000000 4.00000000000000 -c 0.00000000000000 0.00000000000000 -c 0.00000000000000 0.00000000000000 -c 8.00000000000000 0.00000000000000 - -c Similarly, the following call to subroutine blpchk - - call blpchk(x,mr,errl2,type) - print *,'errl2 = ',errl2 - print *,'type = ',type - print *,'mr = ' - print 10, (mr(i),i=1,ik1+3*ik2) - -c should yield the following values for the returned parameters: - -c errl2 = 0. -c type = 2 -c mr = -c 0.00000000000000 2.00000000000000 -c 1.00000000000000 0.00000000000000 -c 1.00000000000000 0.00000000000000 -c 2.00000000000000 0.00000000000000 -c 1.00000000000000 4.00000000000000 -c 2.00000000000000 - - print *,'** VERIFICATION TEST 1 completed **' - print *,' Check output above against values in' - print *,' program verify.' - endif - -c ******************* -c VERIFICATION TEST 2 -c ******************* - -c Subroutines tested: -c blpmkm, blpmx -c ******************* - -c Note: the typical user will never need to call these -c routines directly - -c Since the transformation matrix M is defined by the -c product H*D*H, where H is a Householder matrix and -c D is a diagonal matrix, we can easily check subroutine -c blpmkm (which constructs M) and blpmx (that forms the -c product of M and a vector x). We do this by forcing -c D to be the identity matrix (by setting lcondm=0) and -c setting x to the ones vector. In this case, M should -c be the identity matrix and y=M*x should equal x. -c We then set y=y-x and form y'y which should equal 0. - -c Note: n can be set to any integer in the range [1,4*maxk2] -c and non0 can be set to any integer in the range -c [1,n] - - n=11 - non0=10 - lcondm=0 - call dload(n,one,x,1) - call blpmkm(lcondm,iseed) - call blpmx(x,y,.false.) - call axpy(n,minus1,x,1,y,1) - error=nrm2(n,y,1) - print *,'** VERIFICATION TEST 2 completed **' - print *,' error = ',error - print *,' error should (approximately) equal 0.' - -c ******************* -c VERIFICATION TEST 3 -c -c Subroutine tested: -c blpax -c ******************* - -c Note: the typical user will never need to call this -c routine directly - -c To test routine blpax we set the all components of x to one. - - call dload(ik1+3*ik2,one,x,1) - -c The following call to subroutine blpax - - call blpax(x,y) - print *,'y = ' - print 10, (y(i),i=1,2*ik1+4*ik2) - -c should return the value of the product of the matrix -c -c ( A 0 ) -c ( 0 B ) -c -c by the vector x in the vector y. -c In this particular case the values returned should -c (approximately) be: - -c y = -c 1.00000000000000 -3.00000000000000 -c 1.00000000000000 1.00000000000000 -c -3.00000000000000 1.00000000000000 -c 1.00000000000000 -3.00000000000000 -c 1.00000000000000 -1.00000000000000 -c 3.00000000000000 -2.00000000000000 -c -3.00000000000000 5.00000000000000 -c -2.00000000000000 1.00000000000000 - - print *,'** VERIFICATION TEST 3 completed **' - print *,' Check output above against values in' - print *,' program verify.' - - 10 format(9x,2f19.14) - stop -c last card of program verify - end #####