/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:55 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "srn2gb.h" #include /* PARAMETER translations */ #define D0INIT 40 #define DINIT 38 #define DTINIT 39 #define DTYPE 16 #define F 10 #define G 28 #define HALF 0.5e0 #define JCN 66 #define JTOL 59 #define MODE 35 #define NEXTV 47 #define NF0 68 #define NF00 81 #define NF1 69 #define NFCALL 6 #define NFCOV 52 #define NFGCAL 7 #define QTR 77 #define RDREQ 57 #define REGD 67 #define RESTOR 9 #define RLIMIT 46 #define RMAT 78 #define TOOBIG 2 #define VNEED 4 #define ZERO 0.e0 /* end of PARAMETER translations */ void /*FUNCTION*/ srn2gb( float b[][2], float d[], float *dr, long iv[], long liv, long lv, long n, long nd, long *n1, long *n2, long p, float r[], float rd[], float v[], float x[]) { #define DR(I_,J_) (*(dr+(I_)*(nd)+(J_))) long int g1, gi, i, iv1, ivmode, jtol1, l, lh, nn, qtr1, rd1, rmat1, y1, yi; float t; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const D = &d[0] - 1; long *const Iv = &iv[0] - 1; float *const R = &r[0] - 1; float *const Rd = &rd[0] - 1; float *const V = &v[0] - 1; float *const X = &x[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. * File: SRN2GB.for Subrs used by the * David Gay & Linda Kaufman nonlinear LS package. * Needed for versions that allow Bounded variables. * SRN2GB is called by SNLAFB, SNLAGB, & SRNSGB. * The other subrs in this file are needed by SRN2GB. * *>> 1998-10-29 SRN2GB Krogh Moved external statement up for mangle. *>> 1996-05-16 SRN2GB Krogh Changes for conversion to C. *>> 1994-11-02 SRN2GB Krogh Changes to use M77CON *>> 1993-03-10 SRN2GB CLL Moved stmt NN = ... to follow IF (IV1 ... *>> 1992-04-27 CLL Removed unreferenced stmt labels. *>> 1991-06-06 CLL Corrected declarations in [D/S]S7DMP *>> 1991-05-21 CLL Changed (1) to (*) in declarations. *>> 1990-06-12 CLL (Revised SRN2GB and SG7ITB from DMG) *>> 1990-03-20 CLL @ JPL *>> 1990-03-14 CLL @ JPL *>> 1990-02-21 CLL @ JPL *>> 1990-02-16 Cll @ JPL *** from netlib, Wed Feb 7 19:41:53 EST 1990 *** *>> 1990-06-12 CLL *>> 1990-04-23 CLL (Recent revision by DMG) *** from netlib, Thu Apr 19 11:58:57 EDT 1990 *** *--S replaces "?": ?RN2GB,?IVSET,?D7TPR,?D7UPD,?G7ITB,?ITSUM,?L7VML *--& ?Q7APL,?Q7RAD,?R7TVM,?V7CPY,?V7SCP,?V2NRM,?RN2G,?A7SST,?F7DHB *--& ?G7QSB,?L7MSB,?L7SQR,?L7TVM,?PARCK,?Q7RSH,?RLDST,?S7DMP,?S7IPR *--& ?S7LUP,?S7LVM,?V2AXY,?V7IPR,?V7VMP,?H2RFG,?H2RFA,?G7QTS,?S7BQN *--& ?D7MLP,?L7MST,?L7ITV,?L7IVM,?R7MDC,?V7SHF,?NLAFB,?NLAGB,?RNSGB * * *** REVISED ITERATION DRIVER FOR NL2SOL WITH SIMPLE BOUNDS *** * */ /* ------------------------ PARAMETER USAGE -------------------------- * * B........ BOUNDS ON X. * D........ SCALE VECTOR. * DR....... DERIVATIVES OF R AT X. * IV....... INTEGER VALUES ARRAY. * LIV...... LENGTH OF IV... LIV MUST BE AT LEAST P + 80. * LV....... LENGTH OF V... LV MUST BE AT LEAST 105 + P*(2*P+16). * N........ TOTAL NUMBER OF RESIDUALS. * ND....... MAX. NO. OF RESIDUALS PASSED ON ONE CALL. * N1....... LOWEST ROW INDEX FOR RESIDUALS SUPPLIED THIS TIME. * N2....... HIGHEST ROW INDEX FOR RESIDUALS SUPPLIED THIS TIME. * P........ NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED. * R........ RESIDUALS. * V........ FLOATING-POINT VALUES ARRAY. * X........ PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS, * OUTPUT = BEST VALUE FOUND). * * *** DISCUSSION *** * * THIS ROUTINE CARRIES OUT ITERATIONS FOR SOLVING NONLINEAR * LEAST SQUARES PROBLEMS. IT IS SIMILAR TO SRN2G, EXCEPT THAT * THIS ROUTINE ENFORCES THE BOUNDS B(1,I) .LE. X(I) .LE. B(2,I), * I = 1(1)P. * * *** GENERAL *** * * CODED BY DAVID M. GAY. * * ++++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++ * * *** EXTERNAL FUNCTIONS AND SUBROUTINES *** * */ /* ------------------------------------------------------------------ * SIVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. * SD7TPR... COMPUTES INNER PRODUCT OF TWO VECTORS. * SD7UPD... UPDATES SCALE VECTOR D. * SG7ITB... PERFORMS BASIC MINIMIZATION ALGORITHM. * SITSUM.... PRINTS ITERATION SUMMARY, INFO ABOUT INITIAL AND FINAL X. * SL7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. * SQ7APL... APPLIES QR TRANSFORMATIONS STORED BY SQ7RAD. * SQ7RAD.... ADDS A NEW BLOCK OF ROWS TO QR DECOMPOSITION. * SR7TVM... MULT. VECTOR BY TRANS. OF UPPER TRIANG. MATRIX FROM QR FACT. * SV7CPY.... COPIES ONE VECTOR TO ANOTHER. * SV7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. * SV2NRM... RETURNS THE 2-NORM OF A VECTOR. * * * *** LOCAL VARIABLES *** * */ /* *** SUBSCRIPTS FOR IV AND V *** * */ /* *** IV SUBSCRIPT VALUES *** * */ /* *** V SUBSCRIPT VALUES *** * */ /* ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ * */ lh = p*(p + 1)/2; if (Iv[1] == 0) sivset( 1, iv, liv, lv, v ); iv1 = Iv[1]; if (iv1 > 2) goto L_10; nn = *n2 - *n1 + 1; Iv[RESTOR] = 0; i = iv1 + 4; if (Iv[TOOBIG] == 0) switch (i) { case 1: goto L_150; case 2: goto L_130; case 3: goto L_150; case 4: goto L_120; case 5: goto L_120; case 6: goto L_150; } if (i != 5) Iv[1] = 2; goto L_40; /* *** FRESH START OR RESTART -- CHECK INPUT INTEGERS *** * */ L_10: if (nd <= 0) goto L_220; if (p <= 0) goto L_220; if (n <= 0) goto L_220; if (iv1 == 14) goto L_30; if (iv1 > 16) goto L_270; if (iv1 < 12) goto L_40; if (iv1 == 12) Iv[1] = 13; if (Iv[1] != 13) goto L_20; Iv[VNEED] += p*(p + 15)/2; L_20: sg7itb( b, d, x, iv, liv, lv, p, p, v, x, x ); if (Iv[1] != 14) goto L_999; /* *** STORAGE ALLOCATION *** * */ Iv[G] = Iv[NEXTV]; Iv[JCN] = Iv[G] + 2*p; Iv[RMAT] = Iv[JCN] + p; Iv[QTR] = Iv[RMAT] + lh; Iv[JTOL] = Iv[QTR] + 2*p; Iv[NEXTV] = Iv[JTOL] + 2*p; /* *** TURN OFF COVARIANCE COMPUTATION *** */ Iv[RDREQ] = 0; if (iv1 == 13) goto L_999; L_30: jtol1 = Iv[JTOL]; if (V[DINIT] >= ZERO) sv7scp( p, d, V[DINIT] ); if (V[DTINIT] > ZERO) sv7scp( p, &V[jtol1], V[DTINIT] ); i = jtol1 + p; if (V[D0INIT] > ZERO) sv7scp( p, &V[i], V[D0INIT] ); Iv[NF0] = 0; Iv[NF1] = 0; if (nd >= n) goto L_40; /* *** SPECIAL CASE HANDLING OF FIRST FUNCTION AND GRADIENT EVALUATION * *** -- ASK FOR BOTH RESIDUAL AND JACOBIAN AT ONCE * */ g1 = Iv[G]; y1 = g1 + p; sg7itb( b, d, &V[g1], iv, liv, lv, p, p, v, x, &V[y1] ); if (Iv[1] != 1) goto L_260; V[F] = ZERO; sv7scp( p, &V[g1], ZERO ); Iv[1] = -1; qtr1 = Iv[QTR]; sv7scp( p, &V[qtr1], ZERO ); Iv[REGD] = 0; rmat1 = Iv[RMAT]; goto L_100; L_40: g1 = Iv[G]; y1 = g1 + p; sg7itb( b, d, &V[g1], iv, liv, lv, p, p, v, x, &V[y1] ); switch (IARITHIF(Iv[1] - 2)) { case -1: goto L_50; case 0: goto L_60; case 1: goto L_260; } L_50: V[F] = ZERO; if (Iv[NF1] == 0) goto L_240; if (Iv[RESTOR] != 2) goto L_240; Iv[NF0] = Iv[NF1]; sv7cpy( n, rd, r ); Iv[REGD] = 0; goto L_240; L_60: sv7scp( p, &V[g1], ZERO ); if (Iv[MODE] > 0) goto L_230; rmat1 = Iv[RMAT]; qtr1 = Iv[QTR]; rd1 = qtr1 + p; sv7scp( p, &V[qtr1], ZERO ); Iv[REGD] = 0; if (nd < n) goto L_90; if (*n1 != 1) goto L_90; if (Iv[MODE] < 0) goto L_100; if (Iv[NF1] == Iv[NFGCAL]) goto L_70; if (Iv[NF0] != Iv[NFGCAL]) goto L_90; sv7cpy( n, r, rd ); goto L_80; L_70: sv7cpy( n, rd, r ); L_80: sq7apl( nd, n, p, dr, rd, 0 ); sr7tvm( nd, p, &V[y1], &V[rd1], dr, rd ); Iv[REGD] = 0; goto L_110; L_90: Iv[1] = -2; if (Iv[MODE] < 0) Iv[1] = -3; L_100: sv7scp( p, &V[y1], ZERO ); L_110: sv7scp( lh, &V[rmat1], ZERO ); goto L_240; /* *** COMPUTE F(X) *** * */ L_120: t = sv2nrm( nn, r ); if (t > V[RLIMIT]) goto L_210; V[F] += HALF*SQ(t); if (*n2 < n) goto L_250; if (*n1 == 1) Iv[NF1] = Iv[NFCALL]; goto L_40; /* *** COMPUTE Y *** * */ L_130: y1 = Iv[G] + p; yi = y1; for (l = 1; l <= p; l++) { V[yi] += sd7tpr( nn, &DR(l - 1,0), r ); yi += 1; } if (*n2 < n) goto L_250; Iv[1] = 2; if (*n1 > 1) Iv[1] = -3; goto L_240; /* *** COMPUTE GRADIENT INFORMATION *** * */ L_150: g1 = Iv[G]; ivmode = Iv[MODE]; if (ivmode < 0) goto L_170; if (ivmode == 0) goto L_180; Iv[1] = 2; /* *** COMPUTE GRADIENT ONLY (FOR USE IN COVARIANCE COMPUTATION) *** * */ gi = g1; for (l = 1; l <= p; l++) { V[gi] += sd7tpr( nn, r, &DR(l - 1,0) ); gi += 1; } goto L_200; /* *** COMPUTE INITIAL FUNCTION VALUE WHEN ND .LT. N *** * */ L_170: if (n <= nd) goto L_180; t = sv2nrm( nn, r ); if (t > V[RLIMIT]) goto L_210; V[F] += HALF*SQ(t); /* *** UPDATE D IF DESIRED *** * */ L_180: if (Iv[DTYPE] > 0) sd7upd( d, dr, iv, liv, lv, n, nd, nn, *n2, p, v ); /* *** COMPUTE RMAT AND QTR *** * */ qtr1 = Iv[QTR]; rmat1 = Iv[RMAT]; sq7rad( nn, nd, p, &V[qtr1], TRUE, &V[rmat1], dr, r ); Iv[NF1] = 0; if (*n1 > 1) goto L_200; if (*n2 < n) goto L_250; /* *** SAVE DIAGONAL OF R FOR COMPUTING Y LATER *** * */ rd1 = qtr1 + p; l = rmat1 - 1; for (i = 1; i <= p; i++) { l += i; V[rd1] = V[l]; rd1 += 1; } L_200: if (*n2 < n) goto L_250; if (ivmode > 0) goto L_40; Iv[NF00] = Iv[NFGCAL]; /* *** COMPUTE G FROM RMAT AND QTR *** * */ sl7vml( p, &V[g1], &V[rmat1], &V[qtr1] ); Iv[1] = 2; if (ivmode == 0) goto L_40; if (n <= nd) goto L_40; /* *** FINISH SPECIAL CASE HANDLING OF FIRST FUNCTION AND GRADIENT * */ y1 = g1 + p; Iv[1] = 1; sg7itb( b, d, &V[g1], iv, liv, lv, p, p, v, x, &V[y1] ); if (Iv[1] != 2) goto L_260; goto L_40; /* *** MISC. DETAILS *** * * *** X IS OUT OF RANGE (OVERSIZE STEP) *** * */ L_210: Iv[TOOBIG] = 1; goto L_40; /* *** BAD N, ND, OR P *** * */ L_220: Iv[1] = 66; goto L_270; /* *** RECORD EXTRA EVALUATIONS FOR FINITE-DIFFERENCE HESSIAN *** * */ L_230: Iv[NFCOV] += 1; Iv[NFCALL] += 1; Iv[NFGCAL] = Iv[NFCALL]; Iv[1] = -1; /* *** RETURN FOR MORE FUNCTION OR GRADIENT INFORMATION *** * */ L_240: *n2 = 0; L_250: *n1 = *n2 + 1; *n2 += nd; if (*n2 > n) *n2 = n; goto L_999; /* *** PRINT SUMMARY OF FINAL ITERATION AND OTHER REQUESTED ITEMS *** * */ L_260: g1 = Iv[G]; L_270: sitsum( d, &V[g1], iv, liv, lv, p, v, x ); L_999: return; /* *** LAST CARD OF SRN2GB FOLLOWS *** */ #undef DR } /* end of function */ /* ================================================================== */ /* PARAMETER translations */ #define CNVCOD 55 #define COSMIN 47 #define COVMAT 26 #define COVREQ 15 #define DGNORM 1 #define DIG 37 #define DSTNRM 2 #define F0 13 #define FDH 74 #define FDIF 11 #define FUZZ 45 #define GTSTEP 4 #define H 56 #define HC 71 #define IERR 75 #define INCFAC 23 #define INITS 25 #define IPIVOT 76 #define IRC 29 #define IVNEED 3 #define KAGQT 33 #define KALM 34 #define LMAT 42 #define LMAX0 35 #define LMAXS 36 #define MODEL 5 #define MXFCAL 17 #define MXITER 18 #define NEGONE (-1.e0) #define NEXTIV 46 #define NGCALL 30 #define NGCOV 53 #define NITER 31 #define NVSAVE 9 #define ONE 1.e0 #define ONEP2 1.2e0 #define P0 48 #define PC 41 #define PERM 58 #define PHMXFC 21 #define PREDUC 7 #define RAD0 9 #define RADFAC 16 #define RADINC 8 #define RADIUS 8 #define RELDX 17 #define S 62 #define SIZE 55 #define STEP 40 #define STGLIM 11 #define STPPAR 5 #define SUSED 64 #define SWITCH_ 12 #define TUNER4 29 #define TUNER5 30 #define VSAVE 60 #define W 65 #define WSCALE 56 #define X0 43 #define XIRC 13 /* end of PARAMETER translations */ void /*FUNCTION*/ sg7itb( float b[][2], float d[], float gg[], long iv[], long liv, long lv, long p, long ps, float v[], float x[], float y[]) { LOGICAL32 havqtr, havrm; long int dig1, g01, h1, hc1, i, i1, ipi, ipiv0, ipiv1, ipiv2, ipn, j, k, l, lmat1, lstgst, p1, p1len, pp1, pp1o2, qtr1, rmat1, s1, step1, stpmod, td1, temp1, temp2, tg1, w1, wlm1, x01; float e, gi, sttsst, t, t1, xi; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const D = &d[0] - 1; float *const Gg = &gg[0] - 1; long *const Iv = &iv[0] - 1; float *const V = &v[0] - 1; float *const X = &x[0] - 1; float *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /*>> 1990-06-12 CLL @ JPL *>> 1990-04-23 CLL (Recent revision by DMG) *** from netlib, Mon Apr 23 20:37:24 EDT 1990 *** *>> 1990-02-20 CLL @ JPL * * *** CARRY OUT NL2SOL-LIKE ITERATIONS FOR GENERALIZED LINEAR *** * *** HAVING SIMPLE BOUNDS ON THE PARAMETERS BEING ESTIMATED. *** * * *** PARAMETER DECLARATIONS *** * */ /* ------------------------- PARAMETER USAGE -------------------------- * * B.... VECTOR OF LOWER AND UPPER BOUNDS ON X. * D.... SCALE VECTOR. * IV... INTEGER VALUE ARRAY. * LIV.. LENGTH OF IV. MUST BE AT LEAST 80. * LH... LENGTH OF H = P*(P+1)/2. * LV... LENGTH OF V. MUST BE AT LEAST P*(3*P + 19)/2 + 7. * GG... GRADIENT AT X (WHEN IV(1) = 2). * HC... GAUSS-NEWTON HESSIAN AT X (WHEN IV(1) = 2). * P.... NUMBER OF PARAMETERS (COMPONENTS IN X). * PS... NUMBER OF NONZERO ROWS AND COLUMNS IN S. * V.... FLOATING-POINT VALUE ARRAY. * X.... PARAMETER VECTOR. * Y.... PART OF YIELD VECTOR (WHEN IV(1)= 2, SCRATCH OTHERWISE). * * *** DISCUSSION *** * * SG7ITB IS SIMILAR TO DG7LIT, EXCEPT FOR THE EXTRA PARAMETER B * -- SG7ITB ENFORCES THE BOUNDS B(1,I) .LE. X(I) .LE. B(2,I), * I = 1(1)P. * SG7ITB PERFORMS NL2SOL-LIKE ITERATIONS FOR A VARIETY OF * REGRESSION PROBLEMS THAT ARE SIMILAR TO NONLINEAR LEAST-SQUARES * IN THAT THE HESSIAN IS THE SUM OF TWO TERMS, A READILY-COMPUTED * FIRST-ORDER TERM AND A SECOND-ORDER TERM. THE CALLER SUPPLIES * THE FIRST-ORDER TERM OF THE HESSIAN IN HC (LOWER TRIANGLE, STORED * COMPACTLY BY ROWS), AND SG7ITB BUILDS AN APPROXIMATION, S, TO THE * SECOND-ORDER TERM. THE CALLER ALSO PROVIDES THE FUNCTION VALUE, * GRADIENT, AND PART OF THE YIELD VECTOR USED IN UPDATING S. * SG7ITB DECIDES DYNAMICALLY WHETHER OR NOT TO USE S WHEN CHOOSING * THE NEXT STEP TO TRY... THE HESSIAN APPROXIMATION USED IS EITHER * HC ALONE (GAUSS-NEWTON MODEL) OR HC + S (AUGMENTED MODEL). * IF PS .LT. P, THEN ROWS AND COLUMNS PS+1...P OF S ARE KEPT * CONSTANT. THEY WILL BE ZERO UNLESS THE CALLER SETS IV(INITS) TO * 1 OR 2 AND SUPPLIES NONZERO VALUES FOR THEM, OR THE CALLER SETS * IV(INITS) TO 3 OR 4 AND THE FINITE-DIFFERENCE INITIAL S THEN * COMPUTED HAS NONZERO VALUES IN THESE ROWS. * * IF IV(INITS) IS 3 OR 4, THEN THE INITIAL S IS COMPUTED BY * FINITE DIFFERENCES. 3 MEANS USE FUNCTION DIFFERENCES, 4 MEANS * USE GRADIENT DIFFERENCES. FINITE DIFFERENCING IS DONE THE SAME * WAY AS IN COMPUTING A COVARIANCE MATRIX (WITH IV(COVREQ) = -1, -2, * 1, OR 2). * * FOR UPDATING S, SG7ITB ASSUMES THAT THE GRADIENT HAS THE FORM * OF A SUM OVER I OF RHO(I,X)*GRAD(R(I,X)), WHERE GRAD DENOTES THE * GRADIENT WITH RESPECT TO X. THE TRUE SECOND-ORDER TERM THEN IS * THE SUM OVER I OF RHO(I,X)*HESSIAN(R(I,X)). IF X = X0 + STEP, * THEN WE WISH TO UPDATE S SO THAT S*STEP IS THE SUM OVER I OF * RHO(I,X)*(GRAD(R(I,X)) - GRAD(R(I,X0))). THE CALLER MUST SUPPLY * PART OF THIS IN Y, NAMELY THE SUM OVER I OF * RHO(I,X)*GRAD(R(I,X0)), WHEN CALLING SG7ITB WITH IV(1) = 2 AND * IV(MODE) = 0 (WHERE MODE = 38). GG THEN CONTANS THE OTHER PART, * SO THAT THE DESIRED YIELD VECTOR IS GG - Y. IF PS .LT. P, THEN * THE ABOVE DISCUSSION APPLIES ONLY TO THE FIRST PS COMPONENTS OF * GRAD(R(I,X)), STEP, AND Y. * * PARAMETERS IV, P, V, AND X ARE THE SAME AS THE CORRESPONDING * ONES TO DN2GB (AND NL2SOL), EXCEPT THAT V CAN BE SHORTER * (SINCE THE PART OF V THAT DN2GB USES FOR STORING D, J, AND R IS * NOT NEEDED). MOREOVER, COMPARED WITH DN2GB (AND NL2SOL), IV(1) * MAY HAVE THE TWO ADDITIONAL OUTPUT VALUES 1 AND 2, WHICH ARE * EXPLAINED BELOW, AS IS THE USE OF IV(TOOBIG) AND IV(NFGCAL). * THE VALUES IV(D), IV(J), AND IV(R), WHICH ARE OUTPUT VALUES FROM * DN2GB (AND DN2FB), ARE NOT REFERENCED BY SG7ITB OR THE * SUBROUTINES IT CALLS. * * WHEN SG7ITB IS FIRST CALLED, I.E., WHEN SG7ITB IS CALLED WITH * IV(1) = 0 OR 12, V(F), GG, AND HC NEED NOT BE INITIALIZED. TO * OBTAIN THESE STARTING VALUES, SG7ITB RETURNS FIRST WITH IV(1) = 1, * THEN WITH IV(1) = 2, WITH IV(MODE) = -1 IN BOTH CASES. ON * SUBSEQUENT RETURNS WITH IV(1) = 2, IV(MODE) = 0 IMPLIES THAT * Y MUST ALSO BE SUPPLIED. (NOTE THAT Y IS USED FOR SCRATCH -- ITS * INPUT CONTENTS ARE LOST. BY CONTRAST, HC IS NEVER CHANGED.) * ONCE CONVERGENCE HAS BEEN OBTAINED, IV(RDREQ) AND IV(COVREQ) MAY * IMPLY THAT A FINITE-DIFFERENCE HESSIAN SHOULD BE COMPUTED FOR USE * IN COMPUTING A COVARIANCE MATRIX. IN THIS CASE SG7ITB WILL MAKE * A NUMBER OF RETURNS WITH IV(1) = 1 OR 2 AND IV(MODE) POSITIVE. * WHEN IV(MODE) IS POSITIVE, Y SHOULD NOT BE CHANGED. * * IV(1) = 1 MEANS THE CALLER SHOULD SET V(F) (I.E., V(10)) TO F(X), THE * FUNCTION VALUE AT X, AND CALL SG7ITB AGAIN, HAVING CHANGED * NONE OF THE OTHER PARAMETERS. AN EXCEPTION OCCURS IF F(X) * CANNOT BE EVALUATED (E.G. IF OVERFLOW WOULD OCCUR), WHICH * MAY HAPPEN BECAUSE OF AN OVERSIZED STEP. IN THIS CASE * THE CALLER SHOULD SET IV(TOOBIG) = IV(2) TO 1, WHICH WILL * CAUSE SG7ITB TO IGNORE V(F) AND TRY A SMALLER STEP. NOTE * THAT THE CURRENT FUNCTION EVALUATION COUNT IS AVAILABLE * IN IV(NFCALL) = IV(6). THIS MAY BE USED TO IDENTIFY * WHICH COPY OF SAVED INFORMATION SHOULD BE USED IN COM- * PUTING GG, HC, AND Y THE NEXT TIME SG7ITB RETURNS WITH * IV(1) = 2. SEE MLPIT FOR AN EXAMPLE OF THIS. * IV(1) = 2 MEANS THE CALLER SHOULD SET GG TO GG(X), THE GRADIENT OF F * AT X. THE CALLER SHOULD ALSO SET HC TO THE GAUSS-NEWTON * HESSIAN AT X. IF IV(MODE) = 0, THEN THE CALLER SHOULD * ALSO COMPUTE THE PART OF THE YIELD VECTOR DESCRIBED ABOVE. * THE CALLER SHOULD THEN CALL SG7ITB AGAIN (WITH IV(1) = 2). * THE CALLER MAY ALSO CHANGE D AT THIS TIME, BUT SHOULD NOT * CHANGE X. NOTE THAT IV(NFGCAL) = IV(7) CONTAINS THE */ /* VALUE THAT IV(NFCALL) HAD DURING THE RETURN WITH * IV(1) = 1 IN WHICH X HAD THE SAME VALUE AS IT NOW HAS. * IV(NFGCAL) IS EITHER IV(NFCALL) OR IV(NFCALL) - 1. MLPIT * IS AN EXAMPLE WHERE THIS INFORMATION IS USED. IF GG OR HC * CANNOT BE EVALUATED AT X, THEN THE CALLER MAY SET * IV(NFGCAL) TO 0, IN WHICH CASE SG7ITB WILL RETURN WITH * IV(1) = 15. * * *** GENERAL *** * * CODED BY DAVID M. GAY. * * (SEE NL2SOL FOR REFERENCES.) * ------------------------------------------------------------------ * References to the function STOPX have been commented out of this * subroutine. If one wishes to be able to terminate this package * gracefully using a keybord "Break" key, one can provide a STOPX * function that returns .true. if the Break key has been pressed * since the last call to STOPX, and otherwise returns .false., and * then uncomment the references to STOPX in this subr. * -- CLL 6/12/90 * Commented out references to RSTRST which was set but not fetched. * -- CLL 6/15/90 * ++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++++ * * *** LOCAL VARIABLES *** * */ /* integer DUMMY, RSTRST */ /* *** CONSTANTS *** * */ /* Fortran intrinsic functions: abs * * *** EXTERNAL FUNCTIONS AND SUBROUTINES *** * * external STOPX * LOGICAL STOPX */ /* ------------------------------------------------------------------ * SA7SST.... ASSESSES CANDIDATE STEP. * SD7TPR... RETURNS INNER PRODUCT OF TWO VECTORS. * SF7DHB... COMPUTE FINITE-DIFFERENCE HESSIAN (FOR INIT. S MATRIX). * SG7QSB... COMPUTES GOLDFELD-QUANDT-TROTTER STEP (AUGMENTED MODEL). * I7COPY.... COPIES ONE INTEGER VECTOR TO ANOTHER. * I7PNVR... INVERTS PERMUTATION ARRAY. * I7SHFT... SHIFTS AN INTEGER VECTOR. * SITSUM.... PRINTS ITERATION SUMMARY AND INFO ON INITIAL AND FINAL X. * SL7MSB... COMPUTES LEVENBERG-MARQUARDT STEP (GAUSS-NEWTON MODEL). * SL7SQR... COMPUTES L * L**T FROM LOWER TRIANGULAR MATRIX L. * SL7TVM... COMPUTES L**T * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. * SL7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. * SPARCK.... CHECK VALIDITY OF IV AND V INPUT COMPONENTS. * SQ7RSH... SHIFTS A QR FACTORIZATION. * SRLDST... COMPUTES V(RELDX) = RELATIVE STEP SIZE. * SS7DMP... MULTIPLIES A SYM. MATRIX FORE AND AFT BY A DIAG. MATRIX. * SS7IPR... APPLIES PERMUTATION TO (LOWER TRIANG. OF) SYM. MATRIX. * SS7LUP... PERFORMS QUASI-NEWTON UPDATE ON COMPACTLY STORED LOWER TRI- * ANGLE OF A SYMMETRIC MATRIX. * SS7LVM... MULTIPLIES COMPACTLY STORED SYM. MATRIX TIMES VECTOR. * STOPX... RETURNS .TRUE. IF THE BREAK KEY HAS BEEN PRESSED. * Call to STOPX commented out. -- CLL 6/12/90 * SV2NRM... RETURNS THE 2-NORM OF A VECTOR. * SV2AXY.... COMPUTES SCALAR TIMES ONE VECTOR PLUS ANOTHER. * SV7CPY.... COPIES ONE VECTOR TO ANOTHER. * SV7IPR... APPLIES A PERMUTATION TO A VECTOR. * SV7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. * SV7VMP... MULTIPLIES (DIVIDES) VECTORS COMPONENTWISE. * * *** SUBSCRIPTS FOR IV AND V *** * */ /* *** IV SUBSCRIPT VALUES *** * * *** (NOTE THAT P0 AND PC ARE STORED IN IV(G0) AND IV(STLSTG) RESP.) * */ /* *** V SUBSCRIPT VALUES *** * */ /* ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ * */ i = Iv[1]; if (i == 1) goto L_50; if (i == 2) goto L_60; if (i < 12) goto L_10; if (i > 13) goto L_10; Iv[VNEED] += p*(3*p + 25)/2 + 7; Iv[IVNEED] += 4*p; L_10: sparck( 1, d, iv, liv, lv, p, v ); i = Iv[1] - 2; if (i > 12) goto L_999; switch (i) { case 1: goto L_360; case 2: goto L_360; case 3: goto L_360; case 4: goto L_360; case 5: goto L_360; case 6: goto L_360; case 7: goto L_240; case 8: goto L_190; case 9: goto L_240; case 10: goto L_20; case 11: goto L_20; case 12: goto L_30; } /* *** STORAGE ALLOCATION *** * */ L_20: pp1o2 = p*(p + 1)/2; Iv[S] = Iv[LMAT] + pp1o2; Iv[X0] = Iv[S] + pp1o2; Iv[STEP] = Iv[X0] + 2*p; Iv[DIG] = Iv[STEP] + 3*p; Iv[W] = Iv[DIG] + 2*p; Iv[H] = Iv[W] + 4*p + 7; Iv[NEXTV] = Iv[H] + pp1o2; Iv[IPIVOT] = Iv[PERM] + 3*p; Iv[NEXTIV] = Iv[IPIVOT] + p; if (Iv[1] != 13) goto L_30; Iv[1] = 14; goto L_999; /* *** INITIALIZATION *** * */ L_30: Iv[NITER] = 0; Iv[NFCALL] = 1; Iv[NGCALL] = 1; Iv[NFGCAL] = 1; Iv[MODE] = -1; Iv[STGLIM] = 2; Iv[TOOBIG] = 0; Iv[CNVCOD] = 0; Iv[COVMAT] = 0; Iv[NFCOV] = 0; Iv[NGCOV] = 0; Iv[RADINC] = 0; Iv[PC] = p; V[RAD0] = ZERO; V[STPPAR] = ZERO; V[RADIUS] = V[LMAX0]/(ONE + V[PHMXFC]); /* *** CHECK CONSISTENCY OF B AND INITIALIZE IP ARRAY *** * */ ipi = Iv[IPIVOT]; for (i = 1; i <= p; i++) { Iv[ipi] = i; ipi += 1; if (b[i - 1][0] > b[i - 1][1]) goto L_680; } /* *** SET INITIAL MODEL AND S MATRIX *** * */ Iv[MODEL] = 1; Iv[1] = 1; if (Iv[S] < 0) goto L_710; if (Iv[INITS] > 1) Iv[MODEL] = 2; s1 = Iv[S]; if (Iv[INITS] == 0 || Iv[INITS] > 2) sv7scp( p*(p + 1)/2, &V[s1], ZERO ); goto L_710; /* *** NEW FUNCTION VALUE *** * */ L_50: if (Iv[MODE] == 0) goto L_360; if (Iv[MODE] > 0) goto L_590; if (Iv[TOOBIG] == 0) goto L_690; Iv[1] = 63; goto L_999; /* *** MAKE SURE GRADIENT COULD BE COMPUTED *** * */ L_60: if (Iv[TOOBIG] == 0) goto L_70; Iv[1] = 65; goto L_999; /* *** NEW GRADIENT *** * */ L_70: Iv[KALM] = -1; Iv[KAGQT] = -1; Iv[FDH] = 0; if (Iv[MODE] > 0) goto L_590; if (Iv[HC] <= 0 && Iv[RMAT] <= 0) goto L_670; /* *** CHOOSE INITIAL PERMUTATION *** * */ ipi = Iv[IPIVOT]; ipn = ipi + p - 1; ipiv2 = Iv[PERM] - 1; k = Iv[PC]; p1 = p; pp1 = p + 1; rmat1 = Iv[RMAT]; havrm = rmat1 > 0; qtr1 = Iv[QTR]; havqtr = qtr1 > 0; /* *** MAKE SURE V(QTR1) IS LEGAL (EVEN WHEN NOT REFERENCED) *** */ w1 = Iv[W]; if (!havqtr) qtr1 = w1 + p; for (i = 1; i <= p; i++) { i1 = Iv[ipn]; ipn -= 1; if (b[i1 - 1][0] >= b[i1 - 1][1]) goto L_80; xi = X[i1]; gi = Gg[i1]; if (xi <= b[i1 - 1][0] && gi > ZERO) goto L_80; if (xi >= b[i1 - 1][1] && gi < ZERO) goto L_80; /* *** DISALLOW CONVERGENCE IF X(I1) HAS JUST BEEN FREED *** */ j = ipiv2 + i1; if (Iv[j] > k) Iv[CNVCOD] = 0; goto L_100; L_80: if (i1 >= p1) goto L_90; i1 = pp1 - i; i7shft( p1, i1, &Iv[ipi] ); if (havrm) sq7rsh( i1, p1, havqtr, &V[qtr1], &V[rmat1], &V[w1] ); L_90: p1 -= 1; L_100: ; } Iv[PC] = p1; /* *** COMPUTE V(DGNORM) (AN OUTPUT VALUE IF WE STOP NOW) *** * */ V[DGNORM] = ZERO; if (p1 <= 0) goto L_110; dig1 = Iv[DIG]; sv7vmp( p, &V[dig1], gg, d, -1 ); sv7ipr( p, &Iv[ipi], &V[dig1] ); V[DGNORM] = sv2nrm( p1, &V[dig1] ); L_110: if (Iv[CNVCOD] != 0) goto L_580; if (Iv[MODE] == 0) goto L_510; Iv[MODE] = 0; V[F0] = V[F]; if (Iv[INITS] <= 2) goto L_170; /* *** ARRANGE FOR FINITE-DIFFERENCE INITIAL S *** * */ Iv[XIRC] = Iv[COVREQ]; Iv[COVREQ] = -1; if (Iv[INITS] > 3) Iv[COVREQ] = 1; Iv[CNVCOD] = 70; goto L_600; /* *** COME TO NEXT STMT AFTER COMPUTING F.D. HESSIAN FOR INIT. S *** * */ L_120: h1 = Iv[FDH]; if (h1 <= 0) goto L_660; Iv[CNVCOD] = 0; Iv[MODE] = 0; Iv[NFCOV] = 0; Iv[NGCOV] = 0; Iv[COVREQ] = Iv[XIRC]; s1 = Iv[S]; pp1o2 = ps*(ps + 1)/2; hc1 = Iv[HC]; if (hc1 <= 0) goto L_130; sv2axy( pp1o2, &V[s1], NEGONE, &V[hc1], &V[h1] ); goto L_140; L_130: rmat1 = Iv[RMAT]; lmat1 = Iv[LMAT]; sl7sqr( p, &V[lmat1], &V[rmat1] ); ipi = Iv[IPIVOT]; ipiv1 = Iv[PERM] + p; i7pnvr( p, &Iv[ipiv1], &Iv[ipi] ); ss7ipr( p, &Iv[ipiv1], &V[lmat1] ); sv2axy( pp1o2, &V[s1], NEGONE, &V[lmat1], &V[h1] ); /* *** ZERO PORTION OF S CORRESPONDING TO FIXED X COMPONENTS *** * */ L_140: for (i = 1; i <= p; i++) { if (b[i - 1][0] < b[i - 1][1]) goto L_160; k = s1 + i*(i - 1)/2; sv7scp( i, &V[k], ZERO ); if (i >= p) goto L_170; k += 2*i - 1; i1 = i + 1; for (j = i1; j <= p; j++) { V[k] = ZERO; k += j; } L_160: ; } L_170: Iv[1] = 2; /* ---------------------------- MAIN LOOP ----------------------------- * * * *** PRINT ITERATION SUMMARY, CHECK ITERATION LIMIT *** * */ L_180: sitsum( d, gg, iv, liv, lv, p, v, x ); L_190: k = Iv[NITER]; if (k < Iv[MXITER]) goto L_200; Iv[1] = 10; goto L_999; L_200: Iv[NITER] = k + 1; /* *** UPDATE RADIUS *** * */ if (k == 0) goto L_220; step1 = Iv[STEP]; for (i = 1; i <= p; i++) { V[step1] *= D[i]; step1 += 1; } step1 = Iv[STEP]; t = V[RADFAC]*sv2nrm( p, &V[step1] ); if (V[RADFAC] < ONE || t > V[RADIUS]) V[RADIUS] = t; /* *** INITIALIZE FOR START OF NEXT ITERATION *** * */ L_220: x01 = Iv[X0]; V[F0] = V[F]; Iv[IRC] = 4; Iv[H] = -labs( Iv[H] ); Iv[SUSED] = Iv[MODEL]; /* *** COPY X TO X0 *** * */ sv7cpy( p, &V[x01], x ); /* *** CHECK STOPX AND FUNCTION EVALUATION LIMIT *** * */ L_230: ; /* if (STOPX(DUMMY)) then * IV(1) = 11 * GO TO 260 * else */ goto L_250; /* endif * * *** COME HERE WHEN RESTARTING AFTER FUNC. EVAL. LIMIT OR STOPX. * */ L_240: if (V[F] >= V[F0]) goto L_250; V[RADFAC] = ONE; k = Iv[NITER]; goto L_200; L_250: if (Iv[NFCALL] < Iv[MXFCAL] + Iv[NFCOV]) goto L_270; Iv[1] = 9; /* 260 continue */ if (V[F] >= V[F0]) goto L_999; /* *** IN CASE OF STOPX OR FUNCTION EVALUATION LIMIT WITH * *** IMPROVED V(F), EVALUATE THE GRADIENT AT X. * */ Iv[CNVCOD] = Iv[1]; goto L_500; /*. . . . . . . . . . . . . COMPUTE CANDIDATE STEP . . . . . . . . . . * */ L_270: step1 = Iv[STEP]; tg1 = Iv[DIG]; td1 = tg1 + p; x01 = Iv[X0]; w1 = Iv[W]; h1 = Iv[H]; p1 = Iv[PC]; ipi = Iv[PERM]; ipiv1 = ipi + p; ipiv2 = ipiv1 + p; ipiv0 = Iv[IPIVOT]; if (Iv[MODEL] == 2) goto L_280; /* *** COMPUTE LEVENBERG-MARQUARDT STEP IF POSSIBLE... * */ rmat1 = Iv[RMAT]; if (rmat1 <= 0) goto L_280; qtr1 = Iv[QTR]; if (qtr1 <= 0) goto L_280; lmat1 = Iv[LMAT]; wlm1 = w1 + p; sl7msb( b, d, gg, Iv[IERR], &Iv[ipiv0], &Iv[ipiv1], &Iv[ipiv2], &Iv[KALM], &V[lmat1], lv, p, &Iv[P0], Iv[PC], &V[qtr1], &V[rmat1], &V[step1], &V[td1], &V[tg1], v, &V[w1], &V[wlm1], x, &V[x01] ); /* *** H IS STORED IN THE END OF W AND HAS JUST BEEN OVERWRITTEN, * *** SO WE MARK IT INVALID... */ Iv[H] = -labs( h1 ); /* *** EVEN IF H WERE STORED ELSEWHERE, IT WOULD BE NECESSARY TO * *** MARK INVALID THE INFORMATION SG7QTS MAY HAVE STORED IN V... */ Iv[KAGQT] = -1; goto L_330; L_280: if (h1 > 0) goto L_320; /* *** SET H TO D**-1 * (HC + T1*S) * D**-1. *** * */ p1len = p1*(p1 + 1)/2; h1 = -h1; Iv[H] = h1; Iv[FDH] = 0; if (p1 <= 0) goto L_320; /* *** MAKE TEMPORARY PERMUTATION ARRAY *** */ i7copy( p, &Iv[ipi], &Iv[ipiv0] ); j = Iv[HC]; if (j > 0) goto L_290; j = h1; rmat1 = Iv[RMAT]; sl7sqr( p1, &V[h1], &V[rmat1] ); goto L_300; L_290: sv7cpy( p*(p + 1)/2, &V[h1], &V[j] ); ss7ipr( p, &Iv[ipi], &V[h1] ); L_300: if (Iv[MODEL] == 1) goto L_310; lmat1 = Iv[LMAT]; s1 = Iv[S]; sv7cpy( p*(p + 1)/2, &V[lmat1], &V[s1] ); ss7ipr( p, &Iv[ipi], &V[lmat1] ); sv2axy( p1len, &V[h1], ONE, &V[lmat1], &V[h1] ); L_310: sv7cpy( p, &V[td1], d ); sv7ipr( p, &Iv[ipi], &V[td1] ); ss7dmp( p1, &V[h1], &V[h1], &V[td1], -1 ); Iv[KAGQT] = -1; /* *** COMPUTE ACTUAL GOLDFELD-QUANDT-TROTTER STEP *** * */ L_320: lmat1 = Iv[LMAT]; sg7qsb( b, d, &V[h1], gg, &Iv[ipi], &Iv[ipiv1], &Iv[ipiv2], &Iv[KAGQT], &V[lmat1], lv, p, &Iv[P0], p1, &V[step1], &V[td1], &V[tg1], v, &V[w1], x, &V[x01] ); if (Iv[KALM] > 0) Iv[KALM] = 0; L_330: if (Iv[IRC] != 6) goto L_340; if (Iv[RESTOR] != 2) goto L_360; /* RSTRST = 2 */ goto L_370; /* *** CHECK WHETHER EVALUATING F(X0 + STEP) LOOKS WORTHWHILE *** * */ L_340: Iv[TOOBIG] = 0; if (V[DSTNRM] <= ZERO) goto L_360; if (Iv[IRC] != 5) goto L_350; if (V[RADFAC] <= ONE) goto L_350; if (V[PREDUC] > ONEP2*V[FDIF]) goto L_350; if (Iv[RESTOR] != 2) goto L_360; /* RSTRST = 0 */ goto L_370; /* *** COMPUTE F(X0 + STEP) *** * */ L_350: x01 = Iv[X0]; step1 = Iv[STEP]; sv2axy( p, x, ONE, &V[step1], &V[x01] ); Iv[NFCALL] += 1; Iv[1] = 1; goto L_710; /*. . . . . . . . . . . . . ASSESS CANDIDATE STEP . . . . . . . . . . . * */ L_360: ; /* RSTRST = 3 */ L_370: x01 = Iv[X0]; V[RELDX] = srldst( p, d, x, &V[x01] ); sa7sst( iv, liv, lv, v ); step1 = Iv[STEP]; lstgst = x01 + p; i = Iv[RESTOR] + 1; switch (i) { case 1: goto L_410; case 2: goto L_380; case 3: goto L_390; case 4: goto L_400; } L_380: sv7cpy( p, x, &V[x01] ); goto L_410; L_390: sv7cpy( p, &V[lstgst], &V[step1] ); goto L_410; L_400: sv7cpy( p, &V[step1], &V[lstgst] ); sv2axy( p, x, ONE, &V[step1], &V[x01] ); V[RELDX] = srldst( p, d, x, &V[x01] ); /* *** IF NECESSARY, SWITCH MODELS *** * */ L_410: if (Iv[SWITCH_] == 0) goto L_420; Iv[H] = -labs( Iv[H] ); Iv[SUSED] += 2; l = Iv[VSAVE]; sv7cpy( NVSAVE, v, &V[l] ); L_420: sv2axy( p, &V[step1], NEGONE, &V[x01], x ); l = Iv[IRC] - 4; stpmod = Iv[MODEL]; if (l > 0) switch (l) { case 1: goto L_440; case 2: goto L_450; case 3: goto L_460; case 4: goto L_460; case 5: goto L_460; case 6: goto L_460; case 7: goto L_460; case 8: goto L_460; case 9: goto L_570; case 10: goto L_510; } /* *** DECIDE WHETHER TO CHANGE MODELS *** * */ e = V[PREDUC] - V[FDIF]; s1 = Iv[S]; ss7lvm( ps, y, &V[s1], &V[step1] ); sttsst = HALF*sd7tpr( ps, &V[step1], y ); if (Iv[MODEL] == 1) sttsst = -sttsst; if (fabsf( e + sttsst )*V[FUZZ] >= fabsf( e )) goto L_430; /* *** SWITCH MODELS *** * */ Iv[MODEL] = 3 - Iv[MODEL]; if (-2 < l) goto L_470; Iv[H] = -labs( Iv[H] ); Iv[SUSED] += 2; l = Iv[VSAVE]; sv7cpy( NVSAVE, &V[l], v ); goto L_230; L_430: if (-3 < l) goto L_470; /* *** RECOMPUTE STEP WITH DIFFERENT RADIUS *** * */ L_440: V[RADIUS] = V[RADFAC]*V[DSTNRM]; goto L_230; /* *** COMPUTE STEP OF LENGTH V(LMAXS) FOR SINGULAR CONVERGENCE TEST * */ L_450: V[RADIUS] = V[LMAXS]; goto L_270; /* *** CONVERGENCE OR FALSE CONVERGENCE *** * */ L_460: Iv[CNVCOD] = l; if (V[F] >= V[F0]) goto L_580; if (Iv[XIRC] == 14) goto L_580; Iv[XIRC] = 14; /*. . . . . . . . . . . . PROCESS ACCEPTABLE STEP . . . . . . . . . . . * */ L_470: Iv[COVMAT] = 0; Iv[REGD] = 0; /* *** SEE WHETHER TO SET V(RADFAC) BY GRADIENT TESTS *** * */ if (Iv[IRC] != 3) goto L_500; step1 = Iv[STEP]; temp1 = step1 + p; temp2 = Iv[X0]; /* *** SET TEMP1 = HESSIAN * STEP FOR USE IN GRADIENT TESTS *** * */ hc1 = Iv[HC]; if (hc1 <= 0) goto L_480; ss7lvm( p, &V[temp1], &V[hc1], &V[step1] ); goto L_490; L_480: rmat1 = Iv[RMAT]; ipiv0 = Iv[IPIVOT]; sv7cpy( p, &V[temp1], &V[step1] ); sv7ipr( p, &Iv[ipiv0], &V[temp1] ); sl7tvm( p, &V[temp1], &V[rmat1], &V[temp1] ); sl7vml( p, &V[temp1], &V[rmat1], &V[temp1] ); ipiv1 = Iv[PERM] + p; i7pnvr( p, &Iv[ipiv1], &Iv[ipiv0] ); sv7ipr( p, &Iv[ipiv1], &V[temp1] ); L_490: if (stpmod == 1) goto L_500; s1 = Iv[S]; ss7lvm( ps, &V[temp2], &V[s1], &V[step1] ); sv2axy( ps, &V[temp1], ONE, &V[temp2], &V[temp1] ); /* *** SAVE OLD GRADIENT AND COMPUTE NEW ONE *** * */ L_500: Iv[NGCALL] += 1; g01 = Iv[W]; sv7cpy( p, &V[g01], gg ); goto L_690; /* *** INITIALIZATIONS -- G0 = GG - G0, ETC. *** * */ L_510: g01 = Iv[W]; sv2axy( p, &V[g01], NEGONE, &V[g01], gg ); step1 = Iv[STEP]; temp1 = step1 + p; temp2 = Iv[X0]; if (Iv[IRC] != 3) goto L_540; /* *** SET V(RADFAC) BY GRADIENT TESTS *** * * *** SET TEMP1 = D**-1 * (HESSIAN * STEP + (GG(X0) - G(X))) *** * */ k = temp1; l = g01; for (i = 1; i <= p; i++) { V[k] = (V[k] - V[l])/D[i]; k += 1; l += 1; } /* *** DO GRADIENT TESTS *** * */ if (sv2nrm( p, &V[temp1] ) <= V[DGNORM]*V[TUNER4]) goto L_530; if (sd7tpr( p, gg, &V[step1] ) >= V[GTSTEP]*V[TUNER5]) goto L_540; L_530: V[RADFAC] = V[INCFAC]; /* *** COMPUTE Y VECTOR NEEDED FOR UPDATING S *** * */ L_540: sv2axy( ps, y, NEGONE, y, gg ); /* *** DETERMINE SIZING FACTOR V(SIZE) *** * * *** SET TEMP1 = S * STEP *** */ s1 = Iv[S]; ss7lvm( ps, &V[temp1], &V[s1], &V[step1] ); t1 = fabsf( sd7tpr( ps, &V[step1], &V[temp1] ) ); t = fabsf( sd7tpr( ps, &V[step1], y ) ); V[SIZE] = ONE; if (t < t1) V[SIZE] = t/t1; /* *** SET G0 TO WCHMTD CHOICE OF FLETCHER AND AL-BAALI *** * */ hc1 = Iv[HC]; if (hc1 <= 0) goto L_550; ss7lvm( ps, &V[g01], &V[hc1], &V[step1] ); goto L_560; L_550: rmat1 = Iv[RMAT]; ipiv0 = Iv[IPIVOT]; sv7cpy( p, &V[g01], &V[step1] ); i = g01 + ps; if (ps < p) sv7scp( p - ps, &V[i], ZERO ); sv7ipr( p, &Iv[ipiv0], &V[g01] ); sl7tvm( p, &V[g01], &V[rmat1], &V[g01] ); sl7vml( p, &V[g01], &V[rmat1], &V[g01] ); ipiv1 = Iv[PERM] + p; i7pnvr( p, &Iv[ipiv1], &Iv[ipiv0] ); sv7ipr( p, &Iv[ipiv1], &V[g01] ); L_560: sv2axy( ps, &V[g01], ONE, y, &V[g01] ); /* *** UPDATE S *** * */ ss7lup( &V[s1], V[COSMIN], ps, V[SIZE], &V[step1], &V[temp1], &V[temp2], &V[g01], &V[WSCALE], y ); Iv[1] = 2; goto L_180; /*. . . . . . . . . . . . . . MISC. DETAILS . . . . . . . . . . . . . . * * *** BAD PARAMETERS TO ASSESS *** * */ L_570: Iv[1] = 64; goto L_999; /* *** CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE *** * */ L_580: if (Iv[RDREQ] == 0) goto L_660; if (Iv[FDH] != 0) goto L_660; if (Iv[CNVCOD] >= 7) goto L_660; if (Iv[REGD] > 0) goto L_660; if (Iv[COVMAT] > 0) goto L_660; if (labs( Iv[COVREQ] ) >= 3) goto L_640; if (Iv[RESTOR] == 0) Iv[RESTOR] = 2; goto L_600; /* *** COMPUTE FINITE-DIFFERENCE HESSIAN FOR COMPUTING COVARIANCE *** * */ L_590: Iv[RESTOR] = 0; L_600: sf7dhb( b, d, gg, &i, iv, liv, lv, p, v, x ); switch (i) { case 1: goto L_610; case 2: goto L_620; case 3: goto L_630; } L_610: Iv[NFCOV] += 1; Iv[NFCALL] += 1; Iv[1] = 1; goto L_710; L_620: Iv[NGCOV] += 1; Iv[NGCALL] += 1; Iv[NFGCAL] = Iv[NFCALL] + Iv[NGCOV]; goto L_690; L_630: if (Iv[CNVCOD] == 70) goto L_120; goto L_660; L_640: h1 = labs( Iv[H] ); Iv[FDH] = h1; Iv[H] = -h1; hc1 = Iv[HC]; if (hc1 <= 0) goto L_650; sv7cpy( p*(p + 1)/2, &V[h1], &V[hc1] ); goto L_660; L_650: rmat1 = Iv[RMAT]; sl7sqr( p, &V[h1], &V[rmat1] ); L_660: Iv[MODE] = 0; Iv[1] = Iv[CNVCOD]; Iv[CNVCOD] = 0; goto L_999; /* *** SPECIAL RETURN FOR MISSING HESSIAN INFORMATION -- BOTH * *** IV(HC) .LE. 0 AND IV(RMAT) .LE. 0 * */ L_670: Iv[1] = 1400; goto L_999; /* *** INCONSISTENT B *** * */ L_680: Iv[1] = 70; goto L_999; /* *** SAVE, THEN INITIALIZE IPIVOT ARRAY BEFORE COMPUTING GG *** * */ L_690: Iv[1] = 2; j = Iv[IPIVOT]; ipi = Iv[PERM]; i7pnvr( p, &Iv[ipi], &Iv[j] ); for (i = 1; i <= p; i++) { Iv[j] = i; j += 1; } /* *** PROJECT X INTO FEASIBLE REGION (PRIOR TO COMPUTING F OR GG) *** * */ L_710: for (i = 1; i <= p; i++) { if (X[i] < b[i - 1][0]) X[i] = b[i - 1][0]; if (X[i] > b[i - 1][1]) X[i] = b[i - 1][1]; } Iv[TOOBIG] = 0; L_999: return; /* *** LAST LINE OF SG7ITB FOLLOWS *** */ } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ sr7tvm( long n, long p, float y[], float d[], float *u, float x[]) { #define U(I_,J_) (*(u+(I_)*(n)+(J_))) long int i, ii, pl, pp1; float t; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const D = &d[0] - 1; float *const X = &x[0] - 1; float *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* *** SET Y TO R*X, WHERE R IS THE UPPER TRIANGULAR MATRIX WHOSE * *** DIAGONAL IS IN D AND WHOSE STRICT UPPER TRIANGLE IS IN U. * * *** X AND Y MAY SHARE STORAGE. * */ /* *** LOCAL VARIABLES *** * */ /* *** BODY *** * */ pl = min( n - 1, p ); pp1 = pl + 1; for (ii = 1; ii <= pl; ii++) { i = pp1 - ii; t = X[i]*D[i]; if (i > 1) t += sd7tpr( i - 1, &U(i - 1,0), x ); Y[i] = t; } /* *** LAST LINE OF SR7TVM FOLLOWS *** */ return; #undef U } /* end of function */ /* PARAMETER translations */ #define DELTA 52 #define DELTA0 44 #define DLTFDC 42 #define FX 53 #define HLIM 0.1e0 #define SAVEI 63 #define TWO 2.e0 #define XMSAVE 51 /* end of PARAMETER translations */ void /*FUNCTION*/ sf7dhb( float b[][2], float d[], float gg[], long *irt, long iv[], long liv, long lv, long p, float v[], float x[]) { LOGICAL32 offsid; long int gsave1, hes, hmi, hpi, hpm, i, k, kind, l, m, mm1, mm1o2, newm1, pp1o2, stp0, stpi, stpm; float del, del0, t, xm, xm1; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const D = &d[0] - 1; float *const Gg = &gg[0] - 1; long *const Iv = &iv[0] - 1; float *const V = &v[0] - 1; float *const X = &x[0] - 1; /* end of OFFSET VECTORS */ /* *** COMPUTE FINITE-DIFFERENCE HESSIAN, STORE IT IN V STARTING * *** AT V(IV(FDH)) = V(-IV(H)). HONOR SIMPLE BOUNDS IN B. * * *** IF IV(COVREQ) .GE. 0 THEN SF7DHB USES GRADIENT DIFFERENCES, * *** OTHERWISE FUNCTION DIFFERENCES. STORAGE IN V IS AS IN DG7LIT. * * IRT VALUES... * 1 = COMPUTE FUNCTION VALUE, I.E., V(F). * 2 = COMPUTE GG. * 3 = DONE. * * * *** PARAMETER DECLARATIONS *** * */ /* *** LOCAL VARIABLES *** * */ /* *** EXTERNAL SUBROUTINES *** * */ /* SV7CPY.... COPY ONE VECTOR TO ANOTHER. * SV7SCP... COPY SCALAR TO ALL COMPONENTS OF A VECTOR. * * *** SUBSCRIPTS FOR IV AND V *** * */ /* ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ * */ *irt = 4; kind = Iv[COVREQ]; m = Iv[MODE]; if (m > 0) goto L_10; hes = labs( Iv[H] ); Iv[H] = -hes; Iv[FDH] = 0; Iv[KAGQT] = -1; V[FX] = V[F]; /* *** SUPPLY ZEROS IN CASE B(1,I) = B(2,I) FOR SOME I *** */ sv7scp( p*(p + 1)/2, &V[hes], ZERO ); L_10: if (m > p) goto L_999; if (kind < 0) goto L_120; /* *** COMPUTE FINITE-DIFFERENCE HESSIAN USING BOTH FUNCTION AND * *** GRADIENT VALUES. * */ gsave1 = Iv[W] + p; if (m > 0) goto L_20; /* *** FIRST CALL ON SF7DHB. SET GSAVE = GG, TAKE FIRST STEP *** */ sv7cpy( p, &V[gsave1], gg ); Iv[SWITCH_] = Iv[NFGCAL]; goto L_80; L_20: del = V[DELTA]; X[m] = V[XMSAVE]; if (Iv[TOOBIG] == 0) goto L_30; /* *** HANDLE OVERSIZE V(DELTA) *** * */ del0 = V[DELTA0]*fmaxf( ONE/D[m], fabsf( X[m] ) ); del *= HALF; if (fabsf( del/del0 ) <= HLIM) goto L_140; L_30: hes = -Iv[H]; /* *** SET GG = (GG - GSAVE)/DEL *** * */ del = ONE/del; for (i = 1; i <= p; i++) { Gg[i] = del*(Gg[i] - V[gsave1]); gsave1 += 1; } /* *** ADD GG AS NEW COL. TO FINITE-DIFF. HESSIAN MATRIX *** * */ k = hes + m*(m - 1)/2; l = k + m - 2; if (m == 1) goto L_60; /* *** SET H(I,M) = 0.5 * (H(I,M) + GG(I)) FOR I = 1 TO M-1 *** * */ mm1 = m - 1; for (i = 1; i <= mm1; i++) { if (b[i - 1][0] < b[i - 1][1]) V[k] = HALF*(V[k] + Gg[i]); k += 1; } /* *** ADD H(I,M) = GG(I) FOR I = M TO P *** * */ L_60: l += 1; for (i = m; i <= p; i++) { if (b[i - 1][0] < b[i - 1][1]) V[l] = Gg[i]; l += i; } L_80: m += 1; Iv[MODE] = m; if (m > p) goto L_340; if (b[m - 1][0] >= b[m - 1][1]) goto L_80; /* *** CHOOSE NEXT FINITE-DIFFERENCE STEP, RETURN TO GET GG THERE *** * */ del = V[DELTA0]*fmaxf( ONE/D[m], fabsf( X[m] ) ); xm = X[m]; if (xm < ZERO) goto L_90; xm1 = xm + del; if (xm1 <= b[m - 1][1]) goto L_110; xm1 = xm - del; if (xm1 >= b[m - 1][0]) goto L_100; goto L_280; L_90: xm1 = xm - del; if (xm1 >= b[m - 1][0]) goto L_100; xm1 = xm + del; if (xm1 <= b[m - 1][1]) goto L_110; goto L_280; L_100: del = -del; L_110: V[XMSAVE] = xm; X[m] = xm1; V[DELTA] = del; *irt = 2; goto L_999; /* *** COMPUTE FINITE-DIFFERENCE HESSIAN USING FUNCTION VALUES ONLY. * */ L_120: stp0 = Iv[W] + p - 1; mm1 = m - 1; mm1o2 = m*mm1/2; hes = -Iv[H]; if (m > 0) goto L_130; /* *** FIRST CALL ON SF7DHB. *** */ Iv[SAVEI] = 0; goto L_240; L_130: if (Iv[TOOBIG] == 0) goto L_150; /* *** PUNT IN THE EVENT OF AN OVERSIZE STEP *** */ L_140: Iv[FDH] = -2; goto L_350; L_150: i = Iv[SAVEI]; if (i > 0) goto L_190; /* *** SAVE F(X + STP(M)*E(M)) IN H(P,M) *** * */ pp1o2 = p*(p - 1)/2; hpm = hes + pp1o2 + mm1; V[hpm] = V[F]; /* *** START COMPUTING ROW M OF THE FINITE-DIFFERENCE HESSIAN H. *** * */ newm1 = 1; goto L_260; L_160: hmi = hes + mm1o2; if (mm1 == 0) goto L_180; hpi = hes + pp1o2; for (i = 1; i <= mm1; i++) { t = ZERO; if (b[i - 1][0] < b[i - 1][1]) t = V[FX] - (V[F] + V[hpi]); V[hmi] = t; hmi += 1; hpi += 1; } L_180: V[hmi] = V[F] - TWO*V[FX]; if (offsid) V[hmi] = V[FX] - TWO*V[F]; /* *** COMPUTE FUNCTION VALUES NEEDED TO COMPLETE ROW M OF H. *** * */ i = 0; goto L_200; L_190: X[i] = V[DELTA]; /* *** FINISH COMPUTING H(M,I) *** * */ stpi = stp0 + i; hmi = hes + mm1o2 + i - 1; stpm = stp0 + m; V[hmi] = (V[hmi] + V[F])/(V[stpi]*V[stpm]); L_200: i += 1; if (i > m) goto L_230; if (b[i - 1][0] < b[i - 1][1]) goto L_210; goto L_200; L_210: Iv[SAVEI] = i; stpi = stp0 + i; V[DELTA] = X[i]; X[i] += V[stpi]; *irt = 1; if (i < m) goto L_999; newm1 = 2; goto L_260; L_220: X[m] = V[XMSAVE] - del; if (offsid) X[m] = V[XMSAVE] + TWO*del; goto L_999; L_230: Iv[SAVEI] = 0; X[m] = V[XMSAVE]; L_240: m += 1; Iv[MODE] = m; if (m > p) goto L_330; if (b[m - 1][0] < b[m - 1][1]) goto L_250; goto L_240; /* *** PREPARE TO COMPUTE ROW M OF THE FINITE-DIFFERENCE HESSIAN H. * *** COMPUTE M-TH STEP SIZE STP(M), THEN RETURN TO OBTAIN * *** F(X + STP(M)*E(M)), WHERE E(M) = M-TH STD. UNIT VECTOR. * */ L_250: V[XMSAVE] = X[m]; newm1 = 3; L_260: xm = V[XMSAVE]; del = V[DLTFDC]*fmaxf( ONE/D[m], fabsf( xm ) ); xm1 = xm + del; offsid = FALSE; if (xm1 <= b[m - 1][1]) goto L_270; offsid = TRUE; xm1 = xm - del; if (xm - TWO*del >= b[m - 1][0]) goto L_300; goto L_280; L_270: if (xm - del >= b[m - 1][0]) goto L_290; offsid = TRUE; if (xm + TWO*del <= b[m - 1][1]) goto L_310; L_280: Iv[FDH] = -2; goto L_350; L_290: if (xm >= ZERO) goto L_310; xm1 = xm - del; L_300: del = -del; L_310: switch (newm1) { case 1: goto L_160; case 2: goto L_220; case 3: goto L_320; } L_320: X[m] = xm1; stpm = stp0 + m; V[stpm] = del; *irt = 1; goto L_999; /* *** HANDLE SPECIAL CASE OF B(1,P) = B(2,P) -- CLEAR SCRATCH VALUES * *** FROM LAST ROW OF FDH... * */ L_330: if (b[p - 1][0] < b[p - 1][1]) goto L_340; i = hes + p*(p - 1)/2; sv7scp( p, &V[i], ZERO ); /* *** RESTORE V(F), ETC. *** * */ L_340: Iv[FDH] = hes; L_350: V[F] = V[FX]; *irt = 3; if (kind < 0) goto L_999; Iv[NFGCAL] = Iv[SWITCH_]; gsave1 = Iv[W] + p; sv7cpy( p, gg, &V[gsave1] ); goto L_999; L_999: return; /* *** LAST LINE OF SF7DHB FOLLOWS *** */ } /* end of function */ /* PARAMETER translations */ #undef ZERO #define ZERO 0.0e0 /* end of PARAMETER translations */ float /*FUNCTION*/ sh2rfg( float a, float b, float *x, float *y, float *z) { float a1, b1, c, sh2rfg_v, t; /* *** DETERMINE X, Y, Z SO I + (1,Z)**T * (X,Y) IS A 2X2 * *** HOUSEHOLDER REFLECTION SENDING (A,B)**T INTO (C,0)**T, * *** WHERE C = -SIGN(A)*SQRT(A**2 + B**2) IS THE VALUE SH2RFG * *** RETURNS. * * ------------------------------------------------------------------ */ /* *** BODY *** * */ if (b != ZERO) goto L_10; *x = ZERO; *y = ZERO; *z = ZERO; sh2rfg_v = a; goto L_999; L_10: t = fabsf( a ) + fabsf( b ); a1 = a/t; b1 = b/t; c = sqrtf( SQ(a1) + SQ(b1) ); if (a1 > ZERO) c = -c; a1 -= c; *z = b1/a1; *x = a1/c; *y = b1/c; sh2rfg_v = t*c; L_999: return( sh2rfg_v ); /* *** LAST LINE OF SH2RFG FOLLOWS *** */ } /* end of function */ void /*FUNCTION*/ sh2rfa( long n, float a[], float b[], float x, float y, float z) { long int i; float t; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const A = &a[0] - 1; float *const B = &b[0] - 1; /* end of OFFSET VECTORS */ /* *** APPLY 2X2 HOUSEHOLDER REFLECTION DETERMINED BY X, Y, Z TO * *** N-VECTORS A, B *** * * ------------------------------------------------------------------ */ for (i = 1; i <= n; i++) { t = A[i]*x + B[i]*y; A[i] += t; B[i] += t*z; } /* *** LAST LINE OF SH2RFA FOLLOWS *** */ return; } /* end of function */ /* PARAMETER translations */ #define DST0 3 #define NREDUC 6 /* end of PARAMETER translations */ void /*FUNCTION*/ sg7qsb( float b[][2], float d[], float dihdi[], float gg[], long ipiv[], long ipiv1[], long ipiv2[], long *ka, float l[], long lv, long p, long *p0, long pc, float *stepx, float td[], float tg[], float v[], float ww[], float x[], float xx0[]) { #define STEPX(I_,J_) (*(stepx+(I_)*(p)+(J_))) long int k, kb, kinit, ns, p1, p10; float ds0, nred, pred, rad; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const D = &d[0] - 1; float *const Dihdi = &dihdi[0] - 1; float *const Gg = &gg[0] - 1; long *const Ipiv = &ipiv[0] - 1; long *const Ipiv1 = &ipiv1[0] - 1; long *const Ipiv2 = &ipiv2[0] - 1; float *const L = &l[0] - 1; float *const Td = &td[0] - 1; float *const Tg = &tg[0] - 1; float *const V = &v[0] - 1; float *const Ww = &ww[0] - 1; float *const X = &x[0] - 1; float *const Xx0 = &xx0[0] - 1; /* end of OFFSET VECTORS */ /* *** COMPUTE HEURISTIC BOUNDED NEWTON STEP *** * */ /* DIMENSION DIHDI(P*(P+1)/2), L(P*(P+1)/2) * */ /* *** LOCAL VARIABLES *** * */ /* *** V SUBSCRIPTS *** * */ /* ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ * */ p1 = pc; if (*ka < 0) goto L_10; nred = V[NREDUC]; ds0 = V[DST0]; goto L_20; L_10: *p0 = 0; *ka = -1; L_20: kinit = -1; if (*p0 == p1) kinit = *ka; sv7cpy( p, x, xx0 ); pred = ZERO; rad = V[RADIUS]; kb = -1; V[DSTNRM] = ZERO; if (p1 > 0) goto L_30; nred = ZERO; ds0 = ZERO; sv7scp( p, stepx, ZERO ); goto L_60; L_30: sv7cpy( p, td, d ); sv7ipr( p, ipiv, td ); sv7vmp( p, tg, gg, d, -1 ); sv7ipr( p, ipiv, tg ); L_40: k = kinit; kinit = -1; V[RADIUS] = rad - V[DSTNRM]; sg7qts( td, tg, dihdi, &k, l, p1, stepx, v, ww ); *p0 = p1; if (*ka >= 0) goto L_50; nred = V[NREDUC]; ds0 = V[DST0]; L_50: *ka = k; V[RADIUS] = rad; p10 = p1; ss7bqn( b, d, &STEPX(1,0), ipiv, ipiv1, ipiv2, &kb, l, lv, &ns, p, &p1, stepx, td, tg, v, ww, x, xx0 ); if (ns > 0) ss7ipr( p10, ipiv1, dihdi ); pred += V[PREDUC]; if (ns != 0) *p0 = 0; if (kb <= 0) goto L_40; L_60: V[DST0] = ds0; V[NREDUC] = nred; V[PREDUC] = pred; V[GTSTEP] = sd7tpr( p, gg, stepx ); /* *** LAST LINE OF SG7QSB FOLLOWS *** */ return; #undef STEPX } /* end of function */ void /*FUNCTION*/ sl7msb( float b[][2], float d[], float gg[], long ierr, long ipiv[], long ipiv1[], long ipiv2[], long *ka, float lmat[], long lv, long p, long *p0, long pc, float qtr[], float rmat[], float *stepx, float td[], float tg[], float v[], float ww[], float wlm[], float x[], float xx0[]) { #define STEPX(I_,J_) (*(stepx+(I_)*(p)+(J_))) long int i, j, k, k0, kb, kinit, l, ns, p1, p10, p11; float ds0, nred, pred, rad; static float one = 1.e0; static float zero = 0.e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const D = &d[0] - 1; float *const Gg = &gg[0] - 1; long *const Ipiv = &ipiv[0] - 1; long *const Ipiv1 = &ipiv1[0] - 1; long *const Ipiv2 = &ipiv2[0] - 1; float *const Lmat = &lmat[0] - 1; float *const Qtr = &qtr[0] - 1; float *const Rmat = &rmat[0] - 1; float *const Td = &td[0] - 1; float *const Tg = &tg[0] - 1; float *const V = &v[0] - 1; float *const Wlm = &wlm[0] - 1; float *const Ww = &ww[0] - 1; float *const X = &x[0] - 1; float *const Xx0 = &xx0[0] - 1; /* end of OFFSET VECTORS */ /* *** COMPUTE HEURISTIC BOUNDED NEWTON STEP *** * */ /* DIMENSION LMAT(P*(P+1)/2), RMAT(P*(P+1)/2), WLM(P*(P+5)/2 + 4) * */ /* *** LOCAL VARIABLES *** * */ /* *** V SUBSCRIPTS *** * */ /* ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ * */ p1 = pc; if (*ka < 0) goto L_10; nred = V[NREDUC]; ds0 = V[DST0]; goto L_20; L_10: *p0 = 0; *ka = -1; L_20: kinit = -1; if (*p0 == p1) kinit = *ka; sv7cpy( p, x, xx0 ); sv7cpy( p, td, d ); /* *** USE STEPX(1,3) AS TEMP. COPY OF QTR *** */ sv7cpy( p, &STEPX(2,0), qtr ); sv7ipr( p, ipiv, td ); pred = zero; rad = V[RADIUS]; kb = -1; V[DSTNRM] = zero; if (p1 > 0) goto L_30; nred = zero; ds0 = zero; sv7scp( p, stepx, zero ); goto L_90; L_30: sv7vmp( p, tg, gg, d, -1 ); sv7ipr( p, ipiv, tg ); p10 = p1; L_40: k = kinit; kinit = -1; V[RADIUS] = rad - V[DSTNRM]; sv7vmp( p1, tg, tg, td, 1 ); for (i = 1; i <= p1; i++) { Ipiv1[i] = i; } k0 = max( 0, k ); sl7mst( td, tg, ierr, ipiv1, &k, p1, &STEPX(2,0), rmat, stepx, v, wlm ); sv7vmp( p1, tg, tg, td, -1 ); *p0 = p1; if (*ka >= 0) goto L_60; nred = V[NREDUC]; ds0 = V[DST0]; L_60: *ka = k; V[RADIUS] = rad; l = p1 + 5; if (k <= k0) sd7mlp( p1, lmat, td, rmat, -1 ); if (k > k0) sd7mlp( p1, lmat, td, &Wlm[l], -1 ); ss7bqn( b, d, &STEPX(1,0), ipiv, ipiv1, ipiv2, &kb, lmat, lv, &ns, p, &p1, stepx, td, tg, v, ww, x, xx0 ); pred += V[PREDUC]; if (ns == 0) goto L_80; *p0 = 0; /* *** UPDATE RMAT AND QTR *** * */ p11 = p1 + 1; l = p10 + p11; for (k = p11; k <= p10; k++) { j = l - k; i = Ipiv2[j]; if (i < j) sq7rsh( i, j, TRUE, qtr, rmat, ww ); } L_80: if (kb > 0) goto L_90; /* *** UPDATE LOCAL COPY OF QTR *** * */ sv7vmp( p10, ww, &STEPX(1,0), td, -1 ); sl7tvm( p10, ww, lmat, ww ); sv2axy( p10, &STEPX(2,0), one, ww, qtr ); goto L_40; L_90: V[DST0] = ds0; V[NREDUC] = nred; V[PREDUC] = pred; V[GTSTEP] = sd7tpr( p, gg, stepx ); /* *** LAST LINE OF SL7MSB FOLLOWS *** */ return; #undef STEPX } /* end of function */ /* PARAMETER translations */ #define PHMNFC 20 /* end of PARAMETER translations */ void /*FUNCTION*/ ss7bqn( float b[][2], float d[], float dst[], long ipiv[], long ipiv1[], long ipiv2[], long *kb, float l[], long lv, long *ns, long p, long *p1, float stepx[], float td[], float tg[], float v[], float ww[], float x[], float xx0[]) { long int i, j, k, p0, p1m1; float alpha, dst0, dst1, dstmax, dstmin, dx, gts, t, t1, ti, xi; static float fudge = 1.0001e0; static float half = 0.5e0; static float meps2 = 0.e0; static float one = 1.0e0; static float two = 2.e0; static float zero = 0.e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const D = &d[0] - 1; float *const Dst = &dst[0] - 1; long *const Ipiv = &ipiv[0] - 1; long *const Ipiv1 = &ipiv1[0] - 1; long *const Ipiv2 = &ipiv2[0] - 1; float *const L = &l[0] - 1; float *const Stepx = &stepx[0] - 1; float *const Td = &td[0] - 1; float *const Tg = &tg[0] - 1; float *const V = &v[0] - 1; float *const Ww = &ww[0] - 1; float *const X = &x[0] - 1; float *const Xx0 = &xx0[0] - 1; /* end of OFFSET VECTORS */ /* *** COMPUTE BOUNDED MODIFIED NEWTON STEP *** * */ /* DIMENSION L(P*(P+1)/2) * */ /* *** LOCAL VARIABLES *** * */ /* *** V SUBSCRIPTS *** * */ /* ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ * */ dstmax = fudge*(one + V[PHMXFC])*V[RADIUS]; dstmin = (one + V[PHMNFC])*V[RADIUS]; dst1 = zero; if (meps2 <= zero) meps2 = two*sr7mdc( 3 ); p0 = *p1; *ns = 0; for (i = 1; i <= p; i++) { Ipiv1[i] = i; Ipiv2[i] = i; } for (i = 1; i <= *p1; i++) { Ww[i] = -Stepx[i]*Td[i]; } alpha = fabsf( V[STPPAR] ); V[PREDUC] = zero; gts = -V[GTSTEP]; if (*kb < 0) sv7scp( p, dst, zero ); *kb = 1; /* *** -WW = D TIMES RESTRICTED NEWTON STEP FROM X + DST/D. * * *** FIND T SUCH THAT X - T*WW IS STILL FEASIBLE. * */ L_30: t = one; k = 0; /* DNSGB (8 of 10) */ for (i = 1; i <= *p1; i++) { j = Ipiv[i]; dx = Ww[i]/D[j]; xi = X[j] - dx; if (xi < b[j - 1][0]) goto L_40; if (xi <= b[j - 1][1]) goto L_60; ti = (X[j] - b[j - 1][1])/dx; k = i; goto L_50; L_40: ti = (X[j] - b[j - 1][0])/dx; k = -i; L_50: if (t <= ti) goto L_60; t = ti; L_60: ; } if (p > *p1) sv7cpy( p - *p1, &Stepx[*p1 + 1], &Dst[*p1 + 1] ); sv2axy( *p1, stepx, -t, ww, dst ); dst0 = dst1; dst1 = sv2nrm( p, stepx ); /* *** CHECK FOR OVERSIZE STEP *** * */ if (dst1 <= dstmax) goto L_80; if (*p1 >= p0) goto L_70; if (dst0 < dstmin) *kb = 0; goto L_110; L_70: k = 0; /* *** UPDATE DST, TG, AND V(PREDUC) *** * */ L_80: V[DSTNRM] = dst1; sv7cpy( *p1, dst, stepx ); t1 = one - t; for (i = 1; i <= *p1; i++) { Tg[i] *= t1; } if (alpha > zero) sv2axy( *p1, tg, t*alpha, ww, tg ); V[PREDUC] += t*((one - half*t)*gts + half*alpha*t*sd7tpr( *p1, ww, ww )); if (k == 0) goto L_110; /* *** PERMUTE L, ETC. IF NECESSARY *** * */ p1m1 = *p1 - 1; j = labs( k ); if (j == *p1) goto L_100; *ns += 1; Ipiv2[*p1] = j; sq7rsh( j, *p1, FALSE, tg, l, ww ); i7shft( *p1, j, ipiv ); i7shft( *p1, j, ipiv1 ); sv7shf( *p1, j, tg ); sv7shf( *p1, j, dst ); L_100: if (k < 0) Ipiv[*p1] = -Ipiv[*p1]; *p1 = p1m1; if (*p1 <= 0) goto L_110; sl7ivm( *p1, ww, l, tg ); gts = sd7tpr( *p1, ww, ww ); sl7itv( *p1, ww, l, ww ); goto L_30; /* *** UNSCALE STEPX *** * */ L_110: for (i = 1; i <= p; i++) { j = labs( Ipiv[i] ); Stepx[j] = Dst[i]/D[j]; } /* *** FUDGE STEPX TO ENSURE THAT IT FORCES APPROPRIATE COMPONENTS * *** TO THEIR BOUNDS *** * */ if (*p1 >= p0) goto L_150; k = *p1 + 1; for (i = k; i <= p0; i++) { j = Ipiv[i]; t = meps2; if (j > 0) goto L_130; t = -t; j = -j; Ipiv[i] = j; L_130: t *= fmaxf( fabsf( X[j] ), fabsf( Xx0[j] ) ); Stepx[j] += t; } L_150: sv2axy( p, x, one, stepx, xx0 ); if (*ns > 0) sv7ipr( p0, ipiv1, td ); /* *** LAST LINE OF SS7BQN FOLLOWS *** */ return; } /* end of function */ void /*FUNCTION*/ sq7rsh( long k, long p, LOGICAL32 havqtr, float qtr1[], float r[], float ww[]) { long int i, i1, j, j1, jm1, jp1, k1, km1, pm1; float a, b, t, wj, x, y, z; static float zero = 0.0e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Qtr1 = &qtr1[0] - 1; float *const R = &r[0] - 1; float *const Ww = &ww[0] - 1; /* end of OFFSET VECTORS */ /* *** PERMUTE COLUMN K OF R TO COLUMN P, MODIFY QTR1 ACCORDINGLY *** * */ /* DIMENSION R(P*(P+1)/2) * */ /* *** LOCAL VARIABLES *** * */ /* ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ * */ if (k >= p) goto L_999; km1 = k - 1; k1 = k*km1/2; sv7cpy( k, ww, &R[k1 + 1] ); wj = Ww[k]; pm1 = p - 1; j1 = k1 + km1; for (j = k; j <= pm1; j++) { jm1 = j - 1; jp1 = j + 1; if (jm1 > 0) sv7cpy( jm1, &R[k1 + 1], &R[j1 + 2] ); j1 += jp1; k1 += j; a = R[j1]; b = R[j1 + 1]; if (b != zero) goto L_10; R[k1] = a; x = zero; z = zero; goto L_40; L_10: R[k1] = sh2rfg( a, b, &x, &y, &z ); if (j == pm1) goto L_30; i1 = j1; for (i = jp1; i <= pm1; i++) { i1 += i; sh2rfa( 1, &R[i1], &R[i1 + 1], x, y, z ); } L_30: if (havqtr) sh2rfa( 1, &Qtr1[j], &Qtr1[jp1], x, y, z ); L_40: t = x*wj; Ww[j] = wj + t; wj = t*z; } Ww[p] = wj; sv7cpy( p, &R[k1 + 1], ww ); L_999: return; } /* end of function */ void /*FUNCTION*/ sv7vmp( long n, float x[], float y[], float z[], long k) { long int i; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const X = &x[0] - 1; float *const Y = &y[0] - 1; float *const Z = &z[0] - 1; /* end of OFFSET VECTORS */ /* *** SET X(I) = Y(I) * Z(I)**K, 1 .LE. I .LE. N (FOR K = 1 OR -1) *** * * ------------------------------------------------------------------ */ if (k >= 0) goto L_20; for (i = 1; i <= n; i++) { X[i] = Y[i]/Z[i]; } goto L_999; L_20: for (i = 1; i <= n; i++) { X[i] = Y[i]*Z[i]; } L_999: return; /* *** LAST CARD OF SV7VMP FOLLOWS *** */ } /* end of function */ void /*FUNCTION*/ sv7ipr( long n, long ip[], float x[]) { long int i, j, k; float t; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Ip = &ip[0] - 1; float *const X = &x[0] - 1; /* end of OFFSET VECTORS */ /* PERMUTE X SO THAT X.OUTPUT(I) = X.INPUT(IP(I)). * IP IS UNCHANGED ON OUTPUT. * * ------------------------------------------------------------------ */ for (i = 1; i <= n; i++) { j = Ip[i]; if (j == i) goto L_30; if (j > 0) goto L_10; Ip[i] = -j; goto L_30; L_10: t = X[i]; k = i; L_20: X[k] = X[j]; k = j; j = Ip[k]; Ip[k] = -j; if (j > i) goto L_20; X[k] = t; L_30: ; } /* *** LAST LINE OF SV7IPR FOLLOWS *** */ return; } /* end of function */ void /*FUNCTION*/ sv7shf( long n, long k, float x[]) { long int i, nm1; float t; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const X = &x[0] - 1; /* end of OFFSET VECTORS */ /* *** SHIFT X(K),...,X(N) LEFT CIRCULARLY ONE POSITION *** * * ------------------------------------------------------------------ */ if (k >= n) goto L_999; nm1 = n - 1; t = X[k]; for (i = k; i <= nm1; i++) { X[i] = X[i + 1]; } X[n] = t; L_999: return; } /* end of function */ void /*FUNCTION*/ ss7ipr( long p, long ip[], float hh[]) { long int i, j, j1, jm, k, k1, kk, km, kmj, l, m; float t; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Hh = &hh[0] - 1; long *const Ip = &ip[0] - 1; /* end of OFFSET VECTORS */ /* APPLY THE PERMUTATION DEFINED BY IP TO THE ROWS AND COLUMNS OF THE * P X P SYMMETRIC MATRIX WHOSE LOWER TRIANGLE IS STORED COMPACTLY IN * HH. THUS H.OUTPUT(I,J) = H.INPUT(IP(I), IP(J)). * * ------------------------------------------------------------------ */ /* *** BODY *** * */ for (i = 1; i <= p; i++) { j = Ip[i]; if (j == i) goto L_90; Ip[i] = labs( j ); if (j < 0) goto L_90; k = i; L_10: j1 = j; k1 = k; if (j <= k) goto L_20; j1 = k; k1 = j; L_20: kmj = k1 - j1; l = j1 - 1; jm = j1*l/2; km = k1*(k1 - 1)/2; if (l <= 0) goto L_40; for (m = 1; m <= l; m++) { jm += 1; t = Hh[jm]; km += 1; Hh[jm] = Hh[km]; Hh[km] = t; } L_40: km += 1; kk = km + kmj; jm += 1; t = Hh[jm]; Hh[jm] = Hh[kk]; Hh[kk] = t; j1 = l; l = kmj - 1; if (l <= 0) goto L_60; for (m = 1; m <= l; m++) { jm += j1 + m; t = Hh[jm]; km += 1; Hh[jm] = Hh[km]; Hh[km] = t; } L_60: if (k1 >= p) goto L_80; l = p - k1; k1 -= 1; km = kk; for (m = 1; m <= l; m++) { km += k1 + m; jm = km - kmj; t = Hh[jm]; Hh[jm] = Hh[km]; Hh[km] = t; } L_80: k = j; j = Ip[k]; Ip[k] = -j; if (j > i) goto L_10; L_90: ; } /* *** LAST LINE OF SS7IPR FOLLOWS *** */ return; } /* end of function */ void /*FUNCTION*/ sd7mlp( long n, float x[], float y[], float z[], long k) { long int i, j, l; float t; static float one = 1.e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const X = &x[0] - 1; float *const Y = &y[0] - 1; float *const Z = &z[0] - 1; /* end of OFFSET VECTORS */ /* *** SET X = DIAG(Y)**K * Z * *** FOR X, Z = LOWER TRIANG. MATRICES STORED COMPACTLY BY ROW * *** K = 1 OR -1. * * ------------------------------------------------------------------ */ l = 1; if (k >= 0) goto L_30; for (i = 1; i <= n; i++) { t = one/Y[i]; for (j = 1; j <= i; j++) { X[l] = t*Z[l]; l += 1; } } goto L_999; L_30: for (i = 1; i <= n; i++) { t = Y[i]; for (j = 1; j <= i; j++) { X[l] = t*Z[l]; l += 1; } } L_999: return; /* *** LAST CARD OF SD7MLP FOLLOWS *** */ } /* end of function */ void /*FUNCTION*/ ss7dmp( long n, float x[], float y[], float z[], long k) { long int i, j, l; float t; static float one = 1.e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const X = &x[0] - 1; float *const Y = &y[0] - 1; float *const Z = &z[0] - 1; /* end of OFFSET VECTORS */ /* *** SET X = DIAG(Z)**K * Y * DIAG(Z)**K * *** FOR X, Y = COMPACTLY STORED LOWER TRIANG. MATRICES * *** K = 1 OR -1. * * ------------------------------------------------------------------ */ l = 1; if (k >= 0) goto L_30; for (i = 1; i <= n; i++) { t = one/Z[i]; for (j = 1; j <= i; j++) { X[l] = t*Y[l]/Z[j]; l += 1; } } goto L_999; L_30: for (i = 1; i <= n; i++) { t = Z[i]; for (j = 1; j <= i; j++) { X[l] = t*Y[l]*Z[j]; l += 1; } } L_999: return; /* *** LAST CARD OF SS7DMP FOLLOWS *** */ } /* end of function */