/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:05 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "divset.h" #include #include #include #include /* File: DIVSET.for 28 subrs used by all main subrs of the * David Gay & Linda Kaufman nonlinear LS package. * Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 2001-06-18 DIVSET Krogh Replaced ". AND." with " .AND." *>> 1998-10-29 DIVSET Krogh Moved external statements up for mangle. *>> 1996-06-18 DIVSET Krogh More work for for C conversion. *>> 1995-11-15 DIVSET Krogh Moved formats up for C conversion. *>> 1994-11-02 DIVSET Krogh Changes to use M77CON *>> 1992-04-27 DIVSET CLL Removed unreferenced stmt labels. *>> 1992-04-13 CLL Change from Hollerith to '...' syntax in formats. *>> 1991-05-21 CLL JPL Changing (1) to (*) in declarations. *>> 1990-06-15 CLL JPL *>> 1990-03-21 CLL JPL *>> 1990-02-16 CLL JPL */ /* PARAMETER translations */ #define ALGSAV 51 #define COVPRT 14 #define COVREQ 15 #define DRADPR 101 #define DTYPE 16 #define HC 71 #define IERR 75 #define INITH 25 #define INITS 25 #define IPIVOT 76 #define IVNEED 3 #define LASTIV 44 #define LASTV 45 #define LMAT 42 #define MXFCAL 17 #define MXITER 18 #define NFCOV 52 #define NGCOV 53 #define NVDFLT 50 #define NVSAVE 9 #define OUTLEV 19 #define PARPRT 20 #define PARSAV 49 #define PERM 58 #define PRUNIT 21 #define QRTYP 80 #define RDREQ 57 #define RMAT 78 #define SOLPRT 22 #define STATPR 23 #define VNEED 4 #define VSAVE 60 #define X0PRT 24 /* end of PARAMETER translations */ void /*FUNCTION*/ divset( long alg, long iv[], long liv, long lv, double v[]) { long int _l0, alg1, miv, mv, _i, _r; static long int miniv[4], minv[4]; static int _aini = 1; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Iv = &iv[0] - 1; long *const Miniv = &miniv[0] - 1; long *const Minv = &minv[0] - 1; double *const V = &v[0] - 1; /* end of OFFSET VECTORS */ if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */ Miniv[1] = 82; Miniv[2] = 59; Miniv[3] = 103; Miniv[4] = 103; Minv[1] = 98; Minv[2] = 71; Minv[3] = 101; Minv[4] = 85; _aini = 0; } /* *** SUPPLY ***SOL (VERSION 2.3) DEFAULT VALUES TO IV AND V *** * * *** ALG = 1 MEANS REGRESSION CONSTANTS. * *** ALG = 2 MEANS GENERAL UNCONSTRAINED OPTIMIZATION CONSTANTS. * ------------------------------------------------------------------ * 1990-06-25 CLL changed default settings for printing so the * default is no printing. This involved setting COVPRT, OUTLEV, * PARPRT, SOLPRT, STATPR, and X0PRT to zero. Previously these were * all set to 1 except that COVPRT = 3. * ------------------------------------------------------------------ * */ /*--D replaces "?": ?IVSET,?V7DFL,?R7MDC,?L7SVN,?D7TPR,?V2AXY,?V2NRM *--& ?Q7APL,?D7UPD,?V7SCP,?Q7RAD,?V7SCL,?PARCK,?V7CPY,?S7LVM,?S7LUP *--& ?L7MST,?L7ITV,?L7IVM,?G7QTS,?L7SRT,?L7SVX,?A7SST,?ITSUM,?L7SQR *--& ?L7TVM,?L7VML,?RLDST * I1MACH... RETURNS MACHINE-DEPENDENT INTEGER CONSTANTS. * DV7DFL.... PROVIDES DEFAULT VALUES TO V. * */ /* *** SUBSCRIPTS FOR IV *** * */ /* *** IV SUBSCRIPT VALUES *** * */ /* ------------------------------ BODY -------------------------------- * * I1MACH(2) returns the unit number for standard output. * */ if (PRUNIT <= liv) Iv[PRUNIT] = 6; if (ALGSAV <= liv) Iv[ALGSAV] = alg; if (alg < 1 || alg > 4) goto L_40; miv = Miniv[alg]; if (liv < miv) goto L_20; mv = Minv[alg]; if (lv < mv) goto L_30; alg1 = ((alg - 1)%2) + 1; dv7dfl( alg1, lv, v ); Iv[1] = 12; if (alg > 2) Iv[DRADPR] = 1; Iv[IVNEED] = 0; Iv[LASTIV] = miv; Iv[LASTV] = mv; Iv[LMAT] = mv + 1; Iv[MXFCAL] = 200; Iv[MXITER] = 150; Iv[OUTLEV] = 0; Iv[PARPRT] = 0; Iv[PERM] = miv + 1; Iv[SOLPRT] = 0; Iv[STATPR] = 0; Iv[VNEED] = 0; Iv[X0PRT] = 0; if (alg1 >= 2) goto L_10; /* *** REGRESSION VALUES * */ Iv[COVPRT] = 0; Iv[COVREQ] = 1; Iv[DTYPE] = 1; Iv[HC] = 0; Iv[IERR] = 0; Iv[INITS] = 0; Iv[IPIVOT] = 0; Iv[NVDFLT] = 32; Iv[VSAVE] = 58; if (alg > 2) Iv[VSAVE] += 3; Iv[PARSAV] = Iv[VSAVE] + NVSAVE; Iv[QRTYP] = 1; Iv[RDREQ] = 3; Iv[RMAT] = 0; goto L_999; /* *** GENERAL OPTIMIZATION VALUES * */ L_10: Iv[DTYPE] = 0; Iv[INITH] = 1; Iv[NFCOV] = 0; Iv[NGCOV] = 0; Iv[NVDFLT] = 25; Iv[PARSAV] = 47; if (alg > 2) Iv[PARSAV] = 61; goto L_999; L_20: Iv[1] = 15; goto L_999; L_30: Iv[1] = 16; goto L_999; L_40: Iv[1] = 67; L_999: return; /* *** LAST CARD OF DIVSET FOLLOWS *** */ } /* end of function */ /* PARAMETER translations */ #define AFCTOL 31 #define BIAS 43 #define COSMIN 47 #define D0INIT 40 #define DECFAC 22 #define DELTA0 44 #define DFAC 41 #define DINIT 38 #define DLTFDC 42 #define DLTFDJ 43 #define DTINIT 39 #define EPSLON 19 #define ETA0 42 #define FUZZ 45 #define HUBERC 48 #define INCFAC 23 #define LMAX0 35 #define LMAXS 36 #define ONE 1.e0 #define PHMNFC 20 #define PHMXFC 21 #define RDFCMN 24 #define RDFCMX 25 #define RFCTOL 32 #define RLIMIT 46 #define RSPTOL 49 #define SCTOL 37 #define SIGMIN 50 #define THREE 3.e0 #define TUNER1 26 #define TUNER2 27 #define TUNER3 28 #define TUNER4 29 #define TUNER5 30 #define XCTOL 33 #define XFTOL 34 /* end of PARAMETER translations */ void /*FUNCTION*/ dv7dfl( long alg, long lv, double v[]) { double machep, mepcrt, sqteps; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const V = &v[0] - 1; /* end of OFFSET VECTORS */ /* *** SUPPLY ***SOL (VERSION 2.3) DEFAULT VALUES TO V *** * * *** ALG = 1 MEANS REGRESSION CONSTANTS. * *** ALG = 2 MEANS GENERAL UNCONSTRAINED OPTIMIZATION CONSTANTS. * */ /* DR7MDC... RETURNS MACHINE-DEPENDENT CONSTANTS * */ /* *** SUBSCRIPTS FOR V *** * */ /* *** V SUBSCRIPT VALUES *** * */ /* ------------------------------ BODY -------------------------------- * */ machep = dr7mdc( 3 ); V[AFCTOL] = 1.e-20; if (machep > 1.e-10) V[AFCTOL] = SQ(machep); V[DECFAC] = 0.5e0; sqteps = dr7mdc( 4 ); V[DFAC] = 0.6e0; V[DTINIT] = 1.e-6; mepcrt = pow(machep,ONE/THREE); V[D0INIT] = 1.e0; V[EPSLON] = 0.1e0; V[INCFAC] = 2.e0; V[LMAX0] = 1.e0; V[LMAXS] = 1.e0; V[PHMNFC] = -0.1e0; V[PHMXFC] = 0.1e0; V[RDFCMN] = 0.1e0; V[RDFCMX] = 4.e0; V[RFCTOL] = fmax( 1.e-10, SQ(mepcrt) ); V[SCTOL] = V[RFCTOL]; V[TUNER1] = 0.1e0; V[TUNER2] = 1.e-4; V[TUNER3] = 0.75e0; V[TUNER4] = 0.5e0; V[TUNER5] = 0.75e0; V[XCTOL] = sqteps; V[XFTOL] = 1.e2*machep; if (alg >= 2) goto L_10; /* *** REGRESSION VALUES * */ V[COSMIN] = fmax( 1.e-6, 1.e2*machep ); V[DINIT] = 0.e0; V[DELTA0] = sqteps; V[DLTFDC] = mepcrt; V[DLTFDJ] = sqteps; V[FUZZ] = 1.5e0; V[HUBERC] = 0.7e0; V[RLIMIT] = dr7mdc( 5 ); V[RSPTOL] = 1.e-3; V[SIGMIN] = 1.e-4; goto L_999; /* *** GENERAL OPTIMIZATION VALUES * */ L_10: V[BIAS] = 0.8e0; V[DINIT] = -1.0e0; V[ETA0] = 1.0e3*machep; L_999: return; /* *** LAST CARD OF DV7DFL FOLLOWS *** */ } /* end of function */ double /*FUNCTION*/ dr7mdc( long k) { long int _l0; double dr7mdc_v; static double big = 0.e0; static double eta = 0.e0; static double machep = 0.e0; static double zero = 0.e0; /* *** RETURN MACHINE DEPENDENT CONSTANTS USED BY NL2SOL *** * */ /* *** THE CONSTANT RETURNED DEPENDS ON K... * * *** K = 1... SMALLEST POS. ETA SUCH THAT -ETA EXISTS. * *** K = 2... SQUARE ROOT OF ETA. * *** K = 3... UNIT ROUNDOFF = SMALLEST POS. NO. MACHEP SUCH * *** THAT 1 + MACHEP .GT. 1 .AND. 1 - MACHEP .LT. 1. * *** K = 4... SQUARE ROOT OF MACHEP. * *** K = 5... SQUARE ROOT OF BIG (SEE K = 6). * *** K = 6... LARGEST MACHINE NO. BIG SUCH THAT -BIG EXISTS. * */ if (big > zero) goto L_1; big = DBL_MAX; eta = DBL_MIN; machep = DBL_EPSILON; L_1: ; /* ------------------------------ BODY -------------------------------- * */ switch (k) { case 1: goto L_10; case 2: goto L_20; case 3: goto L_30; case 4: goto L_40; case 5: goto L_50; case 6: goto L_60; } L_10: dr7mdc_v = eta; goto L_999; L_20: dr7mdc_v = sqrt( 256.e0*eta )/16.e0; goto L_999; L_30: dr7mdc_v = machep; goto L_999; L_40: dr7mdc_v = sqrt( machep ); goto L_999; L_50: dr7mdc_v = sqrt( big/256.e0 )*16.e0; goto L_999; L_60: dr7mdc_v = big; L_999: return( dr7mdc_v ); /* *** LAST CARD OF DR7MDC FOLLOWS *** */ } /* end of function */ /* PARAMETER translations */ #define HALF 0.5e0 #define R9973 9973.e0 #define ZERO 0.e0 /* end of PARAMETER translations */ double /*FUNCTION*/ dl7svn( long p, double l[], double x[], double y[]) { long int i, ii, ix, j, j0, ji, jj, jjj, jm1, pm1; double b, dl7svn_v, sminus, splus, t, xminus, xplus; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const L = &l[0] - 1; double *const X = &x[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* *** ESTIMATE SMALLEST SING. VALUE OF PACKED LOWER TRIANG. MATRIX L * * *** PARAMETER DECLARATIONS *** * */ /* DIMENSION L(P*(P+1)/2) * * ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * * *** PURPOSE *** * * THIS FUNCTION RETURNS A GOOD OVER-ESTIMATE OF THE SMALLEST * SINGULAR VALUE OF THE PACKED LOWER TRIANGULAR MATRIX L. * * *** PARAMETER DESCRIPTION *** * * P (IN) = THE ORDER OF L. L IS A P X P LOWER TRIANGULAR MATRIX. * L (IN) = ARRAY HOLDING THE ELEMENTS OF L IN ROW ORDER, I.E. * L(1,1), L(2,1), L(2,2), L(3,1), L(3,2), L(3,3), ETC. * X (OUT) IF DL7SVN RETURNS A POSITIVE VALUE, THEN X IS A NORMALIZED * APPROXIMATE LEFT SINGULAR VECTOR CORRESPONDING TO THE * SMALLEST SINGULAR VALUE. THIS APPROXIMATION MAY BE VERY * CRUDE. IF DL7SVN RETURNS ZERO, THEN SOME COMPONENTS OF X * ARE ZERO AND THE REST RETAIN THEIR INPUT VALUES. * Y (OUT) IF DL7SVN RETURNS A POSITIVE VALUE, THEN Y = (L**-1)*X IS AN * UNNORMALIZED APPROXIMATE RIGHT SINGULAR VECTOR CORRESPOND- * ING TO THE SMALLEST SINGULAR VALUE. THIS APPROXIMATION * MAY BE CRUDE. IF DL7SVN RETURNS ZERO, THEN Y RETAINS ITS * INPUT VALUE. THE CALLER MAY PASS THE SAME VECTOR FOR X * AND Y (NONSTANDARD FORTRAN USAGE), IN WHICH CASE Y OVER- * WRITES X (FOR NONZERO DL7SVN RETURNS). * * *** ALGORITHM NOTES *** * * THE ALGORITHM IS BASED ON (1), WITH THE ADDITIONAL PROVISION THAT * DL7SVN = 0 IS RETURNED IF THE SMALLEST DIAGONAL ELEMENT OF L * (IN MAGNITUDE) IS NOT MORE THAN THE UNIT ROUNDOFF TIMES THE * LARGEST. THE ALGORITHM USES A RANDOM NUMBER GENERATOR PROPOSED * IN (4), WHICH PASSES THE SPECTRAL TEST WITH FLYING COLORS -- SEE * (2) AND (3). * * *** SUBROUTINES AND FUNCTIONS CALLED *** * * DV2NRM - FUNCTION, RETURNS THE 2-NORM OF A VECTOR. * * *** REFERENCES *** * * (1) CLINE, A., MOLER, C., STEWART, G., AND WILKINSON, J.H.(1977), * AN ESTIMATE FOR THE CONDITION NUMBER OF A MATRIX, REPORT * TM-310, APPLIED MATH. DIV., ARGONNE NATIONAL LABORATORY. * * (2) HOAGLIN, D.C. (1976), THEORETICAL PROPERTIES OF CONGRUENTIAL * RANDOM-NUMBER GENERATORS -- AN EMPIRICAL VIEW, * MEMORANDUM NS-340, DEPT. OF STATISTICS, HARVARD UNIV. * * (3) KNUTH, D.E. (1969), THE ART OF COMPUTER PROGRAMMING, VOL. 2 * (SEMINUMERICAL ALGORITHMS), ADDISON-WESLEY, READING, MASS. * * (4) SMITH, C.S. (1971), MULTIPLICATIVE PSEUDO-RANDOM NUMBER * GENERATORS WITH PRIME MODULUS, J. ASSOC. COMPUT. MACH. 18, * PP. 586-593. * * *** HISTORY *** * * DESIGNED AND CODED BY DAVID M. GAY (WINTER 1977/SUMMER 1978). * * *** GENERAL *** * * THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH * SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS * MCS-7600324, DCR75-10143, 76-14311DSS, AND MCS76-11989. * * ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * * *** LOCAL VARIABLES *** * */ /* *** CONSTANTS *** * */ /* *** EXTERNAL FUNCTIONS AND SUBROUTINES *** * */ /* *** BODY *** * */ ix = 2; pm1 = p - 1; /* *** FIRST CHECK WHETHER TO RETURN DL7SVN = 0 AND INITIALIZE X *** * */ ii = 0; j0 = p*pm1/2; jj = j0 + p; if (L[jj] == ZERO) goto L_110; ix = (3432*ix)%9973; b = HALF*(ONE + (double)( ix )/R9973); xplus = b/L[jj]; X[p] = xplus; if (p <= 1) goto L_60; for (i = 1; i <= pm1; i++) { ii += i; if (L[ii] == ZERO) goto L_110; ji = j0 + i; X[i] = xplus*L[ji]; } /* *** SOLVE (L**T)*X = B, WHERE THE COMPONENTS OF B HAVE RANDOMLY * *** CHOSEN MAGNITUDES IN (.5,1) WITH SIGNS CHOSEN TO MAKE X LARGE. * * DO J = P-1 TO 1 BY -1... */ for (jjj = 1; jjj <= pm1; jjj++) { j = p - jjj; /* *** DETERMINE X(J) IN THIS ITERATION. NOTE FOR I = 1,2,...,J * *** THAT X(I) HOLDS THE CURRENT PARTIAL SUM FOR ROW I. */ ix = (3432*ix)%9973; b = HALF*(ONE + (double)( ix )/R9973); xplus = b - X[j]; xminus = -b - X[j]; splus = fabs( xplus ); sminus = fabs( xminus ); jm1 = j - 1; j0 = j*jm1/2; jj = j0 + j; xplus /= L[jj]; xminus /= L[jj]; if (jm1 == 0) goto L_30; for (i = 1; i <= jm1; i++) { ji = j0 + i; splus += fabs( X[i] + L[ji]*xplus ); sminus += fabs( X[i] + L[ji]*xminus ); } L_30: if (sminus > splus) xplus = xminus; X[j] = xplus; /* *** UPDATE PARTIAL SUMS *** */ if (jm1 > 0) dv2axy( jm1, x, xplus, &L[j0 + 1], x ); } /* *** NORMALIZE X *** * */ L_60: t = ONE/dv2nrm( p, x ); for (i = 1; i <= p; i++) { X[i] *= t; } /* *** SOLVE L*Y = X AND RETURN DL7SVN = 1/TWONORM(Y) *** * */ for (j = 1; j <= p; j++) { jm1 = j - 1; j0 = j*jm1/2; jj = j0 + j; t = ZERO; if (jm1 > 0) t = dd7tpr( jm1, &L[j0 + 1], y ); Y[j] = (X[j] - t)/L[jj]; } dl7svn_v = ONE/dv2nrm( p, y ); goto L_999; L_110: dl7svn_v = ZERO; L_999: return( dl7svn_v ); /* *** LAST CARD OF DL7SVN FOLLOWS *** */ } /* end of function */ void /*FUNCTION*/ dq7apl( long nn, long n, long p, double *j, double r[], long ierr) { #define J(I_,J_) (*(j+(I_)*(nn)+(J_))) long int k, l, nl1; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const R = &r[0] - 1; /* end of OFFSET VECTORS */ /* *****PARAMETERS. */ /* .................................................................. * .................................................................. * * *****PURPOSE. * THIS SUBROUTINE APPLIES TO R THE ORTHOGONAL TRANSFORMATIONS * STORED IN J BY QRFACT * * *****PARAMETER DESCRIPTION. * ON INPUT. * * NN IS THE ROW DIMENSION OF THE MATRIX J AS DECLARED IN * THE CALLING PROGRAM DIMENSION STATEMENT * * N IS THE NUMBER OF ROWS OF J AND THE SIZE OF THE VECTOR R * * P IS THE NUMBER OF COLUMNS OF J AND THE SIZE OF SIGMA * * J CONTAINS ON AND BELOW ITS DIAGONAL THE COLUMN VECTORS * U WHICH DETERMINE THE HOUSEHOLDER TRANSFORMATIONS * IDENT - U*U.TRANSPOSE * * R IS THE RIGHT HAND SIDE VECTOR TO WHICH THE ORTHOGONAL * TRANSFORMATIONS WILL BE APPLIED * * IERR IF NON-ZERO INDICATES THAT NOT ALL THE TRANSFORMATIONS * WERE SUCCESSFULLY DETERMINED AND ONLY THE FIRST * ABS(IERR) - 1 TRANSFORMATIONS WILL BE USED * * ON OUTPUT. * * R HAS BEEN OVERWRITTEN BY ITS TRANSFORMED IMAGE * * *****APPLICATION AND USAGE RESTRICTIONS. * NONE * * *****ALGORITHM NOTES. * THE VECTORS U WHICH DETERMINE THE HOUSEHOLDER TRANSFORMATIONS * ARE NORMALIZED SO THAT THEIR 2-NORM SQUARED IS 2. THE USE OF * THESE TRANSFORMATIONS HERE IS IN THE SPIRIT OF (1). * * *****SUBROUTINES AND FUNCTIONS CALLED. * * DD7TPR - FUNCTION, RETURNS THE INNER PRODUCT OF VECTORS * * *****REFERENCES. * (1) BUSINGER, P. A., AND GOLUB, G. H. (1965), LINEAR LEAST SQUARES * SOLUTIONS BY HOUSEHOLDER TRANSFORMATIONS, NUMER. MATH. 7, * PP. 269-276. * * *****HISTORY. * DESIGNED BY DAVID M. GAY, CODED BY STEPHEN C. PETERS (WINTER 1977) * CALL ON DV2AXY SUBSTITUTED FOR DO LOOP, FALL 1983. * * *****GENERAL. * * THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH * SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS * MCS-7600324, DCR75-10143, 76-14311DSS, AND MCS76-11989. * * .................................................................. * .................................................................. * * *****LOCAL VARIABLES. */ /* *****FUNCTIONS. */ /* ------------------------------------------------------------------ * * *** BODY *** * */ k = p; if (ierr != 0) k = labs( ierr ) - 1; if (k == 0) goto L_999; for (l = 1; l <= k; l++) { nl1 = n - l + 1; dv2axy( nl1, &R[l], -dd7tpr( nl1, &J(l - 1,l - 1), &R[l] ), &J(l - 1,l - 1), &R[l] ); } L_999: return; /* *** LAST LINE OF DQ7APL FOLLOWS *** */ #undef J } /* end of function */ double /*FUNCTION*/ dv2nrm( long p, double x[]) { long int i, j; double dv2nrm_v, r, scale, t, xi; static double sqteta = 0.e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const X = &x[0] - 1; /* end of OFFSET VECTORS */ /* *** RETURN THE 2-NORM OF THE P-VECTOR X, TAKING *** * *** CARE TO AVOID THE MOST LIKELY UNDERFLOWS. *** * */ /* ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ if (p > 0) goto L_10; dv2nrm_v = ZERO; goto L_999; L_10: for (i = 1; i <= p; i++) { if (X[i] != ZERO) goto L_30; } dv2nrm_v = ZERO; goto L_999; L_30: scale = fabs( X[i] ); if (i < p) goto L_40; dv2nrm_v = scale; goto L_999; L_40: t = ONE; if (sqteta == ZERO) sqteta = dr7mdc( 2 ); /* *** SQTETA IS (SLIGHTLY LARGER THAN) THE SQUARE ROOT OF THE * *** SMALLEST POSITIVE FLOATING POINT NUMBER ON THE MACHINE. * *** THE TESTS INVOLVING SQTETA ARE DONE TO PREVENT UNDERFLOWS. * */ j = i + 1; for (i = j; i <= p; i++) { xi = fabs( X[i] ); if (xi > scale) goto L_50; r = xi/scale; if (r > sqteta) t += r*r; goto L_60; L_50: r = scale/xi; if (r <= sqteta) r = ZERO; t = ONE + t*r*r; scale = xi; L_60: ; } dv2nrm_v = scale*sqrt( t ); L_999: return( dv2nrm_v ); /* *** LAST LINE OF DV2NRM FOLLOWS *** */ } /* end of function */ double /*FUNCTION*/ dd7tpr( long p, double x[], double y[]) { long int i; double dd7tpr_v, t; static double sqteta = 0.e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const X = &x[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* *** RETURN THE INNER PRODUCT OF THE P-VECTORS X AND Y. *** * */ /* *** DR7MDC(2) RETURNS A MACHINE-DEPENDENT CONSTANT, SQTETA, WHICH * *** IS SLIGHTLY LARGER THAN THE SMALLEST POSITIVE NUMBER THAT * *** CAN BE SQUARED WITHOUT UNDERFLOWING. * */ dd7tpr_v = ZERO; if (p <= 0) goto L_999; if (sqteta == ZERO) sqteta = dr7mdc( 2 ); for (i = 1; i <= p; i++) { t = fmax( fabs( X[i] ), fabs( Y[i] ) ); if (t > ONE) goto L_10; if (t < sqteta) goto L_20; t = (X[i]/sqteta)*Y[i]; if (fabs( t ) < sqteta) goto L_20; L_10: dd7tpr_v += X[i]*Y[i]; L_20: ; } L_999: return( dd7tpr_v ); /* *** LAST LINE OF DD7TPR FOLLOWS *** */ } /* end of function */ /* PARAMETER translations */ #define JCN 66 #define JTOL 59 #define NITER 31 #define S 62 /* end of PARAMETER translations */ void /*FUNCTION*/ dd7upd( double d[], double *dr, long iv[], long liv, long lv, long n, long nd, long nn, long n2, long p, double v[]) { #define DR(I_,J_) (*(dr+(I_)*(nd)+(J_))) long int d0, i, jcn0, jcn1, jcni, jtol0, jtoli, k, sii; double t, vdfac; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const D = &d[0] - 1; long *const Iv = &iv[0] - 1; double *const V = &v[0] - 1; /* end of OFFSET VECTORS */ /* *** UPDATE SCALE VECTOR D FOR NL2IT *** * * *** PARAMETER DECLARATIONS *** * */ /* *** LOCAL VARIABLES *** * */ /* *** CONSTANTS *** * */ /* *** EXTERNAL SUBROUTINE *** * */ /* DV7SCP... SETS ALL COMPONENTS OF A VECTOR TO A SCALAR. * * *** SUBSCRIPTS FOR IV AND V *** * */ /* ------------------------------ BODY -------------------------------- * */ if (Iv[DTYPE] != 1 && Iv[NITER] > 0) goto L_999; jcn1 = Iv[JCN]; jcn0 = labs( jcn1 ) - 1; if (jcn1 < 0) goto L_10; Iv[JCN] = -jcn1; dv7scp( p, &V[jcn1], ZERO ); L_10: for (i = 1; i <= p; i++) { jcni = jcn0 + i; t = V[jcni]; for (k = 1; k <= nn; k++) { t = fmax( t, fabs( DR(i - 1,k - 1) ) ); } V[jcni] = t; } if (n2 < n) goto L_999; vdfac = V[DFAC]; jtol0 = Iv[JTOL] - 1; d0 = jtol0 + p; sii = Iv[S] - 1; for (i = 1; i <= p; i++) { sii += i; jcni = jcn0 + i; t = V[jcni]; if (V[sii] > ZERO) t = fmax( sqrt( V[sii] ), t ); jtoli = jtol0 + i; d0 += 1; if (t < V[jtoli]) t = fmax( V[d0], V[jtoli] ); D[i] = fmax( vdfac*D[i], t ); } L_999: return; /* *** LAST CARD OF DD7UPD FOLLOWS *** */ #undef DR } /* end of function */ void /*FUNCTION*/ dq7rad( long n, long nn, long p, double qtr[], LOGICAL32 qtrset, double rmat[], double *w, double y[]) { #define W(I_,J_) (*(w+(I_)*(nn)+(J_))) long int i, ii, ij, ip1, j, k, nk; double ari, qri, ri, ss, t, wi; static double big = -1.e0; static double bigrt = -1.e0; static double one = 1.e0; static double tiny = 0.e0; static double tinyrt = 0.e0; static double zero = 0.e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Qtr = &qtr[0] - 1; double *const Rmat = &rmat[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* *** ADD ROWS W TO QR FACTORIZATION WITH R MATRIX RMAT AND * *** Q**T * RESIDUAL = QTR. Y = NEW COMPONENTS OF RESIDUAL * *** CORRESPONDING TO W. QTR, Y REFERENCED ONLY IF QTRSET = .TRUE. * */ /* DIMENSION RMAT(P*(P+1)/2) */ /* *** LOCAL VARIABLES *** * */ /* ----------------------------- BODY ----------------------------------- * */ if (tiny > zero) goto L_10; tiny = dr7mdc( 1 ); big = dr7mdc( 6 ); if (tiny*big < one) tiny = one/big; L_10: k = 1; nk = n; ii = 0; for (i = 1; i <= p; i++) { ii += i; ip1 = i + 1; ij = ii + i; if (nk <= 1) t = fabs( W(i - 1,k - 1) ); if (nk > 1) t = dv2nrm( nk, &W(i - 1,k - 1) ); if (t < tiny) goto L_180; ri = Rmat[ii]; if (ri != zero) goto L_100; if (nk > 1) goto L_30; ij = ii; for (j = i; j <= p; j++) { Rmat[ij] = W(j - 1,k - 1); ij += j; } if (qtrset) Qtr[i] = Y[k]; W(i - 1,k - 1) = zero; goto L_999; L_30: wi = W(i - 1,k - 1); if (bigrt > zero) goto L_40; bigrt = dr7mdc( 5 ); tinyrt = dr7mdc( 2 ); L_40: if (t <= tinyrt) goto L_50; if (t >= bigrt) goto L_50; if (wi < zero) t = -t; wi += t; ss = sqrt( t*wi ); goto L_70; L_50: ss = sqrt( t ); if (wi < zero) goto L_60; wi += t; ss *= sqrt( wi ); goto L_70; L_60: t = -t; wi += t; ss *= sqrt( -wi ); L_70: W(i - 1,k - 1) = wi; dv7scl( nk, &W(i - 1,k - 1), one/ss, &W(i - 1,k - 1) ); Rmat[ii] = -t; if (!qtrset) goto L_80; dv2axy( nk, &Y[k], -dd7tpr( nk, &Y[k], &W(i - 1,k - 1) ), &W(i - 1,k - 1), &Y[k] ); Qtr[i] = Y[k]; L_80: if (ip1 > p) goto L_999; for (j = ip1; j <= p; j++) { dv2axy( nk, &W(j - 1,k - 1), -dd7tpr( nk, &W(j - 1,k - 1), &W(i - 1,k - 1) ), &W(i - 1,k - 1), &W(j - 1,k - 1) ); Rmat[ij] = W(j - 1,k - 1); ij += j; } if (nk <= 1) goto L_999; k += 1; nk -= 1; goto L_180; L_100: ari = fabs( ri ); if (ari > t) goto L_110; t *= sqrt( one + powi(ari/t,2) ); goto L_120; L_110: t = ari*sqrt( one + powi(t/ari,2) ); L_120: if (ri < zero) t = -t; ri += t; Rmat[ii] = -t; ss = -ri/t; if (nk <= 1) goto L_150; dv7scl( nk, &W(i - 1,k - 1), one/ri, &W(i - 1,k - 1) ); if (!qtrset) goto L_130; qri = Qtr[i]; t = ss*(qri + dd7tpr( nk, &Y[k], &W(i - 1,k - 1) )); Qtr[i] = qri + t; L_130: if (ip1 > p) goto L_999; if (qtrset) dv2axy( nk, &Y[k], t, &W(i - 1,k - 1), &Y[k] ); for (j = ip1; j <= p; j++) { ri = Rmat[ij]; t = ss*(ri + dd7tpr( nk, &W(j - 1,k - 1), &W(i - 1,k - 1) )); dv2axy( nk, &W(j - 1,k - 1), t, &W(i - 1,k - 1), &W(j - 1,k - 1) ); Rmat[ij] = ri + t; ij += j; } goto L_180; L_150: wi = W(i - 1,k - 1)/ri; W(i - 1,k - 1) = wi; if (!qtrset) goto L_160; qri = Qtr[i]; t = ss*(qri + Y[k]*wi); Qtr[i] = qri + t; L_160: if (ip1 > p) goto L_999; if (qtrset) Y[k] += t*wi; for (j = ip1; j <= p; j++) { ri = Rmat[ij]; t = ss*(ri + W(j - 1,k - 1)*wi); W(j - 1,k - 1) += t*wi; Rmat[ij] = ri + t; ij += j; } L_180: ; } L_999: return; /* *** LAST LINE OF DQ7RAD FOLLOWS *** */ #undef W } /* end of function */ #include /* PARAMETER translations */ #define DTYPE0 54 #define NEXTIV 46 #define NEXTV 47 #define OLDN 38 /* end of PARAMETER translations */ void /*FUNCTION*/ dparck( long alg, double d[], long iv[], long liv, long lv, long n, double v[]) { char which[3][5]; static char cngd[3][5], dflt[3][5], vn[34][2][5]; static byte sh[2], varnm[2]; long int _n, alg1, i, ii, iv1, j, k, l, m, miv1, miv2, ndfalt, parsv1, pu, _i, _r; static long int jlim[4], miniv[4], ndflt[4]; double vk; static double vm[34], vx[34]; static double big = 0.e0; static double machep = -1.e0; static double tiny = 1.e0; static double zero = 0.e0; static long ijmp = 33; static int _aini = 1; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const D = &d[0] - 1; long *const Iv = &iv[0] - 1; long *const Jlim = &jlim[0] - 1; long *const Miniv = &miniv[0] - 1; long *const Ndflt = &ndflt[0] - 1; byte *const Sh = &sh[0] - 1; double *const V = &v[0] - 1; byte *const Varnm = &varnm[0] - 1; double *const Vm = &vm[0] - 1; double *const Vx = &vx[0] - 1; /* end of OFFSET VECTORS */ if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */ strcpy( vn[0][0], "EPSL" ); strcpy( vn[0][1], "ON.." ); strcpy( vn[1][0], "PHMN" ); strcpy( vn[1][1], "FC.." ); strcpy( vn[2][0], "PHMX" ); strcpy( vn[2][1], "FC.." ); strcpy( vn[3][0], "DECF" ); strcpy( vn[3][1], "AC.." ); strcpy( vn[4][0], "INCF" ); strcpy( vn[4][1], "AC.." ); strcpy( vn[5][0], "RDFC" ); strcpy( vn[5][1], "MN.." ); strcpy( vn[6][0], "RDFC" ); strcpy( vn[6][1], "MX.." ); strcpy( vn[7][0], "TUNE" ); strcpy( vn[7][1], "R1.." ); strcpy( vn[8][0], "TUNE" ); strcpy( vn[8][1], "R2.." ); strcpy( vn[9][0], "TUNE" ); strcpy( vn[9][1], "R3.." ); strcpy( vn[10][0], "TUNE" ); strcpy( vn[10][1], "R4.." ); strcpy( vn[11][0], "TUNE" ); strcpy( vn[11][1], "R5.." ); strcpy( vn[12][0], "AFCT" ); strcpy( vn[12][1], "OL.." ); strcpy( vn[13][0], "RFCT" ); strcpy( vn[13][1], "OL.." ); strcpy( vn[14][0], "XCTO" ); strcpy( vn[14][1], "L..." ); strcpy( vn[15][0], "XFTO" ); strcpy( vn[15][1], "L..." ); strcpy( vn[16][0], "LMAX" ); strcpy( vn[16][1], "0..." ); strcpy( vn[17][0], "LMAX" ); strcpy( vn[17][1], "S..." ); strcpy( vn[18][0], "SCTO" ); strcpy( vn[18][1], "L..." ); strcpy( vn[19][0], "DINI" ); strcpy( vn[19][1], "T..." ); strcpy( vn[20][0], "DTIN" ); strcpy( vn[20][1], "IT.." ); strcpy( vn[21][0], "D0IN" ); strcpy( vn[21][1], "IT.." ); strcpy( vn[22][0], "DFAC" ); strcpy( vn[22][1], "...." ); strcpy( vn[23][0], "DLTF" ); strcpy( vn[23][1], "DC.." ); strcpy( vn[24][0], "DLTF" ); strcpy( vn[24][1], "DJ.." ); strcpy( vn[25][0], "DELT" ); strcpy( vn[25][1], "A0.." ); strcpy( vn[26][0], "FUZZ" ); strcpy( vn[26][1], "...." ); strcpy( vn[27][0], "RLIM" ); strcpy( vn[27][1], "IT.." ); strcpy( vn[28][0], "COSM" ); strcpy( vn[28][1], "IN.." ); strcpy( vn[29][0], "HUBE" ); strcpy( vn[29][1], "RC.." ); strcpy( vn[30][0], "RSPT" ); strcpy( vn[30][1], "OL.." ); strcpy( vn[31][0], "SIGM" ); strcpy( vn[31][1], "IN.." ); strcpy( vn[32][0], "ETA0" ); strcpy( vn[32][1], "...." ); strcpy( vn[33][0], "BIAS" ); strcpy( vn[33][1], "...." ); Vm[1] = 1.0e-3; Vm[2] = -0.99e0; Vm[3] = 1.0e-3; Vm[4] = 1.0e-2; Vm[5] = 1.2e0; Vm[6] = 1.e-2; Vm[7] = 1.2e0; Vm[8] = 0.e0; Vm[9] = 0.e0; Vm[10] = 1.e-3; Vm[11] = -1.e0; Vm[13] = 0.e0; Vm[15] = 0.e0; Vm[16] = 0.e0; Vm[19] = 0.e0; Vm[20] = -10.e0; Vm[21] = 0.e0; Vm[22] = 0.e0; Vm[23] = 0.e0; Vm[27] = 1.01e0; Vm[28] = 1.e10; Vm[30] = 0.e0; Vm[31] = 0.e0; Vm[32] = 0.e0; Vm[34] = 0.e0; Vx[1] = 0.9e0; Vx[2] = -1.e-3; Vx[3] = 1.e1; Vx[4] = 0.8e0; Vx[5] = 1.e2; Vx[6] = 0.8e0; Vx[7] = 1.e2; Vx[8] = 0.5e0; Vx[9] = 0.5e0; Vx[10] = 1.e0; Vx[11] = 1.e0; Vx[14] = 0.1e0; Vx[15] = 1.e0; Vx[16] = 1.e0; Vx[19] = 1.e0; Vx[23] = 1.e0; Vx[24] = 1.e0; Vx[25] = 1.e0; Vx[26] = 1.e0; Vx[27] = 1.e10; Vx[29] = 1.e0; Vx[31] = 1.e0; Vx[32] = 1.e0; Vx[33] = 1.e0; Vx[34] = 1.e0; Varnm[1] = 'P'; Varnm[2] = 'P'; Sh[1] = 'S'; Sh[2] = 'H'; strcpy( cngd[0], "---C" ); strcpy( cngd[1], "HANG" ); strcpy( cngd[2], "ED V" ); strcpy( dflt[0], "NOND" ); strcpy( dflt[1], "EFAU" ); strcpy( dflt[2], "LT V" ); Jlim[1] = 0; Jlim[2] = 24; Jlim[3] = 0; Jlim[4] = 24; Ndflt[1] = 32; Ndflt[2] = 25; Ndflt[3] = 32; Ndflt[4] = 25; Miniv[1] = 82; Miniv[2] = 59; Miniv[3] = 103; Miniv[4] = 103; _aini = 0; } /* *** CHECK ***SOL (VERSION 2.3) PARAMETERS, PRINT CHANGED VALUES *** * * *** ALG = 1 FOR REGRESSION, ALG = 2 FOR GENERAL UNCONSTRAINED OPT. * */ /* DIVSET -- SUPPLIES DEFAULT VALUES TO BOTH IV AND V. * DR7MDC -- RETURNS MACHINE-DEPENDENT CONSTANTS. * DV7CPY -- COPIES ONE VECTOR TO ANOTHER. * DV7DFL -- SUPPLIES DEFAULT PARAMETER VALUES TO V ALONE. * * *** LOCAL VARIABLES *** * */ /* *** IV AND V SUBSCRIPTS *** * */ /*............................... BODY ................................ * */ /*++(~.C.) Default UNITNO='(PU,' *++(.C.) Default UNITNO='(*,' *++ Replace "(*, " = UNITNO */ pu = 0; if (PRUNIT <= liv) pu = Iv[PRUNIT]; if (ALGSAV > liv) goto L_20; if (alg == Iv[ALGSAV]) goto L_20; if (pu != 0) { printf("\n THE FIRST PARAMETER TO DIVSET SHOULD BE%3ld RATHER THAN%3ld\n", alg, Iv[ALGSAV]); } Iv[1] = 67; goto L_999; L_20: if (alg < 1 || alg > 4) goto L_340; miv1 = Miniv[alg]; if (Iv[1] == 15) goto L_360; alg1 = ((alg - 1)%2) + 1; if (Iv[1] == 0) divset( alg, iv, liv, lv, v ); iv1 = Iv[1]; if (iv1 != 13 && iv1 != 12) goto L_30; if (PERM <= liv) miv1 = max( miv1, Iv[PERM] - 1 ); if (IVNEED <= liv) miv2 = miv1 + max( Iv[IVNEED], 0 ); if (LASTIV <= liv) Iv[LASTIV] = miv2; if (liv < miv1) goto L_300; Iv[IVNEED] = 0; Iv[LASTV] = max( Iv[VNEED], 0 ) + Iv[LMAT] - 1; Iv[VNEED] = 0; if (liv < miv2) goto L_300; if (lv < Iv[LASTV]) goto L_320; L_30: if (iv1 < 12 || iv1 > 14) goto L_60; if (n >= 1) goto L_50; Iv[1] = 81; if (pu == 0) goto L_999; printf("\n /// BAD%c =%5ld\n", Varnm[alg1], n); goto L_999; L_50: if (iv1 != 14) Iv[NEXTIV] = Iv[PERM]; if (iv1 != 14) Iv[NEXTV] = Iv[LMAT]; if (iv1 == 13) goto L_999; k = Iv[PARSAV] - EPSLON; dv7dfl( alg1, lv - k, &V[k + 1] ); Iv[DTYPE0] = 2 - alg1; Iv[OLDN] = n; strcpy( which[0], dflt[0] ); strcpy( which[1], dflt[1] ); strcpy( which[2], dflt[2] ); goto L_110; L_60: if (n == Iv[OLDN]) goto L_80; Iv[1] = 17; if (pu == 0) goto L_999; printf("\n /// %c CHANGED FROM %5ld TO %5ld\n", Varnm[alg1], Iv[OLDN], n); goto L_999; L_80: if (iv1 <= 11 && iv1 >= 1) goto L_100; Iv[1] = 80; if (pu != 0) { printf("\n /// IV(1) =%5ld SHOULD BE BETWEEN 0 AND 14.\n", iv1); } goto L_999; L_100: strcpy( which[0], cngd[0] ); strcpy( which[1], cngd[1] ); strcpy( which[2], cngd[2] ); L_110: if (iv1 == 14) iv1 = 12; if (big > tiny) goto L_120; tiny = dr7mdc( 1 ); machep = dr7mdc( 3 ); big = dr7mdc( 6 ); Vm[12] = machep; Vx[12] = big; Vx[13] = big; Vm[14] = machep; Vm[17] = tiny; Vx[17] = big; Vm[18] = tiny; Vx[18] = big; Vx[20] = big; Vx[21] = big; Vx[22] = big; Vm[24] = machep; Vm[25] = machep; Vm[26] = machep; Vx[28] = dr7mdc( 5 ); Vm[29] = machep; Vx[30] = big; Vm[33] = machep; L_120: m = 0; i = 1; j = Jlim[alg1]; k = EPSLON; ndfalt = Ndflt[alg1]; for (l = 1; l <= ndfalt; l++) { vk = V[k]; if (vk >= Vm[i] && vk <= Vx[i]) goto L_140; m = k; if (pu != 0) { printf("\n /// %4.4s%4.4s.. V(%2ld) =%11.3g should be between%11.3g and%11.3g\n", vn[i - 1][0], vn[i - 1][1], k, vk, Vm[i], Vx[i]); } L_140: k += 1; i += 1; if (i == j) i = ijmp; } if (Iv[NVDFLT] == ndfalt) goto L_170; Iv[1] = 51; if (pu == 0) goto L_999; printf("\n IV(NVDFLT) =%5ld RATHER THAN %5ld\n", Iv[NVDFLT], ndfalt); goto L_999; L_170: if ((Iv[DTYPE] > 0 || V[DINIT] > zero) && iv1 == 12) goto L_200; for (i = 1; i <= n; i++) { if (D[i] > zero) goto L_190; m = 18; if (pu != 0) { printf("\n /// D(%3ld) =%11.3g SHOULD BE POSITIVE\n", i, D[i]); } L_190: ; } L_200: if (m == 0) goto L_210; Iv[1] = m; goto L_999; L_210: if (pu == 0 || Iv[PARPRT] == 0) goto L_999; if (iv1 != 12 || Iv[INITS] == alg1 - 1) goto L_230; m = 1; printf("\n NONDEFAULT VALUES....\n INIT%c..... IV(25) =%3ld\n", Sh[alg1], Iv[INITS]); L_230: if (Iv[DTYPE] == Iv[DTYPE0]) goto L_250; if (m == 0) { printf("\n "); for(_n=0L; _n < 3; _n++) printf("%4.4s", (char*)which+_n*5); printf("\n"); } m = 1; printf(" DTYPE..... IV(16) =%3ld\n", Iv[DTYPE]); L_250: i = 1; j = Jlim[alg1]; k = EPSLON; l = Iv[PARSAV]; ndfalt = Ndflt[alg1]; for (ii = 1; ii <= ndfalt; ii++) { if (V[k] == V[l]) goto L_280; if (m == 0) { printf("\n "); for(_n=0L; _n < 3; _n++) printf("%4.4s", (char*)which+_n*5); printf("\n"); } m = 1; printf(" %4.4s%4.4s.. V(%2ld) =%15.7g\n", vn[i - 1][0] , vn[i - 1][1], k, V[k]); L_280: k += 1; l += 1; i += 1; if (i == j) i = ijmp; } Iv[DTYPE0] = Iv[DTYPE]; parsv1 = Iv[PARSAV]; dv7cpy( Iv[NVDFLT], &V[parsv1], &V[EPSLON] ); goto L_999; L_300: Iv[1] = 15; if (pu == 0) goto L_999; printf("\n /// LIV =%5ld MUST BE AT LEAST%5ld\n", liv, miv2); if (liv < miv1) goto L_999; if (lv < Iv[LASTV]) goto L_320; goto L_999; L_320: Iv[1] = 16; if (pu != 0) { printf("\n /// LV =%5ld MUST BE AT LEAST%5ld\n", lv, Iv[LASTV]); } goto L_999; L_340: Iv[1] = 67; if (pu != 0) { printf("\n /// ALG =%5ld MUST BE 1 2, 3, OR 4\n", alg); } goto L_999; L_360: if (pu != 0) { printf("\n /// LIV =%5ld MUST BE AT LEAST%5ld TO COMPUTE TRUE MIN. LIV AND MIN. LV\n", liv, miv1); } if (LASTIV <= liv) Iv[LASTIV] = miv1; if (LASTV <= liv) Iv[LASTV] = 0; L_999: return; /* *** LAST LINE OF DPARCK FOLLOWS *** */ } /* end of function */ void /*FUNCTION*/ ds7lvm( long p, double y[], double ss[], double x[]) { long int i, im1, j, k; double xi; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Ss = &ss[0] - 1; double *const X = &x[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* *** SET Y = SS * X, SS = P X P SYMMETRIC MATRIX. *** * *** LOWER TRIANGLE OF SS STORED ROWWISE. *** * * *** PARAMETER DECLARATIONS *** * */ /* DIMENSION SS(P*(P+1)/2) * * *** LOCAL VARIABLES *** * */ /* *** EXTERNAL FUNCTION *** * */ /* ---------------------------------------------------------------------- * */ j = 1; for (i = 1; i <= p; i++) { Y[i] = dd7tpr( i, &Ss[j], x ); j += i; } if (p <= 1) goto L_999; j = 1; for (i = 2; i <= p; i++) { xi = X[i]; im1 = i - 1; j += 1; for (k = 1; k <= im1; k++) { Y[k] += Ss[j]*xi; j += 1; } } L_999: return; /* *** LAST CARD OF DS7LVM FOLLOWS *** */ } /* end of function */ void /*FUNCTION*/ ds7lup( double a[], double cosmin, long p, double size, double step[], double u[], double w[], double wchmtd[], double *wscale, double y[]) { long int i, j, k; double denmin, sdotwm, t, ui, wi; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const A = &a[0] - 1; double *const Step = &step[0] - 1; double *const U = &u[0] - 1; double *const W = &w[0] - 1; double *const Wchmtd = &wchmtd[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* *** UPDATE SYMMETRIC A SO THAT A * STEP = Y *** * *** (LOWER TRIANGLE OF A STORED ROWWISE *** * * *** PARAMETER DECLARATIONS *** * */ /* DIMENSION A(P*(P+1)/2) * * *** LOCAL VARIABLES *** * */ /* *** CONSTANTS *** */ /* *** EXTERNAL FUNCTIONS AND SUBROUTINES *** * */ /* ------------------------------------------------------------------ * */ sdotwm = dd7tpr( p, step, wchmtd ); denmin = cosmin*dv2nrm( p, step )*dv2nrm( p, wchmtd ); *wscale = ONE; if (denmin != ZERO) *wscale = fmin( ONE, fabs( sdotwm/denmin ) ); t = ZERO; if (sdotwm != ZERO) t = *wscale/sdotwm; for (i = 1; i <= p; i++) { W[i] = t*Wchmtd[i]; } ds7lvm( p, u, a, step ); t = HALF*(size*dd7tpr( p, step, u ) - dd7tpr( p, step, y )); for (i = 1; i <= p; i++) { U[i] = t*W[i] + Y[i] - size*U[i]; } /* *** SET A = A + U*(W**T) + W*(U**T) *** * */ k = 1; for (i = 1; i <= p; i++) { ui = U[i]; wi = W[i]; for (j = 1; j <= i; j++) { A[k] = size*A[k] + ui*W[j] + wi*U[j]; k += 1; } } /* *** LAST CARD OF DS7LUP FOLLOWS *** */ return; } /* end of function */ /* PARAMETER translations */ #undef DFAC #define DFAC 256.e0 #define DGNORM 1 #define DST0 3 #define DSTNRM 2 #define EIGHT 8.e0 #define GTSTEP 4 #define NEGONE (-1.e0) #define NREDUC 6 #define P001 1.e-3 #define PREDUC 7 #define RAD0 9 #define RADIUS 8 #define STPPAR 5 #define TTOL 2.5e0 /* end of PARAMETER translations */ void /*FUNCTION*/ dl7mst( double d[], double g[], long ierr, long ipivot[], long *ka, long p, double qtr[], double r[], double step[], double v[], double w[]) { long int dstsav, i, i1, ip1, j1, k, kalim, l, lk0, phipin, pp1o2, res, res0, rmat, rmat0, uk0; double a, adi, alphak, b, d1, d2, dfacsq, dst, dtol, lk, oldphi, phi, phimax, phimin, psifac, rad, si, sj, sqrtak, t, twopsi, uk, wl; static double big = 0.e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const D = &d[0] - 1; double *const G = &g[0] - 1; long *const Ipivot = &ipivot[0] - 1; double *const Qtr = &qtr[0] - 1; double *const R = &r[0] - 1; double *const Step = &step[0] - 1; double *const V = &v[0] - 1; double *const W = &w[0] - 1; /* end of OFFSET VECTORS */ /* *** COMPUTE LEVENBERG-MARQUARDT STEP USING MORE-HEBDEN TECHNIQUE ** * *** NL2SOL VERSION 2.2. *** * * *** PARAMETER DECLARATIONS *** * */ /* DIMENSION R(*), W(P*(P+5)/2 + 4) * * ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * * *** PURPOSE *** * * GIVEN THE R MATRIX FROM THE QR DECOMPOSITION OF A JACOBIAN * MATRIX, J, AS WELL AS Q-TRANSPOSE TIMES THE CORRESPONDING * RESIDUAL VECTOR, RESID, THIS SUBROUTINE COMPUTES A LEVENBERG- * MARQUARDT STEP OF APPROXIMATE LENGTH V(RADIUS) BY THE MORE- * TECHNIQUE. * * *** PARAMETER DESCRIPTION *** * * D (IN) = THE SCALE VECTOR. * G (IN) = THE GRADIENT VECTOR (J**T)*R. * IERR (I/O) = RETURN CODE FROM QRFACT OR DQ7RGS -- 0 MEANS R HAS * FULL RANK. * IPIVOT (I/O) = PERMUTATION ARRAY FROM QRFACT OR DQ7RGS, WHICH COMPUTE * QR DECOMPOSITIONS WITH COLUMN PIVOTING. * KA (I/O). KA .LT. 0 ON INPUT MEANS THIS IS THE FIRST CALL ON * DL7MST FOR THE CURRENT R AND QTR. ON OUTPUT KA CON- * TAINS THE NUMBER OF HEBDEN ITERATIONS NEEDED TO DETERMINE * STEP. KA = 0 MEANS A GAUSS-NEWTON STEP. * P (IN) = NUMBER OF PARAMETERS. * QTR (IN) = (Q**T)*RESID = Q-TRANSPOSE TIMES THE RESIDUAL VECTOR. * R (IN) = THE R MATRIX, STORED COMPACTLY BY COLUMNS. * STEP (OUT) = THE LEVENBERG-MARQUARDT STEP COMPUTED. * V (I/O) CONTAINS VARIOUS CONSTANTS AND VARIABLES DESCRIBED BELOW. * W (I/O) = WORKSPACE OF LENGTH P*(P+5)/2 + 4. * * *** ENTRIES IN V *** * * V(DGNORM) (I/O) = 2-NORM OF (D**-1)*G. * V(DSTNRM) (I/O) = 2-NORM OF D*STEP. * V(DST0) (I/O) = 2-NORM OF GAUSS-NEWTON STEP (FOR NONSING. J). * V(EPSLON) (IN) = MAX. REL. ERROR ALLOWED IN TWONORM(R)**2 MINUS * TWONORM(R - J*STEP)**2. (SEE ALGORITHM NOTES BELOW.) * V(GTSTEP) (OUT) = INNER PRODUCT BETWEEN G AND STEP. * V(NREDUC) (OUT) = HALF THE REDUCTION IN THE SUM OF SQUARES PREDICTED * FOR A GAUSS-NEWTON STEP. * V(PHMNFC) (IN) = TOL. (TOGETHER WITH V(PHMXFC)) FOR ACCEPTING STEP * (MORE*S SIGMA). THE ERROR V(DSTNRM) - V(RADIUS) MUST LIE * BETWEEN V(PHMNFC)*V(RADIUS) AND V(PHMXFC)*V(RADIUS). * V(PHMXFC) (IN) (SEE V(PHMNFC).) * V(PREDUC) (OUT) = HALF THE REDUCTION IN THE SUM OF SQUARES PREDICTED * BY THE STEP RETURNED. * V(RADIUS) (IN) = RADIUS OF CURRENT (SCALED) TRUST REGION. * V(RAD0) (I/O) = VALUE OF V(RADIUS) FROM PREVIOUS CALL. * V(STPPAR) (I/O) = MARQUARDT PARAMETER (OR ITS NEGATIVE IF THE SPECIAL * CASE MENTIONED BELOW IN THE ALGORITHM NOTES OCCURS). * * NOTE -- SEE DATA STATEMENT BELOW FOR VALUES OF ABOVE SUBSCRIPTS. * * *** USAGE NOTES *** * * IF IT IS DESIRED TO RECOMPUTE STEP USING A DIFFERENT VALUE OF * V(RADIUS), THEN THIS ROUTINE MAY BE RESTARTED BY CALLING IT * WITH ALL PARAMETERS UNCHANGED EXCEPT V(RADIUS). (THIS EXPLAINS * WHY MANY PARAMETERS ARE LISTED AS I/O). ON AN INTIIAL CALL (ONE * WITH KA = -1), THE CALLER NEED ONLY HAVE INITIALIZED D, G, KA, P, * QTR, R, V(EPSLON), V(PHMNFC), V(PHMXFC), V(RADIUS), AND V(RAD0). * * *** APPLICATION AND USAGE RESTRICTIONS *** * * THIS ROUTINE IS CALLED AS PART OF THE NL2SOL (NONLINEAR LEAST- * SQUARES) PACKAGE (REF. 1). * * *** ALGORITHM NOTES *** * * THIS CODE IMPLEMENTS THE STEP COMPUTATION SCHEME DESCRIBED IN * REFS. 2 AND 4. FAST GIVENS TRANSFORMATIONS (SEE REF. 3, PP. 60- * 62) ARE USED TO COMPUTE STEP WITH A NONZERO MARQUARDT PARAMETER. * A SPECIAL CASE OCCURS IF J IS (NEARLY) SINGULAR AND V(RADIUS) * IS SUFFICIENTLY LARGE. IN THIS CASE THE STEP RETURNED IS SUCH * THAT TWONORM(R)**2 - TWONORM(R - J*STEP)**2 DIFFERS FROM ITS * OPTIMAL VALUE BY LESS THAN V(EPSLON) TIMES THIS OPTIMAL VALUE, * WHERE J AND R DENOTE THE ORIGINAL JACOBIAN AND RESIDUAL. (SEE * REF. 2 FOR MORE DETAILS.) * * *** FUNCTIONS AND SUBROUTINES CALLED *** * * DD7TPR - RETURNS INNER PRODUCT OF TWO VECTORS. * DL7ITV - APPLY INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX. * DL7IVM - APPLY INVERSE OF COMPACT LOWER TRIANG. MATRIX. * DV7CPY - COPIES ONE VECTOR TO ANOTHER. * DV2NRM - RETURNS 2-NORM OF A VECTOR. * * *** REFERENCES *** * * 1. DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1981), AN ADAPTIVE * NONLINEAR LEAST-SQUARES ALGORITHM, ACM TRANS. MATH. * SOFTWARE, VOL. 7, NO. 3. * 2. GAY, D.M. (1981), COMPUTING OPTIMAL LOCALLY CONSTRAINED STEPS, * SIAM J. SCI. STATIST. COMPUTING, VOL. 2, NO. 2, PP. * 186-197. * 3. LAWSON, C.L., AND HANSON, R.J. (1974), SOLVING LEAST SQUARES * PROBLEMS, PRENTICE-HALL, ENGLEWOOD CLIFFS, N.J. * 4. MORE, J.J. (1978), THE LEVENBERG-MARQUARDT ALGORITHM, IMPLEMEN- * TATION AND THEORY, PP.105-116 OF SPRINGER LECTURE NOTES */ /* IN MATHEMATICS NO. 630, EDITED BY G.A. WATSON, SPRINGER- * VERLAG, BERLIN AND NEW YORK. * * *** GENERAL *** * * CODED BY DAVID M. GAY. * THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH * SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS * MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND * MCS-7906671. * * ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * * *** LOCAL VARIABLES *** * */ /* *** CONSTANTS *** */ /* *** EXTERNAL FUNCTIONS AND SUBROUTINES *** * */ /* *** SUBSCRIPTS FOR V *** * */ /* *** BODY *** * * *** FOR USE IN RECOMPUTING STEP, THE FINAL VALUES OF LK AND UK, * *** THE INVERSE DERIVATIVE OF MORE*S PHI AT 0 (FOR NONSING. J) * *** AND THE VALUE RETURNED AS V(DSTNRM) ARE STORED AT W(LK0), * *** W(UK0), W(PHIPIN), AND W(DSTSAV) RESPECTIVELY. */ lk0 = p + 1; phipin = lk0 + 1; uk0 = phipin + 1; dstsav = uk0 + 1; rmat0 = dstsav; /* *** A COPY OF THE R-MATRIX FROM THE QR DECOMPOSITION OF J IS * *** STORED IN W STARTING AT W(RMAT), AND A COPY OF THE RESIDUAL * *** VECTOR IS STORED IN W STARTING AT W(RES). THE LOOPS BELOW * *** THAT UPDATE THE QR DECOMP. FOR A NONZERO MARQUARDT PARAMETER * *** WORK ON THESE COPIES. */ rmat = rmat0 + 1; pp1o2 = p*(p + 1)/2; res0 = pp1o2 + rmat0; res = res0 + 1; rad = V[RADIUS]; if (rad > ZERO) psifac = V[EPSLON]/((EIGHT*(V[PHMNFC] + ONE) + THREE)*SQ(rad)); if (big <= ZERO) big = dr7mdc( 6 ); phimax = V[PHMXFC]*rad; phimin = V[PHMNFC]*rad; /* *** DTOL, DFAC, AND DFACSQ ARE USED IN RESCALING THE FAST GIVENS * *** REPRESENTATION OF THE UPDATED QR DECOMPOSITION. */ dtol = ONE/DFAC; dfacsq = DFAC*DFAC; /* *** OLDPHI IS USED TO DETECT LIMITS OF NUMERICAL ACCURACY. IF * *** WE RECOMPUTE STEP AND IT DOES NOT CHANGE, THEN WE ACCEPT IT. */ oldphi = ZERO; lk = ZERO; uk = ZERO; kalim = *ka + 12; /* *** START OR RESTART, DEPENDING ON KA *** * */ switch (IARITHIF(*ka)) { case -1: goto L_10; case 0: goto L_20; case 1: goto L_370; } /* *** FRESH START -- COMPUTE V(NREDUC) *** * */ L_10: *ka = 0; kalim = 12; k = p; if (ierr != 0) k = labs( ierr ) - 1; V[NREDUC] = HALF*dd7tpr( k, qtr, qtr ); /* *** SET UP TO TRY INITIAL GAUSS-NEWTON STEP *** * */ L_20: V[DST0] = NEGONE; if (ierr != 0) goto L_90; t = dl7svn( p, r, step, &W[res] ); if (t >= ONE) goto L_30; if (dv2nrm( p, qtr ) >= big*t) goto L_90; /* *** COMPUTE GAUSS-NEWTON STEP *** * * *** NOTE -- THE R-MATRIX IS STORED COMPACTLY BY COLUMNS IN * *** R(1), R(2), R(3), ... IT IS THE TRANSPOSE OF A * *** LOWER TRIANGULAR MATRIX STORED COMPACTLY BY ROWS, AND WE * *** TREAT IT AS SUCH WHEN USING DL7ITV AND DL7IVM. */ L_30: dl7itv( p, w, r, qtr ); /* *** TEMPORARILY STORE PERMUTED -D*STEP IN STEP. */ for (i = 1; i <= p; i++) { j1 = Ipivot[i]; Step[i] = D[j1]*W[i]; } dst = dv2nrm( p, step ); V[DST0] = dst; phi = dst - rad; if (phi <= phimax) goto L_410; /* *** IF THIS IS A RESTART, GO TO 110 *** */ if (*ka > 0) goto L_110; /* *** GAUSS-NEWTON STEP WAS UNACCEPTABLE. COMPUTE L0 *** * */ for (i = 1; i <= p; i++) { j1 = Ipivot[i]; Step[i] = D[j1]*(Step[i]/dst); } dl7ivm( p, step, r, step ); t = ONE/dv2nrm( p, step ); W[phipin] = (t/rad)*t; lk = phi*W[phipin]; /* *** COMPUTE U0 *** * */ L_90: for (i = 1; i <= p; i++) { W[i] = G[i]/D[i]; } V[DGNORM] = dv2nrm( p, w ); uk = V[DGNORM]/rad; if (uk <= ZERO) goto L_390; /* *** ALPHAK WILL BE USED AS THE CURRENT MARQUARDT PARAMETER. WE * *** USE MORE*S SCHEME FOR INITIALIZING IT. * */ alphak = fabs( V[STPPAR] )*V[RAD0]/rad; alphak = fmin( uk, fmax( alphak, lk ) ); /* *** TOP OF LOOP -- INCREMENT KA, COPY R TO RMAT, QTR TO RES *** * */ L_110: *ka += 1; dv7cpy( pp1o2, &W[rmat], r ); dv7cpy( p, &W[res], qtr ); /* *** SAFEGUARD ALPHAK AND INITIALIZE FAST GIVENS SCALE VECTOR. *** * */ if ((alphak <= ZERO || alphak < lk) || alphak >= uk) alphak = uk*fmax( P001, sqrt( lk/uk ) ); if (alphak <= ZERO) alphak = HALF*uk; sqrtak = sqrt( alphak ); for (i = 1; i <= p; i++) { W[i] = ONE; } /* *** ADD ALPHAK*D AND UPDATE QR DECOMP. USING FAST GIVENS TRANS. *** * */ for (i = 1; i <= p; i++) { /* *** GENERATE, APPLY 1ST GIVENS TRANS. FOR ROW I OF ALPHAK*D. * *** (USE STEP TO STORE TEMPORARY ROW) *** */ l = i*(i + 1)/2 + rmat0; wl = W[l]; d2 = ONE; d1 = W[i]; j1 = Ipivot[i]; adi = sqrtak*D[j1]; if (adi >= fabs( wl )) goto L_150; L_130: a = adi/wl; b = d2*a/d1; t = a*b + ONE; if (t > TTOL) goto L_150; W[i] = d1/t; d2 /= t; W[l] = t*wl; a = -a; for (j1 = i; j1 <= p; j1++) { l += j1; Step[j1] = a*W[l]; } goto L_170; L_150: b = wl/adi; a = d1*b/d2; t = a*b + ONE; if (t > TTOL) goto L_130; W[i] = d2/t; d2 = d1/t; W[l] = t*adi; for (j1 = i; j1 <= p; j1++) { l += j1; wl = W[l]; Step[j1] = -wl; W[l] = a*wl; } L_170: if (i == p) goto L_280; /* *** NOW USE GIVENS TRANS. TO ZERO ELEMENTS OF TEMP. ROW *** * */ ip1 = i + 1; for (i1 = ip1; i1 <= p; i1++) { l = i1*(i1 + 1)/2 + rmat0; wl = W[l]; si = Step[i1 - 1]; d1 = W[i1]; /* *** RESCALE ROW I1 IF NECESSARY *** * */ if (d1 >= dtol) goto L_190; d1 *= dfacsq; wl /= DFAC; k = l; for (j1 = i1; j1 <= p; j1++) { k += j1; W[k] /= DFAC; } /* *** USE GIVENS TRANS. TO ZERO NEXT ELEMENT OF TEMP. ROW * */ L_190: if (fabs( si ) > fabs( wl )) goto L_220; if (si == ZERO) goto L_260; L_200: a = si/wl; b = d2*a/d1; t = a*b + ONE; if (t > TTOL) goto L_220; W[l] = t*wl; W[i1] = d1/t; d2 /= t; for (j1 = i1; j1 <= p; j1++) { l += j1; wl = W[l]; sj = Step[j1]; W[l] = wl + b*sj; Step[j1] = sj - a*wl; } goto L_240; L_220: b = wl/si; a = d1*b/d2; t = a*b + ONE; if (t > TTOL) goto L_200; W[i1] = d2/t; d2 = d1/t; W[l] = t*si; for (j1 = i1; j1 <= p; j1++) { l += j1; wl = W[l]; sj = Step[j1]; W[l] = a*wl + sj; Step[j1] = b*sj - wl; } /* *** RESCALE TEMP. ROW IF NECESSARY *** * */ L_240: if (d2 >= dtol) goto L_260; d2 *= dfacsq; for (k = i1; k <= p; k++) { Step[k] /= DFAC; } L_260: ; } } /* *** COMPUTE STEP *** * */ L_280: dl7itv( p, &W[res], &W[rmat], &W[res] ); /* *** RECOVER STEP AND STORE PERMUTED -D*STEP AT W(RES) *** */ for (i = 1; i <= p; i++) { j1 = Ipivot[i]; k = res0 + i; t = W[k]; Step[j1] = -t; W[k] = t*D[j1]; } dst = dv2nrm( p, &W[res] ); phi = dst - rad; if (phi <= phimax && phi >= phimin) goto L_430; if (oldphi == phi) goto L_430; oldphi = phi; /* *** CHECK FOR (AND HANDLE) SPECIAL CASE *** * */ if (phi > ZERO) goto L_310; if (*ka >= kalim) goto L_430; twopsi = alphak*dst*dst - dd7tpr( p, step, g ); if (alphak >= twopsi*psifac) goto L_310; V[STPPAR] = -alphak; goto L_440; /* *** UNACCEPTABLE STEP -- UPDATE LK, UK, ALPHAK, AND TRY AGAIN *** * */ L_300: if (phi < ZERO) uk = fmin( uk, alphak ); goto L_320; L_310: if (phi < ZERO) uk = alphak; L_320: for (i = 1; i <= p; i++) { j1 = Ipivot[i]; k = res0 + i; Step[i] = D[j1]*(W[k]/dst); } dl7ivm( p, step, &W[rmat], step ); for (i = 1; i <= p; i++) { Step[i] /= sqrt( W[i] ); } t = ONE/dv2nrm( p, step ); alphak += t*phi*t/rad; lk = fmax( lk, alphak ); alphak = lk; goto L_110; /* *** RESTART *** * */ L_370: lk = W[lk0]; uk = W[uk0]; if (V[DST0] > ZERO && V[DST0] - rad <= phimax) goto L_20; alphak = fabs( V[STPPAR] ); dst = W[dstsav]; phi = dst - rad; t = V[DGNORM]/rad; if (rad > V[RAD0]) goto L_380; /* *** SMALLER RADIUS *** */ uk = t; if (alphak <= ZERO) lk = ZERO; if (V[DST0] > ZERO) lk = fmax( lk, (V[DST0] - rad)*W[phipin] ); goto L_300; /* *** BIGGER RADIUS *** */ L_380: if (alphak <= ZERO || uk > t) uk = t; lk = ZERO; if (V[DST0] > ZERO) lk = fmax( lk, (V[DST0] - rad)*W[phipin] ); goto L_300; /* *** SPECIAL CASE -- RAD .LE. 0 OR (G = 0 AND J IS SINGULAR) *** * */ L_390: V[STPPAR] = ZERO; dst = ZERO; lk = ZERO; uk = ZERO; V[GTSTEP] = ZERO; V[PREDUC] = ZERO; for (i = 1; i <= p; i++) { Step[i] = ZERO; } goto L_450; /* *** ACCEPTABLE GAUSS-NEWTON STEP -- RECOVER STEP FROM W *** * */ L_410: alphak = ZERO; for (i = 1; i <= p; i++) { j1 = Ipivot[i]; Step[j1] = -W[i]; } /* *** SAVE VALUES FOR USE IN A POSSIBLE RESTART *** * */ L_430: V[STPPAR] = alphak; L_440: V[GTSTEP] = fmin( dd7tpr( p, step, g ), ZERO ); V[PREDUC] = HALF*(alphak*dst*dst - V[GTSTEP]); L_450: V[DSTNRM] = dst; W[dstsav] = dst; W[lk0] = lk; W[uk0] = uk; V[RAD0] = rad; /* *** LAST CARD OF DL7MST FOLLOWS *** */ return; } /* end of function */ /* PARAMETER translations */ #define EPSFAC 50.0e0 #define FOUR 4.0e0 #define KAPPA 2.0e0 #undef NEGONE #define NEGONE (-1.0e0) #undef ONE #define ONE 1.0e0 #undef P001 #define P001 1.0e-3 #define SIX 6.0e0 #undef THREE #define THREE 3.0e0 #define TWO 2.0e0 #undef ZERO #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ dg7qts( double d[], double dig[], double dihdi[], long *ka, double l[], long p, double step[], double v[], double w[]) { LOGICAL32 restrt; long int dggdmx, diag, diag0, dstsav, emax, emin, i, im1, inc, irc, j, k, k1, kalim, kamin, lk0, phipin, q, q0, uk0, x; double aki, akk, alphak, delta, dst, eps, gtsta, lk, oldphi, phi, phimax, phimin, psifac, rad, radsq, root, si, sk, sw, t, t1, t2, twopsi, uk, wi; static double big = 0.e0; static double dgxfac = 0.e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const D = &d[0] - 1; double *const Dig = &dig[0] - 1; double *const Dihdi = &dihdi[0] - 1; double *const L = &l[0] - 1; double *const Step = &step[0] - 1; double *const V = &v[0] - 1; double *const W = &w[0] - 1; /* end of OFFSET VECTORS */ /* *** COMPUTE GOLDFELD-QUANDT-TROTTER STEP BY MORE-HEBDEN TECHNIQUE *** * *** (NL2SOL VERSION 2.2), MODIFIED A LA MORE AND SORENSEN *** * * *** PARAMETER DECLARATIONS *** * */ /* DIMENSION DIHDI(P*(P+1)/2), L(P*(P+1)/2), W(4*P+7) * * ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * * *** PURPOSE *** * * GIVEN THE (COMPACTLY STORED) LOWER TRIANGLE OF A SCALED * HESSIAN (APPROXIMATION) AND A NONZERO SCALED GRADIENT VECTOR, * THIS SUBROUTINE COMPUTES A GOLDFELD-QUANDT-TROTTER STEP OF * APPROXIMATE LENGTH V(RADIUS) BY THE MORE-HEBDEN TECHNIQUE. IN * OTHER WORDS, STEP IS COMPUTED TO (APPROXIMATELY) MINIMIZE * PSI(STEP) = (G**T)*STEP + 0.5*(STEP**T)*H*STEP SUCH THAT THE * 2-NORM OF D*STEP IS AT MOST (APPROXIMATELY) V(RADIUS), WHERE * G IS THE GRADIENT, H IS THE HESSIAN, AND D IS A DIAGONAL * SCALE MATRIX WHOSE DIAGONAL IS STORED IN THE PARAMETER D. * (DG7QTS ASSUMES DIG = D**-1 * G AND DIHDI = D**-1 * H * D**-1.) * * *** PARAMETER DESCRIPTION *** * * D (IN) = THE SCALE VECTOR, I.E. THE DIAGONAL OF THE SCALE * MATRIX D MENTIONED ABOVE UNDER PURPOSE. * DIG (IN) = THE SCALED GRADIENT VECTOR, D**-1 * G. IF G = 0, THEN * STEP = 0 AND V(STPPAR) = 0 ARE RETURNED. * DIHDI (IN) = LOWER TRIANGLE OF THE SCALED HESSIAN (APPROXIMATION), * I.E., D**-1 * H * D**-1, STORED COMPACTLY BY ROWS., I.E., * IN THE ORDER (1,1), (2,1), (2,2), (3,1), (3,2), ETC. * KA (I/O) = THE NUMBER OF HEBDEN ITERATIONS (SO FAR) TAKEN TO DETER- * MINE STEP. KA .LT. 0 ON INPUT MEANS THIS IS THE FIRST * ATTEMPT TO DETERMINE STEP (FOR THE PRESENT DIG AND DIHDI) * -- KA IS INITIALIZED TO 0 IN THIS CASE. OUTPUT WITH * KA = 0 (OR V(STPPAR) = 0) MEANS STEP = -(H**-1)*G. * L (I/O) = WORKSPACE OF LENGTH P*(P+1)/2 FOR CHOLESKY FACTORS. * P (IN) = NUMBER OF PARAMETERS -- THE HESSIAN IS A P X P MATRIX. * STEP (I/O) = THE STEP COMPUTED. * V (I/O) CONTAINS VARIOUS CONSTANTS AND VARIABLES DESCRIBED BELOW. * W (I/O) = WORKSPACE OF LENGTH 4*P + 6. * * *** ENTRIES IN V *** * * V(DGNORM) (I/O) = 2-NORM OF (D**-1)*G. * V(DSTNRM) (OUTPUT) = 2-NORM OF D*STEP. * V(DST0) (I/O) = 2-NORM OF D*(H**-1)*G (FOR POS. DEF. H ONLY), OR * OVERESTIMATE OF SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1). * V(EPSLON) (IN) = MAX. REL. ERROR ALLOWED FOR PSI(STEP). FOR THE * STEP RETURNED, PSI(STEP) WILL EXCEED ITS OPTIMAL VALUE * BY LESS THAN -V(EPSLON)*PSI(STEP). SUGGESTED VALUE = 0.1. * V(GTSTEP) (OUT) = INNER PRODUCT BETWEEN G AND STEP. * V(NREDUC) (OUT) = PSI(-(H**-1)*G) = PSI(NEWTON STEP) (FOR POS. DEF. * H ONLY -- V(NREDUC) IS SET TO ZERO OTHERWISE). * V(PHMNFC) (IN) = TOL. (TOGETHER WITH V(PHMXFC)) FOR ACCEPTING STEP * (MORE*S SIGMA). THE ERROR V(DSTNRM) - V(RADIUS) MUST LIE * BETWEEN V(PHMNFC)*V(RADIUS) AND V(PHMXFC)*V(RADIUS). * V(PHMXFC) (IN) (SEE V(PHMNFC).) * SUGGESTED VALUES -- V(PHMNFC) = -0.25, V(PHMXFC) = 0.5. * V(PREDUC) (OUT) = PSI(STEP) = PREDICTED OBJ. FUNC. REDUCTION FOR STEP. * V(RADIUS) (IN) = RADIUS OF CURRENT (SCALED) TRUST REGION. * V(RAD0) (I/O) = VALUE OF V(RADIUS) FROM PREVIOUS CALL. * V(STPPAR) (I/O) IS NORMALLY THE MARQUARDT PARAMETER, I.E. THE ALPHA * DESCRIBED BELOW UNDER ALGORITHM NOTES. IF H + ALPHA*D**2 * (SEE ALGORITHM NOTES) IS (NEARLY) SINGULAR, HOWEVER, * THEN V(STPPAR) = -ALPHA. * * *** USAGE NOTES *** * * IF IT IS DESIRED TO RECOMPUTE STEP USING A DIFFERENT VALUE OF * V(RADIUS), THEN THIS ROUTINE MAY BE RESTARTED BY CALLING IT * WITH ALL PARAMETERS UNCHANGED EXCEPT V(RADIUS). (THIS EXPLAINS * WHY STEP AND W ARE LISTED AS I/O). ON AN INITIAL CALL (ONE WITH * KA .LT. 0), STEP AND W NEED NOT BE INITIALIZED AND ONLY COMPO- * NENTS V(EPSLON), V(STPPAR), V(PHMNFC), V(PHMXFC), V(RADIUS), AND * V(RAD0) OF V MUST BE INITIALIZED. * * *** ALGORITHM NOTES *** * * THE DESIRED G-Q-T STEP (REF. 2, 3, 4, 6) SATISFIES * (H + ALPHA*D**2)*STEP = -G FOR SOME NONNEGATIVE ALPHA SUCH THAT * H + ALPHA*D**2 IS POSITIVE SEMIDEFINITE. ALPHA AND STEP ARE * COMPUTED BY A SCHEME ANALOGOUS TO THE ONE DESCRIBED IN REF. 5. * ESTIMATES OF THE SMALLEST AND LARGEST EIGENVALUES OF THE HESSIAN * ARE OBTAINED FROM THE GERSCHGORIN CIRCLE THEOREM ENHANCED BY A * SIMPLE FORM OF THE SCALING DESCRIBED IN REF. 7. CASES IN WHICH * H + ALPHA*D**2 IS NEARLY (OR EXACTLY) SINGULAR ARE HANDLED BY * THE TECHNIQUE DISCUSSED IN REF. 2. IN THESE CASES, A STEP OF * (EXACT) LENGTH V(RADIUS) IS RETURNED FOR WHICH PSI(STEP) EXCEEDS * ITS OPTIMAL VALUE BY LESS THAN -V(EPSLON)*PSI(STEP). THE TEST * SUGGESTED IN REF. 6 FOR DETECTING THE SPECIAL CASE IS PERFORMED * ONCE TWO MATRIX FACTORIZATIONS HAVE BEEN DONE -- DOING SO SOONER * SEEMS TO DEGRADE THE PERFORMANCE OF OPTIMIZATION ROUTINES THAT * CALL THIS ROUTINE. * * *** FUNCTIONS AND SUBROUTINES CALLED *** * * DD7TPR - RETURNS INNER PRODUCT OF TWO VECTORS. * DL7ITV - APPLIES INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX. * DL7IVM - APPLIES INVERSE OF COMPACT LOWER TRIANG. MATRIX. * DL7SRT - FINDS CHOLESKY FACTOR (OF COMPACTLY STORED LOWER TRIANG.). * DL7SVN - RETURNS APPROX. TO MIN. SING. VALUE OF LOWER TRIANG. MATRIX. * DR7MDC - RETURNS MACHINE-DEPENDENT CONSTANTS. * DV2NRM - RETURNS 2-NORM OF A VECTOR. * */ /* *** REFERENCES *** * * 1. DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1981), AN ADAPTIVE * NONLINEAR LEAST-SQUARES ALGORITHM, ACM TRANS. MATH. * SOFTWARE, VOL. 7, NO. 3. * 2. GAY, D.M. (1981), COMPUTING OPTIMAL LOCALLY CONSTRAINED STEPS, * SIAM J. SCI. STATIST. COMPUTING, VOL. 2, NO. 2, PP. * 186-197. * 3. GOLDFELD, S.M., QUANDT, R.E., AND TROTTER, H.F. (1966), * MAXIMIZATION BY QUADRATIC HILL-CLIMBING, ECONOMETRICA 34, * PP. 541-551. * 4. HEBDEN, M.D. (1973), AN ALGORITHM FOR MINIMIZATION USING EXACT * SECOND DERIVATIVES, REPORT T.P. 515, THEORETICAL PHYSICS * DIV., A.E.R.E. HARWELL, OXON., ENGLAND. * 5. MORE, J.J. (1978), THE LEVENBERG-MARQUARDT ALGORITHM, IMPLEMEN- * TATION AND THEORY, PP.105-116 OF SPRINGER LECTURE NOTES * IN MATHEMATICS NO. 630, EDITED BY G.A. WATSON, SPRINGER- * VERLAG, BERLIN AND NEW YORK. * 6. MORE, J.J., AND SORENSEN, D.C. (1981), COMPUTING A TRUST REGION * STEP, TECHNICAL REPORT ANL-81-83, ARGONNE NATIONAL LAB. * 7. VARGA, R.S. (1965), MINIMAL GERSCHGORIN SETS, PACIFIC J. MATH. 15, * PP. 719-729. * * *** GENERAL *** * * CODED BY DAVID M. GAY. * THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH * SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS * MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND * MCS-7906671. * * ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * * *** LOCAL VARIABLES *** * */ /* *** CONSTANTS *** */ /* *** EXTERNAL FUNCTIONS AND SUBROUTINES *** * */ /* *** SUBSCRIPTS FOR V *** * */ /* *** BODY *** * */ if (big <= ZERO) big = dr7mdc( 6 ); /* *** STORE LARGEST ABS. ENTRY IN (D**-1)*H*(D**-1) AT W(DGGDMX). */ dggdmx = p + 1; /* *** STORE GERSCHGORIN OVER- AND UNDERESTIMATES OF THE LARGEST * *** AND SMALLEST EIGENVALUES OF (D**-1)*H*(D**-1) AT W(EMAX) * *** AND W(EMIN) RESPECTIVELY. */ emax = dggdmx + 1; emin = emax + 1; /* *** FOR USE IN RECOMPUTING STEP, THE FINAL VALUES OF LK, UK, DST, * *** AND THE INVERSE DERIVATIVE OF MORE*S PHI AT 0 (FOR POS. DEF. * *** H) ARE STORED IN W(LK0), W(UK0), W(DSTSAV), AND W(PHIPIN) * *** RESPECTIVELY. */ lk0 = emin + 1; phipin = lk0 + 1; uk0 = phipin + 1; dstsav = uk0 + 1; /* *** STORE DIAG OF (D**-1)*H*(D**-1) IN W(DIAG),...,W(DIAG0+P). */ diag0 = dstsav; diag = diag0 + 1; /* *** STORE -D*STEP IN W(Q),...,W(Q0+P). */ q0 = diag0 + p; q = q0 + 1; /* *** ALLOCATE STORAGE FOR SCRATCH VECTOR X *** */ x = q + p; rad = V[RADIUS]; radsq = SQ(rad); /* *** PHITOL = MAX. ERROR ALLOWED IN DST = V(DSTNRM) = 2-NORM OF * *** D*STEP. */ phimax = V[PHMXFC]*rad; phimin = V[PHMNFC]*rad; psifac = big; t1 = TWO*V[EPSLON]/(THREE*(FOUR*(V[PHMNFC] + ONE)*(KAPPA + ONE) + KAPPA + TWO)*rad); if (t1 < big*fmin( rad, ONE )) psifac = t1/rad; /* *** OLDPHI IS USED TO DETECT LIMITS OF NUMERICAL ACCURACY. IF * *** WE RECOMPUTE STEP AND IT DOES NOT CHANGE, THEN WE ACCEPT IT. */ oldphi = ZERO; eps = V[EPSLON]; irc = 0; restrt = FALSE; kalim = *ka + 50; /* *** START OR RESTART, DEPENDING ON KA *** * */ if (*ka >= 0) goto L_290; /* *** FRESH START *** * */ k = 0; uk = NEGONE; *ka = 0; kalim = 50; V[DGNORM] = dv2nrm( p, dig ); V[NREDUC] = ZERO; V[DST0] = ZERO; kamin = 3; if (V[DGNORM] == ZERO) kamin = 0; /* *** STORE DIAG(DIHDI) IN W(DIAG0+1),...,W(DIAG0+P) *** * */ j = 0; for (i = 1; i <= p; i++) { j += i; k1 = diag0 + i; W[k1] = Dihdi[j]; } /* *** DETERMINE W(DGGDMX), THE LARGEST ELEMENT OF DIHDI *** * */ t1 = ZERO; j = p*(p + 1)/2; for (i = 1; i <= j; i++) { t = fabs( Dihdi[i] ); if (t1 < t) t1 = t; } W[dggdmx] = t1; /* *** TRY ALPHA = 0 *** * */ L_30: dl7srt( 1, p, l, dihdi, &irc ); if (irc == 0) goto L_50; /* *** INDEF. H -- UNDERESTIMATE SMALLEST EIGENVALUE, USE THIS * *** ESTIMATE TO INITIALIZE LOWER BOUND LK ON ALPHA. */ j = irc*(irc + 1)/2; t = L[j]; L[j] = ONE; for (i = 1; i <= irc; i++) { W[i] = ZERO; } W[irc] = ONE; dl7itv( irc, w, l, w ); t1 = dv2nrm( irc, w ); lk = -t/t1/t1; V[DST0] = -lk; if (restrt) goto L_210; goto L_70; /* *** POSITIVE DEFINITE H -- COMPUTE UNMODIFIED NEWTON STEP. *** */ L_50: lk = ZERO; t = dl7svn( p, l, &W[q], &W[q] ); if (t >= ONE) goto L_60; if (V[DGNORM] >= t*t*big) goto L_70; L_60: dl7ivm( p, &W[q], l, dig ); gtsta = dd7tpr( p, &W[q], &W[q] ); V[NREDUC] = HALF*gtsta; dl7itv( p, &W[q], l, &W[q] ); dst = dv2nrm( p, &W[q] ); V[DST0] = dst; phi = dst - rad; if (phi <= phimax) goto L_260; if (restrt) goto L_210; /* *** PREPARE TO COMPUTE GERSCHGORIN ESTIMATES OF LARGEST (AND * *** SMALLEST) EIGENVALUES. *** * */ L_70: k = 0; for (i = 1; i <= p; i++) { wi = ZERO; if (i == 1) goto L_90; im1 = i - 1; for (j = 1; j <= im1; j++) { k += 1; t = fabs( Dihdi[k] ); wi += t; W[j] += t; } L_90: W[i] = wi; k += 1; } /* *** (UNDER-)ESTIMATE SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1) *** * */ k = 1; t1 = W[diag] - W[1]; if (p <= 1) goto L_120; for (i = 2; i <= p; i++) { j = diag0 + i; t = W[j] - W[i]; if (t >= t1) goto L_110; t1 = t; k = i; L_110: ; } L_120: sk = W[k]; j = diag0 + k; akk = W[j]; k1 = k*(k - 1)/2 + 1; inc = 1; t = ZERO; for (i = 1; i <= p; i++) { if (i == k) goto L_130; aki = fabs( Dihdi[k1] ); si = W[i]; j = diag0 + i; t1 = HALF*(akk - W[j] + si - aki); t1 += sqrt( t1*t1 + sk*aki ); if (t < t1) t = t1; if (i < k) goto L_140; L_130: inc = i; L_140: k1 += inc; } W[emin] = akk - t; uk = V[DGNORM]/rad - W[emin]; if (V[DGNORM] == ZERO) uk += P001 + P001*uk; if (uk <= ZERO) uk = P001; /* *** COMPUTE GERSCHGORIN (OVER-)ESTIMATE OF LARGEST EIGENVALUE *** * */ k = 1; t1 = W[diag] + W[1]; if (p <= 1) goto L_170; for (i = 2; i <= p; i++) { j = diag0 + i; t = W[j] + W[i]; if (t <= t1) goto L_160; t1 = t; k = i; L_160: ; } L_170: sk = W[k]; j = diag0 + k; akk = W[j]; k1 = k*(k - 1)/2 + 1; inc = 1; t = ZERO; for (i = 1; i <= p; i++) { if (i == k) goto L_180; aki = fabs( Dihdi[k1] ); si = W[i]; j = diag0 + i; t1 = HALF*(W[j] + si - aki - akk); t1 += sqrt( t1*t1 + sk*aki ); if (t < t1) t = t1; if (i < k) goto L_190; L_180: inc = i; L_190: k1 += inc; } W[emax] = akk + t; lk = fmax( lk, V[DGNORM]/rad - W[emax] ); /* *** ALPHAK = CURRENT VALUE OF ALPHA (SEE ALG. NOTES ABOVE). WE * *** USE MORE*S SCHEME FOR INITIALIZING IT. */ alphak = fabs( V[STPPAR] )*V[RAD0]/rad; alphak = fmin( uk, fmax( alphak, lk ) ); if (irc != 0) goto L_210; /* *** COMPUTE L0 FOR POSITIVE DEFINITE H *** * */ dl7ivm( p, w, l, &W[q] ); t = dv2nrm( p, w ); W[phipin] = rad/t/t; lk = fmax( lk, phi*W[phipin] ); /* *** SAFEGUARD ALPHAK AND ADD ALPHAK*I TO (D**-1)*H*(D**-1) *** * */ L_210: *ka += 1; if ((-V[DST0] >= alphak || alphak < lk) || alphak >= uk) alphak = uk*fmax( P001, sqrt( lk/uk ) ); if (alphak <= ZERO) alphak = HALF*uk; if (alphak <= ZERO) alphak = uk; k = 0; for (i = 1; i <= p; i++) { k += i; j = diag0 + i; Dihdi[k] = W[j] + alphak; } /* *** TRY COMPUTING CHOLESKY DECOMPOSITION *** * */ dl7srt( 1, p, l, dihdi, &irc ); if (irc == 0) goto L_240; /* *** (D**-1)*H*(D**-1) + ALPHAK*I IS INDEFINITE -- OVERESTIMATE * *** SMALLEST EIGENVALUE FOR USE IN UPDATING LK *** * */ j = (irc*(irc + 1))/2; t = L[j]; L[j] = ONE; for (i = 1; i <= irc; i++) { W[i] = ZERO; } W[irc] = ONE; dl7itv( irc, w, l, w ); t1 = dv2nrm( irc, w ); lk = alphak - t/t1/t1; V[DST0] = -lk; if (alphak < lk) goto L_210; /* *** NASTY CASE -- EXACT GERSCHGORIN BOUNDS. FUDGE LK, UK... * */ t = P001*alphak; if (t <= ZERO) t = P001; lk = alphak + t; if (uk <= lk) uk = lk + t; goto L_210; /* *** ALPHAK MAKES (D**-1)*H*(D**-1) POSITIVE DEFINITE. * *** COMPUTE Q = -D*STEP, CHECK FOR CONVERGENCE. *** * */ L_240: dl7ivm( p, &W[q], l, dig ); gtsta = dd7tpr( p, &W[q], &W[q] ); dl7itv( p, &W[q], l, &W[q] ); dst = dv2nrm( p, &W[q] ); phi = dst - rad; if (phi <= phimax && phi >= phimin) goto L_270; if (phi == oldphi) goto L_270; oldphi = phi; if (phi < ZERO) goto L_330; /* *** UNACCEPTABLE ALPHAK -- UPDATE LK, UK, ALPHAK *** * */ L_250: if (*ka >= kalim) goto L_270; /* *** THE FOLLOWING min IS NECESSARY BECAUSE OF RESTARTS *** */ if (phi < ZERO) uk = fmin( uk, alphak ); /* *** KAMIN = 0 ONLY IFF THE GRADIENT VANISHES *** */ if (kamin == 0) goto L_210; dl7ivm( p, w, l, &W[q] ); /* *** THE FOLLOWING, COMMENTED CALCULATION OF ALPHAK IS SOMETIMES * *** SAFER BUT WORSE IN PERFORMANCE... * T1 = DST / DV2NRM(P, W) * ALPHAK = ALPHAK + T1 * (PHI/RAD) * T1 */ t1 = dv2nrm( p, w ); alphak += (phi/t1)*(dst/t1)*(dst/rad); lk = fmax( lk, alphak ); alphak = lk; goto L_210; /* *** ACCEPTABLE STEP ON FIRST TRY *** * */ L_260: alphak = ZERO; /* *** SUCCESSFUL STEP IN GENERAL. COMPUTE STEP = -(D**-1)*Q *** * */ L_270: for (i = 1; i <= p; i++) { j = q0 + i; Step[i] = -W[j]/D[i]; } V[GTSTEP] = -gtsta; V[PREDUC] = HALF*(fabs( alphak )*dst*dst + gtsta); goto L_410; /* *** RESTART WITH NEW RADIUS *** * */ L_290: if (V[DST0] <= ZERO || V[DST0] - rad > phimax) goto L_310; /* *** PREPARE TO RETURN NEWTON STEP *** * */ restrt = TRUE; *ka += 1; k = 0; for (i = 1; i <= p; i++) { k += i; j = diag0 + i; Dihdi[k] = W[j]; } uk = NEGONE; goto L_30; L_310: kamin = *ka + 3; if (V[DGNORM] == ZERO) kamin = 0; if (*ka == 0) goto L_50; dst = W[dstsav]; alphak = fabs( V[STPPAR] ); phi = dst - rad; t = V[DGNORM]/rad; uk = t - W[emin]; if (V[DGNORM] == ZERO) uk += P001 + P001*uk; if (uk <= ZERO) uk = P001; if (rad > V[RAD0]) goto L_320; /* *** SMALLER RADIUS *** */ lk = ZERO; if (alphak > ZERO) lk = W[lk0]; lk = fmax( lk, t - W[emax] ); if (V[DST0] > ZERO) lk = fmax( lk, (V[DST0] - rad)*W[phipin] ); goto L_250; /* *** BIGGER RADIUS *** */ L_320: if (alphak > ZERO) uk = fmin( uk, W[uk0] ); lk = fmax( fmax( ZERO, -V[DST0] ), t - W[emax] ); if (V[DST0] > ZERO) lk = fmax( lk, (V[DST0] - rad)*W[phipin] ); goto L_250; /* *** DECIDE WHETHER TO CHECK FOR SPECIAL CASE... IN PRACTICE (FROM * *** THE STANDPOINT OF THE CALLING OPTIMIZATION CODE) IT SEEMS BEST * *** NOT TO CHECK UNTIL A FEW ITERATIONS HAVE FAILED -- HENCE THE * *** TEST ON KAMIN BELOW. * */ L_330: delta = alphak + fmin( ZERO, V[DST0] ); twopsi = alphak*dst*dst + gtsta; if (*ka >= kamin) goto L_340; /* *** IF THE TEST IN REF. 2 IS SATISFIED, FALL THROUGH TO HANDLE * *** THE SPECIAL CASE (AS SOON AS THE MORE-SORENSEN TEST DETECTS * *** IT). */ if (psifac >= big) goto L_340; if (delta >= psifac*twopsi) goto L_370; /* *** CHECK FOR THE SPECIAL CASE OF H + ALPHA*D**2 (NEARLY) * *** SINGULAR. USE ONE STEP OF INVERSE POWER METHOD WITH START * *** FROM DL7SVN TO OBTAIN APPROXIMATE EIGENVECTOR CORRESPONDING * *** TO SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1). DL7SVN RETURNS * *** X AND W WITH L*W = X. * */ L_340: t = dl7svn( p, l, &W[x], w ); /* *** NORMALIZE W *** */ for (i = 1; i <= p; i++) { W[i] *= t; } /* *** COMPLETE CURRENT INV. POWER ITER. -- REPLACE W BY (L**-T)*W. */ dl7itv( p, w, l, w ); t2 = ONE/dv2nrm( p, w ); for (i = 1; i <= p; i++) { W[i] *= t2; } t *= t2; /* *** NOW W IS THE DESIRED APPROXIMATE (UNIT) EIGENVECTOR AND * *** T*X = ((D**-1)*H*(D**-1) + ALPHAK*I)*W. * */ sw = dd7tpr( p, &W[q], w ); t1 = (rad + dst)*(rad - dst); root = sqrt( sw*sw + t1 ); if (sw < ZERO) root = -root; si = t1/(sw + root); /* *** THE ACTUAL TEST FOR THE SPECIAL CASE... * */ if (powi(t2*si,2) <= eps*(SQ(dst) + alphak*radsq)) goto L_380; /* *** UPDATE UPPER BOUND ON SMALLEST EIGENVALUE (WHEN NOT POSITIVE) * *** (AS RECOMMENDED BY MORE AND SORENSEN) AND CONTINUE... * */ if (V[DST0] <= ZERO) V[DST0] = fmin( V[DST0], SQ(t2) - alphak ); lk = fmax( lk, -V[DST0] ); /* *** CHECK WHETHER WE CAN HOPE TO DETECT THE SPECIAL CASE IN * *** THE AVAILABLE ARITHMETIC. ACCEPT STEP AS IT IS IF NOT. * * *** IF NOT YET AVAILABLE, OBTAIN MACHINE DEPENDENT VALUE DGXFAC. */ L_370: if (dgxfac == ZERO) dgxfac = EPSFAC*dr7mdc( 3 ); if (delta > dgxfac*W[dggdmx]) goto L_250; goto L_270; /* *** SPECIAL CASE DETECTED... NEGATE ALPHAK TO INDICATE SPECIAL CASE * */ L_380: alphak = -alphak; V[PREDUC] = HALF*twopsi; /* *** ACCEPT CURRENT STEP IF ADDING SI*W WOULD LEAD TO A * *** FURTHER RELATIVE REDUCTION IN PSI OF LESS THAN V(EPSLON)/3. * */ t1 = ZERO; t = si*(alphak*sw - HALF*si*(alphak + t*dd7tpr( p, &W[x], w ))); if (t < eps*twopsi/SIX) goto L_390; V[PREDUC] += t; dst = rad; t1 = -si; L_390: for (i = 1; i <= p; i++) { j = q0 + i; W[j] = t1*W[i] - W[j]; Step[i] = W[j]/D[i]; } V[GTSTEP] = dd7tpr( p, dig, &W[q] ); /* *** SAVE VALUES FOR USE IN A POSSIBLE RESTART *** * */ L_410: V[DSTNRM] = dst; V[STPPAR] = alphak; W[lk0] = lk; W[uk0] = uk; V[RAD0] = rad; W[dstsav] = dst; /* *** RESTORE DIAGONAL OF DIHDI *** * */ j = 0; for (i = 1; i <= p; i++) { j += i; k = diag0 + i; Dihdi[j] = W[k]; } /* *** LAST CARD OF DG7QTS FOLLOWS *** */ return; } /* end of function */ /* PARAMETER translations */ #undef ZERO #define ZERO 0.e0 /* end of PARAMETER translations */ void /*FUNCTION*/ dl7ivm( long n, double x[], double l[], double y[]) { long int i, j, k; double t; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const L = &l[0] - 1; double *const X = &x[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* *** SOLVE L*X = Y, WHERE L IS AN N X N LOWER TRIANGULAR * *** MATRIX STORED COMPACTLY BY ROWS. X AND Y MAY OCCUPY THE SAME * *** STORAGE. *** * */ for (k = 1; k <= n; k++) { if (Y[k] != ZERO) goto L_20; X[k] = ZERO; } goto L_999; L_20: j = k*(k + 1)/2; X[k] = Y[k]/L[j]; if (k >= n) goto L_999; k += 1; for (i = k; i <= n; i++) { t = dd7tpr( i - 1, &L[j + 1], x ); j += i; X[i] = (Y[i] - t)/L[j]; } L_999: return; /* *** LAST CARD OF DL7IVM FOLLOWS *** */ } /* end of function */ /* PARAMETER translations */ #undef ONE #define ONE 1.e0 /* end of PARAMETER translations */ double /*FUNCTION*/ dl7svx( long p, double l[], double x[], double y[]) { long int i, ix, j, j0, ji, jj, jjj, jm1, pm1, pplus1; double b, blji, dl7svx_v, sminus, splus, t, yi; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const L = &l[0] - 1; double *const X = &x[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* *** ESTIMATE LARGEST SING. VALUE OF PACKED LOWER TRIANG. MATRIX L * * *** PARAMETER DECLARATIONS *** * */ /* DIMENSION L(P*(P+1)/2) * * ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * * *** PURPOSE *** * * THIS FUNCTION RETURNS A GOOD UNDER-ESTIMATE OF THE LARGEST * SINGULAR VALUE OF THE PACKED LOWER TRIANGULAR MATRIX L. * * *** PARAMETER DESCRIPTION *** * * P (IN) = THE ORDER OF L. L IS A P X P LOWER TRIANGULAR MATRIX. * L (IN) = ARRAY HOLDING THE ELEMENTS OF L IN ROW ORDER, I.E. * L(1,1), L(2,1), L(2,2), L(3,1), L(3,2), L(3,3), ETC. * X (OUT) IF DL7SVX RETURNS A POSITIVE VALUE, THEN X = (L**T)*Y IS AN * (UNNORMALIZED) APPROXIMATE RIGHT SINGULAR VECTOR * CORRESPONDING TO THE LARGEST SINGULAR VALUE. THIS * APPROXIMATION MAY BE CRUDE. * Y (OUT) IF DL7SVX RETURNS A POSITIVE VALUE, THEN Y = L*X IS A * NORMALIZED APPROXIMATE LEFT SINGULAR VECTOR CORRESPOND- * ING TO THE LARGEST SINGULAR VALUE. THIS APPROXIMATION * MAY BE VERY CRUDE. THE CALLER MAY PASS THE SAME VECTOR * FOR X AND Y (NONSTANDARD FORTRAN USAGE), IN WHICH CASE X * OVER-WRITES Y. * * *** ALGORITHM NOTES *** * * THE ALGORITHM IS BASED ON ANALOGY WITH (1). IT USES A * RANDOM NUMBER GENERATOR PROPOSED IN (4), WHICH PASSES THE * SPECTRAL TEST WITH FLYING COLORS -- SEE (2) AND (3). * * *** SUBROUTINES AND FUNCTIONS CALLED *** * * DV2NRM - FUNCTION, RETURNS THE 2-NORM OF A VECTOR. * * *** REFERENCES *** * * (1) CLINE, A., MOLER, C., STEWART, G., AND WILKINSON, J.H.(1977), * AN ESTIMATE FOR THE CONDITION NUMBER OF A MATRIX, REPORT * TM-310, APPLIED MATH. DIV., ARGONNE NATIONAL LABORATORY. * * (2) HOAGLIN, D.C. (1976), THEORETICAL PROPERTIES OF CONGRUENTIAL * RANDOM-NUMBER GENERATORS -- AN EMPIRICAL VIEW, * MEMORANDUM NS-340, DEPT. OF STATISTICS, HARVARD UNIV. * * (3) KNUTH, D.E. (1969), THE ART OF COMPUTER PROGRAMMING, VOL. 2 * (SEMINUMERICAL ALGORITHMS), ADDISON-WESLEY, READING, MASS. * * (4) SMITH, C.S. (1971), MULTIPLICATIVE PSEUDO-RANDOM NUMBER * GENERATORS WITH PRIME MODULUS, J. ASSOC. COMPUT. MACH. 18, * PP. 586-593. * * *** HISTORY *** * * DESIGNED AND CODED BY DAVID M. GAY (WINTER 1977/SUMMER 1978). * * *** GENERAL *** * * THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH * SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS * MCS-7600324, DCR75-10143, 76-14311DSS, AND MCS76-11989. * * ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * * *** LOCAL VARIABLES *** * */ /* *** CONSTANTS *** * */ /* *** EXTERNAL FUNCTIONS AND SUBROUTINES *** * */ /* *** BODY *** * */ ix = 2; pplus1 = p + 1; pm1 = p - 1; /* *** FIRST INITIALIZE X TO PARTIAL SUMS *** * */ j0 = p*pm1/2; jj = j0 + p; ix = (3432*ix)%9973; b = HALF*(ONE + (double)( ix )/R9973); X[p] = b*L[jj]; if (p <= 1) goto L_40; for (i = 1; i <= pm1; i++) { ji = j0 + i; X[i] = b*L[ji]; } /* *** COMPUTE X = (L**T)*B, WHERE THE COMPONENTS OF B HAVE RANDOMLY * *** CHOSEN MAGNITUDES IN (.5,1) WITH SIGNS CHOSEN TO MAKE X LARGE. * * DO J = P-1 TO 1 BY -1... */ for (jjj = 1; jjj <= pm1; jjj++) { j = p - jjj; /* *** DETERMINE X(J) IN THIS ITERATION. NOTE FOR I = 1,2,...,J * *** THAT X(I) HOLDS THE CURRENT PARTIAL SUM FOR ROW I. */ ix = (3432*ix)%9973; b = HALF*(ONE + (double)( ix )/R9973); jm1 = j - 1; j0 = j*jm1/2; splus = ZERO; sminus = ZERO; for (i = 1; i <= j; i++) { ji = j0 + i; blji = b*L[ji]; splus += fabs( blji + X[i] ); sminus += fabs( blji - X[i] ); } if (sminus > splus) b = -b; X[j] = ZERO; /* *** UPDATE PARTIAL SUMS *** */ dv2axy( j, x, b, &L[j0 + 1], x ); } /* *** NORMALIZE X *** * */ L_40: t = dv2nrm( p, x ); if (t <= ZERO) goto L_80; t = ONE/t; for (i = 1; i <= p; i++) { X[i] *= t; } /* *** COMPUTE L*X = Y AND RETURN SVMAX = TWONORM(Y) *** * */ for (jjj = 1; jjj <= p; jjj++) { j = pplus1 - jjj; ji = j*(j - 1)/2 + 1; Y[j] = dd7tpr( j, &L[ji], x ); } /* *** NORMALIZE Y AND SET X = (L**T)*Y *** * */ t = ONE/dv2nrm( p, y ); ji = 1; for (i = 1; i <= p; i++) { yi = t*Y[i]; X[i] = ZERO; dv2axy( i, x, yi, &L[ji], x ); ji += i; } dl7svx_v = dv2nrm( p, x ); goto L_999; L_80: dl7svx_v = ZERO; L_999: return( dl7svx_v ); /* *** LAST CARD OF DL7SVX FOLLOWS *** */ } /* end of function */ /* ================================================================== */ /* PARAMETER translations */ #define DSTSAV 18 #define F 10 #define F0 13 #define FDIF 11 #define FLSTGD 12 #define GTSLST 14 #define IRC 29 #define MLSTGD 32 #define MODEL 5 #define NFCALL 6 #define NFGCAL 7 #define ONEP2 1.2e0 #define PLSTGD 15 #define RADFAC 16 #define RADINC 8 #define RELDX 17 #define RESTOR 9 #define STAGE 10 #define STGLIM 11 #define SWITCH_ 12 #define TOOBIG 2 #undef TWO #define TWO 2.e0 #define XIRC 13 /* end of PARAMETER translations */ void /*FUNCTION*/ da7sst( long iv[], long liv, long lv, double v[]) { LOGICAL32 goodx; long int i, nfc; double emax, emaxs, gts, rfac1, xmax; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Iv = &iv[0] - 1; double *const V = &v[0] - 1; /* end of OFFSET VECTORS */ /* *** ASSESS CANDIDATE STEP (***SOL VERSION 2.3) *** * * ------------------------------------------------------------------ */ /* *** PURPOSE *** * * THIS SUBROUTINE IS CALLED BY AN UNCONSTRAINED MINIMIZATION * ROUTINE TO ASSESS THE NEXT CANDIDATE STEP. IT MAY RECOMMEND ONE * OF SEVERAL COURSES OF ACTION, SUCH AS ACCEPTING THE STEP, RECOM- * PUTING IT USING THE SAME OR A NEW QUADRATIC MODEL, OR HALTING DUE * TO CONVERGENCE OR FALSE CONVERGENCE. SEE THE RETURN CODE LISTING * BELOW. * * ------------------------- PARAMETER USAGE -------------------------- * * IV (I/O) INTEGER PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTION * BELOW OF IV VALUES REFERENCED. * LIV (IN) LENGTH OF IV ARRAY. * LV (IN) LENGTH OF V ARRAY. * V (I/O) REAL PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTION * BELOW OF V VALUES REFERENCED. * * *** IV VALUES REFERENCED *** * * IV(IRC) (I/O) ON INPUT FOR THE FIRST STEP TRIED IN A NEW ITERATION, * IV(IRC) SHOULD BE SET TO 3 OR 4 (THE VALUE TO WHICH IT IS * SET WHEN STEP IS DEFINITELY TO BE ACCEPTED). ON INPUT * AFTER STEP HAS BEEN RECOMPUTED, IV(IRC) SHOULD BE * UNCHANGED SINCE THE PREVIOUS RETURN OF DA7SST. * ON OUTPUT, IV(IRC) IS A RETURN CODE HAVING ONE OF THE * FOLLOWING VALUES... * 1 = SWITCH MODELS OR TRY SMALLER STEP. * 2 = SWITCH MODELS OR ACCEPT STEP. * 3 = ACCEPT STEP AND DETERMINE V(RADFAC) BY GRADIENT * TESTS. * 4 = ACCEPT STEP, V(RADFAC) HAS BEEN DETERMINED. * 5 = RECOMPUTE STEP (USING THE SAME MODEL). * 6 = RECOMPUTE STEP WITH RADIUS = V(LMAXS) BUT DO NOT * EVAULATE THE OBJECTIVE FUNCTION. * 7 = X-CONVERGENCE (SEE V(XCTOL)). * 8 = RELATIVE FUNCTION CONVERGENCE (SEE V(RFCTOL)). * 9 = BOTH X- AND RELATIVE FUNCTION CONVERGENCE. * 10 = ABSOLUTE FUNCTION CONVERGENCE (SEE V(AFCTOL)). * 11 = SINGULAR CONVERGENCE (SEE V(LMAXS)). * 12 = FALSE CONVERGENCE (SEE V(XFTOL)). * 13 = IV(IRC) WAS OUT OF RANGE ON INPUT. * RETURN CODE I HAS PRECDENCE OVER I+1 FOR I = 9, 10, 11. * IV(MLSTGD) (I/O) SAVED VALUE OF IV(MODEL). * IV(MODEL) (I/O) ON INPUT, IV(MODEL) SHOULD BE AN INTEGER IDENTIFYING * THE CURRENT QUADRATIC MODEL OF THE OBJECTIVE FUNCTION. * IF A PREVIOUS STEP YIELDED A BETTER FUNCTION REDUCTION, * THEN IV(MODEL) WILL BE SET TO IV(MLSTGD) ON OUTPUT. * IV(NFCALL) (IN) INVOCATION COUNT FOR THE OBJECTIVE FUNCTION. * IV(NFGCAL) (I/O) VALUE OF IV(NFCALL) AT STEP THAT GAVE THE BIGGEST * FUNCTION REDUCTION THIS ITERATION. IV(NFGCAL) REMAINS * UNCHANGED UNTIL A FUNCTION REDUCTION IS OBTAINED. * IV(RADINC) (I/O) THE NUMBER OF RADIUS INCREASES (OR MINUS THE NUMBER * OF DECREASES) SO FAR THIS ITERATION. * IV(RESTOR) (OUT) SET TO 1 IF V(F) HAS BEEN RESTORED AND X SHOULD BE * RESTORED TO ITS INITIAL VALUE, TO 2 IF X SHOULD BE SAVED, * TO 3 IF X SHOULD BE RESTORED FROM THE SAVED VALUE, AND TO * 0 OTHERWISE. * IV(STAGE) (I/O) COUNT OF THE NUMBER OF MODELS TRIED SO FAR IN THE * CURRENT ITERATION. * IV(STGLIM) (IN) MAXIMUM NUMBER OF MODELS TO CONSIDER. * IV(SWITCH) (OUT) SET TO 0 UNLESS A NEW MODEL IS BEING TRIED AND IT * GIVES A SMALLER FUNCTION VALUE THAN THE PREVIOUS MODEL, * IN WHICH CASE DA7SST SETS IV(SWITCH) = 1. * IV(TOOBIG) (IN) IS NONZERO IF STEP WAS TOO BIG (E.G. IF IT CAUSED * OVERFLOW). * IV(XIRC) (I/O) VALUE THAT IV(IRC) WOULD HAVE IN THE ABSENCE OF * CONVERGENCE, FALSE CONVERGENCE, AND OVERSIZED STEPS. * * *** V VALUES REFERENCED *** * * V(AFCTOL) (IN) ABSOLUTE FUNCTION CONVERGENCE TOLERANCE. IF THE * ABSOLUTE VALUE OF THE CURRENT FUNCTION VALUE V(F) IS LESS * THAN V(AFCTOL), THEN DA7SST RETURNS WITH IV(IRC) = 10. * V(DECFAC) (IN) FACTOR BY WHICH TO DECREASE RADIUS WHEN IV(TOOBIG) IS * NONZERO. * V(DSTNRM) (IN) THE 2-NORM OF D*STEP. * V(DSTSAV) (I/O) VALUE OF V(DSTNRM) ON SAVED STEP. * V(DST0) (IN) THE 2-NORM OF D TIMES THE NEWTON STEP (WHEN DEFINED, * I.E., FOR V(NREDUC) .GE. 0). * V(F) (I/O) ON BOTH INPUT AND OUTPUT, V(F) IS THE OBJECTIVE FUNC- * TION VALUE AT X. IF X IS RESTORED TO A PREVIOUS VALUE, * THEN V(F) IS RESTORED TO THE CORRESPONDING VALUE. * V(FDIF) (OUT) THE FUNCTION REDUCTION V(F0) - V(F) (FOR THE OUTPUT * VALUE OF V(F) IF AN EARLIER STEP GAVE A BIGGER FUNCTION * DECREASE, AND FOR THE INPUT VALUE OF V(F) OTHERWISE). * V(FLSTGD) (I/O) SAVED VALUE OF V(F). * V(F0) (IN) OBJECTIVE FUNCTION VALUE AT START OF ITERATION. * V(GTSLST) (I/O) VALUE OF V(GTSTEP) ON SAVED STEP. * V(GTSTEP) (IN) INNER PRODUCT BETWEEN STEP AND GRADIENT. * V(INCFAC) (IN) MINIMUM FACTOR BY WHICH TO INCREASE RADIUS. * V(LMAXS) (IN) MAXIMUM REASONABLE STEP SIZE (AND INITIAL STEP BOUND). * IF THE ACTUAL FUNCTION DECREASE IS NO MORE THAN TWICE * WHAT WAS PREDICTED, IF A RETURN WITH IV(IRC) = 7, 8, 9, * OR 10 DOES NOT OCCUR, IF V(DSTNRM) .GT. V(LMAXS), AND IF * V(PREDUC) .LE. V(SCTOL) * ABS(V(F0)), THEN DA7SST RE- * TURNS WITH IV(IRC) = 11. IF SO DOING APPEARS WORTHWHILE, * THEN DA7SST REPEATS THIS TEST WITH V(PREDUC) COMPUTED FOR * A STEP OF LENGTH V(LMAXS) (BY A RETURN WITH IV(IRC) = 6). */ /* V(NREDUC) (I/O) FUNCTION REDUCTION PREDICTED BY QUADRATIC MODEL FOR * NEWTON STEP. IF DA7SST IS CALLED WITH IV(IRC) = 6, I.E., * IF V(PREDUC) HAS BEEN COMPUTED WITH RADIUS = V(LMAXS) FOR * USE IN THE SINGULAR CONVERVENCE TEST, THEN V(NREDUC) IS * SET TO -V(PREDUC) BEFORE THE LATTER IS RESTORED. * V(PLSTGD) (I/O) VALUE OF V(PREDUC) ON SAVED STEP. * V(PREDUC) (I/O) FUNCTION REDUCTION PREDICTED BY QUADRATIC MODEL FOR * CURRENT STEP. * V(RADFAC) (OUT) FACTOR TO BE USED IN DETERMINING THE NEW RADIUS, * WHICH SHOULD BE V(RADFAC)*DST, WHERE DST IS EITHER THE * OUTPUT VALUE OF V(DSTNRM) OR THE 2-NORM OF * DIAG(NEWD)*STEP FOR THE OUTPUT VALUE OF STEP AND THE * UPDATED VERSION, NEWD, OF THE SCALE VECTOR D. FOR * IV(IRC) = 3, V(RADFAC) = 1.0 IS RETURNED. * V(RDFCMN) (IN) MINIMUM VALUE FOR V(RADFAC) IN TERMS OF THE INPUT * VALUE OF V(DSTNRM) -- SUGGESTED VALUE = 0.1. * V(RDFCMX) (IN) MAXIMUM VALUE FOR V(RADFAC) -- SUGGESTED VALUE = 4.0. * V(RELDX) (IN) SCALED RELATIVE CHANGE IN X CAUSED BY STEP, COMPUTED * (E.G.) BY FUNCTION DRLDST AS * MAX (D(I)*ABS(X(I)-X0(I)), 1 .LE. I .LE. P) / * MAX (D(I)*(ABS(X(I))+ABS(X0(I))), 1 .LE. I .LE. P). * V(RFCTOL) (IN) RELATIVE FUNCTION CONVERGENCE TOLERANCE. IF THE * ACTUAL FUNCTION REDUCTION IS AT MOST TWICE WHAT WAS PRE- * DICTED AND V(NREDUC) .LE. V(RFCTOL)*ABS(V(F0)), THEN * DA7SST RETURNS WITH IV(IRC) = 8 OR 9. * V(STPPAR) (IN) MARQUARDT PARAMETER -- 0 MEANS FULL NEWTON STEP. * V(TUNER1) (IN) TUNING CONSTANT USED TO DECIDE IF THE FUNCTION * REDUCTION WAS MUCH LESS THAN EXPECTED. SUGGESTED * VALUE = 0.1. * V(TUNER2) (IN) TUNING CONSTANT USED TO DECIDE IF THE FUNCTION * REDUCTION WAS LARGE ENOUGH TO ACCEPT STEP. SUGGESTED * VALUE = 10**-4. * V(TUNER3) (IN) TUNING CONSTANT USED TO DECIDE IF THE RADIUS * SHOULD BE INCREASED. SUGGESTED VALUE = 0.75. * V(XCTOL) (IN) X-CONVERGENCE CRITERION. IF STEP IS A NEWTON STEP * (V(STPPAR) = 0) HAVING V(RELDX) .LE. V(XCTOL) AND GIVING * AT MOST TWICE THE PREDICTED FUNCTION DECREASE, THEN * DA7SST RETURNS IV(IRC) = 7 OR 9. * V(XFTOL) (IN) FALSE CONVERGENCE TOLERANCE. IF STEP GAVE NO OR ONLY * A SMALL FUNCTION DECREASE AND V(RELDX) .LE. V(XFTOL), * THEN DA7SST RETURNS WITH IV(IRC) = 12. * * ------------------------------ NOTES ------------------------------- * * *** APPLICATION AND USAGE RESTRICTIONS *** * * THIS ROUTINE IS CALLED AS PART OF THE NL2SOL (NONLINEAR * LEAST-SQUARES) PACKAGE. IT MAY BE USED IN ANY UNCONSTRAINED * MINIMIZATION SOLVER THAT USES DOGLEG, GOLDFELD-QUANDT-TROTTER, * OR LEVENBERG-MARQUARDT STEPS. * * *** ALGORITHM NOTES *** * * SEE (1) FOR FURTHER DISCUSSION OF THE ASSESSING AND MODEL * SWITCHING STRATEGIES. WHILE NL2SOL CONSIDERS ONLY TWO MODELS, * DA7SST IS DESIGNED TO HANDLE ANY NUMBER OF MODELS. * * *** USAGE NOTES *** * * ON THE FIRST CALL OF AN ITERATION, ONLY THE I/O VARIABLES * STEP, X, IV(IRC), IV(MODEL), V(F), V(DSTNRM), V(GTSTEP), AND * V(PREDUC) NEED HAVE BEEN INITIALIZED. BETWEEN CALLS, NO I/O * VALUES EXECPT STEP, X, IV(MODEL), V(F) AND THE STOPPING TOLER- * ANCES SHOULD BE CHANGED. * AFTER A RETURN FOR CONVERGENCE OR FALSE CONVERGENCE, ONE CAN * CHANGE THE STOPPING TOLERANCES AND CALL DA7SST AGAIN, IN WHICH * CASE THE STOPPING TESTS WILL BE REPEATED. * * *** REFERENCES *** * * (1) DENNIS, J.E., JR., GAY, D.M., AND WELSCH, R.E. (1981), * AN ADAPTIVE NONLINEAR LEAST-SQUARES ALGORITHM, * ACM TRANS. MATH. SOFTWARE, VOL. 7, NO. 3. * * (2) POWELL, M.J.D. (1970) A FORTRAN SUBROUTINE FOR SOLVING * SYSTEMS OF NONLINEAR ALGEBRAIC EQUATIONS, IN NUMERICAL * METHODS FOR NONLINEAR ALGEBRAIC EQUATIONS, EDITED BY * P. RABINOWITZ, GORDON AND BREACH, LONDON. * * *** HISTORY *** * * JOHN DENNIS DESIGNED MUCH OF THIS ROUTINE, STARTING WITH * IDEAS IN (2). ROY WELSCH SUGGESTED THE MODEL SWITCHING STRATEGY. * DAVID GAY AND STEPHEN PETERS CAST THIS SUBROUTINE INTO A MORE * PORTABLE FORM (WINTER 1977), AND DAVID GAY CAST IT INTO ITS * PRESENT FORM (FALL 1978). * * *** GENERAL *** * * THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH * SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS * MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND * MCS-7906671. * * ----------------------- EXTERNAL QUANTITIES ------------------------ * * *** NO EXTERNAL FUNCTIONS AND SUBROUTINES *** * * ------------------------- LOCAL VARIABLES -------------------------- * */ /* *** SUBSCRIPTS FOR IV AND V *** * */ /* *** DATA INITIALIZATIONS *** * */ /* ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ * */ nfc = Iv[NFCALL]; Iv[SWITCH_] = 0; Iv[RESTOR] = 0; rfac1 = ONE; goodx = TRUE; i = Iv[IRC]; if (i >= 1 && i <= 12) switch (i) { case 1: goto L_20; case 2: goto L_30; case 3: goto L_10; case 4: goto L_10; case 5: goto L_40; case 6: goto L_280; case 7: goto L_220; case 8: goto L_220; case 9: goto L_220; case 10: goto L_220; case 11: goto L_220; case 12: goto L_170; } Iv[IRC] = 13; goto L_999; /* *** INITIALIZE FOR NEW ITERATION *** * */ L_10: Iv[STAGE] = 1; Iv[RADINC] = 0; V[FLSTGD] = V[F0]; if (Iv[TOOBIG] == 0) goto L_110; Iv[STAGE] = -1; Iv[XIRC] = i; goto L_60; /* *** STEP WAS RECOMPUTED WITH NEW MODEL OR SMALLER RADIUS *** * *** FIRST DECIDE WHICH *** * */ L_20: if (Iv[MODEL] != Iv[MLSTGD]) goto L_30; /* *** OLD MODEL RETAINED, SMALLER RADIUS TRIED *** * *** DO NOT CONSIDER ANY MORE NEW MODELS THIS ITERATION *** */ Iv[STAGE] = Iv[STGLIM]; Iv[RADINC] = -1; goto L_110; /* *** A NEW MODEL IS BEING TRIED. DECIDE WHETHER TO KEEP IT. *** * */ L_30: Iv[STAGE] += 1; /* *** NOW WE ADD THE POSSIBILTIY THAT STEP WAS RECOMPUTED WITH *** * *** THE SAME MODEL, PERHAPS BECAUSE OF AN OVERSIZED STEP. *** * */ L_40: if (Iv[STAGE] > 0) goto L_50; /* *** STEP WAS RECOMPUTED BECAUSE IT WAS TOO BIG. *** * */ if (Iv[TOOBIG] != 0) goto L_60; /* *** RESTORE IV(STAGE) AND PICK UP WHERE WE LEFT OFF. *** * */ Iv[STAGE] = -Iv[STAGE]; i = Iv[XIRC]; switch (i) { case 1: goto L_20; case 2: goto L_30; case 3: goto L_110; case 4: goto L_110; case 5: goto L_70; } L_50: if (Iv[TOOBIG] == 0) goto L_70; /* *** HANDLE OVERSIZE STEP *** * */ if (Iv[RADINC] > 0) goto L_80; Iv[STAGE] = -Iv[STAGE]; Iv[XIRC] = Iv[IRC]; L_60: V[RADFAC] = V[DECFAC]; Iv[RADINC] -= 1; Iv[IRC] = 5; Iv[RESTOR] = 1; goto L_999; L_70: if (V[F] < V[FLSTGD]) goto L_110; /* *** THE NEW STEP IS A LOSER. RESTORE OLD MODEL. *** * */ if (Iv[MODEL] == Iv[MLSTGD]) goto L_80; Iv[MODEL] = Iv[MLSTGD]; Iv[SWITCH_] = 1; /* *** RESTORE STEP, ETC. ONLY IF A PREVIOUS STEP DECREASED V(F). * */ L_80: if (V[FLSTGD] >= V[F0]) goto L_110; Iv[RESTOR] = 1; V[F] = V[FLSTGD]; V[PREDUC] = V[PLSTGD]; V[GTSTEP] = V[GTSLST]; if (Iv[SWITCH_] == 0) rfac1 = V[DSTNRM]/V[DSTSAV]; V[DSTNRM] = V[DSTSAV]; nfc = Iv[NFGCAL]; goodx = FALSE; L_110: V[FDIF] = V[F0] - V[F]; if (V[FDIF] > V[TUNER2]*V[PREDUC]) goto L_140; if (Iv[RADINC] > 0) goto L_140; /* *** NO (OR ONLY A TRIVIAL) FUNCTION DECREASE * *** -- SO TRY NEW MODEL OR SMALLER RADIUS * */ if (V[F] < V[F0]) goto L_120; Iv[MLSTGD] = Iv[MODEL]; V[FLSTGD] = V[F]; V[F] = V[F0]; Iv[RESTOR] = 1; goto L_130; L_120: Iv[NFGCAL] = nfc; L_130: Iv[IRC] = 1; if (Iv[STAGE] < Iv[STGLIM]) goto L_160; Iv[IRC] = 5; Iv[RADINC] -= 1; goto L_160; /* *** NONTRIVIAL FUNCTION DECREASE ACHIEVED *** * */ L_140: Iv[NFGCAL] = nfc; rfac1 = ONE; V[DSTSAV] = V[DSTNRM]; if (V[FDIF] > V[PREDUC]*V[TUNER1]) goto L_190; /* *** DECREASE WAS MUCH LESS THAN PREDICTED -- EITHER CHANGE MODELS * *** OR ACCEPT STEP WITH DECREASED RADIUS. * */ if (Iv[STAGE] >= Iv[STGLIM]) goto L_150; /* *** CONSIDER SWITCHING MODELS *** */ Iv[IRC] = 2; goto L_160; /* *** ACCEPT STEP WITH DECREASED RADIUS *** * */ L_150: Iv[IRC] = 4; /* *** SET V(RADFAC) TO FLETCHER*S DECREASE FACTOR *** * */ L_160: Iv[XIRC] = Iv[IRC]; emax = V[GTSTEP] + V[FDIF]; V[RADFAC] = HALF*rfac1; if (emax < V[GTSTEP]) V[RADFAC] = rfac1*fmax( V[RDFCMN], HALF*V[GTSTEP]/emax ); /* *** DO FALSE CONVERGENCE TEST *** * */ L_170: if (V[RELDX] <= V[XFTOL]) goto L_180; Iv[IRC] = Iv[XIRC]; if (V[F] < V[F0]) goto L_200; goto L_230; L_180: Iv[IRC] = 12; goto L_240; /* *** HANDLE GOOD FUNCTION DECREASE *** * */ L_190: if (V[FDIF] < (-V[TUNER3]*V[GTSTEP])) goto L_210; /* *** INCREASING RADIUS LOOKS WORTHWHILE. SEE IF WE JUST * *** RECOMPUTED STEP WITH A DECREASED RADIUS OR RESTORED STEP * *** AFTER RECOMPUTING IT WITH A LARGER RADIUS. * */ if (Iv[RADINC] < 0) goto L_210; if (Iv[RESTOR] == 1) goto L_210; /* *** WE DID NOT. TRY A LONGER STEP UNLESS THIS WAS A NEWTON * *** STEP. * */ V[RADFAC] = V[RDFCMX]; gts = V[GTSTEP]; if (V[FDIF] < (HALF/V[RADFAC] - ONE)*gts) V[RADFAC] = fmax( V[INCFAC], HALF*gts/(gts + V[FDIF]) ); Iv[IRC] = 4; if (V[STPPAR] == ZERO) goto L_230; if (V[DST0] >= ZERO && (V[DST0] < TWO*V[DSTNRM] || V[NREDUC] < ONEP2*V[FDIF])) goto L_230; /* *** STEP WAS NOT A NEWTON STEP. RECOMPUTE IT WITH * *** A LARGER RADIUS. */ Iv[IRC] = 5; Iv[RADINC] += 1; /* *** SAVE VALUES CORRESPONDING TO GOOD STEP *** * */ L_200: V[FLSTGD] = V[F]; Iv[MLSTGD] = Iv[MODEL]; if (Iv[RESTOR] != 1) Iv[RESTOR] = 2; V[DSTSAV] = V[DSTNRM]; Iv[NFGCAL] = nfc; V[PLSTGD] = V[PREDUC]; V[GTSLST] = V[GTSTEP]; goto L_230; /* *** ACCEPT STEP WITH RADIUS UNCHANGED *** * */ L_210: V[RADFAC] = ONE; Iv[IRC] = 3; goto L_230; /* *** COME HERE FOR A RESTART AFTER CONVERGENCE *** * */ L_220: Iv[IRC] = Iv[XIRC]; if (V[DSTSAV] >= ZERO) goto L_240; Iv[IRC] = 12; goto L_240; /* *** PERFORM CONVERGENCE TESTS *** * */ L_230: Iv[XIRC] = Iv[IRC]; L_240: if (Iv[RESTOR] == 1 && V[FLSTGD] < V[F0]) Iv[RESTOR] = 3; if (fabs( V[F] ) < V[AFCTOL]) Iv[IRC] = 10; if (HALF*V[FDIF] > V[PREDUC]) goto L_999; emax = V[RFCTOL]*fabs( V[F0] ); emaxs = V[SCTOL]*fabs( V[F0] ); if (V[PREDUC] <= emaxs && (V[DSTNRM] > V[LMAXS] || V[STPPAR] == ZERO)) Iv[IRC] = 11; if (V[DST0] < ZERO) goto L_250; i = 0; if ((V[NREDUC] > ZERO && V[NREDUC] <= emax) || (V[NREDUC] == ZERO && V[PREDUC] == ZERO)) i = 2; if ((V[STPPAR] == ZERO && V[RELDX] <= V[XCTOL]) && goodx) i += 1; if (i > 0) Iv[IRC] = i + 6; /* *** CONSIDER RECOMPUTING STEP OF LENGTH V(LMAXS) FOR SINGULAR * *** CONVERGENCE TEST. * */ L_250: if (Iv[IRC] > 5 && Iv[IRC] != 12) goto L_999; if (V[STPPAR] == ZERO) goto L_999; if (V[DSTNRM] > V[LMAXS]) goto L_260; if (V[PREDUC] >= emaxs) goto L_999; if (V[DST0] <= ZERO) goto L_270; if (HALF*V[DST0] <= V[LMAXS]) goto L_999; goto L_270; L_260: if (HALF*V[DSTNRM] <= V[LMAXS]) goto L_999; xmax = V[LMAXS]/V[DSTNRM]; if (xmax*(TWO - xmax)*V[PREDUC] >= emaxs) goto L_999; L_270: if (V[NREDUC] < ZERO) goto L_290; /* *** RECOMPUTE V(PREDUC) FOR USE IN SINGULAR CONVERGENCE TEST *** * */ V[GTSLST] = V[GTSTEP]; V[DSTSAV] = V[DSTNRM]; if (Iv[IRC] == 12) V[DSTSAV] = -V[DSTSAV]; V[PLSTGD] = V[PREDUC]; i = Iv[RESTOR]; Iv[RESTOR] = 2; if (i == 3) Iv[RESTOR] = 0; Iv[IRC] = 6; goto L_999; /* *** PERFORM SINGULAR CONVERGENCE TEST WITH RECOMPUTED V(PREDUC) *** * */ L_280: V[GTSTEP] = V[GTSLST]; V[DSTNRM] = fabs( V[DSTSAV] ); Iv[IRC] = Iv[XIRC]; if (V[DSTSAV] <= ZERO) Iv[IRC] = 12; V[NREDUC] = -V[PREDUC]; V[PREDUC] = V[PLSTGD]; Iv[RESTOR] = 3; L_290: if (-V[NREDUC] <= V[RFCTOL]*fabs( V[F0] )) Iv[IRC] = 11; L_999: return; /* *** LAST CARD OF DA7SST FOLLOWS *** */ } /* end of function */ /* ================================================================== */ /* PARAMETER translations */ #define NEEDHD 36 #define NGCALL 30 #define PRNTIT 39 #define SUSED 64 /* end of PARAMETER translations */ void /*FUNCTION*/ ditsum( double d[], double g[], long iv[], long liv, long lv, long p, double v[], double x[]) { long int alg, i, iv1, m, nf, ng, ol, pu; double nreldf, oldf, preldf, reldf; static char model1[6][5]={" "," "," "," "," G "," S "}; static char model2[6][5]={" G "," S ","G-S ","S-G ","-S-G","-G-S"}; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const D = &d[0] - 1; double *const G = &g[0] - 1; long *const Iv = &iv[0] - 1; double *const V = &v[0] - 1; double *const X = &x[0] - 1; /* end of OFFSET VECTORS */ /* *** PRINT ITERATION SUMMARY FOR ***SOL (VERSION 2.3) *** * * 6/28/90 CLL Added test before the GO TO (...), IV1 * to avoid printing the termination message when IV1 = 1 to 4 and * no other printing has been requested. * *** PARAMETER DECLARATIONS *** * */ /* ------------------------------------------------------------------ * *** LOCAL VARIABLES *** * */ /* *** NO EXTERNAL FUNCTIONS OR SUBROUTINES *** * * *** SUBSCRIPTS FOR IV AND V *** * */ /* *** IV SUBSCRIPT VALUES *** * */ /* *** V SUBSCRIPT VALUES *** * */ /* ------------------------------ BODY -------------------------------- * */ /*340 FORMAT(/' ***** BAD PARAMETERS TO ASSESS *****') */ pu = Iv[PRUNIT]; if (pu == 0) goto L_999; iv1 = Iv[1]; if (iv1 > 62) iv1 -= 51; ol = Iv[OUTLEV]; alg = ((Iv[ALGSAV] - 1)%2) + 1; if (iv1 < 2 || iv1 > 15) goto L_370; if (iv1 >= 12) goto L_120; if (iv1 == 2 && Iv[NITER] == 0) goto L_390; if (ol == 0) goto L_120; if (iv1 >= 10 && Iv[PRNTIT] == 0) goto L_120; if (iv1 > 2) goto L_10; Iv[PRNTIT] += 1; if (Iv[PRNTIT] < labs( ol )) goto L_999; L_10: nf = Iv[NFCALL] - labs( Iv[NFCOV] ); Iv[PRNTIT] = 0; reldf = ZERO; preldf = ZERO; oldf = fmax( fabs( V[F0] ), fabs( V[F] ) ); if (oldf <= ZERO) goto L_20; reldf = V[FDIF]/oldf; preldf = V[PREDUC]/oldf; L_20: if (ol > 0) goto L_60; /* *** PRINT SHORT SUMMARY LINE *** * */ if (Iv[NEEDHD] == 1 && alg == 1) { printf("\n IT NF F RELDF PRELDF RELDX MODEL STPPAR\n"); } if (Iv[NEEDHD] == 1 && alg == 2) { printf("\n IT NF F RELDF PRELDF RELDX STPPAR\n"); } Iv[NEEDHD] = 0; if (alg == 2) goto L_50; m = Iv[SUSED]; printf("%6ld%5ld%10.3g%9.2g%9.2g%8.1g%3.3s%4.4s%8.1g\n", Iv[NITER], nf, V[F], reldf, preldf, V[RELDX], model1[m - 1], model2[m - 1], V[STPPAR]); goto L_120; L_50: printf("%6ld%5ld%11.3g%10.2g%10.2g%9.1g%9.1g\n", Iv[NITER], nf, V[F], reldf, preldf, V[RELDX], V[STPPAR]); goto L_120; /* *** PRINT LONG SUMMARY LINE *** * */ L_60: if (Iv[NEEDHD] == 1 && alg == 1) { printf("\n IT NF F RELDF PRELDF RELDX MODEL STPPAR D*STEP NPRELDF\n"); } if (Iv[NEEDHD] == 1 && alg == 2) { printf("\n IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF\n"); } Iv[NEEDHD] = 0; nreldf = ZERO; if (oldf > ZERO) nreldf = V[NREDUC]/oldf; if (alg == 2) goto L_90; m = Iv[SUSED]; printf("%6ld%5ld%10.3g%9.2g%9.2g%8.1g%3.3s%4.4s%8.1g%8.1g%9.2g\n", Iv[NITER], nf, V[F], reldf, preldf, V[RELDX], model1[m - 1], model2[m - 1], V[STPPAR], V[DSTNRM], nreldf); goto L_120; L_90: printf("%6ld%5ld%11.3g%10.2g%10.2g%9.1g%9.1g%9.1g%10.2g\n", Iv[NITER], nf, V[F], reldf, preldf, V[RELDX], V[STPPAR], V[DSTNRM], nreldf); L_120: if (iv1 <= 2) goto L_999; i = Iv[STATPR]; if (i == (-1)) goto L_460; if (i + iv1 < 0) goto L_460; /* The following test skips printing the termination message * when convergence is satisfactory and no other printing has been * requested. * */ if (((((((iv1 >= 3 && iv1 <= 6) && Iv[COVPRT] == 0) && Iv[OUTLEV] == 0) && Iv[PARPRT] == 0) && Iv[SOLPRT] == 0) && Iv[STATPR] == 0) && Iv[X0PRT] == 0) goto L_430; switch (iv1) { case 1: goto L_999; case 2: goto L_999; case 3: goto L_130; case 4: goto L_150; case 5: goto L_170; case 6: goto L_190; case 7: goto L_210; case 8: goto L_230; case 9: goto L_250; case 10: goto L_270; case 11: goto L_290; case 12: goto L_310; case 13: goto L_330; case 14: goto L_350; case 15: goto L_500; } L_130: printf("\n ***** COEFFICIENT CONVERGENCE *****\n"); goto L_430; L_150: printf("\n ***** RELATIVE FUNCTION CONVERGENCE *****\n"); goto L_430; L_170: printf(" \n ***** COEFFICIENT AND RELATIVE FUNCTION CONVERGENCE *****\n"); goto L_430; L_190: printf("\n ***** ABSOLUTE FUNCTION CONVERGENCE *****\n"); goto L_430; L_210: printf("\n ***** SINGULAR CONVERGENCE *****\n"); goto L_430; L_230: printf("\n ***** FALSE CONVERGENCE *****\n"); goto L_430; L_250: printf("\n ***** FUNCTION EVALUATION LIMIT *****\n"); goto L_430; L_270: printf("\n ***** ITERATION LIMIT *****\n"); goto L_430; L_290: printf("\n ***** STOPX *****\n"); goto L_430; L_310: printf("\n ***** INITIAL F(X) CANNOT BE COMPUTED *****\n"); goto L_390; L_330: printf("\n ***** BAD PARAMETERS (Internal error) *****\n"); goto L_999; L_350: printf("\n ***** GRADIENT COULD NOT BE COMPUTED *****\n"); if (Iv[NITER] > 0) goto L_460; goto L_390; L_370: printf("\n ***** IV(1) =%5ld *****\n", Iv[1]); goto L_999; /* *** INITIAL CALL ON DITSUM *** * */ L_390: if (Iv[X0PRT] != 0) { printf("\n I INITIAL X(I) D(I)\n\n"); for (i = 1; i <= p; i++) { printf(" %5ld%17.6g%14.3g\n", i, X[i], D[i]); } } /* *** THE FOLLOWING ARE TO AVOID UNDEFINED VARIABLES WHEN THE * *** FUNCTION EVALUATION LIMIT IS 1... */ V[DSTNRM] = ZERO; V[FDIF] = ZERO; V[NREDUC] = ZERO; V[PREDUC] = ZERO; V[RELDX] = ZERO; if (iv1 >= 12) goto L_999; Iv[NEEDHD] = 0; Iv[PRNTIT] = 0; if (ol == 0) goto L_999; if (ol < 0 && alg == 1) { printf("\n IT NF F RELDF PRELDF RELDX MODEL STPPAR\n"); } if (ol < 0 && alg == 2) { printf("\n IT NF F RELDF PRELDF RELDX STPPAR\n"); } if (ol > 0 && alg == 1) { printf("\n IT NF F RELDF PRELDF RELDX MODEL STPPAR D*STEP NPRELDF\n"); } if (ol > 0 && alg == 2) { printf("\n IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF\n"); } if (alg == 1) { printf("\n 0%5ld%10.3g\n", Iv[NFCALL], V[F]); } if (alg == 2) { printf("\n 0%5ld%11.3g\n", Iv[NFCALL], V[F]); } goto L_999; /* *** PRINT VARIOUS INFORMATION REQUESTED ON SOLUTION *** * */ L_430: Iv[NEEDHD] = 1; if (Iv[STATPR] <= 0) goto L_460; oldf = fmax( fabs( V[F0] ), fabs( V[F] ) ); preldf = ZERO; nreldf = ZERO; if (oldf <= ZERO) goto L_440; preldf = V[PREDUC]/oldf; nreldf = V[NREDUC]/oldf; L_440: nf = Iv[NFCALL] - Iv[NFCOV]; ng = Iv[NGCALL] - Iv[NGCOV]; printf("\n FUNCTION%17.6g RELDX%17.3g\n FUNC. EVALS%8ld GRAD. EVALS%8ld\n PRELDF%16.3g NPRELDF%15.3g\n", V[F], V[RELDX], nf, ng, preldf, nreldf); L_460: if (Iv[SOLPRT] == 0) goto L_999; Iv[NEEDHD] = 1; if (Iv[ALGSAV] > 2) goto L_999; printf("\n I FINAL X(I) D(I) G(I)\n\n"); for (i = 1; i <= p; i++) { printf(" %5ld%16.6g%14.3g%14.3g\n", i, X[i], D[i], G[i]); } goto L_999; L_500: printf("\n INCONSISTENT DIMENSIONS\n"); L_999: return; /* *** LAST CARD OF DITSUM FOLLOWS *** */ } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dl7itv( long n, double x[], double l[], double y[]) { long int i, i0, ii, ij, im1, j, np1; double xi; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const L = &l[0] - 1; double *const X = &x[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* *** SOLVE (L**T)*X = Y, WHERE L IS AN N X N LOWER TRIANGULAR * *** MATRIX STORED COMPACTLY BY ROWS. X AND Y MAY OCCUPY THE SAME * *** STORAGE. *** * ------------------------------------------------------------------ */ for (i = 1; i <= n; i++) { X[i] = Y[i]; } np1 = n + 1; i0 = n*(n + 1)/2; for (ii = 1; ii <= n; ii++) { i = np1 - ii; xi = X[i]/L[i0]; X[i] = xi; if (i <= 1) goto L_999; i0 -= i; if (xi == ZERO) goto L_30; im1 = i - 1; for (j = 1; j <= im1; j++) { ij = i0 + j; X[j] += -xi*L[ij]; } L_30: ; } L_999: return; /* *** LAST CARD OF DL7ITV FOLLOWS *** */ } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dl7sqr( long n, double a[], double l[]) { long int i, i0, ii, ij, ik, ip1, j, j0, jj, jk, k, np1; double t; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const A = &a[0] - 1; double *const L = &l[0] - 1; /* end of OFFSET VECTORS */ /* *** COMPUTE A = LOWER TRIANGLE OF L*(L**T), WITH BOTH * *** L AND A STORED COMPACTLY BY ROWS. (BOTH MAY OCCUPY THE * *** SAME STORAGE. * * *** PARAMETERS *** * * ------------------------------------------------------------------ */ /* DIMENSION A(N*(N+1)/2), L(N*(N+1)/2) * * *** LOCAL VARIABLES *** * */ np1 = n + 1; i0 = n*(n + 1)/2; for (ii = 1; ii <= n; ii++) { i = np1 - ii; ip1 = i + 1; i0 -= i; j0 = i*(i + 1)/2; for (jj = 1; jj <= i; jj++) { j = ip1 - jj; j0 -= j; t = 0.0e0; for (k = 1; k <= j; k++) { ik = i0 + k; jk = j0 + k; t += L[ik]*L[jk]; } ij = i0 + j; A[ij] = t; } } return; } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dl7srt( long n1, long n, double l[], double a[], long *irc) { long int i, i0, ij, ik, im1, j, j0, jk, jm1, k; double t, td; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const A = &a[0] - 1; double *const L = &l[0] - 1; /* end of OFFSET VECTORS */ /* *** COMPUTE ROWS N1 THROUGH N OF THE CHOLESKY FACTOR L OF * *** A = L*(L**T), WHERE L AND THE LOWER TRIANGLE OF A ARE BOTH * *** STORED COMPACTLY BY ROWS (AND MAY OCCUPY THE SAME STORAGE). * *** IRC = 0 MEANS ALL WENT WELL. IRC = J MEANS THE LEADING * *** PRINCIPAL J X J SUBMATRIX OF A IS NOT POSITIVE DEFINITE -- * *** AND L(J*(J+1)/2) CONTAINS THE (NONPOS.) REDUCED J-TH DIAGONAL. * * *** PARAMETERS *** * * ------------------------------------------------------------------ */ /* DIMENSION L(N*(N+1)/2), A(N*(N+1)/2) * * *** LOCAL VARIABLES *** * */ /* *** BODY *** * */ i0 = n1*(n1 - 1)/2; for (i = n1; i <= n; i++) { td = ZERO; if (i == 1) goto L_40; j0 = 0; im1 = i - 1; for (j = 1; j <= im1; j++) { t = ZERO; if (j == 1) goto L_20; jm1 = j - 1; for (k = 1; k <= jm1; k++) { ik = i0 + k; jk = j0 + k; t += L[ik]*L[jk]; } L_20: ij = i0 + j; j0 += j; t = (A[ij] - t)/L[j0]; L[ij] = t; td += t*t; } L_40: i0 += i; t = A[i0] - td; if (t <= ZERO) goto L_60; L[i0] = sqrt( t ); /* DNSGB (9 of 10) */ } *irc = 0; goto L_999; L_60: L[i0] = t; *irc = i; L_999: return; /* *** LAST CARD OF DL7SRT *** */ } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dl7tvm( long n, double x[], double l[], double y[]) { long int i, i0, ij, j; double yi; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const L = &l[0] - 1; double *const X = &x[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* *** COMPUTE X = (L**T)*Y, WHERE L IS AN N X N LOWER * *** TRIANGULAR MATRIX STORED COMPACTLY BY ROWS. X AND Y MAY * *** OCCUPY THE SAME STORAGE. *** * ------------------------------------------------------------------ */ /* DIMENSION L(N*(N+1)/2) */ i0 = 0; for (i = 1; i <= n; i++) { yi = Y[i]; X[i] = ZERO; for (j = 1; j <= i; j++) { ij = i0 + j; X[j] += yi*L[ij]; } i0 += i; } /* *** LAST CARD OF DL7TVM FOLLOWS *** */ return; } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dl7vml( long n, double x[], double l[], double y[]) { long int i, i0, ii, ij, j, np1; double t; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const L = &l[0] - 1; double *const X = &x[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* *** COMPUTE X = L*Y, WHERE L IS AN N X N LOWER TRIANGULAR * *** MATRIX STORED COMPACTLY BY ROWS. X AND Y MAY OCCUPY THE SAME * *** STORAGE. *** * * ------------------------------------------------------------------ */ /* DIMENSION L(N*(N+1)/2) */ np1 = n + 1; i0 = n*(n + 1)/2; for (ii = 1; ii <= n; ii++) { i = np1 - ii; i0 -= i; t = ZERO; for (j = 1; j <= i; j++) { ij = i0 + j; t += L[ij]*Y[j]; } X[i] = t; } /* *** LAST CARD OF DL7VML FOLLOWS *** */ return; } /* end of function */ /* ================================================================== */ double /*FUNCTION*/ drldst( long p, double d[], double x[], double x0[]) { long int i; double drldst_v, emax, t, xmax; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const D = &d[0] - 1; double *const X = &x[0] - 1; double *const X0 = &x0[0] - 1; /* end of OFFSET VECTORS */ /* *** COMPUTE AND RETURN RELATIVE DIFFERENCE BETWEEN X AND X0 *** * *** NL2SOL VERSION 2.2 *** * ------------------------------------------------------------------ */ /* *** BODY *** * */ emax = ZERO; xmax = ZERO; for (i = 1; i <= p; i++) { t = fabs( D[i]*(X[i] - X0[i]) ); if (emax < t) emax = t; t = D[i]*(fabs( X[i] ) + fabs( X0[i] )); if (xmax < t) xmax = t; } drldst_v = ZERO; if (xmax > ZERO) drldst_v = emax/xmax; /* *** LAST CARD OF DRLDST FOLLOWS *** */ return( drldst_v ); } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dv2axy( long p, double w[], double a, double x[], double y[]) { long int i; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const W = &w[0] - 1; double *const X = &x[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* *** SET W = A*X + Y -- W, X, Y = P-VECTORS, A = SCALAR *** * * ------------------------------------------------------------------ */ for (i = 1; i <= p; i++) { W[i] = a*X[i] + Y[i]; } return; } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dv7cpy( long p, double y[], double x[]) { long int i; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const X = &x[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* *** SET Y = X, WHERE X AND Y ARE P-VECTORS *** * * ------------------------------------------------------------------ */ for (i = 1; i <= p; i++) { Y[i] = X[i]; } return; } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dv7scl( long n, double x[], double a, double y[]) { long int i; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const X = &x[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* *** SET X(I) = A*Y(I), I = 1(1)N *** * * ------------------------------------------------------------------ */ for (i = 1; i <= n; i++) { X[i] = a*Y[i]; } /* *** LAST LINE OF DV7SCL FOLLOWS *** */ return; } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dv7scp( long p, double y[], double ss) { long int i; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* *** SET P-VECTOR Y TO SCALAR SS *** * * ------------------------------------------------------------------ */ for (i = 1; i <= p; i++) { Y[i] = ss; } return; } /* end of function */