#!/bin/sh # to unpack, sh this message PATH=/bin:/usr/bin cat > makefile <<'CUT HERE............' CUT HERE............ cat > README <<'CUT HERE............' This is the README file for the STENMIN package: STENMIN: A Software Package for Large, Sparse Unconstrained Optimization Using Tensor Methods STENMIN minimizes an unconstrained nonlinear function in n unknowns where the Hessian is large and sparse using a new class of methods called tensor methods. The software allows the user to select between a tensor method and a standard method based upon a quadratic model. The tensor method models the objective function by a fourth-order model, where the third- and fourth-order terms are chosen such that the extra cost of forming and solving the model is small. The test results indicate that, in general, the tensor method is often more efficient and more reliable than the standard Newton method for solving large, sparse unconstrained optimization problems. It is especially useful on problems where the Hessian at the solution is singular. Two interfaces are provided with STENMIN. If the user wishes to use all the defaults options provided by the package, then he (or she) should call STUMSD (STUMSS if single-precision is used). The other interface, STUMCD (STUMCS if single-precision is used), requires the user to supply all parameters. The user may specify selected parameters only by first invoking the subroutine DFAULT, which sets all parameters to their default values, and then overriding only the desired values. The two calling sequences are as follows: 1. CALL STUMSD(N, X0, NZ, IRN, LIRN, ICN, LICN, FCN, TYPX, MSG, * XPLS, FPLS, GPLS, HESS, WRK, LWRK, IWRK, LIWRK, HTV, TERMCD) 2. CALL STDFLT(N, TYPX, FSCALE, GRADTL, STEPTL, ILIM, STEPMX, * IPR, METHOD, GRDFLG, HSNFLG, NDIGIT, INFORM, MSG) C USER OVERRIDES SPECIFIC DEFAULT VALUES PARAMETERS, E.G. GRADTL = 1.0D-6 ILIM = 1000 GRDFLG = 1 HSNFLG = 1 CALL STUMCD(N, X0, NZ, IRN, LIRN, ICN, LICN, FCN, UGR, * USH, TYPX, FSCALE, GRADTL, STEPTL, ILIM, STEPMX, IPR, * METHOD, GRDFLG, HSNFLG, NDIGIT, MSG, XPLS, FPLS, GPLS, * HESS, WRK, LWRK, IWRK, LIWRK, TERMCD, HTV, INFORM) You should first check that the following files are in your directory: 1. README 2. todble.f 3. tosngl.f 4. makefile 5. driver1.f 6. driver2.f 7. broyden.f 8. odc.f 9. output1 10. output2 11. stenmin.f 12. colmor.f Note: driver1.f is a sample driver for solving a simple unconstrained minimization problem whereas driver2.f is a sample driver for a more elaborate application of the STENMIN package. To obtain the executable "todble" that converts STENMIN to a double precision version type: f77 -O -o todble todble.f To obtain the executable "tosngl" that converts STENMIN to a single precision version type: f77 -O -o tosngl tosngl.f To obtain a double or a single precision version of STENMIN, compile, and link with the Broyden tridiagonal test problem "broyden.f" type: make broydend (double precision) or make broydens (single precision) To execute test problem "broyden.f" and redirect the output to file "out1" type: stenmin > out1 To obtain a double or a single precision version of STENMIN, compile, and link with the optimal design with composite materials problem "odc.f" type: make odcd (double precision) or make odcs (single precision) To execute test problem "odc.f" and redirect the output to file "out2" type: stenmin > out2 For answers to technical questions about STENMIN, contact: Dr. Ali Bouaricha MCS division Argonne National Laboratory Argonne, Illinois 60439 e-mail: bouarich@mcs.anl.gov IMPORTANT: STENMIN uses the MA27 package. Due to licencing restrictions the MA27 is NOT included as part of STENMIN. The user needs to obtain a copy of the MA27 package directly from HARWELL, UK. The contact person is Dr. Iain Duff (e-mail: duff@cerfacs.fr). CUT HERE............ cat > todble.f <<'CUT HERE............' PROGRAM TODBLE * SPECIALIZES SOURCE FILES TO DOUBLE PRECISION * (SUN/UNIX VERSION) * PH.L. TOINT, APRIL 1991. INTEGER INDEV, OUTDEV, EOL, K CHARACTER*1 LINE(80) CHARACTER*5 HEAD PARAMETER ( INDEV = 5, OUTDEV = 6 ) * LOOP ON THE LINES 10 CONTINUE READ( INDEV, '(80A1)', END = 99 ) ( LINE(K), K = 1, 80 ) * DETECT END OF LINE DO 20 EOL = 80, 1, -1 IF( LINE(EOL) .NE. ' ' ) GO TO 30 20 CONTINUE * ACTIVATE THE STATEMENTS SUITABLE FOR THE DESIRED ENVIRONMENT 30 CONTINUE IF( LINE(1) .EQ. 'C' ) THEN DO 25 K = 1, 5 HEAD(K:K) = LINE(K) 25 CONTINUE IF( HEAD .EQ. 'CD ' ) GO TO 50 IF( HEAD .EQ. 'CSUN ' ) GO TO 50 IF( HEAD .EQ. 'CUNIX' ) GO TO 50 IF( HEAD .EQ. 'CIEEE' ) GO TO 50 END IF GO TO 40 * REPLACE HEAD BY 5 BLANKS 50 CONTINUE DO 60 K = 1, 5 LINE(K) = ' ' 60 CONTINUE * PRINT LINE 40 CONTINUE WRITE( OUTDEV, '(80A1)' ) ( LINE(K), K = 1, EOL ) GO TO 10 * END OF FILE 99 CONTINUE STOP END CUT HERE............ cat > tosngl.f <<'CUT HERE............' PROGRAM TOSNGL * SPECIALIZES SOURCE FILES TO SINGLE PRECISION * (SUN/UNIX VERSION) * PH.L. TOINT, APRIL 1991, FOR CGT PRODUCTIONS. INTEGER INDEV, OUTDEV, EOL, K CHARACTER*1 LINE(80) CHARACTER*5 HEAD PARAMETER ( INDEV = 5, OUTDEV = 6 ) * LOOP ON THE LINES 10 CONTINUE READ( INDEV, '(80A1)', END = 99 ) ( LINE(K), K = 1, 80 ) * DETECT END OF LINE DO 20 EOL = 80, 1, -1 IF( LINE(EOL) .NE. ' ' ) GO TO 30 20 CONTINUE * ACTIVATE THE STATEMENTS SUITABLE FOR THE DESIRED ENVIRONMENT 30 CONTINUE IF( LINE(1) .EQ. 'C' ) THEN DO 25 K = 1, 5 HEAD(K:K) = LINE(K) 25 CONTINUE IF( HEAD .EQ. 'CS ' ) GO TO 50 IF( HEAD .EQ. 'CUNIX' ) GO TO 50 IF( HEAD .EQ. 'CSUN ' ) GO TO 50 IF( HEAD .EQ. 'CIEEE' ) GO TO 50 END IF GO TO 40 * REPLACE HEAD BY 5 BLANKS 50 CONTINUE DO 60 K = 1, 5 LINE(K) = ' ' 60 CONTINUE * PRINT LINE 40 CONTINUE WRITE( OUTDEV, '(80A1)' ) ( LINE(K), K = 1, EOL ) GO TO 10 * END OF FILE 99 CONTINUE STOP END CUT HERE............ cat > makefile <<'CUT HERE............' # # This makefile uses the double precision version of STENMIN # FFLAGS = -O # Files for STENMIN # # The program runs sparse unconstrained optimization problems # using tensor methods MAIND1 = maind1.o MAINS1 = mains1.o MAIND2 = maind2.o MAINS2 = mains2.o OBJD = stenmind.o colmord.o OBJS = stenmins.o colmors.o BROYDEND = broydend.o BROYDENS = broydens.o ODCD = odcd.o ODCS = odcs.o broydend : $(MAIND1) $(OBJD) $(BROYDEND) f77 $(FFLAGS) $(MAIND1) $(OBJD) $(BROYDEND) -o stenmin broydens : $(MAINS1) $(OBJS) $(BROYDENS) f77 $(FFLAGS) $(MAINS1) $(OBJS) $(BROYDENS) -o stenmin odcd : $(MAIND2) $(OBJD) $(ODCD) f77 $(FFLAGS) $(MAIND2) $(OBJD) $(ODCD) -o stenmin odcs : $(MAINS2) $(OBJS) $(ODCS) f77 $(FFLAGS) $(MAINS2) $(OBJS) $(ODCS) -o stenmin broydend.f: broyden.f /bin/rm -f broydend.f; todble < broyden.f > broydend.f broydens.f: broyden.f /bin/rm -f broydens.f; tosngl < broyden.f > broydens.f odcd.f: odc.f /bin/rm -f odcd.f; todble < odc.f > odcd.f odcs.f: odc.f /bin/rm -f odcs.f; tosngl < odc.f > odcs.f maind1.f: driver1.f /bin/rm -f maind1.f; todble < driver1.f > maind1.f mains1.f: driver1.f /bin/rm -f mains1.f; tosngl < driver1.f > mains1.f maind2.f: driver2.f /bin/rm -f maind2.f; todble < driver2.f > maind2.f mains2.f: driver2.f /bin/rm -f mains2.f; tosngl < driver2.f > mains2.f colmord.f: colmor.f /bin/rm -f colmord.f; todble < colmor.f > colmord.f colmors.f: colmor.f /bin/rm -f colmors.f; tosngl < colmor.f > colmors.f stenmind.f: stenmin.f /bin/rm -f stenmind.f; todble < stenmin.f > stenmind.f stenmins.f: stenmin.f /bin/rm -f stenmins.f; tosngl < stenmin.f > stenmins.f CUT HERE............ cat > driver1.f <<'CUT HERE............' C C STENMIN MINIMIZES AN UNCONSTRAINED NONLINEAR FUNCTION IN N C UNKNOWNS WHERE THE HESSIAN IS LARGE AND SPARSE, USING TENSOR C METHODS. C C EXAMPLE OF USE FOR STENMIN. THE TEST PROBLEM IS C THE BROYDEN TRIDIAGONAL (SOURCE: BUCKLEY#78 (P. 42)). C C ALI BOUARICHA, OCTOBER 1994. C MCS DIVISION, ARGONNE NATIONAL LAB. C INTEGER NMAX, N, NZ, LIRN, LICN, ILIM, IPR, METHOD INTEGER GRDFLG, HSNFLG, NDIGIT, MSG, LWRK, LIWRK INTEGER TERMCD, INFORM, I CD DOUBLE PRECISION FSCALE, GRADTL, STEPTL, FPLS, STEPMX, ONE CS REAL FSCALE, GRADTL, STEPTL, FPLS, STEPMX, ONE PARAMETER ( NMAX = 10000, LIRN = 50000, LICN = 500000 ) PARAMETER ( LIWRK = 2 * LIRN + 12 * NMAX + 2 ) PARAMETER ( LWRK = 7 * NMAX ) INTEGER IRN ( LIRN ), ICN ( LICN ) INTEGER IWRK( LIWRK ) CD DOUBLE PRECISION X ( NMAX ), TYPX( NMAX ), XPLS( NMAX ) CD DOUBLE PRECISION GPLS( NMAX ), HESS( LICN ), WRK ( LWRK ) CD DOUBLE PRECISION HTV ( NMAX ) CS REAL X ( NMAX ), TYPX( NMAX ), XPLS( NMAX ) CS REAL GPLS( NMAX ), HESS( LICN ), WRK ( LWRK ) CS REAL HTV ( NMAX ) EXTERNAL FCN, UGRAD, UHESS CD DATA ONE / 1.0D0 / CS DATA ONE / 1.0E0 / C RUN BROYDEN PROBLEM WITH N = 10000. N = 10000 C COMPUTE THE STANDARD STARTING POINT. DO 10 I = 1, N X(I) = -ONE 10 CONTINUE C SET THE DEFAULT VALUES. CALL STDFLT(N, TYPX, FSCALE, GRADTL, STEPTL, ILIM, STEPMX, * IPR, METHOD, GRDFLG, HSNFLG, NDIGIT, INFORM, MSG) CD GRADTL = 1.0D-5 CS GRADTL = 1.0E-3 GRDFLG = 2 HSNFLG = 2 C CALL THE SPARSE OPTIMIZER. CD CALL STUMCD(N,X,NZ,IRN,LIRN,ICN,LICN,FCN,UGRAD, CS CALL STUMCS(N,X,NZ,IRN,LIRN,ICN,LICN,FCN,UGRAD, * UHESS,TYPX,FSCALE,GRADTL,STEPTL,ILIM,STEPMX,IPR, * METHOD,GRDFLG,HSNFLG,NDIGIT,MSG,XPLS,FPLS,GPLS, * HESS,WRK,LWRK,IWRK,LIWRK,TERMCD,HTV,INFORM) STOP END CUT HERE............ cat > broyden.f <<'CUT HERE............' C C THE FOLLOWING IS A SUBROUTINE FOR THE BROYDEN TRIDIAGONAL C PROBLEM (SOURCE: BUCKLEY#78 (P. 42)) C SUBROUTINE FCN(N, X, F) INTEGER N CD DOUBLE PRECISION X(N), F CS REAL X(N), F C C LOCAL VARIABLES C INTEGER I CD DOUBLE PRECISION ONE, TWO, THREE CS REAL ONE, TWO, THREE CD DATA ONE, TWO, THREE / 1.0D0, 2.0D0, 3.0D0 / CS DATA ONE, TWO, THREE / 1.0E0, 2.0E0, 3.0E0 / C F = ((THREE - TWO * X(1)) * X(1) - TWO * X(2) + ONE) ** 2 DO 10 I = 2, N-1 F = F + ((THREE - TWO * X(I)) * X(I) - X(I-1) - * TWO * X(I+1) + ONE) ** 2 10 CONTINUE F = F + ((THREE - TWO * X(N)) * X(N) - X(N-1) + ONE) ** 2 RETURN END C C THE FOLLOWING IS A SUBROUTINE FOR THE GRADIENT OF THE BROYDEN C TRIDIAGONAL PROBLEM C SUBROUTINE UGRAD(N, X, G) INTEGER N CD DOUBLE PRECISION X(N), G(N) CS REAL X(N), G(N) C C LOCAL VARIABLES C INTEGER I CD DOUBLE PRECISION RL, RM, RR, ONE, TWO, THREE, FOUR CS REAL RL, RM, RR, ONE, TWO, THREE, FOUR CD DATA ONE, TWO, THREE, FOUR/1.0D0, 2.0D0, 3.0D0, 4.0D0/ CS DATA ONE, TWO, THREE, FOUR/ 1.0E0, 2.0E0, 3.0E0, 4.0E0/ C RL = (THREE - TWO * X(1)) * X(1) - TWO * X(2) + ONE RR = (THREE - TWO * X(2)) * X(2) - X(1) - TWO * X(3) + ONE G(1) = TWO * (RL * (THREE - FOUR * X(1)) - RR) DO 10 I = 2, N-1 IF(I .NE. 2) THEN RL = (THREE - TWO * X(I-1)) * X(I-1) - X(I-2) - * TWO * X(I) + ONE ENDIF RM = (THREE - TWO * X(I)) * X(I) - X(I-1) - * TWO * X(I+1) + ONE IF(I .EQ. N-1) THEN RR = (THREE - TWO * X(N)) * X(N) - X(N-1) + ONE ELSE RR = (THREE - TWO * X(I+1)) * X(I+1) - X(I) - * TWO * X(I+2) + ONE ENDIF G(I) = -TWO * (TWO * RL - RM * (THREE - FOUR * X(I)) + RR) 10 CONTINUE G(N) = -TWO * (TWO * RM - RR * (THREE - FOUR * X(N))) RETURN END C C THE FOLLOWING IS A SUBROUTINE FOR THE HESSIAN OF THE BROYDEN C TRIDIAGONAL PROBLEM C SUBROUTINE UHESS(N,X,NZ,LICN,HESS,IRN,ICN) INTEGER N, NZ, LICN INTEGER IRN(*), ICN(LICN) CD DOUBLE PRECISION X(N), HESS(LICN) CS REAL X(N), HESS(LICN) C C LOCAL VARIABLES C INTEGER I CD DOUBLE PRECISION RL,RM,RR CD DOUBLE PRECISION ONE,TWO,THREE,FOUR CS REAL RL,RM,RR CS REAL ONE,TWO,THREE,FOUR CD DATA ONE, TWO, THREE, FOUR/1.0D0, 2.0D0, 3.0D0, 4.0D0/ CS DATA ONE, TWO, THREE, FOUR/1.0E0, 2.0E0, 3.0E0, 4.0E0/ C NZ = 1 RL = (THREE - TWO * X(1)) * X(1) - TWO * X(2) + ONE HESS(NZ) = TWO * ((THREE - FOUR * X(1))**2 - * FOUR * RL + ONE) IRN(NZ) = 1 ICN(NZ) = 1 DO 10 I = 2, N-1 IF(I .NE. 2) THEN NZ = NZ + 1 HESS(NZ) = FOUR IRN(NZ) = I ICN(NZ) = I-2 ENDIF NZ = NZ + 1 HESS(NZ) = -TWO * (TWO * (THREE - FOUR * X(I-1)) + * ONE * (THREE - FOUR * X(I))) IRN(NZ) = I ICN(NZ) = I-1 RM = (THREE - TWO * X(I)) * X(I) - X(I-1) - * TWO * X(I+1) + ONE NZ = NZ + 1 HESS(NZ) = -TWO * (-FOUR - (THREE - FOUR * X(I))**2 + * FOUR * RM - ONE) IRN(NZ) = I ICN(NZ) = I 10 CONTINUE RR = (THREE - TWO * X(N)) * X(N) - X(N-1) + ONE NZ = NZ + 1 HESS(NZ) = FOUR IRN(NZ) = N ICN(NZ) = N-2 NZ = NZ + 1 HESS(NZ) = -TWO * (TWO * (THREE - FOUR * X(N-1)) + * THREE - FOUR * X(N)) IRN(NZ) = N ICN(NZ) = N-1 NZ = NZ + 1 HESS(NZ) = TWO * (FOUR + (THREE - FOUR * X(N))**2 - FOUR * RR) IRN(NZ) = N ICN(NZ) = N RETURN END CUT HERE............ cat > output1 <<'CUT HERE............' STDRUO GRADIENT FLAG = 2 STDRUO HESSIAN FLAG = 2 STDRUO METHOD = 1 STDRUO ITERATION LIMIT = 500 STDRUO MACHINE EPSILON = 0.2220446049250E-15 STDRUO STEP TOLERANCE = 0.3666852862501E-10 STDRUO GRADIENT TOLERANCE = 0.1000000000000E-04 STDRUO MAXIMUM STEP SIZE = 0.1000000000000E+06 --------------------------------------------- STRSLT ITERATION K = 0 STRSLT FUNCTION AT X(K) STRSLT 0.1001100000000E+05 STRSLT SCALED GRADIENT AT X(K) STRSLT 0.3800000000000E+02 --------------------------------------------- STCHKS RELATIVE GRADIENT CLOSE TO ZERO STCHKS CURRENT ITERATE IS PROBABLY SOLUTION --------------------------------------------- STRSLT ITERATION K = 4 STRSLT FUNCTION AT X(K) STRSLT 0.1884575867777E-13 STRSLT SCALED GRADIENT AT X(K) STRSLT 0.1113397081739E-05 --------------------------------------------- STRSLT NUMBER OF FUNCTION EVALUATIONS 5 STRSLT NUMBER OF GRADIENT EVALUATIONS 5 STRSLT NUMBER OF HESSIAN EVALUATIONS 4 CUT HERE............ cat > driver2.f <<'CUT HERE............' C C STENMIN MINIMIZES AN UNCONSTRAINED NONLINEAR FUNCTION IN N C UNKNOWNS WHERE THE HESSIAN IS LARGE AND SPARSE, USING TENSOR C METHODS. C C*****WARNING***** C C THE SINGLE-PRECISION FINITE DIFFERENCE HESSIAN APPROXIMATION C MAY NOT BE ACCURATE ENOUGH. THIS MAY LEAD TO A CONVERGENCE C DETERIORATION. THE USER IS ADVISED TO PROVIDE AN ANALYTIC C HESSIAN WHEN THE SINGLE-PRECISION VERSION OF STENMIN IS USED. C C***************** C C EXAMPLE OF USE FOR STENMIN. THE TEST PROBLEM IS THE C OPTIMAL DESIGN WITH COMPOSITE MATERIALS PROBLEM FROM C THE MINPACK-2 TEST PROBLEM COLLECTION. C C ALI BOUARICHA, OCTOBER 1994. C MCS DIVISION, ARGONNE NATIONAL LAB. C INTEGER NMAX, N, NZ, LIRN, LICN, ILIM, IPR, METHOD INTEGER GRDFLG, HSNFLG, NDIGIT, MSG, LWRK, LIWRK INTEGER TERMCD, INFORM, I, J, K, NX, NY CD DOUBLE PRECISION FSCALE, GRADTL, STEPTL, FPLS, STEPMX CD DOUBLE PRECISION LAMBDA, HX, HY, TEMP, ONE CS REAL FSCALE, GRADTL, STEPTL, FPLS, STEPMX CS REAL LAMBDA, HX, HY, TEMP, ONE PARAMETER ( NMAX = 10000, LIRN = 50000, LICN = 500000 ) PARAMETER ( LIWRK = 2 * LIRN + 12 * NMAX + 2 ) PARAMETER ( LWRK = 7 * NMAX ) INTEGER IRN ( LIRN ), ICN ( LICN ) INTEGER IWRK( LIWRK ) CD DOUBLE PRECISION X ( NMAX ), TYPX( NMAX ), XPLS( NMAX ) CD DOUBLE PRECISION GPLS( NMAX ), HESS( LICN ), WRK ( LWRK ) CD DOUBLE PRECISION HTV ( NMAX ) CS REAL X ( NMAX ), TYPX( NMAX ), XPLS( NMAX ) CS REAL GPLS( NMAX ), HESS( LICN ), WRK ( LWRK ) CS REAL HTV ( NMAX ) COMMON / PARAM / NX, NY COMMON / OTHER / LAMBDA EXTERNAL DODCF, DODCG, STDUSH CD INTRINSIC DBLE, MIN CS INTRINSIC FLOAT, MIN CD DATA ONE / 1.0D0 / CS DATA ONE / 1.0E0 / C RUN OPTIMAL DESIGN PROBLEM WITH NX = 50, NY = 50, AND LAMBDA = 0.008 NX = 50 NY = 50 CD LAMBDA = 0.008D0 CS LAMBDA = 0.008E0 N = NX * NY C COMPUTE THE STANDARD STARTING POINT. CD HX = ONE/DBLE(NX+1) CD HY = ONE/DBLE(NY+1) CS HX = ONE/FLOAT(NX+1) CS HY = ONE/FLOAT(NY+1) DO 20 J = 1, NY CD TEMP = DBLE(MIN(J,NY-J+1))*HY CS TEMP = FLOAT(MIN(J,NY-J+1))*HY DO 10 I = 1, NX K = NX*(J-1) + I CD X(K) = -(MIN(DBLE(MIN(I,NX-I+1))*HX,TEMP))**2 CS X(K) = -(MIN(FLOAT(MIN(I,NX-I+1))*HX,TEMP))**2 10 CONTINUE 20 CONTINUE C DEFINE THE SPARSITY STRUCTURE OF THE HESSIAN. CALL DODCSP(NX,NY,NZ,IRN,ICN) C SET THE DEFAULT VALUES OF THE PACKAGE. CALL STDFLT(N, TYPX, FSCALE, GRADTL, STEPTL, ILIM, STEPMX, * IPR, METHOD, GRDFLG, HSNFLG, NDIGIT, INFORM, MSG) CS GRADTL = 1.0E-3 CD GRADTL = 1.0D-5 GRDFLG = 2 C CALL THE SPARSE OPTIMIZER. CD CALL STUMCD(N,X,NZ,IRN,LIRN,ICN,LICN,DODCF,DODCG, CS CALL STUMCS(N,X,NZ,IRN,LIRN,ICN,LICN,DODCF,DODCG, * STDUSH,TYPX,FSCALE,GRADTL,STEPTL,ILIM,STEPMX,IPR, * METHOD,GRDFLG,HSNFLG,NDIGIT,MSG,XPLS,FPLS,GPLS, * HESS,WRK,LWRK,IWRK,LIWRK,TERMCD,HTV,INFORM) STOP END CUT HERE............ cat > odc.f <<'CUT HERE............' SUBROUTINE DODCF(N,X,F) C ********** C C PURPOSE C ------- C C THIS SUBROUTINE COMPUTES THE FUNCTION OF THE C OPTIMAL DESIGN WITH COMPOSITE MATERIALS PROBLEM. C C PARAMETERS C ---------- C C N IS THE DIMENSION OF THE PROBLEM C C NX IS AN INTEGER VARIABLE. C ON ENTRY NX IS THE NUMBER OF GRID POINTS IN THE FIRST C COORDINATE DIRECTION. C ON EXIT NX IS UNCHANGED. C C NY IS AN INTEGER VARIABLE. C ON ENTRY NY IS THE NUMBER OF GRID POINTS IN THE SECOND C COORDINATE DIRECTION. C ON EXIT NY IS UNCHANGED. C C X IS A DOUBLE PRECISION (REAL) ARRAY OF DIMENSION N = NX*NY. C C F IS A DOUBLE PRECISION (REAL) VARIABLE. C ON EXIT F IS SET TO THE FUNCTION EVALUATED AT X. C C LAMBDA IS A DOUBLE PRECISION (REAL) VARIABLE. C ON ENTRY LAMBDA IS THE LAGRANGE MULTIPLIER. C ON EXIT LAMBDA IS UNCHANGED. C C SUBPROGRAMS CALLED C C MINPACK-SUPPLIED ... DODCPS C C MINPACK-2 PROJECT. NOVEMBER 1993. C ARGONNE NATIONAL LABORATORY AND UNIVERSITY OF MINNESOTA. C BRETT M. AVERICK. MODIFIED BY ALI BOUARICHA ON OCTOBER 1994. C C ********** C INTEGER N, NX, NY CD DOUBLE PRECISION F, LAMBDA CD DOUBLE PRECISION X(N) CS REAL F, LAMBDA CS REAL X(N) COMMON/PARAM/NX,NY COMMON/OTHER/LAMBDA C C LOCAL VARIABLES C INTEGER I, J, K CD DOUBLE PRECISION MU1, MU2, ONE, P5, TWO, ZERO CD DOUBLE PRECISION AREA, DPSI, DVDX, DVDY, GRADV, HX, HXHY CD DOUBLE PRECISION HY, TEMP, T1, T2, V, VB, VL, VR, VT CS REAL MU1, MU2, ONE, P5, TWO, ZERO CS REAL AREA, DPSI, DVDX, DVDY, GRADV, HX, HXHY CS REAL HY, TEMP, T1, T2, V, VB, VL, VR, VT CD PARAMETER (ZERO=0.0D0,P5=0.5D0,ONE=1.0D0,TWO=2.0D0) CS PARAMETER (ZERO=0.0E0,P5=0.5E0,ONE=1.0E0,TWO=2.0E0) PARAMETER (MU1=ONE,MU2=TWO) CD INTRINSIC DBLE, SQRT CS INTRINSIC FLOAT, SQRT EXTERNAL DODCPS C C INITIALIZATION. C CD HX = ONE/DBLE(NX+1) CD HY = ONE/DBLE(NY+1) CS HX = ONE/FLOAT(NX+1) CS HY = ONE/FLOAT(NY+1) HXHY = HX*HY AREA = P5*HXHY C C COMPUTE THE BREAK POINTS. C T1 = SQRT(TWO*LAMBDA*MU1/MU2) T2 = SQRT(TWO*LAMBDA*MU2/MU1) F = ZERO C C COMPUTATION OF THE FUNCTION OVER THE LOWER C TRIANGULAR ELEMENTS. C DO 50 J = 0, NY DO 40 I = 0, NX K = NX*(J-1) + I V = ZERO VR = ZERO VT = ZERO IF (J .GE. 1 .AND. I .GE. 1) V = X(K) IF (I .LT. NX .AND. J .GT. 0) VR = X(K+1) IF (I .GT. 0 .AND. J .LT. NY) VT = X(K+NX) DVDX = (VR-V)/HX DVDY = (VT-V)/HY GRADV = DVDX**2 + DVDY**2 CALL DODCPS(GRADV,MU1,MU2,T1,T2,DPSI,0,LAMBDA) F = F + DPSI 40 CONTINUE 50 CONTINUE C C COMPUTATION OF THE FUNCTION OVER THE UPPER C TRIANGULAR ELEMENTS. C DO 70 J = 1, NY + 1 DO 60 I = 1, NX + 1 K = NX*(J-1) + I VB = ZERO VL = ZERO V = ZERO IF (I .LE. NX .AND. J .GT. 1) VB = X(K-NX) IF (I .GT. 1 .AND. J .LE. NY) VL = X(K-1) IF (I .LE. NX .AND. J .LE. NY) V = X(K) DVDX = (V-VL)/HX DVDY = (V-VB)/HY GRADV = DVDX**2 + DVDY**2 CALL DODCPS(GRADV,MU1,MU2,T1,T2,DPSI,0,LAMBDA) F = F + DPSI 60 CONTINUE 70 CONTINUE C C SCALE THE FUNCTION. C F = AREA*F C C INTEGRATE V OVER THE DOMAIN. C TEMP = ZERO DO 80 K = 1, NX*NY TEMP = TEMP + X(K) 80 CONTINUE F = F + HXHY*TEMP END C SUBROUTINE DODCG(N,X,FGRAD) C ********** C C PURPOSE C ------- C C THIS SUBROUTINE COMPUTES THE GRADIENT OF THE C OPTIMAL DESIGN WITH COMPOSITE MATERIALS PROBLEM. C C PARAMETERS C ---------- C C N IS THE DIMENSION OF THE PROBLEM C C NX IS AN INTEGER VARIABLE. C ON ENTRY NX IS THE NUMBER OF GRID POINTS IN THE FIRST C COORDINATE DIRECTION. C ON EXIT NX IS UNCHANGED. C C NY IS AN INTEGER VARIABLE. C ON ENTRY NY IS THE NUMBER OF GRID POINTS IN THE SECOND C COORDINATE DIRECTION. C ON EXIT NY IS UNCHANGED. C C X IS A DOUBLE PRECISION (REAL) ARRAY OF DIMENSION N = NX*NY. C C FGRAD IS A DOUBLE PRECISION (REAL) ARRAY OF DIMENSION N = NX*NY. C ON ENTRY FGRAD NEED NOT BE SPECIFIED. C ON EXIT FGRAD CONTAINS THE GRADIENT EVALUATED AT X. C C LAMBDA IS A DOUBLE PRECISION (REAL) VARIABLE. C ON ENTRY LAMBDA IS THE LAGRANGE MULTIPLIER. C ON EXIT LAMBDA IS UNCHANGED. C C SUBPROGRAMS CALLED C C MINPACK-SUPPLIED ... DODCPS C C MINPACK-2 PROJECT. NOVEMBER 1993. C ARGONNE NATIONAL LABORATORY AND UNIVERSITY OF MINNESOTA. C BRETT M. AVERICK. MODIFIED BY ALI BOUARICHA ON OCTOBER 1994. C C ********** C INTEGER N, NX, NY CD DOUBLE PRECISION X(N), FGRAD(N), LAMBDA CS REAL X(N), FGRAD(N), LAMBDA COMMON/PARAM/NX,NY COMMON/OTHER/LAMBDA C C LOCAL VARIABLES C INTEGER I, J, K CD DOUBLE PRECISION MU1, MU2, ONE, P5, TWO, ZERO CD DOUBLE PRECISION AREA, DPSIP, DVDX, DVDY, GRADV, HX, HXHY CD DOUBLE PRECISION HY, T1, T2, V, VB, VL, VR, VT CS REAL MU1, MU2, ONE, P5, TWO, ZERO CS REAL AREA, DPSIP, DVDX, DVDY, GRADV, HX, HXHY CS REAL HY, T1, T2, V, VB, VL, VR, VT CD PARAMETER (ZERO=0.0D0,P5=0.5D0,ONE=1.0D0,TWO=2.0D0) CS PARAMETER (ZERO=0.0E0,P5=0.5E0,ONE=1.0E0,TWO=2.0E0) PARAMETER (MU1=ONE,MU2=TWO) CD INTRINSIC DBLE, SQRT CS INTRINSIC FLOAT, SQRT EXTERNAL DODCPS C INITIALIZATION. CD HX = ONE/DBLE(NX+1) CD HY = ONE/DBLE(NY+1) CS HX = ONE/FLOAT(NX+1) CS HY = ONE/FLOAT(NY+1) HXHY = HX*HY AREA = P5*HXHY C C COMPUTE THE BREAK POINTS. C T1 = SQRT(TWO*LAMBDA*MU1/MU2) T2 = SQRT(TWO*LAMBDA*MU2/MU1) DO 30 K = 1, NX*NY FGRAD(K) = ZERO 30 CONTINUE C C COMPUTATION OF THE THE GRADIENT OVER THE LOWER C TRIANGULAR ELEMENTS. C DO 50 J = 0, NY DO 40 I = 0, NX K = NX*(J-1) + I V = ZERO VR = ZERO VT = ZERO IF (J .GE. 1 .AND. I .GE. 1) V = X(K) IF (I .LT. NX .AND. J .GT. 0) VR = X(K+1) IF (I .GT. 0 .AND. J .LT. NY) VT = X(K+NX) DVDX = (VR-V)/HX DVDY = (VT-V)/HY GRADV = DVDX**2 + DVDY**2 CALL DODCPS(GRADV,MU1,MU2,T1,T2,DPSIP,1,LAMBDA) IF (I .GE. 1 .AND. J .GE. 1) + FGRAD(K) = FGRAD(K) - TWO*(DVDX/HX+DVDY/HY)*DPSIP IF (I .LT. NX .AND. J .GT. 0) + FGRAD(K+1) = FGRAD(K+1) + TWO*(DVDX/HX)*DPSIP IF (I .GT. 0 .AND. J .LT. NY) + FGRAD(K+NX) = FGRAD(K+NX) + TWO*(DVDY/HY)*DPSIP 40 CONTINUE 50 CONTINUE C C COMPUTATION OF THE GRADIENT OVER THE UPPER C TRIANGULAR ELEMENTS. C DO 70 J = 1, NY + 1 DO 60 I = 1, NX + 1 K = NX*(J-1) + I VB = ZERO VL = ZERO V = ZERO IF (I .LE. NX .AND. J .GT. 1) VB = X(K-NX) IF (I .GT. 1 .AND. J .LE. NY) VL = X(K-1) IF (I .LE. NX .AND. J .LE. NY) V = X(K) DVDX = (V-VL)/HX DVDY = (V-VB)/HY GRADV = DVDX**2 + DVDY**2 CALL DODCPS(GRADV,MU1,MU2,T1,T2,DPSIP,1,LAMBDA) IF (I .LE. NX .AND. J .GT. 1) + FGRAD(K-NX) = FGRAD(K-NX) - TWO*(DVDY/HY)*DPSIP IF (I .GT. 1 .AND. J .LE. NY) + FGRAD(K-1) = FGRAD(K-1) - TWO*(DVDX/HX)*DPSIP IF (I .LE. NX .AND. J .LE. NY) + FGRAD(K) = FGRAD(K) + TWO*(DVDX/HX+DVDY/HY)*DPSIP 60 CONTINUE 70 CONTINUE C C INTEGRATE V OVER THE DOMAIN. C DO 90 K = 1, NX*NY FGRAD(K) = AREA*FGRAD(K) + HXHY 90 CONTINUE END C SUBROUTINE DODCPS(T,MU1,MU2,T1,T2,RESULT,OPTION,LAMBDA) C C ********** C C THIS SUBROUTINE COMPUTES THE FUNCTION PSI(T) AND THE SCALED C FUNCTIONS PSI'(T)/T AND PSI''(T)/T FOR THE OPTIMAL DESIGN C WITH COMPOSITE MATERIALS PROBLEM. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE DODCPS(T,MU1,MU2,T1,T2,RESULT,OPTION,LAMBDA) C C WHERE C C T IS A DOUBLE PRECISION (REAL) VARIABLE. C ON ENTRY T IS THE VARIABLE T C ON EXIT T IS UNCHANGED C C MU1 IS A DOUBLE PRECISION (REAL) VARIABLE. C ON ENTRY MU1 IS THE RECIPROCAL SHEAR MODULUS OF MATERIAL 1. C ON EXIT MU1 IS UNCHANGED. C C MU2 IS A DOUBLE PRECISION (REAL) VARIABLE. C ON ENTRY MU2 IS THE RECIPROCAL SHEAR MODULUS OF MATERIAL 2. C ON EXIT MU2 IS UNCHANGED. C C T1 IS A DOUBLE PRECISION (REAL) VARIABLE. C ON ENTRY T1 IS THE FIRST BREAKPOINT. C ON EXIT T1 IS UNCHANGED. C C T2 IS A DOUBLE PRECISION (REAL) VARIABLE. C ON ENTRY T2 IS THE SECOND BREAKPOINT. C ON EXIT T2 IS UNCHANGED. C C RESULT IS A DOUBLE PRECISION (REAL) VARIABLE. C ON ENTRY RESULT NEED NOT BE SPECIFIED. C ON EXIT RESULT IS SET ACCORDING TO TASK. C C OPTION IS AN INTEGER VARIABLE. C ON ENTRY OPTION SPECIFIES THE ACTION OF THE SUBROUTINE: C C IF OPTION = 0 THEN EVALUATE THE FUNCTION PSI(T). C IF OPTION = 1 THEN EVALUATE THE SCALED FUNCTION PSI'(T)/T. C IF OPTION = 2 THEN EVALUATE THE SCALED FUNCTION PSI''(T)/T. C C ON OPTION TASK IS UNCHANGED. C C LAMBDA IS A DOUBLE PRECISION (REAL) VARIABLE C ON ENTRY LAMBDA IS THE LAGRANGE MULTIPLIER. C ON EXIT LAMBDA IS UNCHANGED. C C MINPACK-2 PROJECT. NOVEMBER 1993. C ARGONNE NATIONAL LABORATORY AND UNIVERSITY OF MINNESOTA. C BRETT M. AVERICK. MODIFIED BY ALI BOUARICHA ON OCTOBER 1994. C C ********** C INTEGER OPTION CD DOUBLE PRECISION T, MU1, MU2, T1, T2, RESULT, LAMBDA CS REAL T, MU1, MU2, T1, T2, RESULT, LAMBDA C C LOCAL VARIABLES C CD DOUBLE PRECISION P25, P5, ZERO CD PARAMETER (ZERO=0.0D0,P25=0.25D0,P5=0.5D0) CS REAL P25, P5, ZERO CS PARAMETER (ZERO=0.0E0,P25=0.25E0,P5=0.5E0) CD DOUBLE PRECISION SQRTT CS REAL SQRTT INTRINSIC SQRT C SQRTT = SQRT(T) IF (OPTION .EQ. 0) THEN IF (SQRTT .LE. T1) THEN RESULT = P5*MU2*T ELSE IF (SQRTT .GT. T1 .AND. SQRTT .LT. T2) THEN RESULT = MU2*T1*SQRTT - LAMBDA*MU1 ELSE IF (SQRTT .GE. T2) THEN RESULT = P5*MU1*T + LAMBDA*(MU2-MU1) END IF ELSE IF (OPTION .EQ. 1) THEN IF (SQRTT .LE. T1) THEN RESULT = P5*MU2 ELSE IF (SQRTT .GT. T1 .AND. SQRTT .LT. T2) THEN RESULT = P5*MU2*T1/SQRTT ELSE IF (SQRTT .GE. T2) THEN RESULT = P5*MU1 END IF ELSE IF (OPTION .EQ. 2) THEN IF (SQRTT .LE. T1) THEN RESULT = ZERO ELSE IF (SQRTT .GT. T1 .AND. SQRTT .LT. T2) THEN RESULT = -P25*MU2*T1/(SQRTT*T) ELSE IF (SQRTT .GE. T2) THEN RESULT = ZERO END IF END IF END C SUBROUTINE DODCSP(NX,NY,NNZ,INDROW,INDCOL) C C ********** C C SUBROUTINE DODCSP C C THIS SUBROUTINE DEFINES THE SPARSITY STRUCTURE OF THE HESSIAN C MATRIX FOR THE OPTIMAL DESIGN WITH COMPOSITES PROBLEM. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE DODCSP(NX,NY,NNZ,INDROW,INDCOL) C C WHERE C C NX IS AN INTEGER VARIABLE. C ON ENTRY NX IS THE NUMBER OF GRID POINTS IN THE FIRST C COORDINATE DIRECTION. C ON EXIT NX IS UNCHANGED. C C NY IS AN INTEGER VARIABLE. C ON ENTRY NY IS THE NUMBER OF GRID POINTS IN THE SECOND C COORDINATE DIRECTION. C ON EXIT NY IS UNCHANGED. C C NNZ IS AN INTEGER VARIABLE. C ON ENTRY NNZ NEED NOT BE SPECIFIED. C ON EXIT NNZ IS SET TO THE NUMBER OF NONZEROS IN THE C LOWER TRIANGLE OF THE HESSIAN MATRIX. C C INDROW IS AN INTEGER ARRAY OF DIMENSION AT LEAST NNZ. C ON ENTRY INDROW NEED NOT BE SPECIFIED. C ON EXIT INDROW CONTAINS THE ROW INDICES OF THE NONZEROS C IN THE LOWER TRIANGLE OF THE HESSIAN MATRIX. C C INDCOL IS AN INTEGER ARRAY OF DIMENSION AT LEAST NNZ. C ON ENTRY INDCOL NEED NOT BE SPECIFIED. C ON EXIT INDCOL CONTAINS THE COLUMN INDICES OF THE NONZEROS C IN THE LOWER TRIANGLE OF THE HESSIAN MATRIX. C C MINPACK-2 PROJECT. NOVEMBER 1993. C ARGONNE NATIONAL LABORATORY AND UNIVERSITY OF MINNESOTA. C BRETT M. AVERICK. MODIFIED BY ALI BOUARICHA ON OCTOBER 1994. C C ********** C INTEGER NX, NY, NNZ INTEGER INDROW(*), INDCOL(*) C C LOCAL VARIABLES C INTEGER I, J C C COMPUTE THE SPARSITY STRUCTURE. C NNZ = 0 DO 20 J = 1, NY DO 10 I = 1, NX NNZ = NNZ + 1 INDROW(NNZ) = (J-1)*NX + I INDCOL(NNZ) = (J-1)*NX + I IF (I .NE. NX) THEN NNZ = NNZ + 1 INDROW(NNZ) = (J-1)*NX + I + 1 INDCOL(NNZ) = (J-1)*NX + I END IF IF (J .NE. NY) THEN NNZ = NNZ + 1 INDROW(NNZ) = (J-1)*NX + I + NX INDCOL(NNZ) = (J-1)*NX + I IF (I .NE. 1) THEN NNZ = NNZ + 1 INDROW(NNZ) = (J-1)*NX + I + NX - 1 INDCOL(NNZ) = (J-1)*NX + I END IF END IF 10 CONTINUE 20 CONTINUE END CUT HERE............ cat > output2 <<'CUT HERE............' STDRUO GRADIENT FLAG = 2 STDRUO HESSIAN FLAG = 0 STDRUO METHOD = 1 STDRUO ITERATION LIMIT = 500 STDRUO MACHINE EPSILON = 0.2220446049250E-15 STDRUO STEP TOLERANCE = 0.3666852862501E-10 STDRUO GRADIENT TOLERANCE = 0.1000000000000E-04 STDRUO MAXIMUM STEP SIZE = 0.3295193424014E+04 --------------------------------------------- STRSLT ITERATION K = 0 STRSLT FUNCTION AT X(K) STRSLT 0.4808237011427E-01 STRSLT SCALED GRADIENT AT X(K) STRSLT 0.3729334871203E-01 --------------------------------------------- STCHKS RELATIVE GRADIENT CLOSE TO ZERO STCHKS CURRENT ITERATE IS PROBABLY SOLUTION --------------------------------------------- STRSLT ITERATION K = 20 STRSLT FUNCTION AT X(K) STRSLT -0.1135947433813E-01 STRSLT SCALED GRADIENT AT X(K) STRSLT 0.2847027636254E-06 --------------------------------------------- STRSLT NUMBER OF FUNCTION EVALUATIONS 78 STRSLT NUMBER OF GRADIENT EVALUATIONS 21 STRSLT NUMBER OF HESSIAN EVALUATIONS 20 CUT HERE............ cat > stenmin.f <<'CUT HERE............' C C STENMIN: A SOFTWARE PACKAGE FOR LARGE, SPARSE UNCONSTRAINED OPTIMIZATION C USING TENSOR METHODS. C C C AUTHOR: C ------ C ALI BOUARICHA C ARGONNE NATIONAL LABORATORY C MCS DIVISION C 9700 SOUTH CASS AVENUE C ARGONNE, IL 60439 C E-MAIL: BOUARICH@MCS.ANL.GOV C C DATE: OCTOBER, 1994 C ---- C C PURPOSE: C ------- C C STENMIN MINIMIZES AN UNCONSTRAINED NONLINEAR FUNCTION IN N UNKNOWNS C WHERE THE HESSIAN IS LARGE AND SPARSE USING A NEW CLASS OF METHODS C CALLED TENSOR METHODS. THE SOFTWARE ALLOWS THE USER TO SELECT BETWEEN C A TENSOR METHOD AND A STANDARD METHOD BASED UPON A QUADRATIC MODEL. C THE TENSOR METHOD MODELS THE OBJECTIVE FUNCTION BY A FOURTH-ORDER C MODEL, WHERE THE THIRD- AND FOURTH-ORDER TERMS ARE CHOSEN SUCH THAT C THE EXTRA COST OF FORMING AND SOLVING THE MODEL IS SMALL. THE TEST C RESULTS INDICATE THAT, IN GENERAL, THE TENSOR METHOD IS OFTEN MORE C EFFICIENT AND MORE RELIABLE THAN THE STANDARD nEWTON METHOD FOR C SOLVING LARGE, SPARSE UNCONSTRAINED OPTIMIZATION PROBLEMS. IT IS C ESPECIALLY USEFUL ON PROBLEMS WHERE THE HESSIAN AT THE SOLUTION IS C SINGULAR. C TWO INTERFACES ARE PROVIDED WITH STENMIN. IF THE USER WISHES TO C USE ALL THE DEFAULTS OPTIONS PROVIDED BY THE PACKAGE, THEN HE (OR SHE) C SHOULD CALL STUMSD (STUMSS IF SINGLE-PRECISION IS USED). THE OTHER C INTERFACE, STUMCD (STUMCS IF SINGLE-PRECISION IS USED), REQUIRES THE C USER TO SUPPLY ALL PARAMETERS. THE USER MAY SPECIFY SELECTED PARAMETERS C ONLY BY FIRST INVOKING THE SUBROUTINE STDFLT, WHICH SETS ALL PARAMETERS C TO THEIR DEFAULT VALUES, AND THEN OVERRIDING ONLY THE DESIRED VALUES. C C LIST OF SUBROUTINE AND FUNCTION NAMES CALLED BY STENMIN: C C STCHKG,STCHKH,STCHKI,STCHKS,STCZ3P,STDFLT,STDRUO,STDUGR,STDUSH, C STFDGR,STFDHS,STFTSM,STHMUV,STLSCH,STMA27,STMNEC,STMSLV,STRSLT, C STSORT,STUMCD,STUMSD. C C PACKAGES CALLED BY STENMIN: C C MA27 (I. S. DUFF AND J. K. REID, "MA27: A SET OF FORTRAN SUBROUTINES FOR C SOLVING SPARSE SYMMETRIC SETS OF LINEAR EQUATIONS", TECH. REP. R-10533, C AERE HARWELL LABORATORY, HARWELL, UK, 1983.) C C DSSM (T. F. COLEMAN, B. S. GARBOW, and J. J. MORE', "SOFTWARE FOR ESTIMATING C SPARSE HESSIAN MATRICES", ACM TRANS. MATH. SOFTW., 11: 363-377,1985.) C C DSYPRC FROM LANCELOT (A. R. CONN, N. I. M. GOULD, AND C PH. L. TOINT, C "LANCELOT", SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, SPRINGER-VERLAG, C 1992.) C C DPMEPS (W. J. Cody, "MACHAR: A SUBROUTINE TO DYNAMICALLY DETERMINE MACHINE C PARAMETERS", ACM TRANS. MATH. SOFTW., 14: 303-311, 1988.) C C BLAS CALLED BY STENMIN: C C LEVEL 1 BLAS: DCOPY,DDOT,DNRM2,DSCAL C SUBROUTINE STCHKG(N,X,FCN,F,G,TYPX,RNF,ANALTL,WRK1,MSG,IPR) C C PURPOSE: C ------- C C THIS ROUTINE CHECKS THE USER'S ANALYTIC GRADIENT AGAINST A C FINITE-DIFFERENCE GRADIENT. C C PARAMETERS: C ---------- C C N --> DIMENSION OF PROBLEM C X --> CURRENT ITERATE C FCN --> NAME OF SUBROUTINE THAT EVALUATES THE C OPTIMIZATION FUNCTION C F --> FUNCTION VALUE AT X C G --> GRADIENT VALUE AT X C TYPX --> TYPICAL SIZE FOR EACH COMPONENT OF X C RNF --> RELATIVE NOISE IN FCN C ANALTL --> TOLERANCE FOR COMPARISON OF ESTIMATED AND ANALYTICAL C GRADIENTS C WRK1 --> WORKSPACE. IT IS USED TO HOLD THE FINITE-DIFFERENCE C APPROXIMATION OF THE GRADIENT AT X C MSG <-- MESSAGE OR ERROR CODE. IF PROBABLE CODING ERROR OF C GRADIENT IS DETECTED THEN MSG = -8 C IPR --> DEVICE TO WHICH TO SEND OUTPUT C INTEGER N,MSG,IPR CD DOUBLE PRECISION X(N),F,G(N),TYPX(N),RNF,ANALTL,WRK1(N) CS REAL X(N),F,G(N),TYPX(N),RNF,ANALTL,WRK1(N) C C LOCAL VARIABLES C INTEGER KER,I CD DOUBLE PRECISION GS,ONE CS REAL GS,ONE EXTERNAL FCN INTRINSIC ABS,MAX CD DATA ONE/1.0D0/ CS DATA ONE/1.0E0/ C C COMPUTE FIRST-ORDER FINITE DIFFERENCE GRADIENT AND COMPARE TO C ANALYTIC GRADIENT C CALL STFDGR(N,X,FCN,F,RNF,WRK1) KER = 0 DO 10 I=1,N GS=MAX(ABS(F),ONE)/MAX(ABS(X(I)),ONE/TYPX(I)) IF(ABS(G(I)-WRK1(I)).GT.MAX(ABS(G(I)),GS)*ANALTL) KER=1 10 CONTINUE IF(KER.EQ.0) GO TO 20 WRITE(IPR,901) WRITE(IPR,902) (I,G(I),WRK1(I),I=1,N) MSG = -8 20 CONTINUE RETURN 901 FORMAT(47H STCHKG PROBABLE ERROR IN CODING OF ANALYTIC, * 19H GRADIENT FUNCTION./ * 16H STCHKG COMP,12X,8HANALYTIC,12X,8HESTIMATE) 902 FORMAT(11H STCHKG ,I5,3X,E20.13,3X,E20.13) END SUBROUTINE STCHKH(N,X,NPAIRS,NZ,IRN,LIRN,ICN,LICN,FCN,UGR, * SCALE,TYPX,GRDFLG,F,G,RNF,ANALTL,FCALL, * CHECKH,LISTP,NGRP,IPNTR,JPNTR,IWA,LIWA, * FHESD,XD,ETA,HESS,IPR,MSG) C C PURPOSE C ------- C C THIS ROUTINE CHECKS THE USER'S ANALYTIC HESSIAN AGAINST A C FINITE-DIFFERENCE HESSIAN. C C PARAMETERS: C ---------- C C N --> NUMBER OF VARIABLES IN THE PROBLEM C X --> CURRENT ITERATE C NPAIRS --> A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF (IRN,ICN) PAIRS USED TO DESCRIBE THE SPARSITY C PATTERN OF THE HESSIAN MATRIX C NZ <-- ACTUAL NUMBER OF NONZEROS STORED IN THE LOWER C HALF OF THE HESSIAN MATRIX C IRN <-->ROW INDEX OF EACH NONZERO STORED IN THE LOWER C OR UPPER HALF OF THE HESSIAN MATRIX ON ENTRY. ON EXIT, C IT HOLDS THE ROW INDEX OF EACH NONZERO STORED IN THE C LOWER HALF OF THE HESSIAN MATRIX C ICN <-->COLUMN INDEX OF THE NONZEROS STORED IN THE LOWER C OR UPPER HALF OF THE HESSIAN MATRIX ON ENTRY. ON EXIT, C IT HOLDS THE COLUMN INDEX OF EACH NONZERO STORED IN THE C LOWER HALF OF THE HESSIAN MATRIX C FCN --> THE NAME OF A USER SUPPLIED SUBROUTINE THAT EVALUATES C THE FUNCTION F AT AN ARBITRARY VECTOR X C UGR --> THE NAME OF A USER SUPPLIED SUBROUTINE THAT RETURNS IN C THE VALUE OF THE GRADIENT C SCALE --> LOGICAL VARIABLE INDICATING WHETHER OR NOT SCALING C OF THE VARIABLES IS PERFORMED C TYPX --> TYPICAL SIZE OF THE COMPONENTS OF X C GRDFLG --> INTEGER FLAG DESIGNATING WHETHER OR NOT ANALYTIC C GRADIENT HAS BEEN SUPPLIED BY THE USER C F --> FUNCTION VALUE AT CURRENT ITERATE C G --> GRADIENT VALUE AT CURRENT ITERATE C RNF --> RELATIVE NOISE IN FUNCTION FCN C ANALTL --> TOLERANCE FOR COMPARISON OF ESTIMATED AND ANALYTICAL C GRADIENTS C FCALL --> A LOGICAL FLAG. IF THE HESSIAN MATRIX IS EVALUATED FOR C THE FIRST TIME, THEN FCALL IS SET TO TRUE; OTHERWISE C IT IS SET TO FALSE C CHECKH --> A LOGICAL FLAG. IF IT IS SET TO TRUE THEN THE USER'S C ANALYTIC HESSIAN IS CHECKED AGAINST THE FINITE C DIFFERENCE APPROXIMATION ONE C LISTP,NGRP,IPNTR,JPNTR,IWA --> WORSPACE FOR THE STFDHS SUBROUTINE C (SEE SUBROUTINE STFDHS FOR MORE DETAIL) C LIWA --> LENGTH OF ARRAY IWA C FHESD,XD,ETA --> WORKSPACE FOR THE STFDHS SUBROUTINE C HESS <--> ANALYTIC HESSIAN ON INPUT. ON OUTPUT, IT ALSO HOLDS C A FINITE-DIFFERENCE APPROXIMATION OF THE HESSIAN C (FROM NZ+1 TO 2*NZ) C IPR --> DEVICE TO WHICH TO SEND OUTPUT C MSG <-- MESSAGE OR ERROR CODE. IF PROBABLE CODING ERROR OF C HESSIAN IS DETECTED THEN MSG = -9 C INTEGER N,NPAIRS,NZ,GRDFLG,LIWA,IPR,MSG,LIRN,LICN INTEGER IRN(LIRN),ICN(LICN),LISTP(N),NGRP(N),IPNTR(N+1) INTEGER JPNTR(N+1),IWA(LIWA) CD DOUBLE PRECISION X(N),TYPX(N),F,G(N),RNF,ANALTL,FHESD(N) CD DOUBLE PRECISION XD(N),ETA(N),HESS(LICN) CS REAL X(N),TYPX(N),F,G(N),RNF,ANALTL,FHESD(N) CS REAL XD(N),ETA(N),HESS(LICN) LOGICAL SCALE,FCALL,CHECKH C C LOCAL VARIABLES C INTEGER IST,IEND,KER,I,IR,J,JP,K,MAXGRP,MINGRP CD DOUBLE PRECISION HS,ONE CS REAL HS,ONE EXTERNAL FCN,UGR INTRINSIC ABS,MAX,MIN CD DATA ONE/ 1.0D+0 / CS DATA ONE/ 1.0E+0 / C C CHECK USER'S INPUT DATA. C C IF NPAIRS IS NOT POSITIVE THEN SET MSG TO -2 AND ABORT. C IF (NPAIRS .LE. 0) THEN WRITE(IPR,902) NPAIRS MSG = -4 RETURN ENDIF C C IF THE K-TH ELEMENT OF IRN OR THE K-TH ELEMENT OF ICN IS C NOT AN INTEGER BETWEEN 1 AND N THEN SET MSG TO -4 AND ABORT. C DO 10 K = 1, NPAIRS IF (IRN(K) .LT. 1 .OR. IRN(K) .GT. N .OR. * ICN(K) .LT. 1 .OR. ICN(K) .GT. N) THEN WRITE(IPR,903) K,IRN(K),K,ICN(K) MSG = -5 RETURN ENDIF 10 CONTINUE C C GENERATE THE SPARSITY PATTERN FOR THE LOWER C TRIANGULAR PART OF HESSIAN. C DO 20 K = 1, NPAIRS I = IRN(K) J = ICN(K) IRN(K) = MAX(I,J) ICN(K) = MIN(I,J) 20 CONTINUE C C SORT THE DATA STRUCTURE BY COLUMNS C CALL STSORT(N,NPAIRS,HESS,IRN,ICN,JPNTR,IWA) C C COMPRESS THE DATA AND DETERMINE THE NUMBER OF NONZERO C ELEMENTS IN THE LOWER TRIANGULAR PART OF HESSIAN. C DO 60 I = 1, N IWA(I) = 0 60 CONTINUE NZ = 0 DO 80 J = 1, N K = NZ DO 70 JP = JPNTR(J), JPNTR(J+1)-1 IR = IRN(JP) IF (IWA(IR) .NE. J) THEN NZ = NZ + 1 IRN(NZ) = IR ICN(NZ) = J IWA(IR) = J ENDIF 70 CONTINUE JPNTR(J) = K + 1 80 CONTINUE JPNTR(N+1) = NZ + 1 IF(NZ .NE. NPAIRS) THEN MSG = -7 WRITE(IPR,904) NPAIRS-NZ RETURN ENDIF C IF(.NOT. CHECKH) RETURN C C COMPUTE FINITE-DIFFERENCE APPROXIMATION OF THE HESSIAN C CALL STFDHS(N,NZ,X,IRN,ICN,G,RNF,GRDFLG, * FCN,UGR,SCALE,TYPX,FCALL,LISTP,NGRP, * IPNTR,JPNTR,IWA,LIWA,FHESD,XD,ETA, * MAXGRP,MINGRP,HESS(NZ+1)) C C COMPARE USER'S ANALYTIC HESSIAN TO THE FINITE-DIFFERENCE C APPROXIMATION C KER = 0 DO 100 J = 1, N HS = MAX(ABS(G(J)),ONE)/MAX(ABS(X(J)),ONE/TYPX(J)) IST = JPNTR(J) IEND = JPNTR(J+1) DO 95 I = IST,IEND - 1 IF(ABS(HESS(I)-HESS(NZ+I)) .GT. MAX(ABS(HESS(I)),HS) * * ANALTL) KER = 1 95 CONTINUE 100 CONTINUE C IF(KER .EQ. 0) GO TO 120 WRITE(IPR,905) DO 110 I = 1, NZ WRITE(IPR,906) IRN(I),ICN(I),HESS(I),HESS(NZ+I) 110 CONTINUE MSG = -9 120 CONTINUE RETURN C 902 FORMAT(' STCHKH ILLEGAL NUMBER OF NONZEROS NZ=',I10) 903 FORMAT(' STCHKH ILLEGAL ROW OR COLUMN INDEX IN ANALYTIC', * ' HESSIAN'/, * ' STCHKH IRN( ',I6,' )=',I6,2X,' ICN( ',I6,' )=',I6) 904 FORMAT(' STCHKH ',I5,' REDUNDANT ENTRIES IN SPARSITY', * ' PATTERN ENCOUNTERED') 905 FORMAT(' STCHKH PROBABLE ERROR IN CODING OF ANALYTIC', + ' HESSIAN FUNCTION'/ + ' STCHKH IRN ICN',14X,'ANALYTIC',14X,'(ESTIMATE)') 906 FORMAT(' STCHKH ',2I5,2X,E20.13,2X,'(',E20.13,')') END SUBROUTINE STCHKI(N,NPAIRS,NZ,X,IRN,LIRN,ICN,LICN,LIWRK,LWRK, * SCALE,TYPX,FSCALE,GRADTL,STEPTL,ILIM,NDIGIT, * EPS,METHOD,GRDFLG,HSNFLG,STEPMX,IWA,LIWA, * JPNTR,MSG,IPR) C C PURPOSE: C ------- C C THIS ROUTINE CHECKS THE INPUT FOR REASONABLENESS. C C PARAMETERS: C ---------- C N --> DIMENSION OF PROBLEM C NPAIRS --> A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF (IRN,ICN) PAIRS USED TO DESCRIBE THE SPARSITY C PATTERN OF THE HESSIAN MATRIX C NZ <-- NUMBER OF NONZEROS IN THE LOWER TRIANGULAR PART C OF THE HESSIAN MATRIX C X --> ON ENTRY, ESTIMATE TO ROOT OF FCN C IRN <--> AN INTEGER ARRAY OF LENGTH NPAIRS. ON INPUT IRN MUST C CONTAIN THE ROW INDICES OF THE NONZERO ELEMENTS IN C THE LOWER TRIANGULAR PART OF THE HESSIAN. ON OUTPUT C IRN IS PERMUTED SO THAT THE CORRESPONDING COLUMN C INDICES ARE IN NON-DECREASING ORDER. THE COLUMN C INDICES CAN BE RECOVERED FROM THE ARRAY JPNTR C LIRN --> LENGTH OF ARRAY IRN C ICN <--> AN INTEGER ARRAY OF LENGTH NPAIRS. ON INPUT ICN MUST C CONTAIN THE COLUMN INDICES OF THE NONZERO ELEMENTS C IN THE LOWER TRIANGULAR PART OF THE HESSIAN. ON C OUTPUT ICN IS PERMUTED SO THAT THE CORRESPONDING C ROW INDICES ARE IN NON-DECREASING ORDER C LICN --> LENGTH OF ARRAY ICN C LIWRK --> LENGTH OF WORKSPACE ARRAY IWRK C LWRK --> LENGTH OF WORKSPACE ARRAY WRK C SCALE <-- LOGICAL VARIABLE INDICATING WHETHER OR NOT SCALING C OF THE X VARIABLES IS PERFORMED C TYPX <--> TYPICAL SIZE OF EACH COMPONENT OF X C FSCALE <--> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION FCN C GRADTL <--> TOLERANCE AT WHICH GRADIENT CONSIDERED CLOSE C ENOUGH TO ZERO TO TERMINATE ALGORITHM C STEPTL <--> TOLERANCE AT WHICH STEP LENGTH CONSIDERED CLOSE C ENOUGH TO ZERO TO TERMINATE ALGORITHM C ILIM <--> MAXIMUM NUMBER OF ALLOWABLE ITERATIONS C NDIGIT <--> NUMBER OF GOOD DIGITS IN OPTIMIZATION FUNCTION FCN C EPS --> MACHINE PRECISION C METHOD <--> METHOD TO BE USED (NEWTON OR TENSOR) C GRDFLG <--> =1 IF ANALYTIC GRADIENT SUPPLIED C HSNFLG <--> =1 IF ANALYTIC HESSIAN SUPPLIED C STEPMX <--> MAXIMUM STEP SIZE C IWA --> WORKSPACE ARRAY C LIWA <-- ACTUAL LENGTH OF ARRAY IWA ON EXIT C JPNTR <-- AN INTEGER OUTPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE ROW INDICES IN IRN. C THE ROW INDICES FOR COLUMN J ARE C IRN(K), K = JPNTR(J),...,JPNTR(J+1)-1. NOTE THAT C JPNTR(N+1)-1 IS THEN THE NUMBER OF NONZERO ELEMENTS C IN THE LOWER TRIANGULAR PART OF THE MATRIX HESSIAN C MSG <--> AN INTEGER VARIABLE THAT THE USER MAY SET ON INPUT TO C INHIBIT CERTAIN AUTOMATIC CHECKS OR OVERRIDE CERTAIN C DEFAULT CHARACTERISTICS OF THE PACKAGE. ON INPUT: C MSG = 0 : NO OUTPUT WILL BE PRODUCED C MSG = 1 : PRINT THE INPUT STATE, THE FINAL RESULTS, C AND THE STOPPING CONDITIONS C MSG = 2 : PRINT THE INTERMEDIATE RESULTS, THAT IS, THE C INPUT STATE, THE VALUES OF THE OBJECTIVE C FUNCTION AND THE SCALED GRADIENT AT EACH C ITERATION, AND THE FINAL RESULTS INCLUDING C THE STOPPING CONDITIONS AND THE NUMBER OF C FUNCTION, GRADIENT, AND HESSIAN C C ON OUTPUT, IF THE PROGRAM HAS TERMINATED BECAUSE OF C ERRONEOUS INPUT, MSG CONTAINS AN ERROR CODE INDICATING C THE REASON: C C MSG = -1 : ILLEGAL DIMENSION N; N <= 0. THE PROGRAM ABORTS C MSG = -2 : ILLEGAL LENGTH OF LIRN OR LICN; LIRN <= 0 OR C LICN <= 0. THE PROGRAM ABORTS C MSG = -3 : ILLEGAL LENGTH OF LIWRK OR LWRK; LIWRK < 2*LIRN C +12*N+2 OR LWRK < 7*N. THE PROGRAM ABORTS C MSG = -4 : ILLEGAL NUMBER OF NONZEROS NZ; NZ <= 0. THE C PROGRAM ABORTS C MSG = -5 : THE K-TH ELEMENT OF IRN OR THE K-TH ELEMENT OF C ICN IS NOT AN INTEGER BETWEEN 1 AND N; (IRN(K) C < 1 OR IRN(K) > N) OR (ICN(K) < 1 OR ICN(K) > N). C THE PROGRAM ABORTS C MSG = -6 : THE K-TH DIAGONAL ELEMENT IS NOT IN THE SPARSITY C PATTERN. THIS IS CHECKED ONLY IF HSNFLG = 0. C THE PROGRAM ABORTS C MSG = -7 : REDUNDANT ENTRIES IN SPARSITY PATTERN WAS C ENCOUNTERED. WHEN HSNFLG = 1 OR HSNFLG = 2, THE C PROGRAM ABORTS. WHEN HSNFLG = 0, THE PROGRAM C ELIMINATES THE REDUNDANT ENTRIES AND CONTINUES C THE EXECUTION C MSG = -8 : PROBABLE CODING ERROR IN THE USER'S ANALYTIC C GRADIENT ROUTINE UGR. THE PROGRAM ABORTS. (THIS C CHECK CAN BE OVERRIDDEN BY SETTING GRDFLG = 2.) C (SEE SUBROUTINE STCHKG) C MSG = -9 : PROBABLE CODING ERROR IN THE USER'S ANALYTIC C HESSIAN ROUTINE USH. THE PROGRAM ABORTS. (THIS C CHECK CAN BE OVERRIDDEN BY SETTING HSNFLG = 2.) C (SEE SUBROUTINE STCHKH) C C IPR --> DEVICE TO WHICH TO SEND OUTPUT C INTEGER N,NPAIRS,NZ,LIRN,LICN,LIWRK,LWRK,ILIM,NDIGIT,METHOD INTEGER GRDFLG,HSNFLG,LIWA,MSG,IPR INTEGER IRN(LIRN),ICN(LICN),IWA(N),JPNTR(N+1) CD DOUBLE PRECISION X(N),TYPX(N),FSCALE,GRADTL,STEPTL,STEPMX CS REAL X(N),TYPX(N),FSCALE,GRADTL,STEPTL,STEPMX LOGICAL SCALE C C LOCAL VARIABLES C INTEGER I,IR,J,JP,K CD DOUBLE PRECISION STPSIZ,TEMP,EPS,ZERO,ONE,THREE,TEN,THOUS CS REAL STPSIZ,TEMP,EPS,ZERO,ONE,THREE,TEN,THOUS INTRINSIC LOG10, MAX, MIN, SQRT CD DATA ZERO, ONE, THREE, TEN, THOUS / 0.0D+0, 1.0D+0, 3.0D+0, CD * 10.0D+0, 1000.0D+0/ CS DATA ZERO, ONE, THREE, TEN, THOUS / 0.0E+0, 1.0E+0, 3.0E+0, CS * 10.0E+0, 1000.0E+0 / C C CHECK THAT PARAMETERS ONLY TAKE ON ACCEPTABLE VALUES. C IF NOT, SET THEM TO DEFAULT VALUES. C C CHECK DIMENSION OF THE PROBLEM C IF(N .LE. 0) THEN WRITE(IPR,901) N MSG = -1 RETURN ENDIF C IF(LIRN .LE. 0 .OR. LICN .LE. 0) THEN WRITE(IPR,902) LIRN,LICN MSG = -2 RETURN ENDIF C C IF LIWRK IS LESS THAN 2*LIRN+12*N+2 OR LWRK LESS THAN 7*N C SET MSG = -3 AND ABORT. C IF (LIWRK .LT. 2*LIRN+12*N+2 .OR. LWRK .LT. 7*N) THEN WRITE(IPR,903) LIWRK,LWRK MSG = -3 RETURN ENDIF C C SET THE LENGTH OF LIWA C LIWA = 2*LIRN + 8*N C IF(HSNFLG .EQ. 0) THEN C C IF NPAIRS IS NOT POSITIVE THEN SET MSG TO -4 AND ABORT. C IF (NPAIRS .LE. 0) THEN WRITE(IPR,904) NPAIRS MSG = -4 RETURN ENDIF C C IF THE K-TH ELEMENT OF IRN OR THE K-TH ELEMENT OF ICN IS C NOT AN INTEGER BETWEEN 1 AND N THEN SET MSG TO -5 AND ABORT. C DO 10 K = 1, NPAIRS IF (IRN(K) .LT. 1 .OR. IRN(K) .GT. N .OR. * ICN(K) .LT. 1 .OR. ICN(K) .GT. N) THEN WRITE(IPR,905) K,IRN(K),K,ICN(K) MSG = -5 RETURN ENDIF 10 CONTINUE C C GENERATE THE SPARSITY PATTERN FOR THE LOWER C TRIANGULAR PART OF HESSIAN. C DO 20 K = 1, NPAIRS I = IRN(K) J = ICN(K) IRN(K) = MAX(I,J) ICN(K) = MIN(I,J) 20 CONTINUE C C IF HSNFLG = 0 AND THE K-TH DIAGONAL ELEMENT IS NOT IN THE C SPARSITY PATTERN THEN ABORT. C I = 0 DO 30 K = 1, N IWA(K) = 0 30 CONTINUE DO 40 K = 1, NPAIRS IF(IRN(K) .EQ. ICN(K)) IWA(IRN(K)) = 1 40 CONTINUE DO 50 K = 1, N IF(IWA(K) .NE. 1) THEN WRITE(IPR,906) K MSG = -6 RETURN ENDIF 50 CONTINUE C C SORT THE DATA STRUCTURE BY COLUMNS. C CALL SRTDAT(N,NPAIRS,IRN,ICN,JPNTR,IWA) C C COMPRESS THE DATA AND DETERMINE THE NUMBER OF NONZERO C ELEMENTS IN THE LOWER TRIANGULAR PART OF HESSIAN. C DO 60 I = 1, N IWA(I) = 0 60 CONTINUE NZ = 0 DO 80 J = 1, N K = NZ DO 70 JP = JPNTR(J), JPNTR(J+1)-1 IR = IRN(JP) IF (IWA(IR) .NE. J) THEN NZ = NZ + 1 IRN(NZ) = IR ICN(NZ) = J IWA(IR) = J ENDIF 70 CONTINUE JPNTR(J) = K + 1 80 CONTINUE JPNTR(N+1) = NZ + 1 C ENDIF C C CHECK METHOD USED C IF(METHOD.NE.0) METHOD = 1 C C CHECK GRADIENT FLAG C IF(GRDFLG.NE.2 .AND. GRDFLG.NE.1) GRDFLG=0 C C CHECK HESSIAN FLAG C IF(HSNFLG.NE.2 .AND. HSNFLG.NE.1) HSNFLG=0 C C CHECK MSG C IF(MSG.NE.2 .AND. MSG.NE.0) MSG=1 C C CHECK SCALE MATRIX C DO 90 I=1,N IF(TYPX(I).EQ.ZERO) TYPX(I)=ONE IF(TYPX(I).LT.ZERO) TYPX(I)=-TYPX(I) 90 CONTINUE C C SET THE LOGICAL VARIABLE SCALE C SCALE = .FALSE. DO 95 I = 1,N IF(TYPX(I) .NE. ONE) THEN SCALE = .TRUE. GO TO 100 ENDIF 95 CONTINUE 100 CONTINUE C C CHECK MAXIMUM STEP SIZE C IF (STEPMX .GT. ZERO) GO TO 110 STPSIZ = ZERO DO 105 I = 1, N STPSIZ = STPSIZ + X(I)*X(I)/TYPX(I)*TYPX(I) 105 CONTINUE STPSIZ = SQRT(STPSIZ) STEPMX = MAX(THOUS*STPSIZ, THOUS) 110 CONTINUE C C CHECK FUNCTION SCALE C IF(FSCALE.EQ.ZERO) FSCALE=ONE IF(FSCALE.LT.ZERO) FSCALE=-FSCALE C C CHECK GRADIENT TOLERANCE C IF(GRADTL.LT.ZERO) GRADTL=EPS**(ONE/THREE) C C CHECK STEP TOLERANCE C IF(STEPTL.LT.ZERO) THEN TEMP=EPS**(ONE/THREE) STEPTL=TEMP*TEMP ENDIF C C CHECK ITERATION LIMIT C IF(ILIM.LE.0) ILIM=150 C C CHECK NUMBER OF DIGITS OF ACCURACY IN FUNCTION FCN C IF(NDIGIT.LE.0) NDIGIT=-LOG10(EPS) IF(TEN**(-NDIGIT).LE.EPS) NDIGIT=-LOG10(EPS) RETURN C C ERROR EXITS AND WARNINGS C 901 FORMAT(' STCHKI ILLEGAL DIMENSION N=',I10) 902 FORMAT(' STCHKI ILLEGAL LENGTH OF LIRN',I10,' OR LICN=',I10) 903 FORMAT(' STCHKI ILLEGAL LENGTH OF INTEGER WORKSPACE'/, * ' STCHKI IWRK=',I10,' OR REAL WORKSPACE WRK=',I10) 904 FORMAT(' STCHKI ILLEGAL NUMBER OF NONZEROS NZ=',I10) 905 FORMAT(' STCHKI ILLEGAL ROW OR COLUMN INDEX'/, * ' STCHKI IRN( ',I6,' )=',I6,2X,' ICN( ',I6,' )=',I5) 906 FORMAT(' STCHKI DIAGONAL ELEMENT',I6,' IS NOT IN THE', * ' SPARSITY PATTERN'/, * ' STCHKI THE FINITE-DIFFERENCE HESSIAN APPROXIMATION', * ' REQUIRE'/, * ' STCHKI THAT DIAGONAL ELEMENTS BE IN THE SPARSITY', * ' PATTERN') END SUBROUTINE STCHKS(N,XPLS,FPLS,GPLS,X,ITNCNT,ICSCMX,ITRMCD, * GRADTL,STEPTL,FSCALE,ILIM,IRETCD,MXTAKE, * RGX,IPR,MSG) C C PURPOSE: C ------- C C THIS ROUTINE CHECKS THE STOPPING CRITERIA AND TERMINATES THE C OPTIMIZATION ALGORITHM IF ANY OF THE FOLLOWING IS SATISFIED: C C 1. THE SCALED GRADIENT IS LESS THAN GRADTL C 2. THE LENGTH OF THE CURRENT STEP IS LESS THAN STEPTL C 3. THE CURRENT GLOBAL STEP FAILED TO LOCATE A POINT LOWER THAN XPLS C 4. THE ITERATION LIMIT HAS BEEN EXCEEDED C 5. FIVE CONSECUTIVE STEPS OF LENGTH STEPMX HAVE BEEN TAKEN C C PARAMETERS: C ---------- C C N --> DIMENSION OF PROBLEM C XPLS --> CURRENT ITERATE C FPLS --> FUNCTION VALUE AT XPLS C GPLS --> GRADIENT AT XPLS C X --> OLD ITERATE C ITNCNT --> CURRENT ITERATION NUMBER C ICSCMX <--> NUMBER CONSECUTIVE STEPS .GE. STEPMX C (RETAIN VALUE BETWEEN SUCCESSIVE CALLS) C ITRMCD <-- TERMINATION CODE C GRADTL --> TOLERANCE AT WHICH RELATIVE GRADIENT CONSIDERED C CLOSE ENOUGH TO ZERO TO TERMINATE ALGORITHM C STEPTL --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES C CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM C FSCALE --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION C ILIM --> MAXIMUM NUMBER OF ALLOWABLE ITERATIONS C IRETCD --> RETURN CODE FROM SUBROUTINE STLSCH C MXTAKE --> BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH C RGX <-- SCALED GRADIENT C IPR --> USED DEVICE TO WHICH TO SEND OUTPUT C MSG --> CONTROL OUTPUT C MSG = 0 MEANS NO INFORMATION IS OUTPUT C MSG >= 1 MEANS INFORMATION ABOUT STOPPING CONDITIONS C IS OUTPUT C INTEGER N,ITNCNT,ICSCMX,ITRMCD,ILIM,IRETCD,IPR,MSG CD DOUBLE PRECISION XPLS(N),FPLS,GPLS(N),X(N),GRADTL,STEPTL CD DOUBLE PRECISION FSCALE,RGX CS REAL XPLS(N),FPLS,GPLS(N),X(N),GRADTL,STEPTL CS REAL FSCALE,RGX C C LOCAL VARIABLES C INTEGER I CD DOUBLE PRECISION D,RELGRD,RSX,RELSTP,ZERO,ONE CS REAL D,RELGRD,RSX,RELSTP,ZERO,ONE LOGICAL MXTAKE INTRINSIC ABS, MAX CD DATA ZERO,ONE/0.0D+0,1.0D+0/ CS DATA ZERO,ONE/0.0E+0,1.0E+0/ C ITRMCD=0 C C LAST GLOBAL STEP FAILED TO LOCATE A POINT LOWER THAN X C IF(IRETCD.EQ.1) THEN ITRMCD=3 IF(MSG .GE. 1) WRITE(IPR,903) RETURN ENDIF C C COMPUTE SCALED GRADIENT AND CHECK WHETHER IT IS WITHIN TOLERANCE C D=MAX(ABS(FPLS),FSCALE) C D=ONE RGX=ZERO DO 100 I=1,N RELGRD=ABS(GPLS(I))*MAX(ABS(XPLS(I)),ONE)/D RGX=MAX(RGX,RELGRD) 100 CONTINUE IF(RGX.LE.GRADTL) THEN ITRMCD=1 IF(MSG .GE. 1) WRITE(IPR,901) RETURN ENDIF C IF(ITNCNT.EQ.0) RETURN C C FIND DIRECTION IN WHICH RELATIVE STEPSIZE MAXIMUM C CHECK WHETHER WITHIN TOLERANCE C RSX=0.0 DO 120 I=1,N RELSTP=ABS(XPLS(I)-X(I))/MAX(ABS(XPLS(I)),ONE) RSX=MAX(RSX,RELSTP) 120 CONTINUE IF(RSX.LE.STEPTL) THEN ITRMCD=2 IF(MSG .GE. 1) WRITE(IPR,902) RETURN ENDIF C C CHECK ITERATION LIMIT C IF(ITNCNT.GE.ILIM) THEN ITRMCD=4 IF(MSG .GE. 1) WRITE(IPR,904) RETURN ENDIF C C CHECK NUMBER OF CONSECUTIVE STEPS \ STEPMX C IF(MXTAKE) THEN ICSCMX=ICSCMX+1 IF(ICSCMX.GE.5) THEN ITRMCD=5 IF(MSG .GE. 1) WRITE(IPR,905) RETURN ENDIF ELSE ICSCMX=0 RETURN ENDIF C 901 FORMAT(' STCHKS RELATIVE GRADIENT CLOSE TO ZERO'/ * ' STCHKS CURRENT ITERATE IS PROBABLY SOLUTION') 902 FORMAT(' STCHKS SUCCESSIVE ITERATES WITHIN TOLERANCE'/ * ' STCHKS CURRENT ITERATE IS PROBABLY SOLUTION') 903 FORMAT(' STCHKS LAST GLOBAL STEP FAILED TO LOCATE A POINT', * ' LOWER THAN X',/ * ' STCHKS EITHER X IS AN APPROXIMATE LOCAL MINIMUM', * ' OF THE FUNCTION',/ * ' STCHKS THE FUNCTION IS TOO NON-LINEAR FOR THIS', * ' ALGORITHM',/ * ' STCHKS OR STEPTL IS TOO LARGE') 904 FORMAT(' STCHKS ITERATION LIMIT EXCEEDED',/ * ' STCHKS ALGORITHM FAILED') 905 FORMAT(' STCHKS MAXIMUM STEP SIZE EXCEEDED 5', * ' CONSECUTIVE TIMES',/ * ' STCHKS EITHER THE FUNCTION IS UNBOUNDED BELOW',/ * ' STCHKS BECOMES ASYMPTOTIC TO A FINITE VALUE', * ' FROM ABOVE IN SOME DIRECTION',/ * ' STCHKS OR STEPMX IS TOO SMALL') END SUBROUTINE STCZ3P(A1,A2,A3,XP,NZEROS) C C PURPOSE: C -------- C C THIS ROUTINE COMPUTES THE ROOTS OF A THIRD DEGREE C POLYNOMIAL OF THE FORM : X**3 + A1*X**2 + A2*X + A3. C C PARAMETERS: C ---------- C C A1, A2, A3 ---> COEFFICIENTS OF THE CUBIC POLYNOMIAL C XP <--- ZEROS OF THE CUBIC POLYNOMIAL C NZEROS <--- NUMBER OF ZEROS FOUND C C INTEGER NZEROS CD DOUBLE PRECISION A1, A2, A3, XP(3) CS REAL A1, A2, A3, XP(3) C C LOCAL VARIABLES C CD DOUBLE PRECISION CST1, CST2, Q, R, D, S, T, THETA, PI, CONST CD DOUBLE PRECISION ZERO, HALF, ONE, TWO, THREE, FOUR, NINE CD DOUBLE PRECISION TSEVEN, FFOUR CS REAL CST1, CST2, Q, R, D, S, T, THETA, PI, CONST CS REAL ZERO, HALF, ONE, TWO, THREE, FOUR, NINE CS REAL TSEVEN, FFOUR INTRINSIC ACOS, COS, SQRT CD PARAMETER ( PI = 3.141592653589793D+0 ) CS PARAMETER ( PI = 3.141592741012573E+0 ) CD DATA ZERO, HALF, ONE, TWO, THREE, FOUR, NINE, TSEVEN, CD * FFOUR / 0.0D+0, 5.0D-1, 1.0D+0, 2.0D+0, 3.0D+0, CD * 4.0D+0, 9.0D+0, 27.0D+0, 54.0D+0 / CS DATA ZERO, HALF, ONE, TWO, THREE, FOUR, NINE, TSEVEN, CS * FFOUR / 0.0E+0, 5.0E-1, 1.0E+0, 2.0E+0, 3.0E+0, CS * 4.0E+0, 9.0E+0, 27.0E+0, 54.0E+0 / C NZEROS = 0 CST1 = ONE / THREE CST2 = CST1 * A1 Q = (THREE * A2 - A1 ** TWO) / NINE R = (NINE * A1 * A2 - TSEVEN * A3 - TWO * A1 ** THREE) / FFOUR D = Q ** THREE + R ** TWO IF(D .GT. ZERO) THEN S = R + SQRT(D) IF(S .LT. ZERO) THEN S = -((-S) ** CST1) ELSE S = S ** CST1 ENDIF T = R - SQRT(D) IF(T .LT. ZERO) THEN T = -((-T) ** CST1) ELSE T = T ** CST1 ENDIF C C COMPUTE THE ONLY REAL ZERO C NZEROS = 1 XP(1) = S + T - CST2 ELSEIF(D .EQ. ZERO) THEN IF(R .LT. ZERO) THEN S = -((-R) ** CST1) ELSE S = R ** CST1 ENDIF T = S C C COMPUTE THE TWO REAL ZEROS C NZEROS = 2 XP(1) = S + T - CST2 XP(2) = -HALF * (S + T) ELSE THETA = ACOS(R / SQRT(-Q ** THREE)) CONST = TWO * SQRT(-Q) THETA = CST1 * THETA C C COMPUTE THE THREE REAL ZEROS C NZEROS = 3 XP(1) = CONST * COS(THETA) - CST2 XP(2) = CONST * COS(THETA + TWO*PI*CST1) - CST2 XP(3) = CONST * COS(THETA + FOUR*PI*CST1) - CST2 ENDIF RETURN END SUBROUTINE STDFLT(N,TYPX,FSCALE,GRADTL,STEPTL,ILIM,STEPMX,IPR, * METHOD,GRDFLG,HSNFLG,NDIGIT,INFORM,MSG) C C PURPOSE: C-------- C C THIS ROUTINE SETS THE DEFAULT VALUES OF THE PACKAGE. C C N --> A POSITIVE INTEGER VARIABLE SPECIFYING THE NUMBER OF C VARIABLES IN THE PROBLEM C TYPX <-- AN ARRAY OF LENGTH N IN WHICH THE TYPICAL SIZE OF THE C COMPONENTS OF X ARE SPECIFIED C FSCALE <-- A POSITIVE REAL NUMBER ESTIMATING THE MAGNITUDE C OF F(X) NEAR THE MINIMIZER C GRADTL <-- POSITIVE SCALAR GIVING THE TOLERANCE AT WHICH THE C SCALED GRADIENT OF F(X) IS CONSIDERED CLOSE ENOUGH TO C ZERO TO TERMINATE THE ALGORITHM C STEPTL <-- A POSITIVE SCALAR PROVIDING THE MINIMUM ALLOWABLE C RELATIVE STEP LENGTH C ILIM <-- POSITIVE INTEGER SPECIFYING THE MAXIMUM ITERATIONS TO BE C PERFORMED BEFORE THE PROGRAM IS TERMINATED C STEPMX <-- A POSITIVE SCALAR PROVIDING THE MAXIMUM ALLOWABLE SCALED C STEP LENGTH C IPR <-- THE UNIT ON WHICH THE ROUTINE OUTPUTS INFORMATION C METHOD <-- AN INTEGER FLAG DESIGNATING WHICH METHOD TO USE C METHOD = 0 : USE NEWTON'S METHOD C METHOD = 1 : USE THE TENSOR METHOD C GRDFLG <-- INTEGER FLAG DESIGNATING WHETHER OR NOT ANALYTIC HESSIAN C HAS BEEN SUPPLIED BY THE USER C GRDFLG = 0 : NO ANALYTIC GRADIENT SUPPLIED C GRDFLG = 1 : ANALYTIC GRADIENT SUPPLIED (WILL BE CHECKED C AGAINST FINITE DIFFERENCE GRADIENT) C GRDFLG = 2 : ANALYTIC GRADIENT SUPPLIED (WILL NOT BE C CHECKED AGAINST FINITE DIFFERENCE GRADIENT) C HSNFLG <-- INTEGER FLAG DESIGNATING WHETHER OR NOT ANALYTIC HESSIAN C HAS BEEN SUPPLIED BY THE USER C HSNFLG = 0 : NO ANALYTIC HESSIAN SUPPLIED C HSNFLG = 1 : ANALYTIC HESSIAN SUPPLIED (WILL BE CHECKED C AGAINST FINITE DIFFERENCE HESSIAN) C HSNFLG = 2 : ANALYTIC HESSIAN SUPPLIED (WILL NOT BE C CHECKED AGAINST FINITE DIFFERENCE HESSIAN) C NDIGIT <-- INTEGER ESTIMATING THE NUMBER OF ACCURATE DIGITS ON THE C OBJECTIVE FUNCTION F(X) C INFORM <-- AN INTEGER VARIABLE. IF IT IS SET TO 1, THE USER MUST C OBTAIN HESS TIMES VECTOR AND RE-ENTER WITH INFORM C UNCHANGED. THE RESULT OF HESS TIMES VECTOR MUST BE C STORED IN VECTOR. THE DEFAULT VALUE OF INFORM IS 0, C MEANING THAT HESS TIMES VECTOR IS COMPUTED BY THE C PACKAGE C MSG <-- AN INTEGER VARIABLE THAT THE USER MAY SET ON INPUT TO C INHIBIT CERTAIN AUTOMATIC CHECKS OR OVERRIDE CERTAIN C DEFAULT CHARACTERISTICS OF THE PACKAGE C INTEGER N, ILIM, IPR, METHOD, GRDFLG, HSNFLG, NDIGIT, INFORM, MSG CD DOUBLE PRECISION TYPX(N), FSCALE, GRADTL, STEPTL, STEPMX CS REAL TYPX(N), FSCALE, GRADTL, STEPTL, STEPMX C C LOCAL VARIABLES C INTEGER I CD DOUBLE PRECISION EPS, TEMP, ZERO, ONE, THREE CD DOUBLE PRECISION DPMEPS CS REAL EPS, TEMP, ZERO, ONE, THREE CS REAL SPMEPS INTRINSIC LOG10 CD DATA ZERO, ONE, THREE / 0.0D+0, 1.0D+0, 3.0D+0/ CS DATA ZERO, ONE, THREE / 0.0E+0, 1.0E+0, 3.0E+0/ C CD EPS = DPMEPS() CS EPS = SPMEPS() METHOD = 1 FSCALE = ONE GRDFLG = 0 HSNFLG = 0 DO 10 I = 1,N TYPX(I) = ONE 10 CONTINUE TEMP = EPS ** (ONE/THREE) GRADTL = TEMP STEPTL = TEMP * TEMP NDIGIT = -LOG10(EPS) C C SET ACTUAL DFAULT VALUE OF STEPMX IN STCHKI C STEPMX = ZERO ILIM = 500 IPR = 6 MSG = 1 INFORM = 0 RETURN END SUBROUTINE STDRUO(N,X,NPAIRS,IRN,LIRN,ICN,LICN,FCN,UGR,USH, * TYPX,FSCALE,GRADTL,STEPTL,ILIM,STEPMX,IPR, * METHOD,GRDFLG,HSNFLG,NDIGIT,MSG,XPLS,FPLS, * GPLS,HESS,G,S,D,DN,E,XD,BV,LISTP,NGRP,IPNTR, * JPNTR,IWA,LIWRK,LWRK,TERMCD,VECTOR,INFORM) C C PURPOSE: C ------- C C THIS ROUTINE IS THE DRIVER FOR SOLVING LARGE, SPARSE UNCONSTRAINED C OPTIMIZATION PROBLEMS USING TENSOR METHODS. C C PARAMETERS: C ---------- C C N --> A POSITIVE INTEGER VARIABLE SPECIFYING THE NUMBER OF C VARIABLES IN THE PROBLEM C X --> AN ARRAY OF LENGTH N THAT CONTAINS AN INITIAL C ESTIMATE OF THE MINIMIZER C NPAIRS --> AN INTEGER VARIABLE THAT MUST BE SET BY THE USER TO C THE NUMBER OF NONZEROS STORED IN THE LOWER OR UPPER C HALF OF THE HESSIAN MATRIX C IRN --> AN INTEGER ARRAY OF LENGTH LIRN. ON ENTRY, IT MUST C HOLD THE ROW INDEX OF EACH NONZERO STORED IN THE C LOWER OR UPPER HALF OF THE HESSIAN MATRIX C LIRN --> AN INTEGER VARIABLE THAT MUST BE SET BY THE USER TO C THE LENGTH OF ARRAY IRN. LIRN NEED NOT BE AS LARGE C AS LICN; NORMALLY IT NEED NOT BE VERY MUCH GREATER C THAN NPAIRS C ICN --> AN INTEGER ARRAY OF LENGTH LICN. ON ENTRY, IT MUST C HOLD THE COLUMN INDEX OF THE NONZEROS STORED IN THE C LOWER OR UPPER HALF OF THE HESSIAN MATRIX C LICN --> AN INTEGER VARIABLE THAT MUST BE SET BY THE USER TO THE C LENGTH OF THE HESSIAN ARRAY HESS AND ICN. LICN SHOULD C ORDINARILY BE 2 TO 4 TIMES AS LARGE AS NPAIRS C FCN --> THE NAME OF A USER SUPPLIED SUBROUTINE THAT EVALUATES C THE FUNCTION F AT AN ARBITRARY VECTOR X. THE SUBROUTINE C MUST BE DECLARED EXTERNAL IN THE USER'S CALLING PROGRAM C AND MUST CONFORM TO C CALL FCN(N, X, F) C UGR --> THE NAME OF A USER SUPPLIED SUBROUTINE THAT RETURNS IN C THE VALUE OF THE GRADIENT. UGR MUST BE DECLARED EXTERNAL C IN THE USER'S CALLING PROGRAM AND MUST CONFORM TO THE C USAGE C CALL UGR(N, X, G) C IF NO ANALYTIC GRADIENT IS SUPPLIED (GRDFLG = 0), THE C USER MUST USE THE DUMMY NAME STDUGR C USH --> THE NAME OF A USER SUPPLIED SUBROUTINE THAT RETURNS C IN HESS THE VALUE OF THE HESSIAN AT THE CURRENT POINT X. C USH MUST BE DECLARED EXTERNAL IN THE USER'S CALLING C PROGRAM AND MUST CONFORM TO THE USAGE C CALL USH(N,X,NPAIRS,LICN,HESS,IRN,ICN) C ONLY THE LOWER OR UPPER TRIANGULAR PART AND THE DIAGONAL C OF HESS SHOULD BE GIVEN, WITH THEIR CORRESPONDING ROW C AND COLUMN INDICES. IF NO ANALYTIC GRADIENT IS SUPPLIED C (HSNFLG = 0), THE USER MUST USE THE DUMMY NAME STDUSH C TYPX --> AN ARRAY OF LENGTH N IN WHICH THE TYPICAL SIZE OF THE C COMPONENTS OF X ARE SPECIFIED C FSCALE --> A POSITIVE REAL NUMBER ESTIMATING THE MAGNITUDE C OF F(X) NEAR THE MINIMIZER C GRADTL --> POSITIVE SCALAR GIVING THE TOLERANCE AT WHICH THE C SCALED GRADIENT OF F(X) IS CONSIDERED CLOSE ENOUGH TO C ZERO TO TERMINATE THE ALGORITHM C STEPTL --> A POSITIVE SCALAR PROVIDING THE MINIMUM ALLOWABLE C RELATIVE STEP LENGTH C ILIM --> POSITIVE INTEGER SPECIFYING THE MAXIMUM ITERATIONS C TO BE PERFORMED BEFORE THE PROGRAM IS TERMINATED C STEPMX --> A POSITIVE SCALAR PROVIDING THE MAXIMUM ALLOWABLE C SCALED STEP LENGTH C IPR --> THE UNIT ON WHICH THE ROUTINE OUTPUTS INFORMATION C METHOD --> AN INTEGER FLAG DESIGNATING WHICH METHOD TO USE C METHOD = 0 : USE NEWTON'S METHOD C METHOD = 1 : USE THE TENSOR METHOD C GRDFLG --> INTEGER FLAG DESIGNATING WHETHER OR NOT ANALYTIC C HESSIAN HAS BEEN SUPPLIED BY THE USER C GRDFLG = 0 : NO ANALYTIC GRADIENT SUPPLIED C GRDFLG = 1 : ANALYTIC GRADIENT SUPPLIED (WILL BE CHECKED C AGAINST FINITE DIFFERENCE GRADIENT) C GRDFLG = 2 : ANALYTIC GRADIENT SUPPLIED (WILL NOT BE C CHECKED AGAINST FINITE DIFFERENCE GRADIENT) C HSNFLG --> INTEGER FLAG DESIGNATING WHETHER OR NOT ANALYTIC C HESSIAN HAS BEEN SUPPLIED BY THE USER C HSNFLG = 0 : NO ANALYTIC HESSIAN SUPPLIED C HSNFLG = 1 : ANALYTIC HESSIAN SUPPLIED (WILL BE CHECKED C AGAINST FINITE DIFFERENCE HESSIAN) C HSNFLG = 2 : ANALYTIC HESSIAN SUPPLIED (WILL NOT BE C CHECKED AGAINST FINITE DIFFERENCE HESSIAN) C NDIGIT --> INTEGER ESTIMATING THE NUMBER OF ACCURATE DIGITS IN THE C OBJECTIVE FUNCTION F(X) C MSG <--> AN INTEGER VARIABLE THAT THE USER MAY SET ON INPUT TO C INHIBIT CERTAIN AUTOMATIC CHECKS OR OVERRIDE CERTAIN C DEFAULT CHARACTERISTICS OF THE PACKAGE C XPLS <-- AN ARRAY OF LENGTH N CONTAINING THE BEST APPROXIMATION C TO THE MINIMIZER UPON RETURN. (IF THE ALGORITHM HAS NOT C CONVERGED, THE LAST ITERATE IS RETURNED) C FPLS <-- A SCALAR VARIABLE THAT CONTAINS THE FUNCTION VALUE AT C THE FINAL ITERATE XPLS C GPLS <-- AN ARRAY OF LENGTH N CONTAINING THE GRADIENT VALUE C AT XPLS C HESS <-- AN ARRAY THAT IS USED TO STORE THE HESSIAN MATRIX AT C EACH ITERATION. IT NEEDS TO BE AT LEAST OF DIMENSION C LICN. ON EXIT, HESS CONTAINS THE HESSIAN VALUE AT THE C MINIMIZER C G,S,D,DN,E,XD,BV,LISTP,NGRP,IPNTR,JPNTR,IWA --> WORKSPACE C LIWRK --> AN INTEGER VARIABLE. IT MUST BE SET BY THE USER TO THE C LENGTH OF ARRAY IWRK AND IS NOT ALTERED BY THE PACKAGE C (SEE SUBROUTINE STUMCD (STUMCS)) C LWRK --> AN INTEGER VARIABLE. IT MUST BE SET BY THE USER TO THE C LENGTH OF ARRAY WRK AND IS NOT ALTERED BY THE PACKAGE C (SEE SUBROUTINE STUMCD (STUMCS)) C TERMCD <-- AN INTEGER THAT SPECIFIES THE REASON WHY THE ALGORITHM C HAS TERMINATED C VECTOR<--> AN ARRAY OF LENGTH N. IT NEED NOT BE SET BY THE USER C ON ENTRY. IF INFORM IS SET TO 1, A RE-ENTRY MUST BE C MADE WITH VECTOR SET TO HESS TIMES VECTOR (SEE INFORM) C INFORM<--> AN INTEGER VARIABLE. IF IT IS SET TO 1, THE USER MUST C OBTAIN HESS TIMES VECTOR AND RE-ENTER WITH INFORM C UNCHANGED. THE RESULT OF HESS TIMES VECTOR MUST BE C STORED IN VECTOR. THE DEFAULT VALUE OF INFORM IS 0, C MEANING THAT HESS TIMES VECTOR IS COMPUTED BY THE C PACKAGE C C BLAS SUBROUTINES: DCOPY,DDOT,DNRM2 C INTEGER N,NPAIRS,LIRN,LICN,ILIM,IPR,METHOD,GRDFLG,HSNFLG INTEGER NDIGIT,MSG,LIWRK,LWRK,TERMCD,INFORM INTEGER IRN(LIRN),ICN(LICN),LISTP(N),NGRP(N),IPNTR(N+1) INTEGER JPNTR(N+1),IWA(*) CD DOUBLE PRECISION X(N),TYPX(N),FSCALE,GRADTL,STEPTL,STEPMX CD DOUBLE PRECISION XPLS(N),FPLS,GPLS(N),HESS(LICN),G(N) CD DOUBLE PRECISION S(N),D(N),DN(N),E(N),XD(N),BV(N),VECTOR(N) CS REAL X(N),TYPX(N),FSCALE,GRADTL,STEPTL,STEPMX CS REAL XPLS(N),FPLS,GPLS(N),HESS(LICN),G(N) CS REAL S(N),D(N),DN(N),E(N),XD(N),BV(N),VECTOR(N) LOGICAL SCALE C C LOCAL VARIABLES C INTEGER ITNNO,ICSCMX,I,MAXGRP,MINGRP,IFCNT,IGCNT,IHCNT INTEGER IFLAG,IRETCD,NZ,LIWA CD DOUBLE PRECISION RNF,ANALTL,F,GNORM,FP,SUM1,SUM2 CD DOUBLE PRECISION ALMBDA,OCENT,ONE,TWO,THREE,TEN CD DOUBLE PRECISION EPS,GAMMA,ALPHA,BETA,FN,RGX CD DOUBLE PRECISION DDOT,DNRM2,DPMEPS CS REAL RNF,ANALTL,F,GNORM,FP,SUM1,SUM2 CS REAL ALMBDA,OCENT,ONE,TWO,THREE,TEN CS REAL EPS,GAMMA,ALPHA,BETA,FN,RGX CS REAL SDOT,SNRM2,SPMEPS LOGICAL MXTAKE,DSCENT,FCALL,CHECKH COMMON/ALMB/ALMBDA COMMON/COUNT/IFCNT,IGCNT,IHCNT EXTERNAL FCN,UGR,USH CD EXTERNAL DCOPY,DNRM2,DDOT CS EXTERNAL SCOPY,SNRM2,SDOT CD DATA OCENT,ONE,TWO,THREE,TEN/1.0D-2,1.0D+0,2.0D+0, CD * 3.0D+0,10.0D+0/ CS DATA OCENT,ONE,TWO,THREE,TEN/1.0E-2,1.0E+0,2.0E+0, CS * 3.0E+0,10.0E+0/ C C IF INFORM = 2 THEN USER MUST SUPPLY VECTOR = HESS * S C IF(INFORM .EQ. 2) GO TO 100 C C COMPUTE MACHINE EPSILON C CD EPS = DPMEPS() CS EPS = SPMEPS() C C CHECK USER'S INPUT DATA C CALL STCHKI(N,NPAIRS,NZ,X,IRN,LIRN,ICN,LICN,LIWRK,LWRK,SCALE, * TYPX,FSCALE,GRADTL,STEPTL,ILIM,NDIGIT,EPS,METHOD, * GRDFLG,HSNFLG,STEPMX,IWA,LIWA,JPNTR,MSG,IPR) IF (MSG.LT.0) RETURN C C PRINT OUT USED PARAMETERS C IF(MSG .GE. 1) THEN WRITE(IPR,901) GRDFLG WRITE(IPR,902) HSNFLG WRITE(IPR,903) METHOD WRITE(IPR,904) ILIM WRITE(IPR,905) EPS WRITE(IPR,906) STEPTL WRITE(IPR,907) GRADTL WRITE(IPR,908) STEPMX ENDIF C C INITIALIZATION C ITNNO = 0 ICSCMX = 0 TERMCD = 0 IFCNT = 0 IGCNT = 0 IHCNT = 0 C C SCALE X C IF(SCALE) THEN DO 10 I = 1, N X(I) = X(I)/TYPX(I) 10 CONTINUE ENDIF C C INITIAL ITERATION C RNF = MAX(TEN ** (-NDIGIT),EPS) ANALTL = MAX(OCENT,SQRT(RNF)) C C UNSCALE X AND COMPUTE F AND G C IF(SCALE) THEN DO 30 I = 1, N VECTOR(I) = X(I)*TYPX(I) 30 CONTINUE ELSE DO 35 I = 1, N VECTOR(I) = X(I) 35 CONTINUE ENDIF IFCNT = IFCNT + 1 CALL FCN(N,VECTOR,F) IGCNT = IGCNT + 1 IF(GRDFLG .EQ. 1) THEN CALL UGR(N,VECTOR,G) CALL STCHKG(N,VECTOR,FCN,F,G,TYPX,RNF, * ANALTL,XD,MSG,IPR) IF(MSG .LT. 0) RETURN ELSEIF(GRDFLG .EQ. 2) THEN CALL UGR(N,VECTOR,G) ELSE CALL STFDGR(N,VECTOR,FCN,F,RNF,G) ENDIF C C SCALE G C IF(SCALE) THEN DO 40 I = 1, N G(I) = G(I) * TYPX(I) 40 CONTINUE ENDIF C CD GNORM = DNRM2(N,G,1) CS GNORM = SNRM2(N,G,1) C C PRINT OUT INITIAL ITERATION C IF(MSG .GE. 1) THEN CALL STRSLT(N,VECTOR,F,G,RGX,ITNNO,TERMCD,IPR) ENDIF C C TEST WHETHER INITIAL GUESS SATISFIES THE STOPPING CRITERIA C IF(GNORM .LE. GRADTL)THEN TERMCD = 1 RETURN ENDIF C C ITERATION 1 C ITNNO = ITNNO + 1 C C COMPUTE HESSIAN C IHCNT = IHCNT + 1 IF(HSNFLG .EQ. 1) THEN CALL USH(N,VECTOR,NPAIRS,LICN,HESS,IRN,ICN) FCALL = .TRUE. CHECKH = .TRUE. CALL STCHKH(N,VECTOR,NPAIRS,NZ,IRN,LIRN,ICN,LICN,FCN, * UGR,SCALE,TYPX,GRDFLG,F,G,RNF,ANALTL,FCALL, * CHECKH,LISTP,NGRP,IPNTR,JPNTR,IWA,LIWA,D, * XD,E,HESS,IPR,MSG) IF(MSG .LT. 0) RETURN FCALL = .FALSE. ELSEIF(HSNFLG .EQ. 2) THEN CALL USH(N,VECTOR,NPAIRS,LICN,HESS,IRN,ICN) CHECKH = .FALSE. CALL STCHKH(N,VECTOR,NPAIRS,NZ,IRN,LIRN,ICN,LICN,FCN, * UGR,SCALE,TYPX,GRDFLG,F,G,RNF,ANALTL,FCALL, * CHECKH,LISTP,NGRP,IPNTR,JPNTR,IWA,LIWA,D, * XD,E,HESS,IPR,MSG) IF(MSG .LT. 0) RETURN ELSE FCALL = .TRUE. CALL STFDHS(N,NZ,VECTOR,IRN,ICN,G,RNF,GRDFLG, * FCN,UGR,SCALE,TYPX,FCALL,LISTP,NGRP, * IPNTR,JPNTR,IWA,LIWA,D,XD,E,MAXGRP,MINGRP,HESS) FCALL = .FALSE. ENDIF C C SCALE HESS C IF(SCALE) THEN DO 45 I = 1, NZ HESS(I) = HESS(I) * TYPX(IRN(I)) * TYPX(ICN(I)) 45 CONTINUE ENDIF C C SOLVE FOR NEWTON STEP D C DO 50 I = 1,N D(I) = -G(I) 50 CONTINUE C C SPARSE CHOLESKY DECOMPOSITION FOR HESS (HESS = L D LT) C CALL STMA27(N,NZ,IRN,ICN,HESS,LICN,IWA,LIWA, * D,XD,1,.TRUE.,EPS,IFLAG) C C COOMPUTE NEWTON STEP USING A LINE SEARCH C CALL STLSCH(N,X,F,G,D,XPLS,FPLS,MXTAKE,IRETCD,STEPMX, * STEPTL,SCALE,TYPX,FCN,VECTOR) C C UPDATE GRADIENT AT XPLS C C UNSCALE XPLS AND COMPUTE GPLS C IF(SCALE) THEN DO 55 I = 1,N VECTOR(I) = XPLS(I) * TYPX(I) 55 CONTINUE ELSE DO 60 I = 1,N VECTOR(I) = XPLS(I) 60 CONTINUE ENDIF IGCNT = IGCNT + 1 IF(GRDFLG .EQ. 0) THEN CALL STFDGR(N,VECTOR,FCN,FPLS,RNF,GPLS) ELSE CALL UGR(N,VECTOR,GPLS) ENDIF C C SCALE GPLS C IF(SCALE) THEN DO 65 I = 1, N GPLS(I) = GPLS(I) * TYPX(I) 65 CONTINUE ENDIF C C CHECK STOPPING CONDITIONS C CALL STCHKS(N,XPLS,FPLS,GPLS,X,ITNNO,ICSCMX, * TERMCD,GRADTL,STEPTL,FSCALE,ILIM,IRETCD, * MXTAKE,RGX,IPR,MSG) C C IF TERMCD > 0 THEN STOPPING CONDITIONS SATISFIED C IF(TERMCD .GT. 0) GO TO 215 C C UPDATE X, F AND S FOR TENSOR MODEL C FP = F F = FPLS DO 70 I = 1,N S(I) = X(I) - XPLS(I) X(I) = XPLS(I) 70 CONTINUE C C IF MSG >= 2 THEN PRINT OUT EACH ITERATION C IF(MSG .GE. 2)THEN CALL STRSLT(N,XPLS,FPLS,GPLS,RGX,ITNNO,TERMCD,IPR) ENDIF C C ITERATION > 1 C C UNSCALE X AND COMPUTE HESS C 75 CONTINUE C IF(SCALE) THEN DO 80 I = 1,N VECTOR(I) = X(I) * TYPX(I) 80 CONTINUE ELSE DO 85 I = 1,N VECTOR(I) = X(I) 85 CONTINUE ENDIF IHCNT = IHCNT + 1 IF(HSNFLG .EQ. 0) THEN CALL STFDHS(N,NZ,VECTOR,IRN,ICN,GPLS,RNF,GRDFLG, * FCN,UGR,SCALE,TYPX,FCALL,LISTP,NGRP, * IPNTR,JPNTR,IWA,LIWA,D,XD,E,MAXGRP,MINGRP,HESS) ELSE CALL USH(N,VECTOR,NZ,LICN,HESS,IRN,ICN) ENDIF C C SCALE HESS C IF(SCALE) THEN DO 90 I = 1, NZ HESS(I) = HESS(I) * TYPX(IRN(I)) * TYPX(ICN(I)) 90 CONTINUE ENDIF C C IF METHOD = 0 THEN USE NEWTON STEP ONLY C IF(METHOD .EQ. 0) THEN C C COMPUTE NEWTON STEP C DO 95 I = 1,N D(I) = - GPLS(I) 95 CONTINUE C C CHOLESKY DECOMPOSITION FOR HESS C CALL STMA27(N,NZ,IRN,ICN,HESS,LICN,IWA,LIWA, * D,XD,1,.TRUE.,EPS,IFLAG) GO TO 130 C ENDIF C IF(INFORM .EQ. 1) THEN INFORM = 2 CD CALL DCOPY(N,S(1),1,VECTOR(1),1) CS CALL SCOPY(N,S(1),1,VECTOR(1),1) RETURN ELSE CALL STHMUV(N,NZ,IRN,ICN,HESS,S,VECTOR) GO TO 110 ENDIF C 100 INFORM = 1 110 CONTINUE C C FORM TENSOR MODEL C CALL STFTSM(N,F,FP,GPLS,G,S,ALPHA,BETA,VECTOR,D) C CD SUM1 = DDOT(N,S(1),1,S(1),1) CS SUM1 = SDOT(N,S(1),1,S(1),1) CD SUM2 = DDOT(N,S(1),1,D(1),1) CS SUM2 = SDOT(N,S(1),1,D(1),1) C DO 120 I = 1,N BV(I) = (THREE*SUM1*D(I)-TWO*S(I)*SUM2)/(THREE*SUM1**3) 120 CONTINUE GAMMA = BETA/SUM1**4 C C SOLVE TENSOR MODEL AND COMPUTE TENSOR STEP C CALL STMSLV(N,D,S,GPLS,E,VECTOR,DN,BV,GAMMA,NZ,IRN, * ICN,HESS,LICN,IWA,LIWA,XD,EPS,DSCENT,IFLAG) C IF(.NOT. DSCENT) THEN CD CALL DCOPY(N,DN(1),1,D(1),1) CS CALL SCOPY(N,DN(1),1,D(1),1) ENDIF C 130 ITNNO = ITNNO + 1 CD CALL DCOPY(N,GPLS(1),1,G(1),1) CS CALL SCOPY(N,GPLS(1),1,G(1),1) C C COMPUTE TENSOR (OR NEWTON) STEP USING A LINE SEARCH C CALL STLSCH(N,X,F,G,D,XPLS,FPLS,MXTAKE,IRETCD,STEPMX, * STEPTL,SCALE,TYPX,FCN,VECTOR) C IF(METHOD .EQ. 0) GO TO 190 IF(ALMBDA .EQ. ONE) GO TO 170 C IF(DSCENT) THEN C C FULL TENSOR STEP IS A DESCENT DIRECTION BUT DOES NOT C PROVIDE ENOUGH DECREASE IN THE OBJECTIVE FUNCTION. C CALCULATE NEWTON STEP ONLY IF HESSIAN IS INDEFINITE; C OTHERWISE NEWTON STEP WAS ALREADY COMPUTED AS PART C OF THE TENSOR STEP COMPUTATION C IF(IFLAG .EQ. 3) THEN DO 140 I = 1, N DN(I) = -E(I) 140 CONTINUE GO TO 160 ENDIF CD CALL DCOPY(N,G(1),1,DN(1),1) CS CALL SCOPY(N,G(1),1,DN(1),1) CALL STMA27(N,NZ,IRN,ICN,HESS,LICN,IWA,LIWA, * DN,XD,3,.TRUE.,EPS,IFLAG) IF(IFLAG .EQ. 30) THEN DO 150 I = 1, N DN(I) = -E(I) 150 CONTINUE ELSE DO 155 I = 1, N DN(I) = -DN(I) 155 CONTINUE ENDIF 160 CONTINUE C C CALCULATE A STEP IN THE NEWTON DIRECTION C CALL STLSCH(N,X,F,G,DN,XD,FN,MXTAKE,IRETCD,STEPMX, * STEPTL,SCALE,TYPX,FCN,VECTOR) C C COMPARE TENSOR STEP TO NEWTON STEP C IF NEWTON STEP IS BETTER, SET NEXT ITERATE TO NEW NEWTON POINT C IF(FN .LT. FPLS)THEN FPLS = FN CD CALL DCOPY(N,DN(1),1,D(1),1) CS CALL SCOPY(N,DN(1),1,D(1),1) CD CALL DCOPY(N,XD(1),1,XPLS(1),1) CS CALL SCOPY(N,XD(1),1,XPLS(1),1) ENDIF ENDIF 170 CONTINUE DO 180 I = 1,N D(I) = XPLS(I) - X(I) 180 CONTINUE C 190 CONTINUE C C UNSCALE XPLS, AND COMPUTE FPLS AND GPLS C IF(SCALE) THEN DO 200 I = 1,N VECTOR(I) = XPLS(I) * TYPX(I) 200 CONTINUE ELSE DO 205 I = 1,N VECTOR(I) = XPLS(I) 205 CONTINUE ENDIF IGCNT = IGCNT + 1 IF(GRDFLG .EQ. 0) THEN CALL STFDGR(N,VECTOR,FCN,FPLS,RNF,GPLS) ELSE CALL UGR(N,VECTOR,GPLS) ENDIF C C SCALE GPLS C IF(SCALE) THEN DO 210 I = 1, N GPLS(I) = GPLS(I) * TYPX(I) 210 CONTINUE ENDIF C C CHECK STOPPING CONDITIONS C CALL STCHKS(N,XPLS,FPLS,GPLS,X,ITNNO,ICSCMX, * TERMCD,GRADTL,STEPTL,FSCALE,ILIM,IRETCD, * MXTAKE,RGX,IPR,MSG) C C IF TERMCD = 0 THEN EXECUTION CONTINUES C IF(TERMCD .EQ. 0) GO TO 230 C 215 CONTINUE C C TRANSFORM XPLS BACK TO ORIGINAL SPACE C IF(SCALE) THEN DO 220 I = 1,N XPLS(I) = XPLS(I) * TYPX(I) 220 CONTINUE ENDIF C C IF MSG >= 1 THEN PRINT OUT FINAL ITERATION C IF(MSG .GE. 1) THEN CALL STRSLT(N,XPLS,FPLS,GPLS,RGX,ITNNO,TERMCD,IPR) ENDIF C C UPDATE THE HESSIAN C IHCNT = IHCNT + 1 IF(HSNFLG .EQ. 0) THEN CALL STFDHS(N,NZ,VECTOR,IRN,ICN,GPLS,RNF,GRDFLG,FCN, * UGR,SCALE,TYPX,FCALL,LISTP,NGRP,IPNTR, * JPNTR,IWA,LIWA,D,XD,E,MAXGRP,MINGRP,HESS) ELSE CALL USH(N,XPLS,NZ,LICN,HESS,IRN,ICN) ENDIF C C SCALE HESS C IF(SCALE) THEN DO 225 I = 1, NZ HESS(I) = HESS(I) * TYPX(IRN(I)) * TYPX(ICN(I)) 225 CONTINUE ENDIF RETURN C C UPDATE INFORMATION AT THE CURRENT POINT C 230 CONTINUE CD CALL DCOPY(N,XPLS(1),1,X(1),1) CS CALL SCOPY(N,XPLS(1),1,X(1),1) DO 240 I = 1,N S(I) = -D(I) 240 CONTINUE C C IF ITERATION LIMIT EXCEEDED THEN RETURN C IF(ITNNO .GT. ILIM) GO TO 215 C C IF MSG >= 2 THEN PRINT OUT EACH ITERATION C IF(MSG .GE. 2)THEN CALL STRSLT(N,XPLS,FPLS,GPLS,RGX,ITNNO,TERMCD,IPR) ENDIF C C UPDATE F C FP = F F = FPLS C C PERFORM NEXT ITERATION C GO TO 75 C C END OF ITERATION > 1 C 901 FORMAT(' STDRUO GRADIENT FLAG = ',I1) 902 FORMAT(' STDRUO HESSIAN FLAG = ',I1) 903 FORMAT(' STDRUO METHOD = ',I1) 904 FORMAT(' STDRUO ITERATION LIMIT =',I5) 905 FORMAT(' STDRUO MACHINE EPSILON = ',E20.13) 906 FORMAT(' STDRUO STEP TOLERANCE = ',E20.13) 907 FORMAT(' STDRUO GRADIENT TOLERANCE = ',E20.13) 908 FORMAT(' STDRUO MAXIMUM STEP SIZE = ',E20.13) C END SUBROUTINE STDUGR(N,X,G) C C PURPOSE: C ------- C C THIS IS A DUMMY ROUTINE TO PREVENT UNSATISFIED EXTERNAL DIAGNOSTIC C WHEN SPECIFIC ANALYTIC GRADIENT IS NOT SUPPLIED. C INTEGER N CD DOUBLE PRECISION X(N),G(N) CS REAL X(N),G(N) RETURN END SUBROUTINE STDUSH(N,X,NZ,LICN,HESS,IRN,ICN) C C PURPOSE: C ------- C C THIS IS A DUMMY ROUTINE TO PREVENT UNSATISFIED EXTERNAL DIAGNOSTIC C WHEN SPECIFIC ANALYTIC HESSIAN IS NOT SUPPLIED. C INTEGER N,NZ,LICN INTEGER IRN(NZ),ICN(LICN) CD DOUBLE PRECISION X(N),HESS(LICN) CS REAL X(N),HESS(LICN) RETURN END SUBROUTINE STFDGR(N,XPLS,FCN,FPLS,RNOISE,FDG) C C PURPOSE C ------- C C THIS ROUTINE FINDS A FORWARD FINITE-DIFFERENCE APPROXIMATION "FDG" C TO THE FIRST DERIVATIVE OF THE FUNCTION DEFINED BY THE SUBPROGRAM C "FCN" EVALUATED AT THE NEW ITERATE "XPLS". C C PARAMETERS: C ---------- C N --> DIMENSION OF PROBLEM C XPLS --> NEW ITERATE C FCN --> NAME OF SUBROUTINE THAT EVALUATES THE C OPTIMIZATION FUNCTION C FPLS --> FUNCTION VALUE AT NEW ITERATE C RNOISE --> RELATIVE NOISE IN FCN C FDG <-- FINITE-DIFFERENCE APPROXIMATION OF GRADIENT AT XPLS C INTEGER N CD DOUBLE PRECISION XPLS(N),FPLS,RNOISE,FDG(N) CS REAL XPLS(N),FPLS,RNOISE,FDG(N) C C LOCAL VARIABLES C INTEGER J CD DOUBLE PRECISION XTMPJ,STEPSZ,FHAT,ONE CS REAL XTMPJ,STEPSZ,FHAT,ONE EXTERNAL FCN INTRINSIC ABS, MAX, SQRT CD DATA ONE/ 1.0D+0 / CS DATA ONE/ 1.0E+0 / C C FIND DERIVATIVE OF FCN WITH RESPECT TO XPLS(J), J=1,N C DO 30 J=1,N XTMPJ=XPLS(J) STEPSZ=SQRT(RNOISE)*MAX(ABS(XPLS(J)),ONE) XPLS(J)=XTMPJ+STEPSZ CALL FCN(N,XPLS,FHAT) XPLS(J)=XTMPJ FDG(J)=(FHAT-FPLS)/STEPSZ 30 CONTINUE RETURN END SUBROUTINE STFDHS(N,NZ,X,IRN,ICN,G,RNF,GRDFLG, * FCN,UGR,SCALE,TYPX,FCALL,LISTP, * NGRP,IPNTR,JPNTR,IWA,LIWA,FHESD,XD, * ETA,MAXGRP,MINGRP,HESS) C C PURPOSE: C ------- C C THIS ROUTINE FINDS A FORWARD FINITE-DIFFERENCE APPROXIMATION "HESS" C TO THE SECOND DERIVATIVE (HESSIAN) OF THE FUNCTION DEFINED BY THE C SUBPROGRAM "FCN" EVALUATED AT THE NEW ITERATE "XPLS". C C PARAMETERS: C ---------- C C N --> DIMENSION OF PROBLEM C NZ --> NUMBER OF NONZEROS IN THE LOWER TRIANGULAR PART C OF THE HESSIAN MATRIX C X --> CURRENT ITERATE C IRN --> ROW INDEX OF THE NONZERO ELEMENTS IN THE LOWER C TRIANGULAR PART OF THE HESSIAN C ICN --> COLUMN INDEX OF THE NONZERO ELEMENTS IN THE LOWER C TRIANGULAR PART OF THE HESSIAN C G --> GRADIENT VALUE AT X C RNF --> RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN C GRDFLG --> =1 IF ANALYTIC GRADIENT SUPPLIED C FCN --> NAME OF A USER SUPPLIED SUBROUTINE THAT EVALUATES C THE FUNCTION F AT AN ARBITRARY VECTOR X C UGR --> NAME OF A USER SUPPLIED SUBROUTINE THAT RETURNS IN C THE VALUE OF THE GRADIENT AT AN ARBITRARY VECTOR X C SCALE --> LOGICAL VARIABLE INDICATING WHETHER OR NOT SCALING C OF THE VARIABLES IS PERFORMED C TYPX --> AN ARRAY OF LENGTH N IN WHICH THE TYPICAL SIZE OF THE C COMPONENTS OF X ARE SPECIFIED C FCALL --> LOGICAL FLAG. IF THE HESSIAN MATRIX IS EVALUATED FOR C THE FIRST TIME, THEN FCALL IS SET TO TRUE; OTHERWISE C IT IS SET TO FALSE C LISTP,NGRP,IPNTR,JPNTR,IWA --> WORKSPACE C (SEE SUBROUTINE DSSM FOR DETAILS) C LIWA --> LENGTH OF IWA ARRAY C FHESD,XD,ETA --> WORKSPACE C MAXGRP,MINGRP<-- SEE SUBROUTINE DSSM FOR DETAILS C HESS <-- FINITE-DIFFERENCE HESSIAN ON EXIT C C INTEGER N,NZ,GRDFLG,LIWA,MAXGRP,MINGRP INTEGER IRN(NZ),ICN(NZ),LISTP(N),NGRP(N),IPNTR(N+1) INTEGER JPNTR(N+1),IWA(LIWA) CD DOUBLE PRECISION X(N),G(N),TYPX(N),RNF,FHESD(N) CD DOUBLE PRECISION XD(N),ETA(N),HESS(NZ) CS REAL X(N),G(N),TYPX(N),RNF,FHESD(N) CS REAL XD(N),ETA(N),HESS(NZ) LOGICAL SCALE C C LOCAL VARIABLES C INTEGER JP,I,INFO,J,METHOD,NUMGRP CD DOUBLE PRECISION FD,ONE CS REAL FD,ONE LOGICAL FCALL EXTERNAL FCN,UGR INTRINSIC ABS, MAX, SQRT CD DATA ONE / 1.0D+0 / CS DATA ONE / 1.0E+0 / C METHOD = 1 C C CALL DSM ONLY IF IT IS THE FIRST TIME C IF(FCALL) THEN CALL DSSM(N,NZ,IRN,ICN,METHOD,LISTP,NGRP,MAXGRP, * MINGRP,INFO,IPNTR,JPNTR,IWA,LIWA) ELSE CALL SETR(N,N,IRN,JPNTR,ICN,IPNTR,IWA) ENDIF C C UNSCALE G C IF(SCALE) THEN DO 20 I = 1, N G(I) = G(I)/TYPX(I) 20 CONTINUE ENDIF C C APPROXIMATE THE HESSIAN MATRIX C DO 60 NUMGRP = 1, MAXGRP DO 40 J = 1, N ETA(J) = SQRT(RNF)*MAX(ABS(X(J)),ONE) XD(J) = X(J) IF (NGRP(J) .EQ. NUMGRP) XD(J) = X(J) + ETA(J) 40 CONTINUE IF(GRDFLG .EQ. 0) THEN CALL FCN(N,XD,FD) CALL STFDGR(N,XD,FCN,FD,RNF,FHESD) ELSE CALL UGR(N,XD,FHESD) ENDIF DO 50 I = 1, N FHESD(I) = FHESD(I) - G(I) 50 CONTINUE CALL FDHS(N,IRN,JPNTR,ICN,IPNTR, * LISTP,NGRP,MAXGRP,NUMGRP,ETA,FHESD,HESS,IWA) 60 CONTINUE DO 80 J = 1, N DO 70 JP = JPNTR(J), JPNTR(J+1)-1 ICN(JP) = J 70 CONTINUE 80 CONTINUE C C SCALE G C IF(SCALE) THEN DO 90 I = 1, N G(I) = G(I)*TYPX(I) 90 CONTINUE ENDIF C END SUBROUTINE STFTSM(N,F,FP,G,GP,S,ALPHA,BETA,SH,A) C C PURPOSE: C ------- C C THIS ROUTINE FORMS THE TENSOR MODEL. C C PARAMETERS: C ---------- C C N --> DIMENSION OF PROBLEM C F --> CURRENT FUNCTION VALUE C FP --> PREVIOUS FUNCTION VALUE C G --> CURRENT GRADIENT C GP --> PREVIOUS GRADIENT C S --> STEP TO PREVIOUS POINT C ALPHA <-- SCALAR TO FORM 3RD ORDER TERM OF TENSOR MODEL C BETA <-- SCALAR TO FORM 4TH ORDER TERM OF TENSOR MODEL C SH <-- HESSIAN MATRIX TIMES S (SEE SUBROUTINE STHMUV) C A <-- A=2*(GP-G-SH-S*BETA/(6*(S-TRANS*S))) C C BLAS SUBROUTINES: DDOT C INTEGER N CD DOUBLE PRECISION F,FP,G(N),GP(N),S(N),ALPHA,BETA,SH(N),A(N) CS REAL F,FP,G(N),GP(N),S(N),ALPHA,BETA,SH(N),A(N) C C LOCAL VARIABLES C INTEGER I CD DOUBLE PRECISION GS, GPS, SHS, B1, B2, STS CD DOUBLE PRECISION HALF, TWO, SIX, TFOUR, STWO CD DOUBLE PRECISION DDOT CS REAL GS, GPS, SHS, B1, B2, STS CS REAL HALF, TWO, SIX, TFOUR, STWO CS REAL SDOT CD EXTERNAL DDOT CS EXTERNAL SDOT CD DATA HALF, TWO, SIX, TFOUR, STWO / 5.0D-1, 2.0D+0, CD * 6.0D+0, 24.0D+0, 72.0D+0/ CS DATA HALF, TWO, SIX, TFOUR, STWO / 5.0E-1, 2.0E+0, CS * 6.0E+0, 24.0E+0, 72.0E+0/ C CD GS = DDOT(N,G(1),1,S(1),1) CS GS = SDOT(N,G(1),1,S(1),1) CD GPS = DDOT(N,GP(1),1,S(1),1) CS GPS = SDOT(N,GP(1),1,S(1),1) CD SHS = DDOT(N,SH(1),1,S(1),1) CS SHS = SDOT(N,SH(1),1,S(1),1) B1 = GPS - GS - SHS B2 = FP - F - GS - HALF * SHS ALPHA = TFOUR * B2 - SIX * B1 BETA = TFOUR * B1 - STWO * B2 C C COMPUTE A C CD STS = DDOT(N,S(1),1,S(1),1) CS STS = SDOT(N,S(1),1,S(1),1) DO 50 I = 1, N A(I) = TWO * (GP(I) - G(I) - SH(I) - S(I)*BETA /(SIX*STS)) 50 CONTINUE RETURN END SUBROUTINE STHMUV(N,NZ,IRN,ICN,HESS,S,SH) C C PURPOSE: C ------- C C THIS ROUTINE COMPUTES SH = HESS * S. C C PARAMETERS: C ---------- C N --> A POSITIVE INTEGER VARIABLE SPECIFYING THE NUMBER OF C VARIABLES IN THE PROBLEM C NZ --> NUMBER OF NONZEROS STORED IN THE LOWER HALF OF C THE HESSIAN MATRIX C IRN --> AN INTEGER ARRAY OF LENGTH LIRN. ON ENTRY, IT MUST C HOLD THE ROW INDEX OF EACH NONZERO STORED IN THE LOWER C HALF OF THE HESSIAN MATRIX C ICN --> AN INTEGER ARRAY OF LENGTH LICN. ON ENTRY, IT MUST C HOLD THE COLUMN INDEX OF THE NONZEROS STORED IN LOWER C HALF OF THE HESSIAN MATRIX C HESS --> HESSIAN MATRIX AT CURRENT ITERATE C S --> STEP TO PREVIOUS POINT C SH <-- HESS * S ON EXIT C C INTEGER N,NZ INTEGER IRN(NZ),ICN(NZ) CD DOUBLE PRECISION HESS(NZ),S(N),SH(N) CS REAL HESS(NZ),S(N),SH(N) C C LOCAL VARIABLES C INTEGER I,II,JJ CD DOUBLE PRECISION ZERO CS REAL ZERO CD DATA ZERO/0.0D0/ CS DATA ZERO/0.0E0/ C DO 110 I = 1,N SH(I) = ZERO 110 CONTINUE DO 120 I = 1,NZ II = IRN(I) JJ = ICN(I) SH(II) = SH(II) + HESS(I) * S(JJ) IF(II .NE. JJ) THEN SH(JJ) = SH(JJ) + HESS(I) * S(II) ENDIF 120 CONTINUE RETURN END SUBROUTINE STLSCH(N,X,F,G,P,XPLS,FPLS,MXTAKE,IRETCD,STEPMX, * STEPTL,SCALE,TYPX,FCN,W2) C C PURPOSE: C ------- C C THIS ROUTINE FIND A NEXT ITERATE BY A LINE SEARCH. C THIS IS THE ALPHA CONDITION ONLY LINE SEARCH. C C PARAMETERS: C ---------- C N --> DIMENSION OF PROBLEM C X --> OLD ITERATE C F --> FUNCTION VALUE AT OLD ITERATE C G --> GRADIENT AT OLD ITERATE C P --> NEWTON OR TENSOR DIRECTION C XPLS <-- NEW ITERATE C FPLS <-- FUNCTION VALUE AT NEW ITERATE C MXTAKE <-- BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED C IRETCD <-- RETURN CODE C STEPMX --> MAXIMUM ALLOWABLE STEP SIZE C STEPTL --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES C CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM C SCALE --> LOGICAL VARIABLE INDICATING WHETHER OR NOT SCALING C OF THE VARIABLES IS PERFORMED C TYPX --> DIAGONAL SCALING MATRIX FOR X C FCN --> NAME OF SUBROUTINE TO EVALUATE OPTIMIZATION FUNCTION C W2 --> WORKING SPACE C C BLAS SUBROUTINES: DDOT,DSCAL C INTEGER N,IRETCD CD DOUBLE PRECISION X(N),F,G(N),P(N),XPLS(N),FPLS,STEPMX CD DOUBLE PRECISION STEPTL,TYPX(N),W2(N) CS REAL X(N),F,G(N),P(N),XPLS(N),FPLS,STEPMX CS REAL STEPTL,TYPX(N),W2(N) LOGICAL SCALE C C LOCAL VARIABLES C INTEGER I,K INTEGER IFCNT,IGCNT,IHCNT CD DOUBLE PRECISION TMP, SLN, SCL, SLP, RLN CD DOUBLE PRECISION TEMP, TEMP1, TEMP2, ALMBMN, ALMBDA CD DOUBLE PRECISION TLMBDA, T1, T2, T3, PLMBDA CD DOUBLE PRECISION A, B, DISC, PFPLS, ZERO, TENTH CD DOUBLE PRECISION HALF, Z99, ONE, TWO, THREE, ALPHA CD DOUBLE PRECISION DDOT CS REAL TMP, SLN, SCL, SLP, RLN CS REAL TEMP, TEMP1, TEMP2, ALMBMN, ALMBDA CS REAL TLMBDA, T1, T2, T3, PLMBDA CS REAL A, B, DISC, PFPLS, ZERO, TENTH CS REAL HALF, Z99, ONE, TWO, THREE, ALPHA CS REAL SDOT LOGICAL MXTAKE COMMON / COUNT / IFCNT, IGCNT, IHCNT COMMON/ALMB/ALMBDA CD EXTERNAL FCN,DSCAL,DDOT CS EXTERNAL FCN,SSCAL,SDOT INTRINSIC ABS,MAX,SQRT CD DATA ZERO, TENTH, HALF, Z99, ONE, TWO, THREE, ALPHA CD */ 0.0D+0, 1.0D-1, 5.0D-1, 0.99D+0, 1.0D+0, 2.0D+0, 3.0D+0, CD * 1.0D-4/ CS DATA ZERO, TENTH, HALF, Z99, ONE, TWO, THREE, ALPHA CS */ 0.0E+0, 1.0E-1, 5.0E-1, 0.99E+0, 1.0E+0, 2.0E+0, 3.0E+0, CS * 1.0E-4/ C MXTAKE=.FALSE. IRETCD=2 C$ WRITE(IPR,954) C$ WRITE(IPR,955) (P(I),I=1,N) TMP=ZERO DO 5 I=1,N TMP=TMP+P(I)*P(I) 5 CONTINUE SLN=SQRT(TMP) IF(SLN.LE.STEPMX) GO TO 10 C C STEP LONGER THAN MAXIMUM ALLOWED C SCL=STEPMX/SLN CD CALL DSCAL(N,SCL,P(1),1) CS CALL SSCAL(N,SCL,P(1),1) SLN=STEPMX C$ WRITE(IPR,954) C$ WRITE(IPR,955) (P(I),I=1,N) 10 CONTINUE CD SLP = DDOT(N,G(1),1,P(1),1) CS SLP = SDOT(N,G(1),1,P(1),1) RLN = ZERO DO 15 I = 1,N TEMP = ONE TEMP1 = ABS(X(I)) TEMP2 = MAX(TEMP1,TEMP) TEMP1 = ABS(P(I)) RLN = MAX(RLN,TEMP1/TEMP2) 15 CONTINUE ALMBMN = STEPTL/RLN ALMBDA = ONE C$ WRITE(IPR,952) SLN,SLP,RMNLMB,STEPMX,STEPTL C C LOOP C CHECK IF NEW ITERATE SATISFACTORY. GENERATE NEW LAMBDA IF NECESSARY. C 100 CONTINUE IF(IRETCD.LT.2) THEN RETURN ENDIF DO 105 I=1,N XPLS(I)=X(I) + ALMBDA*P(I) 105 CONTINUE IFCNT = IFCNT + 1 IF(SCALE) THEN DO 101 K = 1,N W2(K) = XPLS(K)*TYPX(K) 101 CONTINUE CALL FCN(N,W2,FPLS) ELSE CALL FCN(N,XPLS,FPLS) ENDIF C$ WRITE(IPR,956) ALMBDA C$ WRITE(IPR,951) C$ WRITE(IPR,955) (XPLS(I),I=1,N) C$ WRITE(IPR,953) FPLS IF(FPLS.GT. F+SLP*ALPHA*ALMBDA) GO TO 130 C IF(FPLS.LE. F+SLP*1.E-4*ALMBDA) C THEN C C SOLUTION FOUND C IRETCD=0 IF(ALMBDA.EQ. TENTH .AND. SLN.GT. Z99*STEPMX) MXTAKE=.TRUE. GO TO 100 C C SOLUTION NOT (YET) FOUND C C ELSE 130 IF(ALMBDA .GE. ALMBMN) GO TO 140 C IF(ALMBDA .LT. ALMBMN) C THEN C C NO SATISFACTORY XPLS FOUND SUFFICIENTLY DISTINCT FROM X C IRETCD=1 GO TO 100 C ELSE C C CALCULATE NEW LAMBDA C 140 IF(ALMBDA.NE.ONE) GO TO 150 C IF(ALMBDA.EQ.1.0) C THEN C C FIRST BACKTRACK: QUADRATIC FIT C TLMBDA=-SLP/(TWO*(FPLS-F-SLP)) GO TO 170 C ELSE C C ALL SUBSEQUENT BACKTRACKS: CUBIC FIT C 150 T1=FPLS-F-ALMBDA*SLP T2=PFPLS-F-PLMBDA*SLP T3=ONE/(ALMBDA-PLMBDA) A=T3*(T1/(ALMBDA*ALMBDA) - T2/(PLMBDA*PLMBDA)) B=T3*(T2*ALMBDA/(PLMBDA*PLMBDA) * - T1*PLMBDA/(ALMBDA*ALMBDA) ) DISC=B*B-THREE*A*SLP IF(DISC.LE. B*B) GO TO 160 C IF(DISC.GT. B*B) C THEN C C ONLY ONE POSITIVE CRITICAL POINT, MUST BE MINIMUM C TLMBDA=(-B+SIGN(ONE,A)*SQRT(DISC))/(THREE*A) GO TO 165 C ELSE C C BOTH CRITICAL POINTS POSITIVE, FIRST IS MINIMUM C 160 TLMBDA=(-B-SIGN(ONE,A)*SQRT(DISC))/(THREE*A) C ENDIF 165 IF(TLMBDA.GT. HALF*ALMBDA) TLMBDA=HALF*ALMBDA C ENDIF 170 PLMBDA=ALMBDA PFPLS=FPLS IF(TLMBDA.GE. ALMBDA*TENTH) GO TO 180 C IF(TLMBDA.LT.ALMBDA/TENTH) C THEN ALMBDA=ALMBDA*TENTH GO TO 190 C ELSE 180 ALMBDA=TLMBDA C ENDIF C ENDIF C ENDIF 190 GO TO 100 C$951 FORMAT(' STLSCH NEW ITERATE (XPLS)') C$952 FORMAT(' STLSCH SLN =',E20.13/ C$ * ' STLSCH SLP =',E20.13/ C$ * ' STLSCH ALMBMN=',E20.13/ C$ * ' STLSCH STEPMX=',E20.13/ C$ * ' STLSCH STEPTL=',E20.13) C$953 FORMAT(' STLSCH F(XPLS)=',E20.13) C$954 FORMAT(' STLSCH COMPUTED STEP (P)') C$955 FORMAT(' STLSCH ',5(E20.13,3X)) C$956 FORMAT(' STLSCH ALMBDA=',E20.13) END SUBROUTINE STMA27(N,NZ,IRN,ICN,A,MAXA,IW,MAXIW, * RHS,W,ISLV,NEWTON,EPS,IFLAG) C C PURPOSE: C -------- C C THIS ROUTINE SOLVES A SPARSE SYMMETRIC SYSTEM OF LINEAR EQUATIONS: C A * X = RHS. C C PARAMETERS: C ---------- C C N --> DIMENSION OF MATRIX A C NZ --> NUMBER OF NONZEROS IN THE LOWER HALF OF MATRIX A C IRN --> ROW INDEX OF A. IRN IS UNALTERED BY MA27 C ICN --> COLUMN INDEX OF A. ICN IS UNALTERED BY MA27 C A <-> A SYMMETRIC MATRIX OF DIMENSION N. ON EXIT, IT HOLDS C THE FACTORS L AND D OF THE L D L-TRANS FACTORIZATION C MAXA --> MAXIMUM LENGTH OF MATRIX A C IW <-- WORKSPACE OF DIMENSION MAXIW. IT IS USED TO C HOLD INTEGER INFORMATION ON THE FACTORS L AND D C RHS <-> RIGHT HAND SIDE VECTOR OF THE EQUATIONS BEING SOLVED. C ON EXIT, IT HOLDS THE SOLUTION VECTOR C W --> WORKSPACE OF DIMENSION N C ISLV --> FLAG WITH THE FOLLOWING MEANINGS: C ISLV = 0 : ANALAYSE AND FACTORIZE MATRIX A C ISLV = 1 : ANALAYSE, FACTORIZE, AND SOLVE SYSTEMS OF C EQUATIONS C ISLV = 2 : PERFORM A SOLVE ONLY (A HAS ALREADY BEEN C FACTORIZED PREVIOUSLY) C ISLV = 3 : CHANGE THE NEGATIVE SIGNS OF THE EIGENVALUES C TO PLUS IF ANY C NEWTON --> LOGICAL PARAMETER. IT IS SET TO TRUE IF NEWTON'S METHOD C IS USED AND TO FALSE IF THE TENSOR METHOD IS USED C EPS --> MACHINE PRECISION C IFLAG <-- A VALUE OF ZERO INDICATES THAT THE SUBROUTINE HAS C PERFORMED SUCCESSFULLY C INTEGER N,NZ,MAXA,MAXIW,ISLV,IFLAG INTEGER IRN(NZ),ICN(MAXA),IW(MAXIW) CD DOUBLE PRECISION A(MAXA),RHS(N),W(N),EPS CS REAL A(MAXA),RHS(N),W(N),EPS LOGICAL NEWTON C C LOCAL VARIABLES C INTEGER NEG1,NEG2,LDIAG,LP,MP,ISTIW,LIW,NSTEPS,MAXFRT CD DOUBLE PRECISION U,TENTH CS REAL U,TENTH COMMON /MA27DD/ U,LP,MP,LDIAG CD DATA TENTH/0.1D+0/ CS DATA TENTH/0.1E+0/ C IF(ISLV .EQ. 2) GO TO 90 IF(ISLV .EQ. 3) GO TO 80 U = TENTH LDIAG = 0 C C SET UP SUBDIVISION OF ARRAYS C ISTIW = 5*N+1 LIW = MAXIW - 5*N IFLAG = 0 C C ANALYSE SPARSITY PATTERN C CD CALL MA27AD(N,NZ,IRN,ICN,IW(ISTIW),LIW,IW, CS CALL MA27A(N,NZ,IRN,ICN,IW(ISTIW),LIW,IW, * IW(3*N+1),NSTEPS,IFLAG) C C FACTORIZE MATRIX C CD CALL MA27BD(N,NZ,IRN,ICN,A,MAXA,IW(ISTIW),LIW, CS CALL MA27B(N,NZ,IRN,ICN,A,MAXA,IW(ISTIW),LIW, * IW,NSTEPS,MAXFRT,IW(3*N+1),IFLAG) C C CHANGE THE NEGATIVE SIGNS OF THE EIGENVALUES TO PLUS C 80 CONTINUE IF(NEWTON .OR. IFLAG .EQ. 3) THEN CALL DSYPRC(MAXA,LIW,A,IW(ISTIW),EPS,NEG1,NEG2) ENDIF C IF(ISLV .EQ. 3 .AND. (NEG1 .EQ. 0 .AND. NEG2 .EQ. 0)) THEN IFLAG = 30 GO TO 100 ENDIF C 90 CONTINUE IF(ISLV .EQ. 0) GO TO 100 C C SOLVE THE EQUATIONS C CD CALL MA27CD(N,A,MAXA,IW(ISTIW),LIW,W, CS CALL MA27C(N,A,MAXA,IW(ISTIW),LIW,W, * MAXFRT,RHS,IW(3*N+1),NSTEPS) 100 RETURN END SUBROUTINE STMSLV(N,DT,S,G,W1,W2,W3,BV,GAMMA,NZ,IRN,ICN, * HESS,LICN,IWA,LIWA,WK,EPS,DSCENT,IFLAG) C C PURPOSE: C ------- C C THIS ROUTINE COMPUTES THE TENSOR STEP. IT ALSO COMPUTES THE C NEWTON STEP IF THE TENSOR STEP IS NOT DESCENT. C C PARAMETERS: C ---------- C C N --> DIMENSION OF PROBLEM C DT <-- TENSOR STEP ON EXIT C S --> STEP TO PREVIOUS POINT C G --> GRADIENT AT CURRENT ITERATE C W1 --> WORKSPACE ARRAY OF DIMENSION N. IT IS USED TO HOLD C THE SOLUTION OF HESS W1 = G C W2 --> WORKSPACE ARRAY OF DIMENSION N. IT IS USED TO HOLD C THE SOLUTION OF HESS W2 = S C W3 <-- NEWTON STEP (IF REQUIRED) C BV --> (3*(S-TRANS*S)*D-2*(S-TRANS*D)*S)/(3*(S-TRANS*S)**3) C (SEE D IN SUBROUTINE STFTSM) C GAMMA --> BETA/(S-TRANS*S)**4 (SEE BETA IN SUBROUTINE STFTSM) C NZ --> NUMBER OF NONZEROS IN THE LOWER TRIANGULAR PART C OF THE HESSIAN MATRIX C IRN --> ROW INDEX OF THE NONZEROS STORED IN THE LOWER C HALF OF THE HESSIAN MATRIX C ICN --> COLUMN INDEX OF THE NONZEROS STORED IN THE LOWER C HALF OF THE HESSIAN MATRIX C HESS --> HESSIAN MATRIX C LICN --> LENGTH OF ARRAY HESS C IWA --> WORKSPACE C LIWA --> LENGTH OF ARRAY IWA C WK --> WORKSPACE OF DIMENSION N C EPS --> MACHINE PRECISION C DSCENT <-- LOGICAL VARIABLE INDICATING WHETHER OR NOT THE TENSOR C STEP IS DESCENT C IFLAG <-- OUTPUT PARAMETER FROM THE MA27 PACKAGE. A VALUE OF C ZERO INDICATES THAT THE FACTORIZATION WAS PERFORMED C SUCCESSFULLY C C BLAS SUBROUTINES: DCOPY,DDOT C INTEGER N,NZ,LICN,LIWA,IFLAG INTEGER IRN(NZ),ICN(LICN),IWA(6*N) CD DOUBLE PRECISION DT(N),S(N),G(N),W1(N),W2(N),W3(N) CD DOUBLE PRECISION BV(N),GAMMA,HESS(LICN),WK(N),EPS CS REAL DT(N),S(N),G(N),W1(N),W2(N),W3(N) CS REAL BV(N),GAMMA,HESS(LICN),WK(N),EPS LOGICAL DSCENT C C LOCAL VARIABLES C INTEGER NZEROS,I CD DOUBLE PRECISION BETA(3), THETA(3), UU, WW, VV, YY, ZZ CD DOUBLE PRECISION AA, A1, A2, A3, AMINB, AMINT,GD CD DOUBLE PRECISION DDOT, ZERO, HALF, ONE, OHALF, TWO, SIX CS REAL BETA(3), THETA(3), UU, WW, VV, YY, ZZ CS REAL AA, A1, A2, A3, AMINB, AMINT,GD CS REAL SDOT, ZERO, HALF, ONE, OHALF, TWO, SIX CD DATA ZERO, HALF, ONE, OHALF, TWO, SIX/ 0.0D0, 0.5D+0, 1.0D+0, CD * 1.5D+0, 2.0D+0, 6.0D+0 / CS DATA ZERO, HALF, ONE, OHALF, TWO, SIX/ 0.0E0, 0.5E+0, 1.0E+0, CS * 1.5E+0, 2.0E+0, 6.0E+0 / C DSCENT = .TRUE. CALL STMA27(N,NZ,IRN,ICN,HESS,LICN,IWA,LIWA, * W1,WK,0,.FALSE.,EPS,IFLAG) C C SOLVE HESS W1 = G FOR W1 C CD CALL DCOPY(N,G(1),1,W1(1),1) CS CALL SCOPY(N,G(1),1,W1(1),1) CALL STMA27(N,NZ,IRN,ICN,HESS,LICN,IWA,LIWA, * W1,WK,2,.FALSE.,EPS,IFLAG) C C SOLVE HESS W2 = S FOR W2 C CD CALL DCOPY(N,S(1),1,W2(1),1) CS CALL SCOPY(N,S(1),1,W2(1),1) CALL STMA27(N,NZ,IRN,ICN,HESS,LICN,IWA,LIWA, * W2,WK,2,.FALSE.,EPS,IFLAG) C C SOLVE HESS W3 = B FOR W3 C CD CALL DCOPY(N,BV(1),1,W3(1),1) CS CALL SCOPY(N,BV(1),1,W3(1),1) CALL STMA27(N,NZ,IRN,ICN,HESS,LICN,IWA,LIWA, * W3,WK,2,.FALSE.,EPS,IFLAG) C C COMPUTE COEFFICIENTS OF THE SYSTEM OF TWO EQUATIONS C OF THE 3RD ORDER IN 2 UNKNOWS C CD UU = DDOT(N,S(1),1,W1(1),1) CS UU = SDOT(N,S(1),1,W1(1),1) CD WW = DDOT(N,S(1),1,W2(1),1) CS WW = SDOT(N,S(1),1,W2(1),1) CD VV = DDOT(N,S(1),1,W3(1),1) CS VV = SDOT(N,S(1),1,W3(1),1) CD YY = DDOT(N,BV(1),1,W1(1),1) CS YY = SDOT(N,BV(1),1,W1(1),1) CD ZZ = DDOT(N,BV(1),1,W3(1),1) CS ZZ = SDOT(N,BV(1),1,W3(1),1) C C CALCULATE COEFFICIENTS OF THIRD DEGREE POLYNOMIAL C AA = HALF*WW*ZZ - (ONE/SIX)*GAMMA*WW - HALF*VV**2 A1 = -OHALF * VV / AA A2 = (YY * WW - UU * VV - ONE) / AA A3 = -UU / AA C C COMPUTE ZEROS OF THIRD DEGREE POLYNOMIAL C CALL STCZ3P(A1,A2,A3,BETA,NZEROS) C DO 10 I = 1,NZEROS THETA(I) = -(UU + BETA(I) + HALF * VV * BETA(I) ** 2 * + (ONE/SIX) * GAMMA * WW * BETA(I) ** 3) / (WW * BETA(I)) 10 CONTINUE C C COMPUTE SMALLEST BETA IN ABSOLUTE VALUE C AMINB = BETA(1) AMINT = THETA(1) IF(NZEROS .EQ. 2) THEN IF(AMINB .GT. ABS(BETA(2))) THEN AMINB = BETA(2) AMINT = THETA(2) ENDIF ELSEIF(NZEROS .EQ. 3) THEN IF(AMINB .GT. ABS(BETA(2))) THEN AMINB = BETA(2) AMINT = THETA(2) ENDIF IF(AMINB .GT. ABS(BETA(3))) THEN AMINB = BETA(3) AMINT = THETA(3) ENDIF ENDIF C C COMPUTE TENSOR STEP C DO 20 I = 1,N DT(I) = -W1(I) - (ONE/TWO) * W3(I) * AMINB ** 2 * -W2(I) * AMINB * AMINT * -(ONE/SIX) * W2(I) * GAMMA * AMINB ** 3 20 CONTINUE C C IF TENSOR STEP IS NOT A DESCENT DIRECTION THEN COMPUTE NEWTON STEP C CD GD = DDOT(N,G(1),1,DT(1),1) CS GD = SDOT(N,G(1),1,DT(1),1) C IF(GD .GT. ZERO) THEN C C COMPUTE THE NEWTON STEP C DSCENT = .FALSE. IF(IFLAG .EQ. 3) THEN DO 30 I = 1, N W3(I) = -W1(I) 30 CONTINUE GO TO 60 ENDIF CD CALL DCOPY(N,G(1),1,W3(1),1) CS CALL SCOPY(N,G(1),1,W3(1),1) CALL STMA27(N,NZ,IRN,ICN,HESS,LICN,IWA,LIWA, * W3,WK,3,.TRUE.,EPS,IFLAG) IF(IFLAG .EQ. 30) THEN DO 40 I = 1,N W3(I) = -W1(I) 40 CONTINUE ELSE DO 50 I = 1,N W3(I) = -W3(I) 50 CONTINUE ENDIF ENDIF C 60 CONTINUE RETURN END SUBROUTINE STRSLT(N,XPLS,FVAL,GPLS,RGX,ITN,TERMCD,IPR) C C PURPOSE: C ------- C C THIS ROUTINE PRINTS INFORMATION. C C PARAMETERS: C ---------- C C N ---> DIMENSION OF PROBLEM C XPLS ---> ITERATE TO BE PRINTED OUT C FVAL ---> FUNCTION VALUE AT XPLS C GPLS ---> GRADIENT AT XPLS C RGX ---> MAXIMUM RELATIVE GRADIENT C ITN ---> ITERATION NUMBER C IPR ---> DEVICE TO WHICH TO SEND OUTPUT C C INTEGER N,ITN,TERMCD,IPR CD DOUBLE PRECISION XPLS(N),FVAL,GPLS(N),RGX CS REAL XPLS(N),FVAL,GPLS(N),RGX C C LOCAL VARIABLES C INTEGER I,IFCNT,IGCNT,IHCNT CD DOUBLE PRECISION D,RELGRD,ZERO,ONE CS REAL D,RELGRD,ZERO,ONE COMMON/COUNT/IFCNT,IGCNT,IHCNT CD DATA ZERO,ONE/0.0D0,1.0D0/ CS DATA ZERO,ONE/0.0E0,1.0E0/ C C FIND DIRECTION IN WHICH RELATIVE GRADIENT MAXIMUM C CHECK WHETHER WITHIN TOLERANCE C IF(ITN .EQ. 0) THEN D = ONE RGX = ZERO DO 100 I = 1,N RELGRD = ABS(GPLS(I))*MAX(ABS(XPLS(I)),ONE)/D RGX = MAX(RGX,RELGRD) 100 CONTINUE ENDIF WRITE(IPR,799) WRITE(IPR,800) WRITE(IPR,801) ITN WRITE(IPR,802) WRITE(IPR,803) FVAL WRITE(IPR,804) WRITE(IPR,805) RGX WRITE(IPR,800) WRITE(IPR,799) IF(TERMCD .GT. 0) THEN WRITE(IPR,806) IFCNT,IGCNT,IHCNT ENDIF C 799 FORMAT(1H ) 800 FORMAT('---------------------------------------------') 801 FORMAT(' STRSLT ITERATION K =',I5) 802 FORMAT(' STRSLT FUNCTION AT X(K)') 803 FORMAT(' STRSLT ',E20.13) 804 FORMAT(' STRSLT SCALED GRADIENT AT X(K)') 805 FORMAT(' STRSLT ',E20.13) 806 FORMAT(' STRSLT NUMBER OF FUNCTION EVALUATIONS',I5/, * ' STRSLT NUMBER OF GRADIENT EVALUATIONS',I5/, * ' STRSLT NUMBER OF HESSIAN EVALUATIONS',I5/) C RETURN END SUBROUTINE STSORT(N,NNZ,A,INDROW,INDCOL,JPNTR,IWA) INTEGER N,NNZ INTEGER INDROW(NNZ),INDCOL(NNZ),JPNTR(N+1),IWA(N) CD DOUBLE PRECISION A(NNZ) CS REAL A(NNZ) C ********** C C SUBROUTINE STSORT C C GIVEN THE NONZERO ELEMENTS OF AN M BY N MATRIX A IN C ARBITRARY ORDER AS SPECIFIED BY THEIR ROW AND COLUMN C INDICES, THIS SUBROUTINE PERMUTES THESE ELEMENTS SO C THAT THEIR COLUMN INDICES ARE IN NON-DECREASING ORDER. C IT ALSO PERMUTES THE CORRESPONDING A ELEMENTS. C C ON INPUT IT IS ASSUMED THAT THE ELEMENTS ARE SPECIFIED IN C C INDROW(K),INDCOL(K), A(K), K = 1,...,NNZ. C C ON OUTPUT THE ELEMENTS ARE PERMUTED SO THAT INDCOL IS C IN NON-DECREASING ORDER. IN ADDITION, THE ARRAY JPNTR C IS SET SO THAT THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE STSORT(N,NNZ,A,INDROW,INDCOL,JPNTR,IWA) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF COLUMNS OF A. C C NNZ IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF NONZERO ELEMENTS OF A. C C A IS AN ARRAY OF LENGTH NNZ. ON INPUT A CONTAINS THE C THE NONZEROS ELEMENTS. ON OUTPUT A IS PERMUTED SO THAT C THE CORRESPONDING COLUMN INDICES OF INDCOL ARE IN C NON-DECREASING ORDER. C C INDROW IS AN INTEGER ARRAY OF LENGTH NNZ. ON INPUT INDROW C MUST CONTAIN THE ROW INDICES OF THE NONZERO ELEMENTS OF A. C ON OUTPUT INDROW IS PERMUTED SO THAT THE CORRESPONDING C COLUMN INDICES OF INDCOL ARE IN NON-DECREASING ORDER. C C INDCOL IS AN INTEGER ARRAY OF LENGTH NNZ. ON INPUT INDCOL C MUST CONTAIN THE COLUMN INDICES OF THE NONZERO ELEMENTS C OF A. ON OUTPUT INDCOL IS PERMUTED SO THAT THESE INDICES C ARE IN NON-DECREASING ORDER. C C JPNTR IS AN INTEGER OUTPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE ROW INDICES IN THE OUTPUT C INDROW. THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT JPNTR(1) IS SET TO 1 AND THAT JPNTR(N+1)-1 C IS THEN NNZ. C C IWA IS AN INTEGER WORK ARRAY OF LENGTH N. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... MAX C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C MODIFIED BY ALI BOUARICHA, OCTOBER 1994. C C ********** INTEGER I,J,K,L CD DOUBLE PRECISION S CS REAL S C C STORE IN ARRAY IWA THE COUNTS OF NONZEROES IN THE COLUMNS. C DO 10 J = 1, N IWA(J) = 0 10 CONTINUE DO 20 K = 1, NNZ IWA(INDCOL(K)) = IWA(INDCOL(K)) + 1 20 CONTINUE C C SET POINTERS TO THE START OF THE COLUMNS IN INDROW. C JPNTR(1) = 1 DO 30 J = 1, N JPNTR(J+1) = JPNTR(J) + IWA(J) IWA(J) = JPNTR(J) 30 CONTINUE K = 1 C C BEGIN IN-PLACE SORT. C 40 CONTINUE J = INDCOL(K) IF (K .GE. JPNTR(J)) THEN C C CURRENT ELEMENT IS IN POSITION. NOW EXAMINE THE C NEXT ELEMENT OR THE FIRST UN-SORTED ELEMENT IN C THE J-TH GROUP. C K = MAX(K+1,IWA(J)) ELSE C C CURRENT ELEMENT IS NOT IN POSITION. PLACE ELEMENT C IN POSITION AND MAKE THE DISPLACED ELEMENT THE C CURRENT ELEMENT. C L = IWA(J) IWA(J) = IWA(J) + 1 I = INDROW(K) S = A(K) INDROW(K) = INDROW(L) INDCOL(K) = INDCOL(L) A(K) = A(L) INDROW(L) = I INDCOL(L) = J A(L) = S END IF IF (K .LE. NNZ) GO TO 40 RETURN C C LAST CARD OF SUBROUTINE SRTDAT. C END CD SUBROUTINE STUMCD(N, X, NPAIRS, IRN, LIRN, ICN, LICN, CS SUBROUTINE STUMCS(N, X, NPAIRS, IRN, LIRN, ICN, LICN, * FCN, UGR, USH, TYPX, FSCALE, GRADTL, STEPTL, * ILIM, STEPMX, IPR, METHOD, GRDFLG, HSNFLG, * NDIGIT, MSG, XPLS, FPLS, GPLS, HESS, WRK, * LWRK, IWRK, LIWRK, TERMCD, VECTOR, INFORM) C C PURPOSE: C ------- C C THIS ROUTINE PROVIDES A COMPLETE INTERFACE TO STENMIN. THE USER C HAS FULL CONTROL OVER THE OPTIONS. C C PARAMETERS: C ---------- C C N --> A POSITIVE INTEGER VARIABLE SPECIFYING THE NUMBER C OF VARIABLES IN THE PROBLEM C X --> AN ARRAY OF LENGTH N THAT CONTAINS AN INITIAL C ESTIMATE OF THE MINIMIZER C NPAIRS --> AN INTEGER VARIABLE THAT MUST BE SET BY THE USER TO C THE NUMBER OF NONZEROS STORED IN THE LOWER OR UPPER C HALF OF THE HESSIAN MATRIX C IRN --> AN INTEGER ARRAY OF LENGTH LIRN. ON ENTRY, IT MUST C HOLD THE ROW INDEX OF EACH NONZERO STORED IN THE C LOWER OR UPPER HALF OF THE HESSIAN MATRIX C LIRN --> AN INTEGER VARIABLE THAT MUST BE SET BY THE USER TO C THE LENGTH OF ARRAY IRN. LIRN NEED NOT BE AS LARGE C AS LICN; NORMALLY IT NEED NOT BE VERY MUCH GREATER C THAN NPAIRS C ICN --> AN INTEGER ARRAY OF LENGTH LICN. ON ENTRY, IT MUST C HOLD THE COLUMN INDEX OF THE NONZEROS STORED IN THE C LOWER OR UPPER HALF OF THE HESSIAN MATRIX C LICN --> AN INTEGER VARIABLE THAT MUST BE SET BY THE USER TO C THE LENGTH OF THE HESSIAN ARRAY HESS AND ICN. LICN C SHOULD ORDINARILY BE 2 TO 4 TIMES AS LARGE AS NPAIRS C FCN --> THE NAME OF A USER SUPPLIED SUBROUTINE THAT EVALUATES C THE FUNCTION F AT AN ARBITRARY VECTOR X. THE SUBROUTINE C MUST BE DECLARED EXTERNAL IN THE USER'S CALLING PROGRAM C AND MUST CONFORM TO C CALL FCN(N, X, F) C UGR --> THE NAME OF A USER SUPPLIED SUBROUTINE THAT RETURNS IN C THE VALUE OF THE GRADIENT. UGR MUST BE DECLARED EXTERNAL C IN THE USER'S CALLING PROGRAM AND MUST CONFORM TO THE C USAGE C CALL UGR(N, X, G) C IF NO ANALYTIC GRADIENT IS SUPPLIED (GRDFLG = 0), THE C USER MUST USE THE DUMMY NAME STDUGR C USH --> THE NAME OF A USER SUPPLIED SUBROUTINE THAT RETURNS C IN HESS THE VALUE OF THE HESSIAN AT THE CURRENT POINT X. C USH MUST BE DECLARED EXTERNAL IN THE USER'S CALLING C PROGRAM AND MUST CONFORM TO THE USAGE C CALL USH(N,X,NPAIRS,LICN,HESS,IRN,ICN) C ONLY THE LOWER OR UPPER TRIANGULAR PART AND THE DIAGONAL C OF HESS SHOULD BE GIVEN, WITH THEIR CORRESPONDING ROW C AND COLUMN INDICES. IF NO ANALYTIC GRADIENT IS SUPPLIED C (HSNFLG = 0), THE USER MUST USE THE DUMMY NAME STDUSH C TYPX --> AN ARRAY OF LENGTH N IN WHICH THE TYPICAL SIZE OF THE C COMPONENTS OF X ARE SPECIFIED C FSCALE --> A POSITIVE REAL NUMBER ESTIMATING THE MAGNITUDE C OF F(X) NEAR THE MINIMIZER C GRADTL --> POSITIVE SCALAR GIVING THE TOLERANCE AT WHICH THE C SCALED GRADIENT OF F(X) IS CONSIDERED CLOSE ENOUGH TO C ZERO TO TERMINATE THE ALGORITHM C STEPTL --> A POSITIVE SCALAR PROVIDING THE MINIMUM ALLOWABLE C RELATIVE STEP LENGTH C ILIM --> POSITIVE INTEGER SPECIFYING THE MAXIMUM ITERATIONS TO BE C PERFORMED BEFORE THE PROGRAM IS TERMINATED C STEPMX --> A POSITIVE SCALAR PROVIDING THE MAXIMUM ALLOWABLE SCALED C STEP LENGTH C IPR --> THE UNIT ON WHICH THE ROUTINE OUTPUTS INFORMATION C METHOD --> AN INTEGER FLAG DESIGNATING WHICH METHOD TO USE C METHOD = 0 : USE NEWTON'S METHOD C METHOD = 1 : USE THE TENSOR METHOD C GRDFLG --> INTEGER FLAG DESIGNATING WHETHER OR NOT ANALYTIC HESSIAN C HAS BEEN SUPPLIED BY THE USER C GRDFLG = 0 : NO ANALYTIC GRADIENT SUPPLIED C GRDFLG = 1 : ANALYTIC GRADIENT SUPPLIED (WILL BE CHECKED C AGAINST FINITE DIFFERENCE GRADIENT) C GRDFLG = 2 : ANALYTIC GRADIENT SUPPLIED (WILL NOT BE C CHECKED AGAINST FINITE DIFFERENCE GRADIENT) C HSNFLG --> INTEGER FLAG DESIGNATING WHETHER OR NOT ANALYTIC HESSIAN C HAS BEEN SUPPLIED BY THE USER C HSNFLG = 0 : NO ANALYTIC HESSIAN SUPPLIED C HSNFLG = 1 : ANALYTIC HESSIAN SUPPLIED (WILL BE CHECKED C AGAINST FINITE DIFFERENCE HESSIAN) C HSNFLG = 2 : ANALYTIC HESSIAN SUPPLIED (WILL NOT BE C CHECKED AGAINST FINITE DIFFERENCE HESSIAN) C NDIGIT --> INTEGER ESTIMATING THE NUMBER OF ACCURATE DIGITS IN THE C OBJECTIVE FUNCTION F(X) C MSG <--> AN INTEGER VARIABLE THAT THE USER MAY SET ON INPUT TO C INHIBIT CERTAIN AUTOMATIC CHECKS OR OVERRIDE CERTAIN C DEFAULT CHARACTERISTICS OF THE PACKAGE C XPLS <-- AN ARRAY OF LENGTH N CONTAINING THE BEST APPROXIMATION C TO THE MINIMIZER UPON RETURN. (IF THE ALGORITHM HAS NOT C CONVERGED, THE LAST ITERATE IS RETURNED) C FPLS <-- A SCALAR VARIABLE THAT CONTAINS THE FUNCTION VALUE AT C THE FINAL ITERATE XPLS C GPLS <-- AN ARRAY OF LENGTH N CONTAINING THE GRADIENT VALUE C AT XPLS C HESS <-- AN ARRAY THAT IS USED TO STORE THE HESSIAN MATRIX AT C EACH ITERATION. IT NEEDS TO BE AT LEAST OF DIMENSION C LICN. ON EXIT, HESS CONTAINS THE HESSIAN VALUE AT THE C MINIMIZER C WRK --> AN ARRAY OF LENGTH LWRK. THIS IS USED AS WORKSPACE C BY THE PACKAGE. ITS LENGTH MUST BE AT LEAST 7*N C LWRK --> AN INTEGER VARIABLE. IT MUST BE SET BY THE USER TO THE C LENGTH OF ARRAY WRK AND IS NOT ALTERED BY THE PACKAGE C IWRK --> AN INTEGER ARRAY OF LENGTH LIWRK. THIS IS USED AS C WORKSPACE BY THE PACKAGE. ITS LENGTH MUST BE AT LEAST C 2*NPAIRS+12*N+2 C LIWRK --> AN INTEGER VARIABLE. IT MUST BE SET BY THE USER TO THE C LENGTH OF ARRAY IWRK AND IS NOT ALTERED BY THE PACKAGE C TERMCD <-- AN INTEGER THAT SPECIFIES THE REASON WHY THE ALGORITHM C HAS TERMINATED C VECTOR<--> AN ARRAY OF LENGTH N. IT NEED NOT BE SET BY THE USER C ON ENTRY. IF INFORM IS SET TO 1, A RE-ENTRY MUST BE C MADE WITH VECTOR SET TO HESS TIMES VECTOR (SEE INFORM) C INFORM<--> AN INTEGER VARIABLE. IF IT IS SET TO 1, THE USER MUST C OBTAIN HESS TIMES VECTOR AND RE-ENTER WITH INFORM C UNCHANGED. THE RESULT OF HESS TIMES VECTOR MUST BE C STORED IN VECTOR. THE DEFAULT VALUE OF INFORM IS 0, C MEANING THAT HESS TIMES VECTOR IS COMPUTED BY THE C PACKAGE C INTEGER N, NPAIRS, LIRN, LICN, ILIM, IPR, METHOD INTEGER GRDFLG, HSNFLG, NDIGIT, MSG, LWRK, LIWRK, TERMCD INTEGER INFORM INTEGER IRN(LIRN), ICN(LICN), IWRK(LIWRK) CD DOUBLE PRECISION X(N), TYPX(N), FSCALE, GRADTL,STEPTL, STEPMX CD DOUBLE PRECISION XPLS(N), FPLS, GPLS(N), HESS(LICN) CD DOUBLE PRECISION WRK(LWRK), VECTOR(N) CS REAL X(N), TYPX(N), FSCALE, GRADTL,STEPTL, STEPMX CS REAL XPLS(N), FPLS, GPLS(N), HESS(LICN), WRK(LWRK) CS REAL VECTOR(N) C EXTERNAL FCN, UGR, USH C C LOCAL VARIABLES C C NONE C C C EQUIVALENCE FOR DOUBLE PRECISION/REAL ARRAYS C C EQUIVALENCE WRK(1) = G(N) C WRK(N+1) = S(N) C WRK(2*N+1) = D(N) C WRK(3*N+1) = DN(N) C WRK(4*N+1) = E(N) C WRK(5*N+1) = XD(N) C WRK(6*N+1) = BV(N) C C EQUIVALENCE FOR INTEGER ARRAYS C C EQUIVALENCE IWRK(1) = LISTP(N) C IWRK(N+1) = NGRP(N) C IWRK(2*N+1) = IPNTR(N+1) C IWRK(3*N+2) = JPNTR(N+1) C IWRK(4*N+3) = IWA(2*NZ+8*N) C CALL STDRUO(N,X,NPAIRS,IRN,LIRN,ICN,LICN,FCN,UGR,USH,TYPX, * FSCALE,GRADTL,STEPTL,ILIM,STEPMX,IPR,METHOD, * GRDFLG,HSNFLG,NDIGIT,MSG,XPLS,FPLS,GPLS,HESS, * WRK(1),WRK(N+1),WRK(2*N+1),WRK(3*N+1),WRK(4*N+1), * WRK(5*N+1),WRK(6*N+1),IWRK(1),IWRK(N+1),IWRK(2*N+1), * IWRK(3*N+2),IWRK(4*N+3),LIWRK,LWRK,TERMCD,VECTOR, * INFORM) RETURN END CD SUBROUTINE STUMSD(N, X, NPAIRS, IRN, LIRN, ICN, LICN, CS SUBROUTINE STUMSS(N, X, NPAIRS, IRN, LIRN, ICN, LICN, * FCN, TYPX, MSG, XPLS, FPLS, GPLS, HESS, * WRK, LWRK, IWRK, LIWRK, TERMCD) C C PURPOSE: C ------- C C THIS ROUTINE PROVIDES A SIMPLE INTERFACE TO STENMIN. THE USER C HAS NO CONTROL OVER THE OPTIONS. C C PARAMETERS: C ---------- C C N --> A POSITIVE INTEGER VARIABLE SPECIFYING THE NUMBER OF C VARIABLES IN THE PROBLEM C X --> AN ARRAY OF LENGTH N THAT CONTAINS AN INITIAL C ESTIMATE OF THE MINIMIZER C NPAIRS --> AN INTEGER VARIABLE THAT MUST BE SET BY THE USER TO C THE NUMBER OF NONZEROS STORED IN THE LOWER OR UPPER C HALF OF THE HESSIAN MATRIX C IRN --> AN INTEGER ARRAY OF LENGTH LIRN. ON ENTRY, IT MUST C HOLD THE ROW INDEX OF EACH NONZERO STORED IN THE LOWER C OR UPPER HALF OF THE HESSIAN MATRIX C LIRN --> AN INTEGER VARIABLE THAT MUST BE SET BY THE USER TO C THE LENGTH OF ARRAY IRN. LIRN NEED NOT BE AS LARGE C AS LICN; NORMALLY IT NEED NOT BE VERY MUCH GREATER C THAN NPAIRS C ICN --> AN INTEGER ARRAY OF LENGTH LICN. ON ENTRY, IT MUST C HOLD THE COLUMN INDEX OF THE NONZEROS STORED IN THE C LOWER OR UPPER HALF OF THE HESSIAN MATRIX. C LICN --> AN INTEGER VARIABLE THAT MUST BE SET BY THE USER TO THE C LENGTH OF THE HESSIAN ARRAY HESS AND ICN. LICN SHOULD C ORDINARILY BE 2 TO 4 TIMES AS LARGE AS NPAIRS C FCN --> THE NAME OF A USER SUPPLIED SUBROUTINE THAT EVALUATES C THE FUNCTION F AT AN ARBITRARY VECTOR X. THE SUBROUTINE C MUST BE DECLARED EXTERNAL IN THE USER'S CALLING PROGRAM C AND MUST CONFORM TO C CALL FCN(N, X, F) C TYPX --> AN ARRAY OF LENGTH N IN WHICH THE TYPICAL SIZE OF THE C COMPONENTS OF X ARE SPECIFIED. IT NEED NOT BE SET ON C ENTRY C MSG <--> AN INTEGER VARIABLE THAT THE USER MAY SET ON INPUT TO C INHIBIT CERTAIN AUTOMATIC CHECKS OR OVERRIDE CERTAIN C DEFAULT CHARACTERISTICS OF THE PACKAGE C XPLS <-- AN ARRAY OF LENGTH N CONTAINING THE BEST APPROXIMATION C TO THE MINIMIZER UPON RETURN. (IF THE ALGORITHM HAS NOT C CONVERGED, THE LAST ITERATE IS RETURNED) C FPLS <-- A SCALAR VARIABLE THAT CONTAINS THE FUNCTION VALUE AT C THE FINAL ITERATE XPLS C GPLS <-- AN ARRAY OF LENGTH N CONTAINING THE GRADIENT VALUE C AT XPLS C HESS <-- AN ARRAY THAT IS USED TO STORE THE HESSIAN MATRIX AT C EACH ITERATION. IT NEEDS TO BE AT LEAST OF DIMENSION C LICN. ON EXIT, HESS CONTAINS THE HESSIAN VALUE AT THE C MINIMIZER C WRK --> AN ARRAY OF LENGTH LWRK. THIS IS USED AS WORKSPACE C BY THE PACKAGE. ITS LENGTH MUST BE AT LEAST 8*N C LWRK --> AN INTEGER VARIABLE. IT MUST BE SET BY THE USER TO THE C LENGTH OF ARRAY WRK AND IS NOT ALTERED BY THE PACKAGE C IWRK --> AN INTEGER ARRAY OF LENGTH LIWRK. THIS IS USED AS C WORKSPACE BY THE PACKAGE. ITS LENGTH MUST BE AT LEAST C 2*NPAIRS+12*N+2 C LIWRK --> AN INTEGER VARIABLE. IT MUST BE SET BY THE USER TO THE C LENGTH OF ARRAY IWRK AND IS NOT ALTERED BY THE PACKAGE C TERMCD <-- AN INTEGER THAT SPECIFIES THE REASON WHY THE ALGORITHM C HAS TERMINATED C C INTEGER N, NPAIRS, LIRN, LICN INTEGER MSG, LWRK, LIWRK, TERMCD INTEGER IRN(LIRN), ICN(LICN), IWRK(LIWRK) CD DOUBLE PRECISION X(N), TYPX(N), XPLS(N), FPLS, GPLS(N) CD DOUBLE PRECISION HESS(LICN), WRK(LWRK) CS REAL X(N), TYPX(N), XPLS(N), FPLS, GPLS(N) CS REAL HESS(LICN), WRK(LWRK) C C LOCAL VARIABLES C INTEGER ILIM, IPR, METHOD, GRDFLG, HSNFLG, NDIGIT, INFORM CD DOUBLE PRECISION FSCALE, GRADTL, STEPTL, STEPMX CS REAL FSCALE, GRADTL, STEPTL, STEPMX C EXTERNAL FCN, STDUGR, STDUSH C C EQUIVALENCE FOR DOUBLE PRECISION/REAL ARRAYS C C EQUIVALENCE WRK(1) = G(N) C WRK(N+1) = S(N) C WRK(2*N+1) = D(N) C WRK(3*N+1) = DN(N) C WRK(4*N+1) = E(N) C WRK(5*N+1) = XD(N) C WRK(6*N+1) = BV(N) C WRK(7*N+1) = VECTOR(N) C C EQUIVALENCE FOR INTEGER ARRAYS C C EQUIVALENCE IWRK(1) = LISTP(N) C IWRK(N+1) = NGRP(N) C IWRK(2*N+1) = IPNTR(N+1) C IWRK(3*N+2) = JPNTR(N+1) C IWRK(4*N+3) = IWA(2*NZ+8*N) C CALL STDFLT(N,TYPX,FSCALE,GRADTL,STEPTL,ILIM,STEPMX,IPR, * METHOD,GRDFLG,HSNFLG,NDIGIT,INFORM,MSG) C CALL STDRUO(N,X,NPAIRS,IRN,LIRN,ICN,LICN,FCN,STDUGR,STDUSH,TYPX, * FSCALE,GRADTL,STEPTL,ILIM,STEPMX,IPR,METHOD,GRDFLG, * HSNFLG,NDIGIT,MSG,XPLS,FPLS,GPLS,HESS,WRK(1),WRK(N+1), * WRK(2*N+1),WRK(3*N+1),WRK(4*N+1),WRK(5*N+1),WRK(6*N+1), * IWRK(1),IWRK(N+1),IWRK(2*N+1),IWRK(3*N+2),IWRK(4*N+3), * LIWRK,LWRK,TERMCD,WRK(7*N+1),INFORM) RETURN END CS SUBROUTINE SCOPY(N,DX,INCX,DY,INCY) CD SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) C C COPIES A VECTOR, X, TO A VECTOR, Y. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C INTEGER I,INCX,INCY,IX,IY,M,MP1,N CD DOUBLE PRECISION DX(*),DY(*) CS REAL DX(*),DY(*) C IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C 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 DY(IY) = DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,7) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DX(I) 30 CONTINUE IF( N .LT. 7 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 DY(I) = DX(I) DY(I + 1) = DX(I + 1) DY(I + 2) = DX(I + 2) DY(I + 3) = DX(I + 3) DY(I + 4) = DX(I + 4) DY(I + 5) = DX(I + 5) DY(I + 6) = DX(I + 6) 50 CONTINUE RETURN END CS REAL FUNCTION SDOT(N,DX,INCX,DY,INCY) CD DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) C 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 INTEGER I,INCX,INCY,IX,IY,M,MP1,N CD DOUBLE PRECISION DX(*),DY(*),DTEMP,ZERO CS REAL DX(*),DY(*),DTEMP,ZERO CS PARAMETER ( ZERO = 0.0E+0 ) CD PARAMETER ( ZERO = 0.0D+0 ) C CS SDOT = ZERO CD DDOT = ZERO DTEMP = ZERO IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C 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 DTEMP = DTEMP + DX(IX)*DY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE CS SDOT = DTEMP CD DDOT = DTEMP RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DTEMP + DX(I)*DY(I) 30 CONTINUE IF( N .LT. 5 ) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,5 DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) + * DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4) 50 CONTINUE 60 CONTINUE CS SDOT = DTEMP CD DDOT = DTEMP RETURN END CS REAL FUNCTION SNRM2 ( N, DX, INCX) CD DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX) INTEGER N, INCX, NEXT INTEGER I, J, NN CD DOUBLE PRECISION DX( * ), CUTLO, CUTHI, HITEST, SUM, CD * XMAX, ZERO, ONE CS REAL DX( * ), CUTLO, CUTHI, HITEST, SUM, CS * XMAX, ZERO, ONE INTRINSIC ABS, SQRT CS PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) CD PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) C C EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() 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 C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE C HOPEFULLY APPLICABLE TO ALL MACHINES. C CUTLO = MAXIMUM OF DSQRT(U/EPS) OVER ALL KNOWN MACHINES. C CUTHI = MINIMUM OF DSQRT(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 CS DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / CD DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C IF ( N .LE. 0) THEN CS SNRM2 = ZERO CD DNRM2 = ZERO ELSE NEXT = 1 SUM = ZERO NN = N * INCX C C BEGIN MAIN LOOP C I = 1 20 CONTINUE GO TO ( 30, 50, 70, 110 ), NEXT 30 CONTINUE IF( ABS( DX( I ) ) .GT. CUTLO ) GO TO 85 NEXT = 2 XMAX = ZERO C C PHASE 1. SUM IS ZERO C 50 CONTINUE IF ( DX( I ) .EQ. ZERO ) GO TO 200 IF ( ABS( DX( I ) ) .GT. CUTLO ) GO TO 85 C C PREPARE FOR PHASE 2. C NEXT = 3 GO TO 105 C C PREPARE FOR PHASE 4. C 100 CONTINUE I = J NEXT = 4 SUM = ( SUM / DX( I ) ) / DX( I ) 105 CONTINUE XMAX = ABS( DX( I ) ) SUM = SUM + ( DX( I ) / XMAX ) ** 2 GO TO 200 C C PHASE 2. SUM IS SMALL. SCALE TO AVOID DESTRUCTIVE UNDERFLOW. C 70 CONTINUE IF ( ABS( DX( I ) ) .GT. CUTLO ) THEN C C PREPARE FOR PHASE 3. C SUM = ( SUM * XMAX) * XMAX GO TO 85 ENDIF C C COMMON CODE FOR PHASES 2 AND 4. C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW. C 110 CONTINUE IF ( ABS( DX( I ) ) .GT. XMAX ) THEN SUM = ONE + SUM * ( XMAX / DX( I ) ) ** 2 XMAX = ABS( DX( I ) ) ELSE SUM = SUM + ( DX( I ) / XMAX ) ** 2 ENDIF 200 CONTINUE I = I + INCX IF ( I .LE. NN ) GO TO 20 C C END OF MAIN LOOP. C C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. C CS SNRM2 = XMAX * SQRT(SUM) CD DNRM2 = XMAX * SQRT(SUM) GO TO 300 C C FOR REAL OR D.P. SET HITEST = CUTHI/N C 85 CONTINUE HITEST = CUTHI/FLOAT( N ) C C PHASE 3. SUM IS MID-RANGE. NO SCALING. C DO 95 J = I, NN, INCX IF( ABS( DX( J ) ) .GE. HITEST ) GO TO 100 SUM = SUM + DX( J ) ** 2 95 CONTINUE CS SNRM2 = SQRT( SUM ) CD DNRM2 = SQRT( SUM ) ENDIF 300 CONTINUE RETURN END CS REAL FUNCTION SPMEPS() CD DOUBLE PRECISION FUNCTION DPMEPS() C ********** C C SUBROUTINE DPMEPS (SPMEPS) C C THIS SUBROUTINE COMPUTES THE MACHINE PRECISION PARAMETER C DPMEPS (SPMEPS) AS THE SMALLEST FLOATING POINT NUMBER SUCH THAT C 1 + DPMEPS DIFFERS FROM 1. C C THIS SUBROUTINE IS BASED ON THE SUBROUTINE MACHAR DESCRIBED IN C C W. J. CODY, MACHAR: A SUBROUTINE TO DYNAMICALLY DETERMINE C MACHINE PARAMETERS, ACM TRANS. MATH. SOFTWARE 14 (1988), 303-311. C C THE FUNCTION STATEMENT IS C C DOUBLE PRECISION DPMEPS() (DOUBLE PRECISION) C REAL SPMEPS() (REAL) C C MINPACK-2 PROJECT. FEBRUARY 1991. C ARGONNE NATIONAL LABORATORY AND UNIVERSITY OF MINNESOTA. C BRETT M. AVERICK. C C ********** INTEGER I, IBETA, IRND, IT, ITEMP, NEGEP CD DOUBLE PRECISION A, B, BETA, BETAIN, BETAH, TEMP, TEMPA, TEMP1 CD DOUBLE PRECISION ZERO, ONE, TWO CS REAL A, B, BETA, BETAIN, BETAH, TEMP, TEMPA, TEMP1 CS REAL ZERO, ONE, TWO CD DATA ZERO, ONE, TWO/0.0D0, 1.0D0, 2.0D0/ CS DATA ZERO, ONE, TWO/0.0E0, 1.0E0, 2.0E0/ C DETERMINE IBETA, BETA ALA MALCOLM. A = ONE B = ONE 10 CONTINUE A = A + A TEMP = A + ONE TEMP1 = TEMP - A IF (TEMP1-ONE .EQ. ZERO) GO TO 10 20 CONTINUE B = B + B TEMP = A + B ITEMP = INT(TEMP-A) IF (ITEMP .EQ. 0) GO TO 20 IBETA = ITEMP BETA = DBLE(IBETA) C DETERMINE IT, IRND. IT = 0 B = ONE 30 CONTINUE IT = IT + 1 B = B*BETA TEMP = B + ONE TEMP1 = TEMP - B IF (TEMP1-ONE .EQ. ZERO) GO TO 30 IRND = 0 BETAH = BETA/TWO TEMP = A + BETAH IF (TEMP-A .NE. ZERO) IRND = 1 TEMPA = A + BETA TEMP = TEMPA + BETAH IF ((IRND .EQ. 0) .AND. (TEMP-TEMPA .NE. ZERO)) IRND = 2 C DETERMINE DPMEPS (SPMEPS). NEGEP = IT + 3 BETAIN = ONE/BETA A = ONE DO 40 I = 1, NEGEP A = A*BETAIN 40 CONTINUE 50 CONTINUE TEMP = ONE + A IF (TEMP-ONE .NE. ZERO) GO TO 60 A = A*BETA GO TO 50 60 CONTINUE CD DPMEPS = A CS SPMEPS = A IF ((IBETA .EQ. 2) .OR. (IRND .EQ. 0)) GO TO 70 A = (A*(ONE+A))/TWO TEMP = ONE + A IF (TEMP-ONE .NE. ZERO) THEN CD DPMEPS = A CS SPMEPS = A ENDIF 70 CONTINUE END CS SUBROUTINE SSCAL(N,DA,DX,INCX) CD SUBROUTINE DSCAL(N,DA,DX,INCX) C C SCALES A VECTOR BY A CONSTANT. C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C INTEGER I,INCX,M,MP1,N,NINCX CD DOUBLE PRECISION DA,DX(1) CS REAL DA,DX(1) C IF(N.LE.0)RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX DX(I) = DA*DX(I) 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DX(I) = DA*DX(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 DX(I) = DA*DX(I) DX(I + 1) = DA*DX(I + 1) DX(I + 2) = DA*DX(I + 2) DX(I + 3) = DA*DX(I + 3) DX(I + 4) = DA*DX(I + 4) 50 CONTINUE RETURN END SUBROUTINE DSYPRC( LA, LIW, A, IW , EPSMCH, NEG1, NEG2 ) INTEGER LA, LIW, NEG1, NEG2 INTEGER IW ( LIW ) CD DOUBLE PRECISION A ( LA ), EPSMCH CS REAL A ( LA ), EPSMCH C C THE GILL-MURRAY-PONCELEON-SAUNDERS CODE FOR MODIFYING THE NEGATIVE C EIGEN-COMPONENTS OBTAINED WHEN FACTORIZING A SYMMETRIC INDEFINITE C MATRIX USING THE HARWELL CODE MA27. (SEE SOL 90-8, P.19-21) C C NICK GOULD, 20TH JULY 1990. C INTEGER ALEN , APOS , IBLK , NBLK , IPOS , NROWS , * NCOLS , J , K CD DOUBLE PRECISION ZERO , ONE , TWO , ALPHA , BETA , CD * GAMMA , TAU, T, C , S , E1, E2 CS REAL ZERO , ONE , TWO , ALPHA , BETA , CS * GAMMA , TAU, T, C , S , E1, E2 LOGICAL SINGLE INTRINSIC ABS , SQRT C C SET DATA. C CD PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 ) CS PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 ) C C NEG1 AND NEG2 ARE THE NUMBER OF NEGATIVE EIGENVALUES WHICH ARISE C FROM NEGATIVE 1x1 AND 2x2 BLOCK PIVOTS. C NEG1 = 0 NEG2 = 0 NBLK = ABS( IW( 1 ) ) IPOS = 2 APOS = 1 C C LOOP OVER ALL THE BLOCK PIVOTS. C DO 100 IBLK = 1, NBLK NCOLS = IW( IPOS ) IF ( NCOLS .LT. 0 ) THEN NROWS = 1 NCOLS = - NCOLS ELSE IPOS = IPOS + 1 NROWS = IW( IPOS ) END IF C C PROCESS THE DIAGONALS IN THIS BLOCK. C ALEN = NCOLS SINGLE = .TRUE. DO 50 K = IPOS + 1, IPOS + NROWS IF ( SINGLE ) THEN ALPHA = A( APOS ) J = IW( K ) SINGLE = J .GT. 0 IF ( SINGLE ) THEN C C NEGATIVE 1x1 BLOCK. C IF ( ALPHA .LT. ZERO ) THEN NEG1 = NEG1 + 1 A( APOS ) = - ALPHA ELSE IF ( ALPHA .GT. ONE / EPSMCH ) THEN NEG1 = NEG1 + 1 A( APOS ) = ONE / EPSMCH END IF END IF ELSE BETA = A( APOS + 1 ) GAMMA = A( APOS + ALEN ) C C 2x2 BLOCK: ( ALPHA BETA ) = ( C S ) ( E1 ) ( C S ) C ( BETA GAMMA ) ( S -C ) ( E2 ) ( S -C ) C IF ( ALPHA * GAMMA .LT. BETA ** 2 ) THEN TAU = ( GAMMA - ALPHA ) / ( TWO * BETA ) T = - ONE / ( ABS( TAU ) + * SQRT( ONE + TAU ** 2 ) ) IF ( TAU .LT. ZERO ) T = - T C = ONE / ( ONE + T ** 2 ) S = T * C E1 = ALPHA + BETA * T E2 = GAMMA - BETA * T C C CHANGE E1 AND E2 TO THEIR ABSOLUTE VALUES AND THEN MULTIPLY THE C THREE 2 * 2 MATRICES TO GET THE MODIFIED ALPHA, BETA AND GAMMA. C IF ( E1 .LT. ZERO ) THEN NEG2 = NEG2 + 1 E1 = - E1 END IF IF ( E2 .LT. ZERO ) THEN NEG2 = NEG2 + 1 E2 = - E2 END IF A( APOS ) = C ** 2 * E1 + S ** 2 * E2 A( APOS + 1 ) = C * S * ( E1 - E2 ) A( APOS + ALEN ) = S ** 2 * E1 + C ** 2 * E2 END IF END IF ELSE SINGLE = .TRUE. END IF APOS = APOS + ALEN ALEN = ALEN - 1 50 CONTINUE IPOS = IPOS + NCOLS + 1 100 CONTINUE RETURN END CUT HERE............ cat > colmor.f <<'CUT HERE............' SUBROUTINE DSSM(N,NPAIRS,INDROW,INDCOL,METHOD,LISTP,NGRP, * MAXGRP,MINGRP,INFO,IPNTR,JPNTR,IWA,LIWA) INTEGER N,NPAIRS,METHOD,MAXGRP,MINGRP,INFO,LIWA INTEGER INDROW(NPAIRS),INDCOL(NPAIRS),LISTP(N),NGRP(N), * IPNTR(N+1),JPNTR(N+1),IWA(LIWA) C ********** C C SUBROUTINE DSSM C C GIVEN THE SPARSITY PATTERN OF A SYMMETRIC MATRIX A OF ORDER N, C THIS SUBROUTINE DETERMINES A SYMMETRIC PERMUTATION OF A AND A C PARTITION OF THE COLUMNS OF A CONSISTENT WITH THE DETERMINATION C OF A BY A LOWER TRIANGULAR SUBSTITUTION METHOD. C C THE SPARSITY PATTERN OF THE MATRIX A IS SPECIFIED BY THE C ARRAYS INDROW AND INDCOL. ON INPUT THE INDICES FOR THE C NON-ZERO ELEMENTS IN THE LOWER TRIANGULAR PART OF A ARE C C (INDROW(K),INDCOL(K)), K = 1,2,...,NPAIRS. C C THE (INDROW(K),INDCOL(K)) PAIRS MAY BE SPECIFIED IN ANY ORDER. C DUPLICATE INPUT PAIRS ARE PERMITTED, BUT THE SUBROUTINE C ELIMINATES THEM. THE SUBROUTINE REQUIRES THAT ALL THE DIAGONAL C ELEMENTS BE PART OF THE SPARSITY PATTERN AND REPLACES ANY PAIR C (INDROW(K),INDCOL(K)) WHERE INDROW(K) IS LESS THAN INDCOL(K) C BY THE PAIR (INDCOL(K),INDROW(K)). C C THE DIRECT METHOD (METHOD = 1) FIRST DETERMINES A PARTITION C OF THE COLUMNS OF A SUCH THAT TWO COLUMNS IN A GROUP HAVE A C NON-ZERO ELEMENT IN ROW K ONLY IF COLUMN K IS IN AN EARLIER C GROUP. USING THIS PARTITION, THE SUBROUTINE THEN COMPUTES A C SYMMETRIC PERMUTATION OF A CONSISTENT WITH THE DETERMINATION C OF A BY A LOWER TRIANGULAR SUBSTITUTION METHOD. C C THE INDIRECT METHOD FIRST COMPUTES A SYMMETRIC PERMUTATION OF A C WHICH MINIMIZES THE MAXIMUM NUMBER OF NON-ZERO ELEMENTS IN ANY C ROW OF L, WHERE L IS THE LOWER TRIANGULAR PART OF THE PERMUTED C MATRIX. THE SUBROUTINE THEN PARTITIONS THE COLUMNS OF L INTO C GROUPS SUCH THAT COLUMNS OF L IN A GROUP DO NOT HAVE A NON-ZERO C IN THE SAME ROW POSITION. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE DSSM(N,NPAIRS,INDROW,INDCOL,METHOD,LISTP,NGRP, C MAXGRP,MINGRP,INFO,IPNTR,JPNTR,IWA,LIWA) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE ORDER OF A. C C NPAIRS IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF (INDROW,INDCOL) PAIRS USED TO DESCRIBE THE SPARSITY C PATTERN OF A. C C INDROW IS AN INTEGER ARRAY OF LENGTH NPAIRS. ON INPUT INDROW C MUST CONTAIN THE ROW INDICES OF THE NON-ZERO ELEMENTS IN C THE LOWER TRIANGULAR PART OF A. ON OUTPUT INDROW IS C PERMUTED SO THAT THE CORRESPONDING COLUMN INDICES ARE IN C NON-DECREASING ORDER. THE COLUMN INDICES CAN BE RECOVERED C FROM THE ARRAY JPNTR. C C INDCOL IS AN INTEGER ARRAY OF LENGTH NPAIRS. ON INPUT INDCOL C MUST CONTAIN THE COLUMN INDICES OF THE NON-ZERO ELEMENTS C IN THE LOWER TRIANGULAR PART OF A. ON OUTPUT INDCOL IS C PERMUTED SO THAT THE CORRESPONDING ROW INDICES ARE IN C NON-DECREASING ORDER. THE ROW INDICES CAN BE RECOVERED C FROM THE ARRAY IPNTR. C C METHOD IS AN INTEGER INPUT VARIABLE. IF METHOD = 1, THE C DIRECT METHOD IS USED TO DETERMINE THE PARTITION AND C SYMMETRIC PERMUTATION. OTHERWISE, THE INDIRECT METHOD IS C USED TO DETERMINE THE SYMMETRIC PERMUTATION AND PARTITION. C C LISTP IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE SYMMETRIC PERMUTATION OF THE MATRIX A. ELEMENT (I,J) C OF A IS THE (LISTP(I),LISTP(J)) ELEMENT OF THE PERMUTED C MATRIX. C C NGRP IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE PARTITION OF THE COLUMNS OF A. COLUMN J BELONGS TO C GROUP NGRP(J). C C MAXGRP IS AN INTEGER OUTPUT VARIABLE WHICH SPECIFIES THE C NUMBER OF GROUPS IN THE PARTITION OF THE COLUMNS OF A. C C MINGRP IS AN INTEGER OUTPUT VARIABLE WHICH SPECIFIES A LOWER C BOUND FOR THE NUMBER OF GROUPS IN ANY PARTITION OF THE C COLUMNS OF A CONSISTENT WITH THE DETERMINATION OF A BY A C LOWER TRIANGULAR SUBSTITUTION METHOD. C C INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS. FOR C NORMAL TERMINATION INFO = 1. C C IPNTR IS AN INTEGER OUTPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL. C THE COLUMN INDICES FOR ROW I ARE C C INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1. C C NOTE THAT IPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS IN THE LOWER TRIANGULAR PART OF THE MATRIX A. C C JPNTR IS AN INTEGER OUTPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW. C THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS IN THE LOWER TRIANGULAR PART OF THE MATRIX A. C C IWA IS AN INTEGER WORK ARRAY OF LENGTH LIWA. C C LIWA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN 6*N. C C SUBPROGRAMS CALLED C C MINPACK-SUPPLIED ... DEGR,IDO,IDOG,NUMSRT,SDPT,SEQ,SETR, C SLO,SLOG,SRTDAT C C FORTRAN-SUPPLIED ... MAX,MIN C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. DECEMBER 1984. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER I,J,JP,MAXID,MAXVD,MAXCLQ,NNZ,NUMGRP C INFO = 1 C C EXTEND THE DATA STRUCTURE TO ROWS. C CALL SETR(N,N,INDROW,JPNTR,INDCOL,IPNTR,IWA) C C DETERMINE THE SMALLEST-LAST ORDERING OF THE VERTICES OF THE C ADJACENCY GRAPH OF A, AND FROM IT DETERMINE A LOWER BOUND C FOR THE NUMBER OF GROUPS. C CALL SLOG(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(1),MAXCLQ, * MAXVD,IWA(N+1),IWA(2*N+1),IWA(3*N+1)) MINGRP = 1 + MAXVD C C USE THE SELECTED METHOD. C IF (METHOD .EQ. 1) THEN C C DIRECT METHOD. DETERMINE A PARTITION OF THE COLUMNS C OF A BY THE POWELL-TOINT METHOD. C CALL SDPT(N,INDROW,JPNTR,INDCOL,IPNTR,NGRP,MAXGRP, * IWA(N+1),IWA(2*N+1)) C C DEFINE A SYMMETRIC PERMUTATION OF A ACCORDING TO THE C ORDERING OF THE COLUMN GROUP NUMBERS IN THE PARTITION. C CALL NUMSRT(N,MAXGRP,NGRP,1,IWA(1),IWA(2*N+1),IWA(N+1)) DO 80 I = 1, N LISTP(IWA(I)) = I 80 CONTINUE ELSE C C INDIRECT METHOD. DETERMINE THE INCIDENCE DEGREE ORDERING C OF THE VERTICES OF THE ADJACENCY GRAPH OF A AND, TOGETHER C WITH THE SMALLEST-LAST ORDERING, DEFINE A SYMMETRIC C PERMUTATION OF A. C CALL IDOG(N,INDROW,JPNTR,INDCOL,IPNTR,LISTP,MAXCLQ, * MAXID,IWA(N+1),IWA(2*N+1),IWA(3*N+1)) IF (MAXID .GT. MAXVD) THEN DO 90 I = 1, N LISTP(I) = IWA(I) 90 CONTINUE END IF C C GENERATE THE SPARSITY PATTERN FOR THE LOWER C TRIANGULAR PART L OF THE PERMUTED MATRIX. C DO 110 J = 1, N DO 100 JP = JPNTR(J), JPNTR(J+1)-1 I = INDROW(JP) INDROW(JP) = MAX(LISTP(I),LISTP(J)) INDCOL(JP) = MIN(LISTP(I),LISTP(J)) 100 CONTINUE 110 CONTINUE C C SORT THE DATA STRUCTURE BY COLUMNS. C CALL SRTDAT(N,NNZ,INDROW,INDCOL,JPNTR,IWA) C C EXTEND THE DATA STRUCTURE TO ROWS. C CALL SETR(N,N,INDROW,JPNTR,INDCOL,IPNTR,IWA) C C DETERMINE THE DEGREE SEQUENCE FOR THE INTERSECTION C GRAPH OF THE COLUMNS OF L. C CALL DEGR(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(5*N+1),IWA(N+1)) C C COLOR THE INTERSECTION GRAPH OF THE COLUMNS OF L C WITH THE SMALLEST-LAST (SL) ORDERING. C CALL SLO(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(5*N+1),IWA(4*N+1), * MAXCLQ,IWA(1),IWA(N+1),IWA(2*N+1),IWA(3*N+1)) CALL SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(4*N+1),IWA(1), * MAXGRP,IWA(N+1)) DO 120 J = 1, N NGRP(J) = IWA(LISTP(J)) 120 CONTINUE C C EXIT IF THE SMALLEST-LAST ORDERING IS OPTIMAL. C IF (MAXGRP .EQ. MAXCLQ) GO TO 140 C C COLOR THE INTERSECTION GRAPH OF THE COLUMNS OF L C WITH THE INCIDENCE DEGREE (ID) ORDERING. C CALL IDO(N,N,INDROW,JPNTR,INDCOL,IPNTR,IWA(5*N+1),IWA(4*N+1), * MAXCLQ,IWA(1),IWA(N+1),IWA(2*N+1),IWA(3*N+1)) CALL SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(4*N+1),IWA(1), * NUMGRP,IWA(N+1)) C C RETAIN THE BETTER OF THE TWO ORDERINGS. C IF (NUMGRP .LT. MAXGRP) THEN MAXGRP = NUMGRP DO 130 J = 1, N NGRP(J) = IWA(LISTP(J)) 130 CONTINUE END IF 140 CONTINUE C C GENERATE THE SPARSITY PATTERN FOR THE LOWER C TRIANGULAR PART OF THE ORIGINAL MATRIX. C DO 150 J = 1, N IWA(LISTP(J)) = J 150 CONTINUE DO 170 J = 1, N DO 160 JP = JPNTR(J), JPNTR(J+1)-1 I = INDROW(JP) INDROW(JP) = MAX(IWA(I),IWA(J)) INDCOL(JP) = MIN(IWA(I),IWA(J)) 160 CONTINUE 170 CONTINUE C C SORT THE DATA STRUCTURE BY COLUMNS. C CALL SRTDAT(N,NNZ,INDROW,INDCOL,JPNTR,IWA) C C EXTEND THE DATA STRUCTURE TO ROWS. C CALL SETR(N,N,INDROW,JPNTR,INDCOL,IPNTR,IWA) END IF RETURN C C LAST CARD OF SUBROUTINE DSSM. C END SUBROUTINE IDOG(N,NGHBRP,NPNTRP,NGHBRS,NPNTRS,LISTP, * MAXCLQ,MAXID,IWA1,IWA2,IWA3) INTEGER N,MAXCLQ,MAXID INTEGER NGHBRP(*),NPNTRP(N+1),NGHBRS(*),NPNTRS(N+1),LISTP(N), * IWA1(0:N-1),IWA2(N),IWA3(N) C ********** C C SUBROUTINE IDOG C C GIVEN A LOOPLESS GRAPH G = (V,E), THIS SUBROUTINE DETERMINES C THE INCIDENCE DEGREE ORDERING OF THE VERTICES OF G. C C THE INCIDENCE DEGREE ORDERING IS DETERMINED RECURSIVELY BY C LETTING LIST(K), K = 1,...,N BE A VERTEX WITH MAXIMAL C INCIDENCE TO THE SUBGRAPH SPANNED BY THE ORDERED VERTICES. C AMONG ALL THE VERTICES OF MAXIMAL INCIDENCE, A VERTEX OF C MAXIMAL DEGREE IS CHOSEN. THIS SUBROUTINE DETERMINES THE C INVERSE OF THE INCIDENCE DEGREE ORDERING, THAT IS, AN ARRAY C LISTP SUCH THAT LISTP(LIST(K)) = K FOR K = 1,2,...,N. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE IDOG(N,NGHBRP,NPNTRP,NGHBRS,NPNTRS,LISTP, C MAXCLQ,MAXID,IWA1,IWA2,IWA3) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF VERTICES OF G. C C NGHBRP IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C PREDECESSOR ADJACENCY LISTS FOR THE GRAPH G. C C NPNTRP IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE PREDECESSOR ADJACENCY C LISTS IN NGHBRP. THE VERTICES PRECEDING AND ADJACENT C TO VERTEX J ARE C C NGHBRP(K), K = NPNTRP(J),...,NPNTRP(J+1)-1. C C NOTE THAT NPNTRP(N+1)-1 IS THEN THE NUMBER OF VERTICES C PLUS EDGES OF THE GRAPH G. C C NGHBRS IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C SUCCESSOR ADJACENCY LISTS FOR THE GRAPH G. C C NPNTRS IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE SUCCESSOR ADJACENCY C LISTS IN NGHBRS. THE VERTICES SUCCEEDING AND ADJACENT C TO VERTEX J ARE C C NGHBRS(K), K = NPNTRS(J),...,NPNTRS(J+1)-1. C C NOTE THAT NPNTRS(N+1)-1 IS THEN THE NUMBER OF VERTICES C PLUS EDGES OF THE GRAPH G. C C LISTP IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE INVERSE OF THE INCIDENCE DEGREE ORDERING OF THE C VERTICES. VERTEX J IS IN POSITION LISTP(J) OF THIS ORDERING. C C MAXCLQ IS AN INTEGER OUTPUT VARIABLE SET TO THE SIZE C OF THE LARGEST CLIQUE FOUND DURING THE ORDERING. C C MAXID IS AN INTEGER OUTPUT VARIABLE SET TO THE MAXIMUM C INCIDENCE DEGREE FOUND DURING THE ORDERING. C C IWA1,IWA2, AND IWA3 ARE INTEGER WORK ARRAYS OF LENGTH N. C C SUBPROGRAMS CALLED C C MINPACK-SUPPLIED ... NUMSRT C C FORTRAN-SUPPLIED ... MAX C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. DECEMBER 1984. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER I,J,K,MAXINC,MAXDEG,MAXLST,NCOMP,NUMDEG,NUMINC,NUMORD C C INITIALIZATION BLOCK. C DO 10 J = 1, N LISTP(J) = (NPNTRP(J+1) - NPNTRP(J) - 1) + * (NPNTRS(J+1) - NPNTRS(J) - 1) 10 CONTINUE MAXLST = (NPNTRP(N+1) + NPNTRS(N+1))/N C C SORT THE DEGREE SEQUENCE. C CALL NUMSRT(N,N-1,LISTP,1,IWA1,IWA2,IWA3) C C CREATE A DOUBLY-LINKED LIST TO ACCESS THE INCIDENCES OF THE C VERTICES. THE POINTERS FOR THE LINKED LIST ARE AS FOLLOWS. C C EACH UN-ORDERED VERTEX I IS IN A LIST (THE INCIDENCE LIST) C OF VERTICES WITH THE SAME INCIDENCE. C C IWA1(NUMINC) IS THE FIRST VERTEX IN THE NUMINC LIST C UNLESS IWA1(NUMINC) = 0. IN THIS CASE THERE ARE C NO VERTICES IN THE NUMINC LIST. C C IWA2(I) IS THE VERTEX BEFORE I IN THE INCIDENCE LIST C UNLESS IWA2(I) = 0. IN THIS CASE I IS THE FIRST C VERTEX IN THIS INCIDENCE LIST. C C IWA3(I) IS THE VERTEX AFTER I IN THE INCIDENCE LIST C UNLESS IWA3(I) = 0. IN THIS CASE I IS THE LAST C VERTEX IN THIS INCIDENCE LIST. C C IF I IS AN UN-ORDERED VERTEX, THEN -LISTP(I) IS THE C INCIDENCE OF I TO THE GRAPH INDUCED BY THE ORDERED C VERTICES. IF J IS AN ORDERED VERTEX, THEN LISTP(J) C IS THE INCIDENCE DEGREE ORDER OF VERTEX J. C MAXINC = 0 DO 20 J = 1, N I = IWA1(J-1) IWA1(J-1) = 0 IWA2(I) = 0 IWA3(I) = IWA1(0) IF (IWA1(0) .GT. 0) IWA2(IWA1(0)) = I IWA1(0) = I LISTP(J) = 0 20 CONTINUE MAXCLQ = 0 MAXID = 0 NUMORD = 1 C C BEGINNING OF ITERATION LOOP. C 30 CONTINUE C C CHOOSE A VERTEX J OF MAXIMAL DEGREE AMONG THE C VERTICES OF MAXIMAL INCIDENCE MAXINC. C 40 CONTINUE K = IWA1(MAXINC) IF (K .GT. 0) GO TO 50 MAXINC = MAXINC - 1 GO TO 40 50 CONTINUE MAXDEG = -1 DO 60 I = 1, MAXLST NUMDEG = (NPNTRP(K+1) - NPNTRP(K) - 1) + * (NPNTRS(K+1) - NPNTRS(K) - 1) IF (NUMDEG .GT. MAXDEG) THEN MAXDEG = NUMDEG J = K END IF K = IWA3(K) IF (K .LE. 0) GO TO 70 60 CONTINUE 70 CONTINUE LISTP(J) = NUMORD MAXID = MAX(MAXID,MAXINC) C C UPDATE THE SIZE OF THE LARGEST CLIQUE C FOUND DURING THE ORDERING. C IF (MAXINC .EQ. 0) NCOMP = 0 NCOMP = NCOMP + 1 IF (MAXINC+1 .EQ. NCOMP) MAXCLQ = MAX(MAXCLQ,NCOMP) C C TERMINATION TEST. C NUMORD = NUMORD + 1 IF (NUMORD .GT. N) GO TO 100 C C DELETE VERTEX J FROM THE MAXINC LIST. C IF (IWA2(J) .EQ. 0) THEN IWA1(MAXINC) = IWA3(J) ELSE IWA3(IWA2(J)) = IWA3(J) END IF IF (IWA3(J) .GT. 0) IWA2(IWA3(J)) = IWA2(J) C C DETERMINE ALL THE NEIGHBORS OF VERTEX J WHICH PRECEDE J C IN THE SUBGRAPH SPANNED BY THE UN-ORDERED VERTICES. C DO 80 K = NPNTRP(J), NPNTRP(J+1)-1 I = NGHBRP(K) C C UPDATE THE POINTERS TO THE CURRENT INCIDENCE LISTS. C NUMINC = -LISTP(I) IF (NUMINC .GE. 0) THEN LISTP(I) = LISTP(I) - 1 MAXINC = MAX(MAXINC,-LISTP(I)) C C DELETE VERTEX I FROM THE NUMINC LIST. C IF (IWA2(I) .EQ. 0) THEN IWA1(NUMINC) = IWA3(I) ELSE IWA3(IWA2(I)) = IWA3(I) END IF IF (IWA3(I) .GT. 0) IWA2(IWA3(I)) = IWA2(I) C C ADD VERTEX I TO THE NUMINC+1 LIST. C IWA2(I) = 0 IWA3(I) = IWA1(NUMINC+1) IF (IWA1(NUMINC+1) .GT. 0) IWA2(IWA1(NUMINC+1)) = I IWA1(NUMINC+1) = I END IF 80 CONTINUE C C DETERMINE ALL THE NEIGHBORS OF VERTEX J WHICH SUCCEED J C IN THE SUBGRAPH SPANNED BY THE UN-ORDERED VERTICES. C DO 90 K = NPNTRS(J), NPNTRS(J+1)-1 I = NGHBRS(K) C C UPDATE THE POINTERS TO THE CURRENT INCIDENCE LISTS. C NUMINC = -LISTP(I) IF (NUMINC .GE. 0) THEN LISTP(I) = LISTP(I) - 1 MAXINC = MAX(MAXINC,-LISTP(I)) C C DELETE VERTEX I FROM THE NUMINC LIST. C IF (IWA2(I) .EQ. 0) THEN IWA1(NUMINC) = IWA3(I) ELSE IWA3(IWA2(I)) = IWA3(I) END IF IF (IWA3(I) .GT. 0) IWA2(IWA3(I)) = IWA2(I) C C ADD VERTEX I TO THE NUMINC+1 LIST. C IWA2(I) = 0 IWA3(I) = IWA1(NUMINC+1) IF (IWA1(NUMINC+1) .GT. 0) IWA2(IWA1(NUMINC+1)) = I IWA1(NUMINC+1) = I END IF 90 CONTINUE C C END OF ITERATION LOOP. C GO TO 30 100 CONTINUE RETURN C C LAST CARD OF SUBROUTINE IDOG. C END SUBROUTINE SDPT(N,NGHBRP,NPNTRP,NGHBRS,NPNTRS,NGRP,MAXGRP, * IWA1,IWA2) INTEGER N,MAXGRP INTEGER NGHBRP(*),NPNTRP(N+1),NGHBRS(*),NPNTRS(N+1),NGRP(N), * IWA1(0:N-1),IWA2(N) C ********** C C SUBROUTINE SDPT C C GIVEN A LOOPLESS GRAPH G = (V,E), THIS SUBROUTINE DETERMINES C A SYMMETRIC COLORING OF G BY THE POWELL-TOINT DIRECT METHOD. C C THE POWELL-TOINT METHOD ASSIGNS THE K-TH COLOR BY EXAMINING C THE UN-COLORED VERTICES U(K) IN ORDER OF NON-INCREASING DEGREE C AND ASSIGNING COLOR K TO VERTEX V IF THERE ARE NO PATHS OF C LENGTH 1 OR 2 (IN THE GRAPH INDUCED BY U(K)) BETWEEN V AND C SOME K-COLORED VERTEX. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE SDPT(N,NGHBRP,NPNTRP,NGHBRS,NPNTRS,NGRP,MAXGRP, C IWA1,IWA2) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF VERTICES OF G. C C NGHBRP IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C PREDECESSOR ADJACENCY LISTS FOR THE GRAPH G. C C NPNTRP IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE PREDECESSOR ADJACENCY C LISTS IN NGHBRP. THE VERTICES PRECEDING AND ADJACENT C TO VERTEX J ARE C C NGHBRP(K), K = NPNTRP(J),...,NPNTRP(J+1)-1. C C NOTE THAT NPNTRP(N+1)-1 IS THEN THE NUMBER OF VERTICES C PLUS EDGES OF THE GRAPH G. C C NGHBRS IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C SUCCESSOR ADJACENCY LISTS FOR THE GRAPH G. C C NPNTRS IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE SUCCESSOR ADJACENCY C LISTS IN NGHBRS. THE VERTICES SUCCEEDING AND ADJACENT C TO VERTEX J ARE C C NGHBRS(K), K = NPNTRS(J),...,NPNTRS(J+1)-1. C C NOTE THAT NPNTRS(N+1)-1 IS THEN THE NUMBER OF VERTICES C PLUS EDGES OF THE GRAPH G. C C NGRP IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE SYMMETRIC COLORING OF G. VERTEX J IS COLORED WITH C COLOR NGRP(J). C C MAXGRP IS AN INTEGER OUTPUT VARIABLE WHICH SPECIFIES THE C NUMBER OF COLORS IN THE SYMMETRIC COLORING OF G. C C IWA1 AND IWA2 ARE INTEGER WORK ARRAYS OF LENGTH N. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... MAX C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. DECEMBER 1984. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER J,JP,K,KP,L,MAXDEG,NUMDEG,NUMV C C INITIALIZATION BLOCK. NUMV IS THE CURRENT NUMBER OF UN-COLORED C VERTICES, MAXDEG IS THE MAXIMUM INDUCED DEGREE OF THESE C VERTICES, AND MAXGRP IS THE CURRENT GROUP NUMBER (COLOR). C NUMV = N MAXDEG = 0 DO 10 J = 1, N NGRP(J) = (NPNTRP(J) - NPNTRP(J+1) + 1) + * (NPNTRS(J) - NPNTRS(J+1) + 1) MAXDEG = MAX(MAXDEG,-NGRP(J)) IWA2(J) = -J 10 CONTINUE MAXGRP = 0 C C BEGINNING OF ITERATION LOOP. C 20 CONTINUE C C SORT THE LIST OF UN-COLORED VERTICES SO THAT THEIR C INDUCED DEGREES ARE IN NON-DECREASING ORDER. C DO 30 NUMDEG = 0, MAXDEG IWA1(NUMDEG) = 0 30 CONTINUE DO 40 L = 1, NUMV NUMDEG = -NGRP(-IWA2(L)) IWA1(NUMDEG) = IWA1(NUMDEG) + 1 40 CONTINUE K = 1 DO 50 NUMDEG = MAXDEG, 0, -1 L = IWA1(NUMDEG) IWA1(NUMDEG) = K K = K + L 50 CONTINUE K = 1 60 CONTINUE J = IWA2(K) IF (J .GT. 0) THEN K = IWA1(-NGRP(J)) ELSE NUMDEG = -NGRP(-J) L = IWA1(NUMDEG) IWA2(K) = IWA2(L) IWA2(L) = -J IWA1(NUMDEG) = IWA1(NUMDEG) + 1 END IF IF (K .LE. NUMV) GO TO 60 MAXGRP = MAXGRP + 1 C C DETERMINE THE VERTICES IN GROUP MAXGRP. C DO 160 L = 1, NUMV J = IWA2(L) C C EXAMINE EACH VERTEX K PRECEDING VERTEX J AND ALL C THE NEIGHBORS OF VERTEX K TO DETERMINE IF VERTEX C J CAN BE CONSIDERED FOR GROUP MAXGRP. C DO 90 JP = NPNTRP(J), NPNTRP(J+1)-1 K = NGHBRP(JP) IF (NGRP(K) .EQ. MAXGRP) GO TO 150 IF (NGRP(K) .LE. 0) THEN DO 70 KP = NPNTRP(K), NPNTRP(K+1)-1 IF (NGRP(NGHBRP(KP)) .EQ. MAXGRP) GO TO 150 70 CONTINUE DO 80 KP = NPNTRS(K), NPNTRS(K+1)-1 IF (NGRP(NGHBRS(KP)) .EQ. MAXGRP) GO TO 150 80 CONTINUE END IF 90 CONTINUE C C EXAMINE EACH VERTEX K SUCCEEDING VERTEX J AND ALL C THE NEIGHBORS OF VERTEX K TO DETERMINE IF VERTEX C J CAN BE ADDED TO GROUP MAXGRP. C DO 120 JP = NPNTRS(J), NPNTRS(J+1)-1 K = NGHBRS(JP) IF (NGRP(K) .EQ. MAXGRP) GO TO 150 IF (NGRP(K) .LE. 0) THEN DO 100 KP = NPNTRP(K), NPNTRP(K+1)-1 IF (NGRP(NGHBRP(KP)) .EQ. MAXGRP) GO TO 150 100 CONTINUE DO 110 KP = NPNTRS(K), NPNTRS(K+1)-1 IF (NGRP(NGHBRS(KP)) .EQ. MAXGRP) GO TO 150 110 CONTINUE END IF 120 CONTINUE C C ADD VERTEX J TO GROUP MAXGRP AND REMOVE VERTEX J C FROM THE LIST OF UN-COLORED VERTICES. C NGRP(J) = MAXGRP IWA2(L) = 0 C C UPDATE THE DEGREES OF THE NEIGHBORS OF VERTEX J. C DO 130 JP = NPNTRP(J), NPNTRP(J+1)-1 K = NGHBRP(JP) IF (NGRP(K) .LT. 0) NGRP(K) = NGRP(K) + 1 130 CONTINUE DO 140 JP = NPNTRS(J), NPNTRS(J+1)-1 K = NGHBRS(JP) IF (NGRP(K) .LT. 0) NGRP(K) = NGRP(K) + 1 140 CONTINUE 150 CONTINUE 160 CONTINUE C C COMPRESS THE UPDATED LIST OF UN-COLORED VERTICES. C RESET NUMV AND RECOMPUTE MAXDEG. C K = 0 MAXDEG = 0 DO 170 L = 1, NUMV IF (IWA2(L) .NE. 0) THEN K = K + 1 IWA2(K) = -IWA2(L) MAXDEG = MAX(MAXDEG,-NGRP(IWA2(L))) END IF 170 CONTINUE NUMV = K C C END OF ITERATION LOOP. C IF (NUMV .GT. 0) GO TO 20 RETURN C C LAST CARD OF SUBROUTINE SDPT. C END SUBROUTINE SLOG(N,NGHBRP,NPNTRP,NGHBRS,NPNTRS,LISTP, * MAXCLQ,MAXVD,IWA1,IWA2,IWA3) INTEGER N,MAXCLQ,MAXVD INTEGER NGHBRP(*),NPNTRP(N+1),NGHBRS(*),NPNTRS(N+1),LISTP(N), * IWA1(0:N-1),IWA2(N),IWA3(N) C ********** C C SUBROUTINE SLOG C C GIVEN A LOOPLESS GRAPH G = (V,E), THIS SUBROUTINE DETERMINES C THE SMALLEST-LAST ORDERING OF THE VERTICES OF G. C C THE SMALLEST-LAST ORDERING IS DETERMINED RECURSIVELY BY C LETTING LIST(K), K = N,...,1 BE A VERTEX WITH LEAST DEGREE C IN THE SUBGRAPH SPANNED BY THE UN-ORDERED VERTICES. C THIS SUBROUTINE DETERMINES THE INVERSE OF THE SMALLEST-LAST C ORDERING, THAT IS, AN ARRAY LISTP SUCH THAT LISTP(LIST(K)) = K C FOR K = 1,2,...,N. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE SLOG(N,NGHBRP,NPNTRP,NGHBRS,NPNTRS,LISTP, C MAXCLQ,MAXVD,IWA1,IWA2,IWA3) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF VERTICES OF G. C C NGHBRP IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C PREDECESSOR ADJACENCY LISTS FOR THE GRAPH G. C C NPNTRP IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE PREDECESSOR ADJACENCY C LISTS IN NGHBRP. THE VERTICES PRECEDING AND ADJACENT C TO VERTEX J ARE C C NGHBRP(K), K = NPNTRP(J),...,NPNTRP(J+1)-1. C C NOTE THAT NPNTRP(N+1)-1 IS THEN THE NUMBER OF VERTICES C PLUS EDGES OF THE GRAPH G. C C NGHBRS IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C SUCCESSOR ADJACENCY LISTS FOR THE GRAPH G. C C NPNTRS IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE SUCCESSOR ADJACENCY C LISTS IN NGHBRS. THE VERTICES SUCCEEDING AND ADJACENT C TO VERTEX J ARE C C NGHBRS(K), K = NPNTRS(J),...,NPNTRS(J+1)-1. C C NOTE THAT NPNTRS(N+1)-1 IS THEN THE NUMBER OF VERTICES C PLUS EDGES OF THE GRAPH G. C C LISTP IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE INVERSE OF THE SMALLEST-LAST ORDERING OF THE VERTICES. C VERTEX J IS IN POSITION LISTP(J) OF THIS ORDERING. C C MAXCLQ IS AN INTEGER OUTPUT VARIABLE SET TO THE SIZE C OF THE LARGEST CLIQUE FOUND DURING THE ORDERING. C C MAXVD IS AN INTEGER OUTPUT VARIABLE SET TO THE MAXIMUM C VERTEX DEGREE FOUND DURING THE ORDERING. C C IWA1,IWA2, AND IWA3 ARE INTEGER WORK ARRAYS OF LENGTH N. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... MAX,MIN C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. DECEMBER 1984. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER I,J,K,MINDEG,NUMDEG,NUMORD C C INITIALIZATION BLOCK. C MINDEG = N DO 10 J = 1, N IWA1(J-1) = 0 LISTP(J) = (NPNTRP(J) - NPNTRP(J+1) + 1) + * (NPNTRS(J) - NPNTRS(J+1) + 1) MINDEG = MIN(MINDEG,-LISTP(J)) 10 CONTINUE C C CREATE A DOUBLY-LINKED LIST TO ACCESS THE DEGREES OF THE C VERTICES. THE POINTERS FOR THE LINKED LIST ARE AS FOLLOWS. C C EACH UN-ORDERED VERTEX I IS IN A LIST (THE DEGREE LIST) C OF VERTICES WITH THE SAME DEGREE. C C IWA1(NUMDEG) IS THE FIRST VERTEX IN THE NUMDEG LIST C UNLESS IWA1(NUMDEG) = 0. IN THIS CASE THERE ARE C NO VERTICES IN THE NUMDEG LIST. C C IWA2(I) IS THE VERTEX BEFORE I IN THE DEGREE LIST C UNLESS IWA2(I) = 0. IN THIS CASE I IS THE FIRST C VERTEX IN THIS DEGREE LIST. C C IWA3(I) IS THE VERTEX AFTER I IN THE DEGREE LIST C UNLESS IWA3(I) = 0. IN THIS CASE I IS THE LAST C VERTEX IN THIS DEGREE LIST. C C IF I IS AN UN-ORDERED VERTEX, THEN -LISTP(I) IS THE C DEGREE OF I IN THE GRAPH INDUCED BY THE UN-ORDERED C VERTICES. IF J IS AN ORDERED VERTEX, THEN LISTP(J) C IS THE SMALLEST-LAST ORDER OF VERTEX J. C DO 20 J = 1, N NUMDEG = -LISTP(J) IWA2(J) = 0 IWA3(J) = IWA1(NUMDEG) IF (IWA1(NUMDEG) .GT. 0) IWA2(IWA1(NUMDEG)) = J IWA1(NUMDEG) = J 20 CONTINUE MAXCLQ = 0 MAXVD = 0 NUMORD = N C C BEGINNING OF ITERATION LOOP. C 30 CONTINUE C C CHOOSE A VERTEX J OF MINIMAL DEGREE MINDEG. C 40 CONTINUE J = IWA1(MINDEG) IF (J .GT. 0) GO TO 50 MINDEG = MINDEG + 1 GO TO 40 50 CONTINUE LISTP(J) = NUMORD MAXVD = MAX(MAXVD,MINDEG) C C MARK THE SIZE OF THE LARGEST CLIQUE C FOUND DURING THE ORDERING. C IF (MINDEG+1 .EQ. NUMORD .AND. MAXCLQ .EQ. 0) * MAXCLQ = NUMORD C C TERMINATION TEST. C NUMORD = NUMORD - 1 IF (NUMORD .EQ. 0) GO TO 80 C C DELETE VERTEX J FROM THE MINDEG LIST. C IWA1(MINDEG) = IWA3(J) IF (IWA3(J) .GT. 0) IWA2(IWA3(J)) = 0 C C DETERMINE ALL THE NEIGHBORS OF VERTEX J WHICH PRECEDE J C IN THE SUBGRAPH SPANNED BY THE UN-ORDERED VERTICES. C DO 60 K = NPNTRP(J), NPNTRP(J+1)-1 I = NGHBRP(K) C C UPDATE THE POINTERS TO THE CURRENT DEGREE LISTS. C NUMDEG = -LISTP(I) IF (NUMDEG .GE. 0) THEN LISTP(I) = LISTP(I) + 1 MINDEG = MIN(MINDEG,-LISTP(I)) C C DELETE VERTEX I FROM THE NUMDEG LIST. C IF (IWA2(I) .EQ. 0) THEN IWA1(NUMDEG) = IWA3(I) ELSE IWA3(IWA2(I)) = IWA3(I) END IF IF (IWA3(I) .GT. 0) IWA2(IWA3(I)) = IWA2(I) C C ADD VERTEX I TO THE NUMDEG-1 LIST. C IWA2(I) = 0 IWA3(I) = IWA1(NUMDEG-1) IF (IWA1(NUMDEG-1) .GT. 0) IWA2(IWA1(NUMDEG-1)) = I IWA1(NUMDEG-1) = I END IF 60 CONTINUE C C DETERMINE ALL THE NEIGHBORS OF VERTEX J WHICH SUCCEED J C IN THE SUBGRAPH SPANNED BY THE UN-ORDERED VERTICES. C DO 70 K = NPNTRS(J), NPNTRS(J+1)-1 I = NGHBRS(K) C C UPDATE THE POINTERS TO THE CURRENT DEGREE LISTS. C NUMDEG = -LISTP(I) IF (NUMDEG .GE. 0) THEN LISTP(I) = LISTP(I) + 1 MINDEG = MIN(MINDEG,-LISTP(I)) C C DELETE VERTEX I FROM THE NUMDEG LIST. C IF (IWA2(I) .EQ. 0) THEN IWA1(NUMDEG) = IWA3(I) ELSE IWA3(IWA2(I)) = IWA3(I) END IF IF (IWA3(I) .GT. 0) IWA2(IWA3(I)) = IWA2(I) C C ADD VERTEX I TO THE NUMDEG-1 LIST. C IWA2(I) = 0 IWA3(I) = IWA1(NUMDEG-1) IF (IWA1(NUMDEG-1) .GT. 0) IWA2(IWA1(NUMDEG-1)) = I IWA1(NUMDEG-1) = I END IF 70 CONTINUE C C END OF ITERATION LOOP. C GO TO 30 80 CONTINUE RETURN C C LAST CARD OF SUBROUTINE SLOG. C END SUBROUTINE FDHS(N,INDROW,JPNTR,INDCOL,IPNTR,LISTP,NGRP, * MAXGRP,NUMGRP,ETA,FHESD,FHES,IWA) INTEGER N,MAXGRP,NUMGRP INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(N+1), * LISTP(N),NGRP(N),IWA(N) CD DOUBLE PRECISION ETA(N),FHESD(N),FHES(*) CS REAL ETA(N),FHESD(N),FHES(*) C ********** C C SUBROUTINE FDHS C C THIS SUBROUTINE COMPUTES AN APPROXIMATION TO THE (SYMMETRIC) C HESSIAN MATRIX OF A FUNCTION BY A SUBSTITUTION METHOD. C THE LOWER TRIANGULAR PART OF THE APPROXIMATION IS STORED C WITH A COLUMN-ORIENTED DEFINITION OF THE SPARSITY PATTERN. C C THIS SUBROUTINE REQUIRES A SYMMETRIC PERMUTATION OF THE C HESSIAN MATRIX AND A PARTITION OF THE COLUMNS OF THE HESSIAN C MATRIX CONSISTENT WITH THE DETERMINATION OF THE HESSIAN C MATRIX BY A LOWER TRIANGULAR SUBSTITUTION METHOD. C THIS INFORMATION CAN BE PROVIDED BY SUBROUTINE DSSM. C C THE SYMMETRIC PERMUTATION OF THE HESSIAN MATRIX IS DEFINED C BY THE ARRAY LISTP. THIS ARRAY IS ONLY USED INTERNALLY. C C THE PARTITION OF THE HESSIAN MATRIX IS DEFINED BY THE ARRAY C NGRP BY SETTING NGRP(J) TO THE GROUP NUMBER OF COLUMN J. C THE USER MUST PROVIDE AN APPROXIMATION TO THE COLUMNS OF C THE HESSIAN MATRIX IN EACH GROUP BY SPECIFYING A DIFFERENCE C PARAMETER VECTOR ETA AND AN APPROXIMATION TO H*D WHERE H IS C THE HESSIAN MATRIX AND THE VECTOR D IS DEFINED BY THE C FOLLOWING SECTION OF CODE. C C DO 10 J = 1, N C D(J) = 0.0 C IF (NGRP(J) .EQ. NUMGRP) D(J) = ETA(J) C 10 CONTINUE C C IN THE ABOVE CODE NUMGRP IS A GROUP NUMBER AND ETA(J) IS THE C DIFFERENCE PARAMETER USED TO APPROXIMATE COLUMN J OF THE C HESSIAN MATRIX. SUITABLE VALUES FOR ETA(J) MUST BE PROVIDED. C C AS MENTIONED ABOVE, AN APPROXIMATION TO H*D MUST BE PROVIDED. C FOR EXAMPLE, IF GRAD(X) IS THE GRADIENT OF THE FUNCTION AT X, C THEN C C GRAD(X+D) - GRAD(X) C C CORRESPONDS TO THE FORWARD DIFFERENCE APPROXIMATION. C C THE LOWER TRIANGULAR SUBSTITUTION METHOD REQUIRES THAT THE C APPROXIMATIONS TO H*D FOR ALL THE GROUPS BE STORED IN SPECIAL C LOCATIONS OF THE ARRAY FHES. THIS IS DONE BY CALLING FDHS C SUCCESSIVELY WITH NUMGRP = 1,2,...,MAXGRP. ON THE CALL WITH C NUMGRP = MAXGRP, THE SUBROUTINE THEN PROCEEDS TO OVERWRITE C FHES WITH THE APPROXIMATION TO THE LOWER TRIANGULAR PART OF C THE HESSIAN MATRIX. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE FDHS(N,INDROW,JPNTR,INDCOL,IPNTR,LISTP,NGRP, C MAXGRP,NUMGRP,ETA,FHESD,FHES,IWA) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE ORDER C OF THE HESSIAN MATRIX. C C INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW C INDICES FOR THE NON-ZEROES IN THE LOWER TRIANGULAR PART C OF THE HESSIAN MATRIX. C C JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW. C THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZEROES C IN THE LOWER TRIANGULAR PART OF THE HESSIAN MATRIX. C C INDCOL IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE COLUMN C INDICES FOR THE NON-ZEROES IN THE LOWER TRIANGULAR PART C OF THE HESSIAN MATRIX. C C IPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL. C THE COLUMN INDICES FOR ROW I ARE C C INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1. C C NOTE THAT IPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZEROES C IN THE LOWER TRIANGULAR PART OF THE HESSIAN MATRIX. C C LISTP IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE SYMMETRIC PERMUTATION OF THE HESSIAN MATRIX. ELEMENT C (I,J) OF THE HESSIAN MATRIX IS THE (LISTP(I),LISTP(J)) C ELEMENT OF THE PERMUTED HESSIAN. C C NGRP IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE PARTITION OF THE COLUMNS OF THE HESSIAN MATRIX. C COLUMN J BELONGS TO GROUP NGRP(J). C C MAXGRP IS A POSITIVE INTEGER INPUT VARIABLE WHICH SPECIFIES C THE NUMBER OF GROUPS IN THE PARTITION OF THE COLUMNS OF C THE HESSIAN MATRIX. C C NUMGRP IS A POSITIVE INTEGER INPUT VARIABLE SET TO A GROUP C NUMBER IN THE PARTITION. C C ETA IS AN INPUT ARRAY OF LENGTH N WHICH CONTAINS THE C DIFFERENCE PARAMETER VECTOR. C C FHESD IS AN INPUT ARRAY OF LENGTH N WHICH CONTAINS AN C APPROXIMATION TO H*D, WHERE H IS THE HESSIAN MATRIX C AND D IS THE DIFFERENCE VECTOR FOR GROUP NUMGRP. C C FHES IS AN OUTPUT ARRAY OF LENGTH NNZ, WHERE NNZ IS THE C NUMBER OF NON-ZERO ELEMENTS IN THE LOWER TRIANGULAR PART C OF THE HESSIAN MATRIX. ON OUTPUT WITH NUMGRP LESS THAN C MAXGRP, THE FHESD ARRAY FOR GROUP NUMGRP HAS BEEN STORED C IN FHES. WHEN NUMGRP = MAXGRP THE SUBROUTINE OVERWRITES C FHES WITH AN APPROXIMATION TO THE LOWER TRIANGULAR PART C OF THE HESSIAN MATRIX. THE APPROXIMATION IS STORED IN C FHES WITH A COLUMN-ORIENTED DEFINITION OF THE SPARSITY C PATTERN. THUS THE ELEMENTS IN COLUMN J OF THE LOWER C TRIANGULAR PART OF THE HESSIAN MATRIX ARE C C FHES(K), K = JPNTR(J),...,JPNTR(J+1)-1, C C AND THE ROW INDICES FOR THESE ELEMENTS ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C IWA IS AN INTEGER WORK ARRAY OF LENGTH N. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... ABS C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. DECEMBER 1984. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER I,IP,IROW,J,JP,K,L,NUMG,NUML CD DOUBLE PRECISION SUM CS REAL SUM C C STORE THE I-TH ELEMENT OF GRADIENT DIFFERENCE FHESD C CORRESPONDING TO GROUP NUMGRP IF THERE IS A POSITION C (I,J) SUCH THAT NGRP(J) = NUMGRP AND (I,J) IS MAPPED C ONTO THE LOWER TRIANGULAR PART OF THE PERMUTED MATRIX. C DO 50 J = 1, N IF (NGRP(J) .EQ. NUMGRP) THEN NUML = LISTP(J) DO 30 IP = IPNTR(J), IPNTR(J+1)-1 I = INDCOL(IP) IF (LISTP(I) .GT. NUML) THEN DO 10 JP = JPNTR(I), JPNTR(I+1)-1 IF (INDROW(JP) .EQ. J) THEN FHES(JP) = FHESD(I) GO TO 20 END IF 10 CONTINUE 20 CONTINUE END IF 30 CONTINUE DO 40 JP = JPNTR(J), JPNTR(J+1)-1 I = INDROW(JP) IF (LISTP(I) .GE. NUML) FHES(JP) = FHESD(I) 40 CONTINUE END IF 50 CONTINUE C C EXIT IF THIS IS NOT THE LAST GROUP. C IF (NUMGRP .LT. MAXGRP) RETURN C C MARK ALL COLUMN INDICES J SUCH THAT (I,J) IS MAPPED ONTO C THE LOWER TRIANGULAR PART OF THE PERMUTED MATRIX. C DO 80 I = 1, N NUML = LISTP(I) DO 60 IP = IPNTR(I), IPNTR(I+1)-1 J = INDCOL(IP) IF (NUML .GE. LISTP(J)) INDCOL(IP) = -INDCOL(IP) 60 CONTINUE DO 70 JP = JPNTR(I), JPNTR(I+1)-1 J = INDROW(JP) IF (NUML .GT. LISTP(J)) INDROW(JP) = -INDROW(JP) 70 CONTINUE 80 CONTINUE C C INVERT THE ARRAY LISTP. C DO 90 J = 1, N IWA(LISTP(J)) = J 90 CONTINUE DO 100 J = 1, N LISTP(J) = IWA(J) 100 CONTINUE C C DETERMINE THE LOWER TRIANGULAR PART OF THE ORIGINAL MATRIX. C DO 220 IROW = N, 1, -1 I = LISTP(IROW) C C FIND THE POSITIONS OF THE ELEMENTS IN THE I-TH ROW OF THE C LOWER TRIANGULAR PART OF THE ORIGINAL MATRIX THAT HAVE C ALREADY BEEN DETERMINED. C DO 130 IP = IPNTR(I), IPNTR(I+1)-1 J = INDCOL(IP) IF (J .GT. 0) THEN DO 110 JP = JPNTR(J), JPNTR(J+1)-1 IF (INDROW(JP) .EQ. I) THEN IWA(J) = JP GO TO 120 END IF 110 CONTINUE 120 CONTINUE END IF 130 CONTINUE C C DETERMINE THE ELEMENTS IN THE I-TH ROW OF THE LOWER C TRIANGULAR PART OF THE ORIGINAL MATRIX WHICH GET MAPPED C ONTO THE LOWER TRIANGULAR PART OF THE PERMUTED MATRIX. C DO 180 K = IPNTR(I), IPNTR(I+1)-1 J = -INDCOL(K) IF (J .GT. 0) THEN INDCOL(K) = J C C DETERMINE THE (I,J) ELEMENT. C NUMG = NGRP(J) SUM = 0.0 DO 140 IP = IPNTR(I), IPNTR(I+1)-1 L = ABS(INDCOL(IP)) IF (NGRP(L) .EQ. NUMG .AND. L .NE. J) * SUM = SUM + FHES(IWA(L))*ETA(L) 140 CONTINUE DO 150 JP = JPNTR(I), JPNTR(I+1)-1 L = ABS(INDROW(JP)) IF (NGRP(L) .EQ. NUMG .AND. L .NE. J) * SUM = SUM + FHES(JP)*ETA(L) 150 CONTINUE C C STORE THE (I,J) ELEMENT. C DO 160 JP = JPNTR(J), JPNTR(J+1)-1 IF (INDROW(JP) .EQ. I) THEN FHES(JP) = (FHES(JP) - SUM)/ETA(J) GO TO 170 END IF 160 CONTINUE 170 CONTINUE END IF 180 CONTINUE C C DETERMINE THE ELEMENTS IN THE I-TH ROW OF THE STRICT UPPER C TRIANGULAR PART OF THE ORIGINAL MATRIX WHICH GET MAPPED C ONTO THE LOWER TRIANGULAR PART OF THE PERMUTED MATRIX. C DO 210 K = JPNTR(I), JPNTR(I+1)-1 J = -INDROW(K) IF (J .GT. 0) THEN INDROW(K) = J C C DETERMINE THE (I,J) ELEMENT. C NUMG = NGRP(J) SUM = 0.0 DO 190 IP = IPNTR(I), IPNTR(I+1)-1 L = ABS(INDCOL(IP)) IF (NGRP(L) .EQ. NUMG) * SUM = SUM + FHES(IWA(L))*ETA(L) 190 CONTINUE DO 200 JP = JPNTR(I), JPNTR(I+1)-1 L = ABS(INDROW(JP)) IF (NGRP(L) .EQ. NUMG .AND. L .NE. J) * SUM = SUM + FHES(JP)*ETA(L) 200 CONTINUE C C STORE THE (I,J) ELEMENT. C FHES(K) = (FHES(K) - SUM)/ETA(J) END IF 210 CONTINUE 220 CONTINUE C C RE-INVERT THE ARRAY LISTP. C DO 230 J = 1, N IWA(LISTP(J)) = J 230 CONTINUE DO 240 J = 1, N LISTP(J) = IWA(J) 240 CONTINUE RETURN C C LAST CARD OF SUBROUTINE FDHS. C END SUBROUTINE DEGR(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,IWA) INTEGER N INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(*),NDEG(N),IWA(N) C ********** C C SUBROUTINE DEGR C C GIVEN THE SPARSITY PATTERN OF AN M BY N MATRIX A, C THIS SUBROUTINE DETERMINES THE DEGREE SEQUENCE FOR C THE INTERSECTION GRAPH OF THE COLUMNS OF A. C C IN GRAPH-THEORY TERMINOLOGY, THE INTERSECTION GRAPH OF C THE COLUMNS OF A IS THE LOOPLESS GRAPH G WITH VERTICES C A(J), J = 1,2,...,N WHERE A(J) IS THE J-TH COLUMN OF A C AND WITH EDGE (A(I),A(J)) IF AND ONLY IF COLUMNS I AND J C HAVE A NON-ZERO IN THE SAME ROW POSITION. C C NOTE THAT THE VALUE OF M IS NOT NEEDED BY DEGR AND IS C THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE DEGR(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,IWA) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF COLUMNS OF A. C C INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW C INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW. C THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS OF THE MATRIX A. C C INDCOL IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C IPNTR IS AN INTEGER INPUT ARRAY OF LENGTH M + 1 WHICH C SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL. C THE COLUMN INDICES FOR ROW I ARE C C INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1. C C NOTE THAT IPNTR(M+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS OF THE MATRIX A. C C NDEG IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH C SPECIFIES THE DEGREE SEQUENCE. THE DEGREE OF THE C J-TH COLUMN OF A IS NDEG(J). C C IWA IS AN INTEGER WORK ARRAY OF LENGTH N. C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER IC,IP,IR,JCOL,JP C C INITIALIZATION BLOCK. C DO 10 JP = 1, N NDEG(JP) = 0 IWA(JP) = 0 10 CONTINUE C C COMPUTE THE DEGREE SEQUENCE BY DETERMINING THE CONTRIBUTIONS C TO THE DEGREES FROM THE CURRENT(JCOL) COLUMN AND FURTHER C COLUMNS WHICH HAVE NOT YET BEEN CONSIDERED. C DO 40 JCOL = 2, N IWA(JCOL) = N C C DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND C TO NON-ZEROES IN THE MATRIX. C DO 30 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1 IR = INDROW(JP) C C FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC) C WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX. C DO 20 IP = IPNTR(IR), IPNTR(IR+1)-1 IC = INDCOL(IP) C C ARRAY IWA MARKS COLUMNS WHICH HAVE CONTRIBUTED TO C THE DEGREE COUNT OF COLUMN JCOL. UPDATE THE DEGREE C COUNTS OF THESE COLUMNS AS WELL AS COLUMN JCOL. C IF (IWA(IC) .LT. JCOL) THEN IWA(IC) = JCOL NDEG(IC) = NDEG(IC) + 1 NDEG(JCOL) = NDEG(JCOL) + 1 END IF 20 CONTINUE 30 CONTINUE 40 CONTINUE RETURN C C LAST CARD OF SUBROUTINE DEGR. C END SUBROUTINE IDO(M,N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST, * MAXCLQ,IWA1,IWA2,IWA3,IWA4) INTEGER M,N,MAXCLQ INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(M+1),NDEG(N), * LIST(N),IWA1(0:N-1),IWA2(N),IWA3(N),IWA4(N) C ********** C C SUBROUTINE IDO C C GIVEN THE SPARSITY PATTERN OF AN M BY N MATRIX A, THIS C SUBROUTINE DETERMINES AN INCIDENCE-DEGREE ORDERING OF THE C COLUMNS OF A. C C THE INCIDENCE-DEGREE ORDERING IS DEFINED FOR THE LOOPLESS C GRAPH G WITH VERTICES A(J), J = 1,2,...,N WHERE A(J) IS THE C J-TH COLUMN OF A AND WITH EDGE (A(I),A(J)) IF AND ONLY IF C COLUMNS I AND J HAVE A NON-ZERO IN THE SAME ROW POSITION. C C THE INCIDENCE-DEGREE ORDERING IS DETERMINED RECURSIVELY BY C LETTING LIST(K), K = 1,...,N BE A COLUMN WITH MAXIMAL C INCIDENCE TO THE SUBGRAPH SPANNED BY THE ORDERED COLUMNS. C AMONG ALL THE COLUMNS OF MAXIMAL INCIDENCE, IDO CHOOSES A C COLUMN OF MAXIMAL DEGREE. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE IDO(M,N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST, C MAXCLQ,IWA1,IWA2,IWA3,IWA4) C C WHERE C C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF ROWS OF A. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF COLUMNS OF A. C C INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW C INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW. C THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS OF THE MATRIX A. C C INDCOL IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C IPNTR IS AN INTEGER INPUT ARRAY OF LENGTH M + 1 WHICH C SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL. C THE COLUMN INDICES FOR ROW I ARE C C INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1. C C NOTE THAT IPNTR(M+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS OF THE MATRIX A. C C NDEG IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE DEGREE SEQUENCE. THE DEGREE OF THE J-TH COLUMN C OF A IS NDEG(J). C C LIST IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE INCIDENCE-DEGREE ORDERING OF THE COLUMNS OF A. THE J-TH C COLUMN IN THIS ORDER IS LIST(J). C C MAXCLQ IS AN INTEGER OUTPUT VARIABLE SET TO THE SIZE C OF THE LARGEST CLIQUE FOUND DURING THE ORDERING. C C IWA1,IWA2,IWA3, AND IWA4 ARE INTEGER WORK ARRAYS OF LENGTH N. C C SUBPROGRAMS CALLED C C MINPACK-SUPPLIED ... NUMSRT C C FORTRAN-SUPPLIED ... MAX C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. AUGUST 1984. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER IC,IP,IR,JCOL,JP, * MAXINC,MAXLST,NCOMP,NUMINC,NUMLST,NUMORD,NUMWGT C C SORT THE DEGREE SEQUENCE. C CALL NUMSRT(N,N-1,NDEG,-1,IWA4,IWA2,IWA3) C C INITIALIZATION BLOCK. C C CREATE A DOUBLY-LINKED LIST TO ACCESS THE INCIDENCES OF THE C COLUMNS. THE POINTERS FOR THE LINKED LIST ARE AS FOLLOWS. C C EACH UN-ORDERED COLUMN IC IS IN A LIST (THE INCIDENCE LIST) C OF COLUMNS WITH THE SAME INCIDENCE. C C IWA1(NUMINC) IS THE FIRST COLUMN IN THE NUMINC LIST C UNLESS IWA1(NUMINC) = 0. IN THIS CASE THERE ARE C NO COLUMNS IN THE NUMINC LIST. C C IWA2(IC) IS THE COLUMN BEFORE IC IN THE INCIDENCE LIST C UNLESS IWA2(IC) = 0. IN THIS CASE IC IS THE FIRST C COLUMN IN THIS INCIDENCE LIST. C C IWA3(IC) IS THE COLUMN AFTER IC IN THE INCIDENCE LIST C UNLESS IWA3(IC) = 0. IN THIS CASE IC IS THE LAST C COLUMN IN THIS INCIDENCE LIST. C C IF IC IS AN UN-ORDERED COLUMN, THEN LIST(IC) IS THE C INCIDENCE OF IC TO THE GRAPH INDUCED BY THE ORDERED C COLUMNS. IF JCOL IS AN ORDERED COLUMN, THEN LIST(JCOL) C IS THE INCIDENCE-DEGREE ORDER OF COLUMN JCOL. C MAXINC = 0 DO 10 JP = N, 1, -1 IC = IWA4(JP) IWA1(N-JP) = 0 IWA2(IC) = 0 IWA3(IC) = IWA1(0) IF (IWA1(0) .GT. 0) IWA2(IWA1(0)) = IC IWA1(0) = IC IWA4(JP) = 0 LIST(JP) = 0 10 CONTINUE C C DETERMINE THE MAXIMAL SEARCH LENGTH FOR THE LIST C OF COLUMNS OF MAXIMAL INCIDENCE. C MAXLST = 0 DO 20 IR = 1, M MAXLST = MAXLST + (IPNTR(IR+1) - IPNTR(IR))**2 20 CONTINUE MAXLST = MAXLST/N MAXCLQ = 0 NUMORD = 1 C C BEGINNING OF ITERATION LOOP. C 30 CONTINUE C C CHOOSE A COLUMN JCOL OF MAXIMAL DEGREE AMONG THE C COLUMNS OF MAXIMAL INCIDENCE MAXINC. C 40 CONTINUE JP = IWA1(MAXINC) IF (JP .GT. 0) GO TO 50 MAXINC = MAXINC - 1 GO TO 40 50 CONTINUE NUMWGT = -1 DO 60 NUMLST = 1, MAXLST IF (NDEG(JP) .GT. NUMWGT) THEN NUMWGT = NDEG(JP) JCOL = JP END IF JP = IWA3(JP) IF (JP .LE. 0) GO TO 70 60 CONTINUE 70 CONTINUE LIST(JCOL) = NUMORD C C UPDATE THE SIZE OF THE LARGEST CLIQUE C FOUND DURING THE ORDERING. C IF (MAXINC .EQ. 0) NCOMP = 0 NCOMP = NCOMP + 1 IF (MAXINC+1 .EQ. NCOMP) MAXCLQ = MAX(MAXCLQ,NCOMP) C C TERMINATION TEST. C NUMORD = NUMORD + 1 IF (NUMORD .GT. N) GO TO 100 C C DELETE COLUMN JCOL FROM THE MAXINC LIST. C IF (IWA2(JCOL) .EQ. 0) THEN IWA1(MAXINC) = IWA3(JCOL) ELSE IWA3(IWA2(JCOL)) = IWA3(JCOL) END IF IF (IWA3(JCOL) .GT. 0) IWA2(IWA3(JCOL)) = IWA2(JCOL) C C FIND ALL COLUMNS ADJACENT TO COLUMN JCOL. C IWA4(JCOL) = N C C DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND C TO NON-ZEROES IN THE MATRIX. C DO 90 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1 IR = INDROW(JP) C C FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC) C WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX. C DO 80 IP = IPNTR(IR), IPNTR(IR+1)-1 IC = INDCOL(IP) C C ARRAY IWA4 MARKS COLUMNS WHICH ARE ADJACENT TO C COLUMN JCOL. C IF (IWA4(IC) .LT. NUMORD) THEN IWA4(IC) = NUMORD C C UPDATE THE POINTERS TO THE CURRENT INCIDENCE LISTS. C NUMINC = LIST(IC) LIST(IC) = LIST(IC) + 1 MAXINC = MAX(MAXINC,LIST(IC)) C C DELETE COLUMN IC FROM THE NUMINC LIST. C IF (IWA2(IC) .EQ. 0) THEN IWA1(NUMINC) = IWA3(IC) ELSE IWA3(IWA2(IC)) = IWA3(IC) END IF IF (IWA3(IC) .GT. 0) IWA2(IWA3(IC)) = IWA2(IC) C C ADD COLUMN IC TO THE NUMINC+1 LIST. C IWA2(IC) = 0 IWA3(IC) = IWA1(NUMINC+1) IF (IWA1(NUMINC+1) .GT. 0) IWA2(IWA1(NUMINC+1)) = IC IWA1(NUMINC+1) = IC END IF 80 CONTINUE 90 CONTINUE C C END OF ITERATION LOOP. C GO TO 30 100 CONTINUE C C INVERT THE ARRAY LIST. C DO 110 JCOL = 1, N IWA2(LIST(JCOL)) = JCOL 110 CONTINUE DO 120 JP = 1, N LIST(JP) = IWA2(JP) 120 CONTINUE RETURN C C LAST CARD OF SUBROUTINE IDO. C END SUBROUTINE NUMSRT(N,NMAX,NUM,MODE,INDEX,LAST,NEXT) INTEGER N,NMAX,MODE INTEGER NUM(N),INDEX(N),LAST(0:NMAX),NEXT(N) C **********. C C SUBROUTINE NUMSRT C C GIVEN A SEQUENCE OF INTEGERS, THIS SUBROUTINE GROUPS C TOGETHER THOSE INDICES WITH THE SAME SEQUENCE VALUE C AND, OPTIONALLY, SORTS THE SEQUENCE INTO EITHER C ASCENDING OR DESCENDING ORDER. C C THE SEQUENCE OF INTEGERS IS DEFINED BY THE ARRAY NUM, C AND IT IS ASSUMED THAT THE INTEGERS ARE EACH FROM THE SET C 0,1,...,NMAX. ON OUTPUT THE INDICES K SUCH THAT NUM(K) = L C FOR ANY L = 0,1,...,NMAX CAN BE OBTAINED FROM THE ARRAYS C LAST AND NEXT AS FOLLOWS. C C K = LAST(L) C WHILE (K .NE. 0) K = NEXT(K) C C OPTIONALLY, THE SUBROUTINE PRODUCES AN ARRAY INDEX SO THAT C THE SEQUENCE NUM(INDEX(I)), I = 1,2,...,N IS SORTED. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE NUMSRT(N,NMAX,NUM,MODE,INDEX,LAST,NEXT) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE. C C NMAX IS A POSITIVE INTEGER INPUT VARIABLE. C C NUM IS AN INPUT ARRAY OF LENGTH N WHICH CONTAINS THE C SEQUENCE OF INTEGERS TO BE GROUPED AND SORTED. IT C IS ASSUMED THAT THE INTEGERS ARE EACH FROM THE SET C 0,1,...,NMAX. C C MODE IS AN INTEGER INPUT VARIABLE. THE SEQUENCE NUM IS C SORTED IN ASCENDING ORDER IF MODE IS POSITIVE AND IN C DESCENDING ORDER IF MODE IS NEGATIVE. IF MODE IS 0, C NO SORTING IS DONE. C C INDEX IS AN INTEGER OUTPUT ARRAY OF LENGTH N SET SO C THAT THE SEQUENCE C C NUM(INDEX(I)), I = 1,2,...,N C C IS SORTED ACCORDING TO THE SETTING OF MODE. IF MODE C IS 0, INDEX IS NOT REFERENCED. C C LAST IS AN INTEGER OUTPUT ARRAY OF LENGTH NMAX + 1. THE C INDEX OF NUM FOR THE LAST OCCURRENCE OF L IS LAST(L) C FOR ANY L = 0,1,...,NMAX UNLESS LAST(L) = 0. IN C THIS CASE L DOES NOT APPEAR IN NUM. C C NEXT IS AN INTEGER OUTPUT ARRAY OF LENGTH N. IF C NUM(K) = L, THEN THE INDEX OF NUM FOR THE PREVIOUS C OCCURRENCE OF L IS NEXT(K) FOR ANY L = 0,1,...,NMAX C UNLESS NEXT(K) = 0. IN THIS CASE THERE IS NO PREVIOUS C OCCURRENCE OF L IN NUM. C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER I,J,JINC,JL,JU,K,L C C DETERMINE THE ARRAYS NEXT AND LAST. C DO 10 I = 0, NMAX LAST(I) = 0 10 CONTINUE DO 20 K = 1, N L = NUM(K) NEXT(K) = LAST(L) LAST(L) = K 20 CONTINUE IF (MODE .EQ. 0) RETURN C C STORE THE POINTERS TO THE SORTED ARRAY IN INDEX. C I = 1 IF (MODE .GT. 0) THEN JL = 0 JU = NMAX JINC = 1 ELSE JL = NMAX JU = 0 JINC = -1 END IF DO 50 J = JL, JU, JINC K = LAST(J) 30 CONTINUE IF (K .EQ. 0) GO TO 40 INDEX(I) = K I = I + 1 K = NEXT(K) GO TO 30 40 CONTINUE 50 CONTINUE RETURN C C LAST CARD OF SUBROUTINE NUMSRT. C END SUBROUTINE SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,LIST,NGRP,MAXGRP, * IWA) INTEGER N,MAXGRP INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(*),LIST(N), * NGRP(N),IWA(N) C ********** C C SUBROUTINE SEQ C C GIVEN THE SPARSITY PATTERN OF AN M BY N MATRIX A, THIS C SUBROUTINE DETERMINES A CONSISTENT PARTITION OF THE C COLUMNS OF A BY A SEQUENTIAL ALGORITHM. C C A CONSISTENT PARTITION IS DEFINED IN TERMS OF THE LOOPLESS C GRAPH G WITH VERTICES A(J), J = 1,2,...,N WHERE A(J) IS THE C J-TH COLUMN OF A AND WITH EDGE (A(I),A(J)) IF AND ONLY IF C COLUMNS I AND J HAVE A NON-ZERO IN THE SAME ROW POSITION. C C A PARTITION OF THE COLUMNS OF A INTO GROUPS IS CONSISTENT C IF THE COLUMNS IN ANY GROUP ARE NOT ADJACENT IN THE GRAPH G. C IN GRAPH-THEORY TERMINOLOGY, A CONSISTENT PARTITION OF THE C COLUMNS OF A CORRESPONDS TO A COLORING OF THE GRAPH G. C C THE SUBROUTINE EXAMINES THE COLUMNS IN THE ORDER SPECIFIED C BY THE ARRAY LIST, AND ASSIGNS THE CURRENT COLUMN TO THE C GROUP WITH THE SMALLEST POSSIBLE NUMBER. C C NOTE THAT THE VALUE OF M IS NOT NEEDED BY SEQ AND IS C THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,LIST,NGRP,MAXGRP, C IWA) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF COLUMNS OF A. C C INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW C INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW. C THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS OF THE MATRIX A. C C INDCOL IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C IPNTR IS AN INTEGER INPUT ARRAY OF LENGTH M + 1 WHICH C SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL. C THE COLUMN INDICES FOR ROW I ARE C C INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1. C C NOTE THAT IPNTR(M+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS OF THE MATRIX A. C C LIST IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE ORDER TO BE USED BY THE SEQUENTIAL ALGORITHM. C THE J-TH COLUMN IN THIS ORDER IS LIST(J). C C NGRP IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE PARTITION OF THE COLUMNS OF A. COLUMN JCOL BELONGS C TO GROUP NGRP(JCOL). C C MAXGRP IS AN INTEGER OUTPUT VARIABLE WHICH SPECIFIES THE C NUMBER OF GROUPS IN THE PARTITION OF THE COLUMNS OF A. C C IWA IS AN INTEGER WORK ARRAY OF LENGTH N. C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER IC,IP,IR,J,JCOL,JP C C INITIALIZATION BLOCK. C MAXGRP = 0 DO 10 JP = 1, N NGRP(JP) = N IWA(JP) = 0 10 CONTINUE C C BEGINNING OF ITERATION LOOP. C DO 60 J = 1, N JCOL = LIST(J) C C FIND ALL COLUMNS ADJACENT TO COLUMN JCOL. C C DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND C TO NON-ZEROES IN THE MATRIX. C DO 30 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1 IR = INDROW(JP) C C FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC) C WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX. C DO 20 IP = IPNTR(IR), IPNTR(IR+1)-1 IC = INDCOL(IP) C C ARRAY IWA MARKS THE GROUP NUMBERS OF THE C COLUMNS WHICH ARE ADJACENT TO COLUMN JCOL. C IWA(NGRP(IC)) = J 20 CONTINUE 30 CONTINUE C C ASSIGN THE SMALLEST UN-MARKED GROUP NUMBER TO JCOL. C DO 40 JP = 1, MAXGRP IF (IWA(JP) .NE. J) GO TO 50 40 CONTINUE MAXGRP = MAXGRP + 1 50 CONTINUE NGRP(JCOL) = JP 60 CONTINUE C C END OF ITERATION LOOP. C RETURN C C LAST CARD OF SUBROUTINE SEQ. C END SUBROUTINE SETR(M,N,INDROW,JPNTR,INDCOL,IPNTR,IWA) INTEGER M,N INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(M+1),IWA(M) C ********** C C SUBROUTINE SETR C C GIVEN A COLUMN-ORIENTED DEFINITION OF THE SPARSITY PATTERN C OF AN M BY N MATRIX A, THIS SUBROUTINE DETERMINES A C ROW-ORIENTED DEFINITION OF THE SPARSITY PATTERN OF A. C C ON INPUT THE COLUMN-ORIENTED DEFINITION IS SPECIFIED BY C THE ARRAYS INDROW AND JPNTR. ON OUTPUT THE ROW-ORIENTED C DEFINITION IS SPECIFIED BY THE ARRAYS INDCOL AND IPNTR. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE SETR(M,N,INDROW,JPNTR,INDCOL,IPNTR,IWA) C C WHERE C C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF ROWS OF A. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF COLUMNS OF A. C C INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW C INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW. C THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS OF THE MATRIX A. C C INDCOL IS AN INTEGER OUTPUT ARRAY WHICH CONTAINS THE C COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C IPNTR IS AN INTEGER OUTPUT ARRAY OF LENGTH M + 1 WHICH C SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL. C THE COLUMN INDICES FOR ROW I ARE C C INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1. C C NOTE THAT IPNTR(1) IS SET TO 1 AND THAT IPNTR(M+1)-1 IS C THEN THE NUMBER OF NON-ZERO ELEMENTS OF THE MATRIX A. C C IWA IS AN INTEGER WORK ARRAY OF LENGTH M. C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER IR,JCOL,JP C C STORE IN ARRAY IWA THE COUNTS OF NON-ZEROES IN THE ROWS. C DO 10 IR = 1, M IWA(IR) = 0 10 CONTINUE DO 20 JP = 1, JPNTR(N+1)-1 IWA(INDROW(JP)) = IWA(INDROW(JP)) + 1 20 CONTINUE C C SET POINTERS TO THE START OF THE ROWS IN INDCOL. C IPNTR(1) = 1 DO 30 IR = 1, M IPNTR(IR+1) = IPNTR(IR) + IWA(IR) IWA(IR) = IPNTR(IR) 30 CONTINUE C C FILL INDCOL. C DO 50 JCOL = 1, N DO 40 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1 IR = INDROW(JP) INDCOL(IWA(IR)) = JCOL IWA(IR) = IWA(IR) + 1 40 CONTINUE 50 CONTINUE RETURN C C LAST CARD OF SUBROUTINE SETR. C END SUBROUTINE SLO(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST, * MAXCLQ,IWA1,IWA2,IWA3,IWA4) INTEGER N,MAXCLQ INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(*),NDEG(N), * LIST(N),IWA1(0:N-1),IWA2(N),IWA3(N),IWA4(N) C ********** C C SUBROUTINE SLO C C GIVEN THE SPARSITY PATTERN OF AN M BY N MATRIX A, THIS C SUBROUTINE DETERMINES THE SMALLEST-LAST ORDERING OF THE C COLUMNS OF A. C C THE SMALLEST-LAST ORDERING IS DEFINED FOR THE LOOPLESS C GRAPH G WITH VERTICES A(J), J = 1,2,...,N WHERE A(J) IS THE C J-TH COLUMN OF A AND WITH EDGE (A(I),A(J)) IF AND ONLY IF C COLUMNS I AND J HAVE A NON-ZERO IN THE SAME ROW POSITION. C C THE SMALLEST-LAST ORDERING IS DETERMINED RECURSIVELY BY C LETTING LIST(K), K = N,...,1 BE A COLUMN WITH LEAST DEGREE C IN THE SUBGRAPH SPANNED BY THE UN-ORDERED COLUMNS. C C NOTE THAT THE VALUE OF M IS NOT NEEDED BY SLO AND IS C THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE SLO(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST, C MAXCLQ,IWA1,IWA2,IWA3,IWA4) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF COLUMNS OF A. C C INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW C INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW. C THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS OF THE MATRIX A. C C INDCOL IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C IPNTR IS AN INTEGER INPUT ARRAY OF LENGTH M + 1 WHICH C SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL. C THE COLUMN INDICES FOR ROW I ARE C C INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1. C C NOTE THAT IPNTR(M+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS OF THE MATRIX A. C C NDEG IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE DEGREE SEQUENCE. THE DEGREE OF THE J-TH COLUMN C OF A IS NDEG(J). C C LIST IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE SMALLEST-LAST ORDERING OF THE COLUMNS OF A. THE J-TH C COLUMN IN THIS ORDER IS LIST(J). C C MAXCLQ IS AN INTEGER OUTPUT VARIABLE SET TO THE SIZE C OF THE LARGEST CLIQUE FOUND DURING THE ORDERING. C C IWA1,IWA2,IWA3, AND IWA4 ARE INTEGER WORK ARRAYS OF LENGTH N. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... MIN C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. AUGUST 1984. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER IC,IP,IR,JCOL,JP,MINDEG,NUMDEG,NUMORD C C INITIALIZATION BLOCK. C MINDEG = N DO 10 JP = 1, N IWA1(JP-1) = 0 IWA4(JP) = N LIST(JP) = NDEG(JP) MINDEG = MIN(MINDEG,NDEG(JP)) 10 CONTINUE C C CREATE A DOUBLY-LINKED LIST TO ACCESS THE DEGREES OF THE C COLUMNS. THE POINTERS FOR THE LINKED LIST ARE AS FOLLOWS. C C EACH UN-ORDERED COLUMN IC IS IN A LIST (THE DEGREE LIST) C OF COLUMNS WITH THE SAME DEGREE. C C IWA1(NUMDEG) IS THE FIRST COLUMN IN THE NUMDEG LIST C UNLESS IWA1(NUMDEG) = 0. IN THIS CASE THERE ARE C NO COLUMNS IN THE NUMDEG LIST. C C IWA2(IC) IS THE COLUMN BEFORE IC IN THE DEGREE LIST C UNLESS IWA2(IC) = 0. IN THIS CASE IC IS THE FIRST C COLUMN IN THIS DEGREE LIST. C C IWA3(IC) IS THE COLUMN AFTER IC IN THE DEGREE LIST C UNLESS IWA3(IC) = 0. IN THIS CASE IC IS THE LAST C COLUMN IN THIS DEGREE LIST. C C IF IC IS AN UN-ORDERED COLUMN, THEN LIST(IC) IS THE C DEGREE OF IC IN THE GRAPH INDUCED BY THE UN-ORDERED C COLUMNS. IF JCOL IS AN ORDERED COLUMN, THEN LIST(JCOL) C IS THE SMALLEST-LAST ORDER OF COLUMN JCOL. C DO 20 JP = 1, N NUMDEG = NDEG(JP) IWA2(JP) = 0 IWA3(JP) = IWA1(NUMDEG) IF (IWA1(NUMDEG) .GT. 0) IWA2(IWA1(NUMDEG)) = JP IWA1(NUMDEG) = JP 20 CONTINUE MAXCLQ = 0 NUMORD = N C C BEGINNING OF ITERATION LOOP. C 30 CONTINUE C C CHOOSE A COLUMN JCOL OF MINIMAL DEGREE MINDEG. C 40 CONTINUE JCOL = IWA1(MINDEG) IF (JCOL .GT. 0) GO TO 50 MINDEG = MINDEG + 1 GO TO 40 50 CONTINUE LIST(JCOL) = NUMORD C C MARK THE SIZE OF THE LARGEST CLIQUE C FOUND DURING THE ORDERING. C IF (MINDEG+1 .EQ. NUMORD .AND. MAXCLQ .EQ. 0) * MAXCLQ = NUMORD C C TERMINATION TEST. C NUMORD = NUMORD - 1 IF (NUMORD .EQ. 0) GO TO 80 C C DELETE COLUMN JCOL FROM THE MINDEG LIST. C IWA1(MINDEG) = IWA3(JCOL) IF (IWA3(JCOL) .GT. 0) IWA2(IWA3(JCOL)) = 0 C C FIND ALL COLUMNS ADJACENT TO COLUMN JCOL. C IWA4(JCOL) = 0 C C DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND C TO NON-ZEROES IN THE MATRIX. C DO 70 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1 IR = INDROW(JP) C C FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC) C WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX. C DO 60 IP = IPNTR(IR), IPNTR(IR+1)-1 IC = INDCOL(IP) C C ARRAY IWA4 MARKS COLUMNS WHICH ARE ADJACENT TO C COLUMN JCOL. C IF (IWA4(IC) .GT. NUMORD) THEN IWA4(IC) = NUMORD C C UPDATE THE POINTERS TO THE CURRENT DEGREE LISTS. C NUMDEG = LIST(IC) LIST(IC) = LIST(IC) - 1 MINDEG = MIN(MINDEG,LIST(IC)) C C DELETE COLUMN IC FROM THE NUMDEG LIST. C IF (IWA2(IC) .EQ. 0) THEN IWA1(NUMDEG) = IWA3(IC) ELSE IWA3(IWA2(IC)) = IWA3(IC) END IF IF (IWA3(IC) .GT. 0) IWA2(IWA3(IC)) = IWA2(IC) C C ADD COLUMN IC TO THE NUMDEG-1 LIST. C IWA2(IC) = 0 IWA3(IC) = IWA1(NUMDEG-1) IF (IWA1(NUMDEG-1) .GT. 0) IWA2(IWA1(NUMDEG-1)) = IC IWA1(NUMDEG-1) = IC END IF 60 CONTINUE 70 CONTINUE C C END OF ITERATION LOOP. C GO TO 30 80 CONTINUE C C INVERT THE ARRAY LIST. C DO 90 JCOL = 1, N IWA2(LIST(JCOL)) = JCOL 90 CONTINUE DO 100 JP = 1, N LIST(JP) = IWA2(JP) 100 CONTINUE RETURN C C LAST CARD OF SUBROUTINE SLO. C END SUBROUTINE SRTDAT(N,NNZ,INDROW,INDCOL,JPNTR,IWA) INTEGER N,NNZ INTEGER INDROW(NNZ),INDCOL(NNZ),JPNTR(N+1),IWA(N) C ********** C C SUBROUTINE SRTDAT C C GIVEN THE NON-ZERO ELEMENTS OF AN M BY N MATRIX A IN C ARBITRARY ORDER AS SPECIFIED BY THEIR ROW AND COLUMN C INDICES, THIS SUBROUTINE PERMUTES THESE ELEMENTS SO C THAT THEIR COLUMN INDICES ARE IN NON-DECREASING ORDER. C C ON INPUT IT IS ASSUMED THAT THE ELEMENTS ARE SPECIFIED IN C C INDROW(K),INDCOL(K), K = 1,...,NNZ. C C ON OUTPUT THE ELEMENTS ARE PERMUTED SO THAT INDCOL IS C IN NON-DECREASING ORDER. IN ADDITION, THE ARRAY JPNTR C IS SET SO THAT THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT THE VALUE OF M IS NOT NEEDED BY SRTDAT AND IS C THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE SRTDAT(N,NNZ,INDROW,INDCOL,JPNTR,IWA) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF COLUMNS OF A. C C NNZ IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF NON-ZERO ELEMENTS OF A. C C INDROW IS AN INTEGER ARRAY OF LENGTH NNZ. ON INPUT INDROW C MUST CONTAIN THE ROW INDICES OF THE NON-ZERO ELEMENTS OF A. C ON OUTPUT INDROW IS PERMUTED SO THAT THE CORRESPONDING C COLUMN INDICES OF INDCOL ARE IN NON-DECREASING ORDER. C C INDCOL IS AN INTEGER ARRAY OF LENGTH NNZ. ON INPUT INDCOL C MUST CONTAIN THE COLUMN INDICES OF THE NON-ZERO ELEMENTS C OF A. ON OUTPUT INDCOL IS PERMUTED SO THAT THESE INDICES C ARE IN NON-DECREASING ORDER. C C JPNTR IS AN INTEGER OUTPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE ROW INDICES IN THE OUTPUT C INDROW. THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT JPNTR(1) IS SET TO 1 AND THAT JPNTR(N+1)-1 C IS THEN NNZ. C C IWA IS AN INTEGER WORK ARRAY OF LENGTH N. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... MAX C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER I,J,K,L C C STORE IN ARRAY IWA THE COUNTS OF NON-ZEROES IN THE COLUMNS. C DO 10 J = 1, N IWA(J) = 0 10 CONTINUE DO 20 K = 1, NNZ IWA(INDCOL(K)) = IWA(INDCOL(K)) + 1 20 CONTINUE C C SET POINTERS TO THE START OF THE COLUMNS IN INDROW. C JPNTR(1) = 1 DO 30 J = 1, N JPNTR(J+1) = JPNTR(J) + IWA(J) IWA(J) = JPNTR(J) 30 CONTINUE K = 1 C C BEGIN IN-PLACE SORT. C 40 CONTINUE J = INDCOL(K) IF (K .GE. JPNTR(J)) THEN C C CURRENT ELEMENT IS IN POSITION. NOW EXAMINE THE C NEXT ELEMENT OR THE FIRST UN-SORTED ELEMENT IN C THE J-TH GROUP. C K = MAX(K+1,IWA(J)) ELSE C C CURRENT ELEMENT IS NOT IN POSITION. PLACE ELEMENT C IN POSITION AND MAKE THE DISPLACED ELEMENT THE C CURRENT ELEMENT. C L = IWA(J) IWA(J) = IWA(J) + 1 I = INDROW(K) INDROW(K) = INDROW(L) INDCOL(K) = INDCOL(L) INDROW(L) = I INDCOL(L) = J END IF IF (K .LE. NNZ) GO TO 40 RETURN C C LAST CARD OF SUBROUTINE SRTDAT. C END