/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:53 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "drn2g.h" #include #include /* PARAMETER translations */ #define CNVCOD 55 #define COVMAT 26 #define COVREQ 15 #define D0INIT 40 #define DINIT 38 #define DTINIT 39 #define DTYPE 16 #define F 10 #define FDH 74 #define G 28 #define H 56 #define HALF 0.5e0 #define IPIVOT 76 #define IVNEED 3 #define JCN 66 #define JTOL 59 #define LMAT 42 #define MODE 35 #define NEXTIV 46 #define NEXTV 47 #define NF0 68 #define NF00 81 #define NF1 69 #define NFCALL 6 #define NFCOV 52 #define NFGCAL 7 #define NGCALL 30 #define NGCOV 53 #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 Y 48 #define ZERO 0.e0 /* end of PARAMETER translations */ void /*FUNCTION*/ drn2g( double d[], double *dr, long iv[], long liv, long lv, long n, long nd, long *n1, long *n2, long p, double r[], double rd[], double v[], double x[]) { #define DR(I_,J_) (*(dr+(I_)*(nd)+(J_))) long int g1, gi, i, iv1, ivmode, jtol1, k, l, lh, nn, qtr1, rmat1, y1, yi; double t; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const D = &d[0] - 1; long *const Iv = &iv[0] - 1; double *const R = &r[0] - 1; double *const Rd = &rd[0] - 1; double *const V = &v[0] - 1; double *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: DRN2G.for Ten subrs used by the * David Gay & Linda Kaufman nonlinear LS package. * Needed for versions that do not allow Bounded variables. * DRN2G is called by DNLAFU, DNLAGU, & DRNSG. * *>> 2015-07-09 DRN2G Krogh Introduced TP to avoid divide by 0, *>> 2000-01-07 DRN2G Krogh Moved COV1 = IV(COMAT) up in DN2CVP. *>> 1998-10-29 DRN2G Krogh Moved external statement up for mangle. *>> 1996-07-09 DRN2G Krogh Changes for conversion to C. *>> 1995-01-26 DRN2G Krogh Moved formats up for C conversion. *>> 1994-11-02 DRN2G Krogh Changes to use M77CON *>> 1993-03-10 DRN2G CLL Moved stmt NN = ... to follow IF (IV1 ... *>> 1992-04-27 CLL Comment out unreferenced stmt labels. *>> 1992-04-13 CLL Change from Hollerith to '...' syntax in formats. *>> 1990-06-29 CLL Changes to formats in DN2CVP. *>> 1990-06-12 CLL Revised DRN2G & DG7LIT from DMG 4/19/90 *>> 1990-03-30 CLL JPL *>> 1990-03-14 CLL JPL *>> 1990-06-12 CLL *>> 1990-04-23 CLL (Recent revision by DMG) *** from netlib, Thu Apr 19 11:58:57 EDT 1990 *** *--D replaces "?": ?RN2G,?C7VFN,?D7TPR,?D7UPD,?G7LIT,?ITSUM,?IVSET *--& ?L7VML,?N2CVP,?N2LRD,?Q7APL,?Q7RAD,?V2NRM,?V7CPY,?V7SCP, ?N2G *--& ?A7SST,?F7HES,?G7QTS,?L7MST,?L7SQR,?L7SRT,?L7SVN,?L7SVX,?L7TVM *--& ?PARCK,?R7MDC,?RLDST,?S7LUP,?S7LVM,?V2AXY,?L7ITV,?L7IVM,?O7PRD *--& ?L7NVR,?L7TSQ,?V7SCL,?N2RDP,?NLAFU,?NLAGU,?RNSG * * *** REVISED ITERATION DRIVER FOR NL2SOL (VERSION 2.3) *** * */ /* ------------------------- PARAMETER USAGE -------------------------- * * D........ SCALE VECTOR. * DR....... DERIVATIVES OF R AT X. * IV....... INTEGER VALUES ARRAY. * LIV...... LENGTH OF IV... LIV MUST BE AT LEAST P + 82. * 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. * RD....... RD(I) = SQRT(G(I)**T * H(I)**-1 * G(I)) ON OUTPUT WHEN * IV(RDREQ) IS NONZERO. DRN2G SETS IV(REGD) = 1 IF RD * IS SUCCESSFULLY COMPUTED, TO 0 IF NO ATTEMPT WAS MADE * TO COMPUTE IT, AND TO -1 IF H (THE FINITE-DIFFERENCE HESSIAN) * WAS INDEFINITE. IF ND .GE. N, THEN RD IS ALSO USED AS * TEMPORARY STORAGE. * V........ FLOATING-POINT VALUES ARRAY. * X........ PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS, * OUTPUT = BEST VALUE FOUND). * * *** DISCUSSION *** * * NOTE... NL2SOL AND NL2ITR (MENTIONED BELOW) ARE DESCRIBED IN * ACM TRANS. MATH. SOFTWARE, VOL. 7, PP. 369-383 (AN ADAPTIVE * NONLINEAR LEAST-SQUARES ALGORITHM, BY J.E. DENNIS, D.M. GAY, * AND R.E. WELSCH). * * THIS ROUTINE CARRIES OUT ITERATIONS FOR SOLVING NONLINEAR * LEAST SQUARES PROBLEMS. WHEN ND = N, IT IS SIMILAR TO NL2ITR * (WITH J = DR), EXCEPT THAT R(X) AND DR(X) NEED NOT BE INITIALIZED * WHEN DRN2G IS CALLED WITH IV(1) = 0 OR 12. DRN2G ALSO ALLOWS * R AND DR TO BE SUPPLIED ROW-WISE -- JUST SET ND = 1 AND CALL * DRN2G ONCE FOR EACH ROW WHEN PROVIDING RESIDUALS AND JACOBIANS. * ANOTHER NEW FEATURE IS THAT CALLING DRN2G WITH IV(1) = 13 * CAUSES STORAGE ALLOCATION ONLY TO BE PERFORMED -- ON RETURN, SUCH * COMPONENTS AS IV(G) (THE FIRST SUBSCRIPT IN G OF THE GRADIENT) * AND IV(S) (THE FIRST SUBSCRIPT IN V OF THE S LOWER TRIANGLE OF * THE S MATRIX) WILL HAVE BEEN SET (UNLESS LIV OR LV IS TOO SMALL), * AND IV(1) WILL HAVE BEEN SET TO 14. CALLING DRN2G WITH IV(1) = 14 * CAUSES EXECUTION OF THE ALGORITHM TO BEGIN UNDER THE ASSUMPTION * THAT STORAGE HAS BEEN ALLOCATED. * * *** SUPPLYING R AND DR *** * * DRN2G USES IV AND V IN THE SAME WAY AS NL2SOL, WITH A SMALL * NUMBER OF OBVIOUS CHANGES. ONE DIFFERENCE BETWEEN DRN2G AND * NL2ITR IS THAT INITIAL FUNCTION AND GRADIENT INFORMATION NEED NOT * BE SUPPLIED IN THE VERY FIRST CALL ON DRN2G, THE ONE WITH * IV(1) = 0 OR 12. ANOTHER DIFFERENCE IS THAT DRN2G RETURNS WITH * IV(1) = -2 WHEN IT WANTS ANOTHER LOOK AT THE OLD JACOBIAN MATRIX * AND THE CURRENT RESIDUAL -- THE ONE CORRESPONDING TO X AND * IV(NFGCAL). IT THEN RETURNS WITH IV(1) = -3 WHEN IT WANTS TO SEE * BOTH THE NEW RESIDUAL AND THE NEW JACOBIAN MATRIX AT ONCE. NOTE * THAT IV(NFGCAL) = IV(7) CONTAINS THE VALUE THAT IV(NFCALL) = IV(6) * HAD WHEN THE CURRENT RESIDUAL WAS EVALUATED. ALSO NOTE THAT THE * VALUE OF X CORRESPONDING TO THE OLD JACOBIAN MATRIX IS STORED IN * V, STARTING AT V(IV(X0)) = V(IV(43)). * ANOTHER NEW RETURN... DRN2G IV(1) = -1 WHEN IT WANTS BOTH THE * RESIDUAL AND THE JACOBIAN TO BE EVALUATED AT X. * A NEW RESIDUAL VECTOR MUST BE SUPPLIED WHEN DRN2G RETURNS WITH * IV(1) = 1 OR -1. THIS TAKES THE FORM OF VALUES OF R(I,X) PASSED * IN R(I-N1+1), I = N1(1)N2. YOU MAY PASS ALL THESE VALUES AT ONCE * (I.E., N1 = 1 AND N2 = N) OR IN PIECES BY MAKING SEVERAL CALLS ON * DRN2G. EACH TIME DRN2G RETURNS WITH IV(1) = 1, N1 WILL HAVE * BEEN SET TO THE INDEX OF THE NEXT RESIDUAL THAT DRN2G EXPECTS TO * SEE, AND N2 WILL BE SET TO THE INDEX OF THE HIGHEST RESIDUAL THAT * COULD BE GIVEN ON THE NEXT CALL, I.E., N2 = N1 + ND - 1. (THUS * WHEN DRN2G FIRST RETURNS WITH IV(1) = 1 FOR A NEW X, IT WILL * HAVE SET N1 TO 1 AND N2 TO MIN(ND,N).) THE CALLER MAY PROVIDE * FEWER THAN N2-N1+1 RESIDUALS ON THE NEXT CALL BY SETTING N2 TO * A SMALLER VALUE. DRN2G ASSUMES IT HAS SEEN ALL THE RESIDUALS * FOR THE CURRENT X WHEN IT IS CALLED WITH N2 .GE. N. * EXAMPLE... SUPPOSE N = 80 AND THAT R IS TO BE PASSED IN 8 * BLOCKS OF SIZE 10. THE FOLLOWING CODE WOULD DO THE JOB. * * N = 80 * ND = 10 * ... * DO 10 K = 1, 8 * *** COMPUTE R(I,X) FOR I = 10*K-9 TO 10*K *** * *** AND STORE THEM IN R(1),...,R(10) *** * CALL DRN2G(..., R, ...) * 10 CONTINUE * * THE SITUATION IS SIMILAR WHEN GRADIENT INFORMATION IS * REQUIRED, I.E., WHEN DRN2G RETURNS WITH IV(1) = 2, -1, OR -2. * NOTE THAT DRN2G OVERWRITES R, BUT THAT IN THE SPECIAL CASE OF * N1 = 1 AND N2 = N ON PREVIOUS CALLS, DRN2G NEVER RETURNS WITH * IV(1) = -2. IT SHOULD BE CLEAR THAT THE PARTIAL DERIVATIVE OF * R(I,X) WITH RESPECT TO X(L) IS TO BE STORED IN DR(I-N1+1,L), * L = 1(1)P, I = N1(1)N2. IT IS ESSENTIAL THAT R(I) AND DR(I,L) * ALL CORRESPOND TO THE SAME RESIDUALS WHEN IV(1) = -1 OR -2. * * *** COVARIANCE MATRIX *** * * IV(RDREQ) = IV(57) TELLS WHETHER TO COMPUTE A COVARIANCE */ /* MATRIX AND/OR REGRESSION DIAGNOSTICS... 0 MEANS NEITHER, * 1 MEANS COVARIANCE MATRIX ONLY, 2 MEANS REG. DIAGNOSTICS ONLY, * 3 MEANS BOTH. AS WITH NL2SOL, IV(COVREQ) = IV(15) TELLS WHAT * HESSIAN APPROXIMATION TO USE IN THIS COMPUTING. * * *** REGRESSION DIAGNOSTICS *** * * SEE THE COMMENTS IN SUBROUTINE DN2G. * * *** GENERAL *** * * CODED BY DAVID M. GAY. * * ++++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++ * * *** EXTERNAL FUNCTIONS AND SUBROUTINES *** * */ /* ------------------------------------------------------------------ * DC7VFN... FINISHES COVARIANCE COMPUTATION. * DIVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. * DD7TPR... COMPUTES INNER PRODUCT OF TWO VECTORS. * DD7UPD... UPDATES SCALE VECTOR D. * DG7LIT.... PERFORMS BASIC MINIMIZATION ALGORITHM. * DITSUM.... PRINTS ITERATION SUMMARY, INFO ABOUT INITIAL AND FINAL X. * DL7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. * DN2CVP... PRINTS COVARIANCE MATRIX. * DN2LRD... COMPUTES REGRESSION DIAGNOSTICS. * DQ7APL... APPLIES QR TRANSFORMATIONS STORED BY DQ7RAD. * DQ7RAD.... ADDS A NEW BLOCK OF ROWS TO QR DECOMPOSITION. * DV7CPY.... COPIES ONE VECTOR TO ANOTHER. * DV7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. * * *** LOCAL VARIABLES *** * */ /* *** SUBSCRIPTS FOR IV AND V *** * */ /* *** IV SUBSCRIPT VALUES *** * */ /* *** V SUBSCRIPT VALUES *** */ /* ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ * */ lh = p*(p + 1)/2; if (Iv[1] == 0) divset( 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_210; if (p <= 0) goto L_210; if (n <= 0) goto L_210; if (iv1 == 14) goto L_30; if (iv1 > 16) goto L_300; if (iv1 < 12) goto L_40; if (iv1 == 12) Iv[1] = 13; if (Iv[1] != 13) goto L_20; Iv[IVNEED] += p; Iv[VNEED] += p*(p + 13)/2; L_20: dg7lit( d, x, iv, liv, lv, p, p, v, x, x ); if (Iv[1] != 14) goto L_999; /* *** STORAGE ALLOCATION *** * */ Iv[IPIVOT] = Iv[NEXTIV]; Iv[NEXTIV] = Iv[IPIVOT] + p; Iv[Y] = Iv[NEXTV]; Iv[G] = Iv[Y] + p; Iv[JCN] = Iv[G] + p; Iv[RMAT] = Iv[JCN] + p; Iv[QTR] = Iv[RMAT] + lh; Iv[JTOL] = Iv[QTR] + p; Iv[NEXTV] = Iv[JTOL] + 2*p; if (iv1 == 13) goto L_999; L_30: jtol1 = Iv[JTOL]; if (V[DINIT] >= ZERO) dv7scp( p, d, V[DINIT] ); if (V[DTINIT] > ZERO) dv7scp( p, &V[jtol1], V[DTINIT] ); i = jtol1 + p; if (V[D0INIT] > ZERO) dv7scp( 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 = Iv[Y]; dg7lit( d, &V[g1], iv, liv, lv, p, p, v, x, &V[y1] ); if (Iv[1] != 1) goto L_220; V[F] = ZERO; dv7scp( p, &V[g1], ZERO ); Iv[1] = -1; qtr1 = Iv[QTR]; dv7scp( p, &V[qtr1], ZERO ); Iv[REGD] = 0; rmat1 = Iv[RMAT]; goto L_100; L_40: g1 = Iv[G]; y1 = Iv[Y]; dg7lit( 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_220; } L_50: V[F] = ZERO; if (Iv[NF1] == 0) goto L_260; if (Iv[RESTOR] != 2) goto L_260; Iv[NF0] = Iv[NF1]; dv7cpy( n, rd, r ); Iv[REGD] = 0; goto L_260; L_60: dv7scp( p, &V[g1], ZERO ); if (Iv[MODE] > 0) goto L_230; rmat1 = Iv[RMAT]; qtr1 = Iv[QTR]; dv7scp( 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; dv7cpy( n, r, rd ); goto L_80; L_70: dv7cpy( n, rd, r ); L_80: dq7apl( nd, n, p, dr, rd, 0 ); dl7vml( p, &V[y1], &V[rmat1], rd ); goto L_110; L_90: Iv[1] = -2; if (Iv[MODE] < 0) Iv[1] = -1; L_100: dv7scp( p, &V[y1], ZERO ); L_110: dv7scp( lh, &V[rmat1], ZERO ); goto L_260; /* *** COMPUTE F(X) *** * */ L_120: t = dv2nrm( nn, r ); if (t > V[RLIMIT]) goto L_200; V[F] += HALF*SQ(t); if (*n2 < n) goto L_270; if (*n1 == 1) Iv[NF1] = Iv[NFCALL]; goto L_40; /* *** COMPUTE Y *** * */ L_130: y1 = Iv[Y]; yi = y1; for (l = 1; l <= p; l++) { V[yi] += dd7tpr( nn, &DR(l - 1,0), r ); yi += 1; } if (*n2 < n) goto L_270; Iv[1] = 2; if (*n1 > 1) Iv[1] = -3; goto L_260; /* *** COMPUTE GRADIENT INFORMATION *** * */ L_150: if (Iv[MODE] > p) goto L_240; 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] += dd7tpr( nn, r, &DR(l - 1,0) ); gi += 1; } goto L_190; /* *** COMPUTE INITIAL FUNCTION VALUE WHEN ND .LT. N *** * */ L_170: if (n <= nd) goto L_180; t = dv2nrm( nn, r ); if (t > V[RLIMIT]) goto L_200; V[F] += HALF*SQ(t); /* *** UPDATE D IF DESIRED *** * */ L_180: if (Iv[DTYPE] > 0) dd7upd( d, dr, iv, liv, lv, n, nd, nn, *n2, p, v ); /* *** COMPUTE RMAT AND QTR *** * */ qtr1 = Iv[QTR]; rmat1 = Iv[RMAT]; dq7rad( nn, nd, p, &V[qtr1], TRUE, &V[rmat1], dr, r ); Iv[NF1] = 0; L_190: if (*n2 < n) goto L_270; if (ivmode > 0) goto L_40; Iv[NF00] = Iv[NFGCAL]; /* *** COMPUTE G FROM RMAT AND QTR *** * */ dl7vml( 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 = Iv[Y]; Iv[1] = 1; dg7lit( d, &V[g1], iv, liv, lv, p, p, v, x, &V[y1] ); if (Iv[1] != 2) goto L_220; goto L_40; /* *** MISC. DETAILS *** * * *** X IS OUT OF RANGE (OVERSIZE STEP) *** * */ L_200: Iv[TOOBIG] = 1; goto L_40; /* *** BAD N, ND, OR P *** * */ L_210: Iv[1] = 66; goto L_300; /* *** CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE *** * */ L_220: if (Iv[COVMAT] != 0) goto L_290; if (Iv[REGD] != 0) goto L_290; /* *** SEE IF CHOLESKY FACTOR OF HESSIAN IS AVAILABLE *** * */ k = Iv[FDH]; if (k <= 0) goto L_280; if (Iv[RDREQ] <= 0) goto L_290; /* *** COMPUTE REGRESSION DIAGNOSTICS AND DEFAULT COVARIANCE IF * DESIRED *** * */ i = 0; if ((Iv[RDREQ]%4) >= 2) i = 1; if ((Iv[RDREQ]%2) == 1 && labs( Iv[COVREQ] ) <= 1) i += 2; if (i == 0) goto L_250; Iv[MODE] = p + i; Iv[NGCALL] += 1; Iv[NGCOV] += 1; Iv[CNVCOD] = Iv[1]; if (i < 2) goto L_230; l = labs( Iv[H] ); dv7scp( lh, &V[l], ZERO ); L_230: Iv[NFCOV] += 1; Iv[NFCALL] += 1; Iv[NFGCAL] = Iv[NFCALL]; Iv[1] = -1; goto L_260; L_240: l = Iv[LMAT]; dn2lrd( dr, iv, &V[l], lh, liv, lv, nd, nn, p, r, rd, v ); if (*n2 < n) goto L_270; if (*n1 > 1) goto L_250; /* *** ENSURE WE CAN RESTART -- AND MAKE RETURN STATE OF DR * *** INDEPENDENT OF WHETHER REGRESSION DIAGNOSTICS ARE COMPUTED. * *** USE STEP VECTOR (ALLOCATED BY DG7LIT) FOR SCRATCH. * */ rmat1 = Iv[RMAT]; dv7scp( lh, &V[rmat1], ZERO ); dq7rad( nn, nd, p, r, FALSE, &V[rmat1], dr, r ); Iv[NF1] = 0; /* *** FINISH COMPUTING COVARIANCE *** * */ L_250: l = Iv[LMAT]; dc7vfn( iv, &V[l], lh, liv, lv, n, p, v ); goto L_290; /* *** RETURN FOR MORE FUNCTION OR GRADIENT INFORMATION *** * */ L_260: *n2 = 0; L_270: *n1 = *n2 + 1; *n2 += nd; if (*n2 > n) *n2 = n; goto L_999; /* *** COME HERE FOR INDEFINITE FINITE-DIFFERENCE HESSIAN *** * */ L_280: Iv[COVMAT] = k; Iv[REGD] = k; /* *** PRINT SUMMARY OF FINAL ITERATION AND OTHER REQUESTED ITEMS *** * */ L_290: g1 = Iv[G]; L_300: ditsum( d, &V[g1], iv, liv, lv, p, v, x ); if (Iv[1] <= 6 && Iv[RDREQ] > 0) dn2cvp( iv, liv, lv, p, v ); L_999: return; /* *** LAST LINE OF DRN2G FOLLOWS *** */ #undef DR } /* end of function */ /* ================================================================== */ /* PARAMETER translations */ #define COSMIN 47 #define DGNORM 1 #define DIG 37 #define DSTNRM 2 #define F0 13 #define FDIF 11 #define FUZZ 45 #define GTSTEP 4 #define HC 71 #define IERR 75 #define INCFAC 23 #define INITS 25 #define IRC 29 #define KAGQT 33 #define KALM 34 #define LMAX0 35 #define LMAXS 36 #define MODEL 5 #define MXFCAL 17 #define MXITER 18 #define NEGONE (-1.e0) #define NITER 31 #define NVSAVE 9 #define ONE 1.e0 #define ONEP2 1.2e0 #define PHMXFC 21 #define PREDUC 7 #define RAD0 9 #define RADFAC 16 #define RADINC 8 #define RADIUS 8 #define RCOND 53 #define RELDX 17 #define S 62 #define SIZE 55 #define STEP 40 #define STGLIM 11 #define STLSTG 41 #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*/ dg7lit( double d[], double gg[], long iv[], long liv, long lv, long p, long ps, double v[], double x[], double yy[]) { long int dig1, g01, h1, hc1, i, ipiv1, j, k, l, lmat1, lstgst, pp1o2, qtr1, rmat1, rstrst, s1, step1, stpmod, temp1, temp2, w1, x01; double e, sttsst, t, t1, tp; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const D = &d[0] - 1; double *const Gg = &gg[0] - 1; long *const Iv = &iv[0] - 1; double *const V = &v[0] - 1; double *const X = &x[0] - 1; double *const Yy = &yy[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 *** * *** REGRESSION PROBLEMS (AND OTHERS OF SIMILAR STRUCTURE) *** * * *** PARAMETER DECLARATIONS *** * */ /* ------------------------- PARAMETER USAGE -------------------------- * * D.... SCALE VECTOR. * IV... INTEGER VALUE ARRAY. * LIV.. LENGTH OF IV. MUST BE AT LEAST 82. * 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). * P.... NUMBER OF PARAMETERS (COMPONENTS IN X). * PS... NUMBER OF NONZERO ROWS AND COLUMNS IN S. * V.... FLOATING-POINT VALUE ARRAY. * X.... PARAMETER VECTOR. * YY... PART OF YIELD VECTOR (WHEN IV(1)= 2, SCRATCH OTHERWISE). * * *** DISCUSSION *** * * DG7LIT 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 IN V, STARTING AT IV(HC)), AND DG7LIT 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. DG7LIT 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,DG7LIT 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 YY, NAMELY THE SUM OVER I OF * RHO(I,X)*GRAD(R(I,X0)), WHEN CALLING DG7LIT 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 - YY. IF PS .LT. P, THEN * THE ABOVE DISCUSSION APPLIES ONLY TO THE FIRST PS COMPONENTS OF * GRAD(R(I,X)), STEP, AND YY. * * PARAMETERS IV, P, V, AND X ARE THE SAME AS THE CORRESPONDING * ONES TO NL2SOL (WHICH SEE), EXCEPT THAT V CAN BE SHORTER * (SINCE THE PART OF V THAT NL2SOL USES FOR STORING D, J, AND R IS * NOT NEEDED). MOREOVER, COMPARED WITH 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 NL2SOL (AND * NL2SNO), ARE NOT REFERENCED BY DG7LIT OR THE SUBROUTINES IT CALLS. * * WHEN DG7LIT IS FIRST CALLED, I.E., WHEN DG7LIT IS CALLED WITH * IV(1) = 0 OR 12, V(F), GG, AND HC NEED NOT BE INITIALIZED. TO * OBTAIN THESE STARTING VALUES,DG7LIT 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 * YY MUST ALSO BE SUPPLIED. (NOTE THAT YY 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 DG7LIT WILL MAKE A * NUMBER OF RETURNS WITH IV(1) = 1 OR 2 AND IV(MODE) POSITIVE. * WHEN IV(MODE) IS POSITIVE, YY 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 DG7LIT 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 DG7LIT 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 YY THE NEXT TIME DG7LIT RETURNS WITH * IV(1) = 2. SEE MLPIT FOR AN EXAMPLE OF THIS. * IV(1) = 2 MEANS THE CALLER SHOULD SET GG TO G(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 DG7LIT 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(TOOBIG) TO 1, IN WHICH CASE DG7LIT WILL RETURN WITH * IV(1) = 15. * * *** GENERAL *** * * CODED BY DAVID M. GAY. * THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH * SUPPORTED IN PART BY D.O.E. GRANT EX-76-A-01-2295 TO MIT/CCREMS. * * (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 * ++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++++ * * *** LOCAL VARIABLES *** * * integer DUMMY */ /* *** CONSTANTS *** * */ /* *** EXTERNAL FUNCTIONS AND SUBROUTINES *** * * external STOPX * LOGICAL STOPX */ /* ------------------------------------------------------------------ * DA7SST.... ASSESSES CANDIDATE STEP. * DD7TPR... RETURNS INNER PRODUCT OF TWO VECTORS. * DF7HES.... COMPUTE FINITE-DIFFERENCE HESSIAN (FOR COVARIANCE). * DG7QTS.... COMPUTES GOLDFELD-QUANDT-TROTTER STEP (AUGMENTED MODEL). * DITSUM.... PRINTS ITERATION SUMMARY AND INFO ON INITIAL AND FINAL X. * DL7MST... COMPUTES LEVENBERG-MARQUARDT STEP (GAUSS-NEWTON MODEL). * DL7SRT.... COMPUTES CHOLESKY FACTOR OF (LOWER TRIANG. OF) SYM. MATRIX. * DL7SQR... COMPUTES L * L**T FROM LOWER TRIANGULAR MATRIX L. * DL7TVM... COMPUTES L**T * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. * DL7SVX... ESTIMATES LARGEST SING. VALUE OF LOWER TRIANG. MATRIX. * DL7SVN... ESTIMATES SMALLEST SING. VALUE OF LOWER TRIANG. MATRIX. * DL7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. * DPARCK.... CHECK VALIDITY OF IV AND V INPUT COMPONENTS. * DRLDST... COMPUTES V(RELDX) = RELATIVE STEP SIZE. * DR7MDC... RETURNS MACHINE-DEPENDENT CONSTANTS. * DS7LUP... PERFORMS QUASI-NEWTON UPDATE ON COMPACTLY STORED LOWER TRI- * ANGLE OF A SYMMETRIC MATRIX. * STOPX.... RETURNS .TRUE. IF THE BREAK KEY HAS BEEN PRESSED. * Call to STOPX commented out. -- CLL 6/12/90 * DV2AXY.... COMPUTES SCALAR TIMES ONE VECTOR PLUS ANOTHER. * DV7CPY.... COPIES ONE VECTOR TO ANOTHER. * DV7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. * DV2NRM... RETURNS THE 2-NORM OF A VECTOR. * * *** SUBSCRIPTS FOR IV AND V *** * */ /* *** IV SUBSCRIPT VALUES *** * */ /* *** V SUBSCRIPT VALUES *** * */ /* ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ * */ i = Iv[1]; if (i == 1) goto L_40; if (i == 2) goto L_50; if (i == 12 || i == 13) Iv[VNEED] += p*(3*p + 19)/2 + 7; dparck( 1, d, iv, liv, lv, p, v ); i = Iv[1] - 2; if (i > 12) goto L_999; switch (i) { case 1: goto L_290; case 2: goto L_290; case 3: goto L_290; case 4: goto L_290; case 5: goto L_290; case 6: goto L_290; case 7: goto L_170; case 8: goto L_120; case 9: goto L_170; case 10: goto L_10; case 11: goto L_10; case 12: goto L_20; } /* *** STORAGE ALLOCATION *** * */ L_10: pp1o2 = p*(p + 1)/2; Iv[S] = Iv[LMAT] + pp1o2; Iv[X0] = Iv[S] + pp1o2; Iv[STEP] = Iv[X0] + p; Iv[STLSTG] = Iv[STEP] + p; Iv[DIG] = Iv[STLSTG] + p; Iv[W] = Iv[DIG] + p; Iv[H] = Iv[W] + 4*p + 7; Iv[NEXTV] = Iv[H] + pp1o2; if (Iv[1] != 13) goto L_20; Iv[1] = 14; goto L_999; /* *** INITIALIZATION *** * */ L_20: 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[RESTOR] = 0; Iv[FDH] = 0; V[RAD0] = ZERO; V[STPPAR] = ZERO; V[RADIUS] = V[LMAX0]/(ONE + V[PHMXFC]); /* *** SET INITIAL MODEL AND S MATRIX *** * */ Iv[MODEL] = 1; if (Iv[S] < 0) goto L_999; if (Iv[INITS] > 1) Iv[MODEL] = 2; s1 = Iv[S]; if (Iv[INITS] == 0 || Iv[INITS] > 2) dv7scp( p*(p + 1)/2, &V[s1], ZERO ); Iv[1] = 1; j = Iv[IPIVOT]; if (j <= 0) goto L_999; for (i = 1; i <= p; i++) { Iv[j] = i; j += 1; } goto L_999; /* *** NEW FUNCTION VALUE *** * */ L_40: if (Iv[MODE] == 0) goto L_290; if (Iv[MODE] > 0) goto L_520; Iv[1] = 2; if (Iv[TOOBIG] == 0) goto L_999; Iv[1] = 63; goto L_999; /* *** NEW GRADIENT *** * */ L_50: Iv[KALM] = -1; Iv[KAGQT] = -1; Iv[FDH] = 0; if (Iv[MODE] > 0) goto L_520; /* *** MAKE SURE GRADIENT COULD BE COMPUTED *** * */ if (Iv[TOOBIG] == 0) goto L_60; Iv[1] = 65; goto L_999; L_60: if (Iv[HC] <= 0 && Iv[RMAT] <= 0) goto L_610; /* *** COMPUTE D**-1 * GRADIENT *** * */ dig1 = Iv[DIG]; k = dig1; for (i = 1; i <= p; i++) { V[k] = Gg[i]/D[i]; k += 1; } V[DGNORM] = dv2nrm( p, &V[dig1] ); if (Iv[CNVCOD] != 0) goto L_510; if (Iv[MODE] == 0) goto L_440; Iv[MODE] = 0; V[F0] = V[F]; if (Iv[INITS] <= 2) goto L_100; /* *** 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_530; /* *** COME TO NEXT STMT AFTER COMPUTING F.D. HESSIAN FOR INIT. S *** * */ L_80: 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_90; dv2axy( pp1o2, &V[s1], NEGONE, &V[hc1], &V[h1] ); goto L_100; L_90: rmat1 = Iv[RMAT]; dl7sqr( ps, &V[s1], &V[rmat1] ); dv2axy( pp1o2, &V[s1], NEGONE, &V[s1], &V[h1] ); L_100: Iv[1] = 2; /* ---------------------------- MAIN LOOP ----------------------------- * * * *** PRINT ITERATION SUMMARY, CHECK ITERATION LIMIT *** * */ L_110: ditsum( d, gg, iv, liv, lv, p, v, x ); L_120: k = Iv[NITER]; if (k < Iv[MXITER]) goto L_130; Iv[1] = 10; goto L_999; L_130: Iv[NITER] = k + 1; /* *** UPDATE RADIUS *** * */ if (k == 0) goto L_150; step1 = Iv[STEP]; for (i = 1; i <= p; i++) { V[step1] *= D[i]; step1 += 1; } step1 = Iv[STEP]; t = V[RADFAC]*dv2nrm( p, &V[step1] ); if (V[RADFAC] < ONE || t > V[RADIUS]) V[RADIUS] = t; /* *** INITIALIZE FOR START OF NEXT ITERATION *** * */ L_150: x01 = Iv[X0]; V[F0] = V[F]; Iv[IRC] = 4; Iv[H] = -labs( Iv[H] ); Iv[SUSED] = Iv[MODEL]; /* *** COPY X TO X0 *** * */ dv7cpy( p, &V[x01], x ); /* *** CHECK STOPX AND FUNCTION EVALUATION LIMIT *** * */ L_160: ; /* if (STOPX(DUMMY)) then * IV(1) = 11 * GO TO 190 * else */ goto L_180; /* endif * * *** COME HERE WHEN RESTARTING AFTER FUNC. EVAL. LIMIT OR STOPX. * */ L_170: if (V[F] >= V[F0]) goto L_180; V[RADFAC] = ONE; k = Iv[NITER]; goto L_130; L_180: if (Iv[NFCALL] < Iv[MXFCAL] + Iv[NFCOV]) goto L_200; Iv[1] = 9; /* 190 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_430; /*. . . . . . . . . . . . . COMPUTE CANDIDATE STEP . . . . . . . . . . * */ L_200: step1 = Iv[STEP]; w1 = Iv[W]; h1 = Iv[H]; t1 = ONE; if (Iv[MODEL] == 2) goto L_210; t1 = ZERO; /* *** COMPUTE LEVENBERG-MARQUARDT STEP IF POSSIBLE... * */ rmat1 = Iv[RMAT]; if (rmat1 <= 0) goto L_210; qtr1 = Iv[QTR]; if (qtr1 <= 0) goto L_210; ipiv1 = Iv[IPIVOT]; dl7mst( d, gg, Iv[IERR], &Iv[ipiv1], &Iv[KALM], p, &V[qtr1], &V[rmat1], &V[step1], v, &V[w1] ); /* *** 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 DG7QTS MAY HAVE STORED IN V... */ Iv[KAGQT] = -1; goto L_260; L_210: if (h1 > 0) goto L_250; /* *** SET H TO D**-1 * (HC + T1*S) * D**-1. *** * */ h1 = -h1; Iv[H] = h1; Iv[FDH] = 0; j = Iv[HC]; if (j > 0) goto L_220; j = h1; rmat1 = Iv[RMAT]; dl7sqr( p, &V[h1], &V[rmat1] ); L_220: s1 = Iv[S]; for (i = 1; i <= p; i++) { t = ONE/D[i]; for (k = 1; k <= i; k++) { V[h1] = t*(V[j] + t1*V[s1])/D[k]; j += 1; h1 += 1; s1 += 1; } } h1 = Iv[H]; Iv[KAGQT] = -1; /* *** COMPUTE ACTUAL GOLDFELD-QUANDT-TROTTER STEP *** * */ L_250: dig1 = Iv[DIG]; lmat1 = Iv[LMAT]; dg7qts( d, &V[dig1], &V[h1], &Iv[KAGQT], &V[lmat1], p, &V[step1], v, &V[w1] ); if (Iv[KALM] > 0) Iv[KALM] = 0; L_260: if (Iv[IRC] != 6) goto L_270; if (Iv[RESTOR] != 2) goto L_290; rstrst = 2; goto L_300; /* *** CHECK WHETHER EVALUATING F(X0 + STEP) LOOKS WORTHWHILE *** * */ L_270: Iv[TOOBIG] = 0; if (V[DSTNRM] <= ZERO) goto L_290; if (Iv[IRC] != 5) goto L_280; if (V[RADFAC] <= ONE) goto L_280; if (V[PREDUC] > ONEP2*V[FDIF]) goto L_280; if (Iv[RESTOR] != 2) goto L_290; rstrst = 0; goto L_300; /* *** COMPUTE F(X0 + STEP) *** * */ L_280: x01 = Iv[X0]; step1 = Iv[STEP]; dv2axy( p, x, ONE, &V[step1], &V[x01] ); Iv[NFCALL] += 1; Iv[1] = 1; goto L_999; /*. . . . . . . . . . . . . ASSESS CANDIDATE STEP . . . . . . . . . . . * */ L_290: rstrst = 3; L_300: x01 = Iv[X0]; V[RELDX] = drldst( p, d, x, &V[x01] ); da7sst( iv, liv, lv, v ); step1 = Iv[STEP]; lstgst = Iv[STLSTG]; i = Iv[RESTOR] + 1; switch (i) { case 1: goto L_340; case 2: goto L_310; case 3: goto L_320; case 4: goto L_330; } L_310: dv7cpy( p, x, &V[x01] ); goto L_340; L_320: dv7cpy( p, &V[lstgst], &V[step1] ); goto L_340; L_330: dv7cpy( p, &V[step1], &V[lstgst] ); dv2axy( p, x, ONE, &V[step1], &V[x01] ); V[RELDX] = drldst( p, d, x, &V[x01] ); Iv[RESTOR] = rstrst; /* *** IF NECESSARY, SWITCH MODELS *** * */ L_340: if (Iv[SWITCH_] == 0) goto L_350; Iv[H] = -labs( Iv[H] ); Iv[SUSED] += 2; l = Iv[VSAVE]; dv7cpy( NVSAVE, v, &V[l] ); L_350: l = Iv[IRC] - 4; stpmod = Iv[MODEL]; if (l > 0) switch (l) { case 1: goto L_370; case 2: goto L_380; case 3: goto L_390; case 4: goto L_390; case 5: goto L_390; case 6: goto L_390; case 7: goto L_390; case 8: goto L_390; case 9: goto L_500; case 10: goto L_440; } /* *** DECIDE WHETHER TO CHANGE MODELS *** * */ e = V[PREDUC] - V[FDIF]; s1 = Iv[S]; ds7lvm( ps, yy, &V[s1], &V[step1] ); sttsst = HALF*dd7tpr( ps, &V[step1], yy ); if (Iv[MODEL] == 1) sttsst = -sttsst; if (fabs( e + sttsst )*V[FUZZ] >= fabs( e )) goto L_360; /* *** SWITCH MODELS *** * */ Iv[MODEL] = 3 - Iv[MODEL]; if (-2 < l) goto L_400; Iv[H] = -labs( Iv[H] ); Iv[SUSED] += 2; l = Iv[VSAVE]; dv7cpy( NVSAVE, &V[l], v ); goto L_160; L_360: if (-3 < l) goto L_400; /* *** RECOMPUTE STEP WITH NEW RADIUS *** * */ L_370: V[RADIUS] = V[RADFAC]*V[DSTNRM]; goto L_160; /* *** COMPUTE STEP OF LENGTH V(LMAXS) FOR SINGULAR CONVERGENCE TEST * */ L_380: V[RADIUS] = V[LMAXS]; goto L_200; /* *** CONVERGENCE OR FALSE CONVERGENCE *** * */ L_390: Iv[CNVCOD] = l; if (V[F] >= V[F0]) goto L_510; if (Iv[XIRC] == 14) goto L_510; Iv[XIRC] = 14; /*. . . . . . . . . . . . PROCESS ACCEPTABLE STEP . . . . . . . . . . . * */ L_400: Iv[COVMAT] = 0; Iv[REGD] = 0; /* *** SEE WHETHER TO SET V(RADFAC) BY GRADIENT TESTS *** * */ if (Iv[IRC] != 3) goto L_430; step1 = Iv[STEP]; temp1 = Iv[STLSTG]; temp2 = Iv[W]; /* *** SET TEMP1 = HESSIAN * STEP FOR USE IN GRADIENT TESTS *** * */ hc1 = Iv[HC]; if (hc1 <= 0) goto L_410; ds7lvm( p, &V[temp1], &V[hc1], &V[step1] ); goto L_420; L_410: rmat1 = Iv[RMAT]; dl7tvm( p, &V[temp1], &V[rmat1], &V[step1] ); dl7vml( p, &V[temp1], &V[rmat1], &V[temp1] ); L_420: if (stpmod == 1) goto L_430; s1 = Iv[S]; ds7lvm( ps, &V[temp2], &V[s1], &V[step1] ); dv2axy( ps, &V[temp1], ONE, &V[temp2], &V[temp1] ); /* *** SAVE OLD GRADIENT AND COMPUTE NEW ONE *** * */ L_430: Iv[NGCALL] += 1; g01 = Iv[W]; dv7cpy( p, &V[g01], gg ); Iv[1] = 2; Iv[TOOBIG] = 0; goto L_999; /* *** INITIALIZATIONS -- G0 = GG - G0, ETC. *** * */ L_440: g01 = Iv[W]; dv2axy( p, &V[g01], NEGONE, &V[g01], gg ); step1 = Iv[STEP]; temp1 = Iv[STLSTG]; temp2 = Iv[W]; if (Iv[IRC] != 3) goto L_470; /* *** SET V(RADFAC) BY GRADIENT TESTS *** * * *** SET TEMP1 = D**-1 * (HESSIAN * STEP + (G(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 (dv2nrm( p, &V[temp1] ) <= V[DGNORM]*V[TUNER4]) goto L_460; if (dd7tpr( p, gg, &V[step1] ) >= V[GTSTEP]*V[TUNER5]) goto L_470; L_460: V[RADFAC] = V[INCFAC]; /* *** COMPUTE YY VECTOR NEEDED FOR UPDATING S *** * */ L_470: dv2axy( ps, yy, NEGONE, yy, gg ); /* *** DETERMINE SIZING FACTOR V(SIZE) *** * * *** SET TEMP1 = S * STEP *** */ s1 = Iv[S]; ds7lvm( ps, &V[temp1], &V[s1], &V[step1] ); t1 = fabs( dd7tpr( ps, &V[step1], &V[temp1] ) ); t = fabs( dd7tpr( ps, &V[step1], yy ) ); 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_480; ds7lvm( ps, &V[g01], &V[hc1], &V[step1] ); goto L_490; L_480: rmat1 = Iv[RMAT]; dl7tvm( ps, &V[g01], &V[rmat1], &V[step1] ); dl7vml( ps, &V[g01], &V[rmat1], &V[g01] ); L_490: dv2axy( ps, &V[g01], ONE, yy, &V[g01] ); /* *** UPDATE S *** * */ ds7lup( &V[s1], V[COSMIN], ps, V[SIZE], &V[step1], &V[temp1], &V[temp2], &V[g01], &V[WSCALE], yy ); Iv[1] = 2; goto L_110; /*. . . . . . . . . . . . . . MISC. DETAILS . . . . . . . . . . . . . . * * *** BAD PARAMETERS TO ASSESS *** * */ L_500: Iv[1] = 64; goto L_999; /* *** CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE *** * */ L_510: if (Iv[RDREQ] == 0) goto L_600; if (Iv[FDH] != 0) goto L_600; if (Iv[CNVCOD] >= 7) goto L_600; if (Iv[REGD] > 0) goto L_600; if (Iv[COVMAT] > 0) goto L_600; if (labs( Iv[COVREQ] ) >= 3) goto L_560; if (Iv[RESTOR] == 0) Iv[RESTOR] = 2; goto L_530; /* *** COMPUTE FINITE-DIFFERENCE HESSIAN FOR COMPUTING COVARIANCE *** * */ L_520: Iv[RESTOR] = 0; L_530: df7hes( d, gg, &i, iv, liv, lv, p, v, x ); switch (i) { case 1: goto L_540; case 2: goto L_550; case 3: goto L_580; } L_540: Iv[NFCOV] += 1; Iv[NFCALL] += 1; Iv[1] = 1; goto L_999; L_550: Iv[NGCOV] += 1; Iv[NGCALL] += 1; Iv[NFGCAL] = Iv[NFCALL] + Iv[NGCOV]; Iv[1] = 2; goto L_999; L_560: h1 = labs( Iv[H] ); Iv[H] = -h1; pp1o2 = p*(p + 1)/2; rmat1 = Iv[RMAT]; if (rmat1 <= 0) goto L_570; lmat1 = Iv[LMAT]; dv7cpy( pp1o2, &V[lmat1], &V[rmat1] ); V[RCOND] = ZERO; goto L_590; L_570: hc1 = Iv[HC]; Iv[FDH] = h1; dv7cpy( p*(p + 1)/2, &V[h1], &V[hc1] ); /* *** COMPUTE CHOLESKY FACTOR OF FINITE-DIFFERENCE HESSIAN * *** FOR USE IN CALLER*S COVARIANCE CALCULATION... * */ L_580: lmat1 = Iv[LMAT]; h1 = Iv[FDH]; if (h1 <= 0) goto L_600; if (Iv[CNVCOD] == 70) goto L_80; dl7srt( 1, p, &V[lmat1], &V[h1], &i ); Iv[FDH] = -1; V[RCOND] = ZERO; if (i != 0) goto L_600; L_590: Iv[FDH] = -1; step1 = Iv[STEP]; t = dl7svn( p, &V[lmat1], &V[step1], &V[step1] ); if (t <= ZERO) goto L_600; tp = dl7svx( p, &V[lmat1], &V[step1], &V[step1] ); if (tp != ZERO) { t /= tp; if (t > dr7mdc( 4 )) Iv[FDH] = h1; V[RCOND] = t; } L_600: 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_610: Iv[1] = 1400; L_999: return; /* *** LAST LINE OF DG7LIT FOLLOWS *** */ } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dn2lrd( double *dr, long iv[], double l[], long lh, long liv, long lv, long nd, long nn, long p, double r[], double rd[], double v[]) { #define DR(I_,J_) (*(dr+(I_)*(nd)+(J_))) long int cov, i, j, m, step1, _i, _r; double a, s, t; static double onev[1]; static int _aini = 1; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Iv = &iv[0] - 1; double *const L = &l[0] - 1; double *const Onev = &onev[0] - 1; double *const R = &r[0] - 1; double *const Rd = &rd[0] - 1; double *const V = &v[0] - 1; /* end of OFFSET VECTORS */ if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */ Onev[1] = 1.e0; _aini = 0; } /* *** COMPUTE REGRESSION DIAGNOSTIC AND DEFAULT COVARIANCE MATRIX FOR * DRN2G *** * * *** PARAMETERS *** * */ /* *** CODED BY DAVID M. GAY (WINTER 1982, FALL 1983) *** * * *** EXTERNAL FUNCTIONS AND SUBROUTINES *** * */ /* *** LOCAL VARIABLES *** * */ /* *** CONSTANTS *** * */ /* *** IV SUBSCRIPTS *** * */ /* +++++++++++++++++++++++++++++++ BODY +++++++++++++++++++++++++++++++ * */ step1 = Iv[STEP]; i = Iv[RDREQ]; if (i <= 0) goto L_999; if ((i%4) < 2) goto L_30; dv7scp( nn, rd, NEGONE ); for (i = 1; i <= nn; i++) { a = SQ(R[i]); m = step1; for (j = 1; j <= p; j++) { V[m] = DR(j - 1,i - 1); m += 1; } dl7ivm( p, &V[step1], l, &V[step1] ); s = dd7tpr( p, &V[step1], &V[step1] ); t = ONE - s; if (t <= ZERO) goto L_20; a = a*s/t; Rd[i] = sqrt( a ); L_20: ; } L_30: if (Iv[MODE] - p < 2) goto L_999; /* *** COMPUTE DEFAULT COVARIANCE MATRIX *** * */ cov = labs( Iv[H] ); for (i = 1; i <= nn; i++) { m = step1; for (j = 1; j <= p; j++) { V[m] = DR(j - 1,i - 1); m += 1; } dl7ivm( p, &V[step1], l, &V[step1] ); dl7itv( p, &V[step1], l, &V[step1] ); do7prd( 1, lh, p, &V[cov], onev, &V[step1], &V[step1] ); } L_999: return; /* *** LAST CARD OF DN2LRD FOLLOWS *** */ #undef DR } /* end of function */ void /*FUNCTION*/ dc7vfn( long iv[], double l[], long lh, long liv, long lv, long n, long p, double v[]) { long int cov, i; static double half = 0.5e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Iv = &iv[0] - 1; double *const L = &l[0] - 1; double *const V = &v[0] - 1; /* end of OFFSET VECTORS */ /* *** FINISH COVARIANCE COMPUTATION FOR DRN2G, DRNSG *** * */ /* *** LOCAL VARIABLES *** * */ /* *** SUBSCRIPTS FOR IV AND V *** * */ /* *** BODY *** * */ Iv[1] = Iv[CNVCOD]; i = Iv[MODE] - p; Iv[MODE] = 0; Iv[CNVCOD] = 0; if (Iv[FDH] <= 0) goto L_999; if (ipow(i - 2,2) == 1) Iv[REGD] = 1; if ((Iv[RDREQ]%2) != 1) goto L_999; /* *** FINISH COMPUTING COVARIANCE MATRIX = INVERSE OF F.D. HESSIAN. * */ cov = labs( Iv[H] ); Iv[FDH] = 0; if (Iv[COVMAT] != 0) goto L_999; if (i >= 2) goto L_10; dl7nvr( p, &V[cov], l ); dl7tsq( p, &V[cov], &V[cov] ); L_10: dv7scl( lh, &V[cov], V[F]/(half*(double)( max( 1, n - p ) )), &V[cov] ); Iv[COVMAT] = cov; L_999: return; /* *** LAST LINE OF DC7VFN FOLLOWS *** */ } /* end of function */ /* PARAMETER translations */ #define DELTA 52 #define DELTA0 44 #define DLTFDC 42 #define FX 53 #define NEGPT5 (-0.5e0) #define SAVEI 63 #define TWO 2.e0 #define XMSAVE 51 /* end of PARAMETER translations */ void /*FUNCTION*/ df7hes( double d[], double gg[], long *irt, long iv[], long liv, long lv, long p, double v[], double x[]) { long int gsave1, hes, hmi, hpi, hpm, i, k, kind, l, m, mm1, mm1o2, pp1o2, stp0, stpi, stpm; double del; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const D = &d[0] - 1; double *const Gg = &gg[0] - 1; long *const Iv = &iv[0] - 1; double *const V = &v[0] - 1; double *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)). * * *** IF IV(COVREQ) .GE. 0 THEN DF7HES 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 G. * 3 = DONE. * * * *** PARAMETER DECLARATIONS *** * */ /* *** LOCAL VARIABLES *** * */ /* *** EXTERNAL SUBROUTINES *** * */ /* DV7CPY.... COPY ONE VECTOR TO ANOTHER. * * *** SUBSCRIPTS FOR IV AND V *** * */ /* ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ * */ *irt = 4; kind = Iv[COVREQ]; m = Iv[MODE]; if (m > 0) goto L_10; Iv[H] = -labs( Iv[H] ); Iv[FDH] = 0; Iv[KAGQT] = -1; V[FX] = V[F]; L_10: if (m > p) goto L_999; if (kind < 0) goto L_110; /* *** COMPUTE FINITE-DIFFERENCE HESSIAN USING BOTH FUNCTION AND * *** GRADIENT VALUES. * */ gsave1 = Iv[W] + p; if (m > 0) goto L_20; /* *** FIRST CALL ON DF7HES. SET GSAVE = G, TAKE FIRST STEP *** */ dv7cpy( p, &V[gsave1], gg ); Iv[SWITCH_] = Iv[NFGCAL]; goto L_90; L_20: del = V[DELTA]; X[m] = V[XMSAVE]; if (Iv[TOOBIG] == 0) goto L_40; /* *** HANDLE OVERSIZE V(DELTA) *** * */ if (del*X[m] > ZERO) goto L_30; /* *** WE ALREADY TRIED SHRINKING V(DELTA), SO QUIT *** */ Iv[FDH] = -2; goto L_220; /* *** TRY SHRINKING V(DELTA) *** */ L_30: del *= NEGPT5; goto L_100; L_40: hes = -Iv[H]; /* *** SET GG = (GG - GSAVE)/DEL *** * */ for (i = 1; i <= p; i++) { Gg[i] = (Gg[i] - V[gsave1])/del; 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_70; /* *** 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++) { V[k] = HALF*(V[k] + Gg[i]); k += 1; } /* *** ADD H(I,M) = GG(I) FOR I = M TO P *** * */ L_70: l += 1; for (i = m; i <= p; i++) { V[l] = Gg[i]; l += i; } L_90: m += 1; Iv[MODE] = m; if (m > p) goto L_210; /* *** CHOOSE NEXT FINITE-DIFFERENCE STEP, RETURN TO GET GG THERE *** * */ del = V[DELTA0]*fmax( ONE/D[m], fabs( X[m] ) ); if (X[m] < ZERO) del = -del; V[XMSAVE] = X[m]; L_100: X[m] += del; V[DELTA] = del; *irt = 2; goto L_999; /* *** COMPUTE FINITE-DIFFERENCE HESSIAN USING FUNCTION VALUES ONLY. * */ L_110: stp0 = Iv[W] + p - 1; mm1 = m - 1; mm1o2 = m*mm1/2; if (m > 0) goto L_120; /* *** FIRST CALL ON DF7HES. *** */ Iv[SAVEI] = 0; goto L_200; L_120: i = Iv[SAVEI]; hes = -Iv[H]; if (i > 0) goto L_180; if (Iv[TOOBIG] == 0) goto L_140; /* *** HANDLE OVERSIZE STEP *** * */ stpm = stp0 + m; del = V[stpm]; if (del*X[XMSAVE] > ZERO) goto L_130; /* *** WE ALREADY TRIED SHRINKING THE STEP, SO QUIT *** */ Iv[FDH] = -2; goto L_220; /* *** TRY SHRINKING THE STEP *** */ L_130: del *= NEGPT5; X[m] = X[XMSAVE] + del; V[stpm] = del; *irt = 1; goto L_999; /* *** SAVE F(X + STP(M)*E(M)) IN H(P,M) *** * */ L_140: pp1o2 = p*(p - 1)/2; hpm = hes + pp1o2 + mm1; V[hpm] = V[F]; /* *** START COMPUTING ROW M OF THE FINITE-DIFFERENCE HESSIAN H. *** * */ hmi = hes + mm1o2; if (mm1 == 0) goto L_160; hpi = hes + pp1o2; for (i = 1; i <= mm1; i++) { V[hmi] = V[FX] - (V[F] + V[hpi]); hmi += 1; hpi += 1; } L_160: V[hmi] = V[F] - TWO*V[FX]; /* *** COMPUTE FUNCTION VALUES NEEDED TO COMPLETE ROW M OF H. *** * */ i = 1; L_170: Iv[SAVEI] = i; stpi = stp0 + i; V[DELTA] = X[i]; X[i] += V[stpi]; if (i == m) X[i] = V[XMSAVE] - V[stpi]; *irt = 1; goto L_999; L_180: X[i] = V[DELTA]; if (Iv[TOOBIG] == 0) goto L_190; /* *** PUNT IN THE EVENT OF AN OVERSIZE STEP *** */ Iv[FDH] = -2; goto L_220; /* *** FINISH COMPUTING H(M,I) *** * */ L_190: stpi = stp0 + i; hmi = hes + mm1o2 + i - 1; stpm = stp0 + m; V[hmi] = (V[hmi] + V[F])/(V[stpi]*V[stpm]); i += 1; if (i <= m) goto L_170; Iv[SAVEI] = 0; X[m] = V[XMSAVE]; L_200: m += 1; Iv[MODE] = m; if (m > p) goto L_210; /* *** 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. * */ del = V[DLTFDC]*fmax( ONE/D[m], fabs( X[m] ) ); if (X[m] < ZERO) del = -del; V[XMSAVE] = X[m]; X[m] += del; stpm = stp0 + m; V[stpm] = del; *irt = 1; goto L_999; /* *** RESTORE V(F), ETC. *** * */ L_210: Iv[FDH] = hes; L_220: V[F] = V[FX]; *irt = 3; if (kind < 0) goto L_999; Iv[NFGCAL] = Iv[SWITCH_]; gsave1 = Iv[W] + p; dv7cpy( p, gg, &V[gsave1] ); goto L_999; L_999: return; /* *** LAST CARD OF DF7HES FOLLOWS *** */ } /* end of function */ /* PARAMETER translations */ #define COVPRT 14 #define NEEDHD 36 #define PRUNIT 21 #define STATPR 23 /* end of PARAMETER translations */ void /*FUNCTION*/ dn2cvp( long iv[], long liv, long lv, long p, double v[]) { long int cov1, i, i1, ii, pu; double t; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Iv = &iv[0] - 1; double *const V = &v[0] - 1; /* end of OFFSET VECTORS */ /* *** PRINT COVARIANCE MATRIX FOR DRN2G *** * 6/27/90 CLL changed 'SCALE' to 'VARFAC' in output labels. * ------------------------------------------------------------------ */ long int j, k; /* INTEGER J */ /* *** LOCAL VARIABLES *** * */ /* *** IV SUBSCRIPTS *** * */ /* *** BODY *** * */ /*++(~.C.) Default UNITNO='(PU,' *++(.C.) Default UNITNO='(*,' *++ Replace "(*, " = UNITNO * */ if (Iv[1] > 8) goto L_999; pu = Iv[PRUNIT]; if (pu == 0) goto L_999; cov1 = Iv[COVMAT]; if (-2 == cov1) { printf("\n ++++++ OVERSIZE STEPS IN COMPUTING COVARIANCE +++++\n"); } if (Iv[STATPR] == 0) goto L_30; if (Iv[NFCOV] > 0) { printf("\n %4ld EXTRA FUNC. EVALS FOR COVARIANCE AND DIAGNOSTICS.\n", Iv[NFCOV]); } if (Iv[NGCOV] > 0) { printf(" %4ld EXTRA GRAD. EVALS FOR COVARIANCE AND DIAGNOSTICS.\n", Iv[NGCOV]); } L_30: if (Iv[COVPRT] <= 0) goto L_999; if (Iv[REGD] <= 0 && cov1 <= 0) goto L_70; Iv[NEEDHD] = 1; t = SQ(V[RCOND]); if (labs( Iv[COVREQ] ) > 2) goto L_50; printf("\n RECIPROCAL CONDITION OF F.D. HESSIAN = AT MOST%10.2g\n", t); goto L_70; L_50: printf("\n RECIPROCAL CONDITION OF (J**T)*J = AT LEAST%10.2g\n", t); L_70: if ((Iv[COVPRT]%2) == 0) goto L_999; Iv[NEEDHD] = 1; switch (IARITHIF(cov1)) { case -1: goto L_80; case 0: goto L_110; case 1: goto L_130; } L_80: if (-1 == cov1) { printf("\n ++++++ INDEFINITE COVARIANCE MATRIX ++++++\n"); } goto L_999; L_110: printf("\n ++++++ COVARIANCE MATRIX NOT COMPUTED ++++++\n"); goto L_999; L_130: i = labs( Iv[COVREQ] ); if (i <= 1) { printf("\n COVARIANCE = VARFAC * H**-1 * (J**T * J) * H**-1\n WHERE H = F.D. HESSIAN\n \n"); } if (i == 2) { printf(" \n COVARIANCE = VARFAC * H**-1, WHERE H = FINITE-DIFFERENCE HESSIAN\n \n"); } if (i > 2) { printf("\n COVARIANCE = VARFAC * J**T * J\n \n"); } ii = cov1 - 1; for (i = 1; i <= p; i++) { i1 = ii + 1; ii += i; printf(" ROW%3ld ", i); for (j = i1; j <= ii; j+=5){ for (k = j; k <= (j <= ii - 5 ? j+4 : ii); k++) printf("%12.3g", V[k]); printf("\n"); if (j <= ii - 5) printf(" ");} } /* WRITE(*, '('' ROW'',I3,2X,5g12.3/(9X,5g12.3))')I,(V(J),J=I1,II) * */ L_999: return; /* *** LAST CARD OF DN2CVP FOLLOWS *** */ } /* end of function */ void /*FUNCTION*/ dn2rdp( long iv[], long liv, long n, double rd[]) { long int pu; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Iv = &iv[0] - 1; double *const Rd = &rd[0] - 1; /* end of OFFSET VECTORS */ /* *** PRINT REGRESSION DIAGNOSTICS FOR MLPSL AND NL2S1 *** * * ------------------------------------------------------------------ *++ Code for .C. is active */ long int j,k; /*++ End */ /* *** IV SUBSCRIPTS *** * */ /* DATA COVPRT/14/, NEEDHD/36/, PRUNIT/21/, RDREQ/57/, REGD/67/ */ /* ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ * */ pu = Iv[PRUNIT]; if (pu == 0) goto L_999; if (Iv[COVPRT] < 2) goto L_999; if (Iv[REGD] <= 0) goto L_999; Iv[NEEDHD] = 1; printf(" REGRESSION DIAGNOSTIC = SQRT(G(I)**T * H(I)**-1 *G(I))...\n\n"); for(j=0; j < n; j+=6){ for (k = j; k < (j < n - 6 ? j+6 : n); k++) printf("%12.3g", rd[k]); printf("\n");} L_999: return; /* WRITE(*, '(6g12.3)') RD * * *** LAST CARD OF DN2RDP FOLLOWS *** */ } /* end of function */ /* PARAMETER translations */ #undef ZERO #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ do7prd( long l, long ls, long pp, double ss[], double ww[], double *yy, double *z) { #define YY(I_,J_) (*(yy+(I_)*(pp)+(J_))) #define Z(I_,J_) (*(z+(I_)*(pp)+(J_))) long int i, j, k, m; double wk, yi; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Ss = &ss[0] - 1; double *const Ww = &ww[0] - 1; /* end of OFFSET VECTORS */ /* *** FOR I = 1..L, SET SS = SS + WW(I)*YY(.,I)*(Z(.,I)**T), I.E., * *** ADD WW(I) TIMES THE OUTER PRODUCT OF Y(.,I) AND Z(.,I). * * ------------------------------------------------------------------ */ /* DIMENSION SS(PP*(PP+1)/2) * */ for (k = 1; k <= l; k++) { wk = Ww[k]; if (wk == ZERO) goto L_30; m = 1; for (i = 1; i <= pp; i++) { yi = wk*YY(k - 1,i - 1); for (j = 1; j <= i; j++) { Ss[m] += yi*Z(k - 1,j - 1); m += 1; } } L_30: ; } return; /* *** LAST CARD OF DO7PRD FOLLOWS *** */ #undef Z #undef YY } /* end of function */ /* PARAMETER translations */ #undef ZERO #define ZERO 0.e0 /* end of PARAMETER translations */ void /*FUNCTION*/ dl7nvr( long n, double lin[], double l[]) { long int i, ii, im1, j0, j1, jj, k, k0, np1; double t; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const L = &l[0] - 1; double *const Lin = &lin[0] - 1; /* end of OFFSET VECTORS */ /* *** COMPUTE LIN = L**-1, BOTH N X N LOWER TRIANG. STORED *** * *** COMPACTLY BY ROWS. LIN AND L MAY SHARE THE SAME STORAGE. *** * * *** PARAMETERS *** * * ------------------------------------------------------------------ */ /* DIMENSION L(N*(N+1)/2), LIN(N*(N+1)/2) * * *** LOCAL VARIABLES *** * */ /* *** BODY *** * */ np1 = n + 1; j0 = n*(np1)/2; for (ii = 1; ii <= n; ii++) { i = np1 - ii; Lin[j0] = ONE/L[j0]; if (i <= 1) goto L_999; j1 = j0; im1 = i - 1; for (jj = 1; jj <= im1; jj++) { t = ZERO; j0 = j1; k0 = j1 - jj; for (k = 1; k <= jj; k++) { t += -L[k0]*Lin[j0]; j0 -= 1; k0 += k - i; } Lin[j0] = t/L[k0]; } j0 -= 1; } L_999: return; /* *** LAST CARD OF DL7NVR FOLLOWS *** */ } /* end of function */ void /*FUNCTION*/ dl7tsq( long n, double a[], double l[]) { long int i, i1, ii, iim1, j, k, m; double lii, lj; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const A = &a[0] - 1; double *const L = &l[0] - 1; /* end of OFFSET VECTORS */ /* *** SET A TO LOWER TRIANGLE OF (L**T) * L *** * * *** L = N X N LOWER TRIANG. MATRIX STORED ROWWISE. *** * *** A IS ALSO STORED ROWWISE AND MAY SHARE STORAGE WITH L. *** * * ------------------------------------------------------------------ */ /* DIMENSION A(N*(N+1)/2), L(N*(N+1)/2) * */ ii = 0; for (i = 1; i <= n; i++) { i1 = ii + 1; ii += i; m = 1; if (i == 1) goto L_30; iim1 = ii - 1; for (j = i1; j <= iim1; j++) { lj = L[j]; for (k = i1; k <= j; k++) { A[m] += lj*L[k]; m += 1; } } L_30: lii = L[ii]; for (j = i1; j <= ii; j++) { A[j] = lii*L[j]; } } return; /* *** LAST CARD OF DL7TSQ FOLLOWS *** */ } /* end of function */