/*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 <math.h>
#include "fcrt.h"
#include "divset.h"
#include <string.h>
#include <float.h>
#include <stdio.h>
#include <stdlib.h>
/* 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 <string.h>
/* 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 */