/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:55 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "srn2g.h"
#include <stdio.h>
#include <stdlib.h>
/* PARAMETER translations */
#define CNVCOD 55
#define COVMAT 26
#define COVREQ 15
#define D0INIT 40
#define DINIT 38
#define DTINIT 39
#define DTYPE 16
#define F 10
#define FDH 74
#define G 28
#define H 56
#define HALF 0.5e0
#define IPIVOT 76
#define IVNEED 3
#define JCN 66
#define JTOL 59
#define LMAT 42
#define MODE 35
#define NEXTIV 46
#define NEXTV 47
#define NF0 68
#define NF00 81
#define NF1 69
#define NFCALL 6
#define NFCOV 52
#define NFGCAL 7
#define NGCALL 30
#define NGCOV 53
#define QTR 77
#define RDREQ 57
#define REGD 67
#define RESTOR 9
#define RLIMIT 46
#define RMAT 78
#define TOOBIG 2
#define VNEED 4
#define Y 48
#define ZERO 0.e0
/* end of PARAMETER translations */
void /*FUNCTION*/ srn2g(
float d[],
float *dr,
long iv[],
long liv,
long lv,
long n,
long nd,
long *n1,
long *n2,
long p,
float r[],
float rd[],
float v[],
float x[])
{
#define DR(I_,J_) (*(dr+(I_)*(nd)+(J_)))
long int g1, gi, i, iv1, ivmode, jtol1, k, l, lh, nn, qtr1, rmat1,
y1, yi;
float t;
/* OFFSET Vectors w/subscript range: 1 to dimension */
float *const D = &d[0] - 1;
long *const Iv = &iv[0] - 1;
float *const R = &r[0] - 1;
float *const Rd = &rd[0] - 1;
float *const V = &v[0] - 1;
float *const X = &x[0] - 1;
/* end of OFFSET VECTORS */
/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
* ALL RIGHTS RESERVED.
* Based on Government Sponsored Research NAS7-03001.
* File: SRN2G.for Ten subrs used by the
* David Gay & Linda Kaufman nonlinear LS package.
* Needed for versions that do not allow Bounded variables.
* SRN2G is called by SNLAFU, SNLAGU, & SRNSG.
*
*>> 2015-07-09 SRN2G Krogh Introduced TP to avoid divide by 0,
*>> 2000-01-07 SRN2G Krogh Moved COV1 = IV(COMAT) up in SN2CVP.
*>> 1998-10-29 SRN2G Krogh Moved external statement up for mangle.
*>> 1996-07-09 SRN2G Krogh Changes for conversion to C.
*>> 1995-01-26 SRN2G Krogh Moved formats up for C conversion.
*>> 1994-11-02 SRN2G Krogh Changes to use M77CON
*>> 1993-03-10 SRN2G CLL Moved stmt NN = ... to follow IF (IV1 ...
*>> 1992-04-27 CLL Comment out unreferenced stmt labels.
*>> 1992-04-13 CLL Change from Hollerith to '...' syntax in formats.
*>> 1990-06-29 CLL Changes to formats in SN2CVP.
*>> 1990-06-12 CLL Revised SRN2G & SG7LIT from DMG 4/19/90
*>> 1990-03-30 CLL JPL
*>> 1990-03-14 CLL JPL
*>> 1990-06-12 CLL
*>> 1990-04-23 CLL (Recent revision by DMG)
*** from netlib, Thu Apr 19 11:58:57 EDT 1990 ***
*--S replaces "?": ?RN2G,?C7VFN,?D7TPR,?D7UPD,?G7LIT,?ITSUM,?IVSET
*--& ?L7VML,?N2CVP,?N2LRD,?Q7APL,?Q7RAD,?V2NRM,?V7CPY,?V7SCP, ?N2G
*--& ?A7SST,?F7HES,?G7QTS,?L7MST,?L7SQR,?L7SRT,?L7SVN,?L7SVX,?L7TVM
*--& ?PARCK,?R7MDC,?RLDST,?S7LUP,?S7LVM,?V2AXY,?L7ITV,?L7IVM,?O7PRD
*--& ?L7NVR,?L7TSQ,?V7SCL,?N2RDP,?NLAFU,?NLAGU,?RNSG
*
* *** REVISED ITERATION DRIVER FOR NL2SOL (VERSION 2.3) ***
* */
/* ------------------------- PARAMETER USAGE --------------------------
*
* D........ SCALE VECTOR.
* DR....... DERIVATIVES OF R AT X.
* IV....... INTEGER VALUES ARRAY.
* LIV...... LENGTH OF IV... LIV MUST BE AT LEAST P + 82.
* LV....... LENGTH OF V... LV MUST BE AT LEAST 105 + P*(2*P+16).
* N........ TOTAL NUMBER OF RESIDUALS.
* ND....... MAX. NO. OF RESIDUALS PASSED ON ONE CALL.
* N1....... LOWEST ROW INDEX FOR RESIDUALS SUPPLIED THIS TIME.
* N2....... HIGHEST ROW INDEX FOR RESIDUALS SUPPLIED THIS TIME.
* P........ NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED.
* R........ RESIDUALS.
* RD....... RD(I) = SQRT(G(I)**T * H(I)**-1 * G(I)) ON OUTPUT WHEN
* IV(RDREQ) IS NONZERO. SRN2G SETS IV(REGD) = 1 IF RD
* IS SUCCESSFULLY COMPUTED, TO 0 IF NO ATTEMPT WAS MADE
* TO COMPUTE IT, AND TO -1 IF H (THE FINITE-DIFFERENCE HESSIAN)
* WAS INDEFINITE. IF ND .GE. N, THEN RD IS ALSO USED AS
* TEMPORARY STORAGE.
* V........ FLOATING-POINT VALUES ARRAY.
* X........ PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS,
* OUTPUT = BEST VALUE FOUND).
*
* *** DISCUSSION ***
*
* NOTE... NL2SOL AND NL2ITR (MENTIONED BELOW) ARE DESCRIBED IN
* ACM TRANS. MATH. SOFTWARE, VOL. 7, PP. 369-383 (AN ADAPTIVE
* NONLINEAR LEAST-SQUARES ALGORITHM, BY J.E. DENNIS, D.M. GAY,
* AND R.E. WELSCH).
*
* THIS ROUTINE CARRIES OUT ITERATIONS FOR SOLVING NONLINEAR
* LEAST SQUARES PROBLEMS. WHEN ND = N, IT IS SIMILAR TO NL2ITR
* (WITH J = DR), EXCEPT THAT R(X) AND DR(X) NEED NOT BE INITIALIZED
* WHEN SRN2G IS CALLED WITH IV(1) = 0 OR 12. SRN2G ALSO ALLOWS
* R AND DR TO BE SUPPLIED ROW-WISE -- JUST SET ND = 1 AND CALL
* SRN2G ONCE FOR EACH ROW WHEN PROVIDING RESIDUALS AND JACOBIANS.
* ANOTHER NEW FEATURE IS THAT CALLING SRN2G WITH IV(1) = 13
* CAUSES STORAGE ALLOCATION ONLY TO BE PERFORMED -- ON RETURN, SUCH
* COMPONENTS AS IV(G) (THE FIRST SUBSCRIPT IN G OF THE GRADIENT)
* AND IV(S) (THE FIRST SUBSCRIPT IN V OF THE S LOWER TRIANGLE OF
* THE S MATRIX) WILL HAVE BEEN SET (UNLESS LIV OR LV IS TOO SMALL),
* AND IV(1) WILL HAVE BEEN SET TO 14. CALLING SRN2G WITH IV(1) = 14
* CAUSES EXECUTION OF THE ALGORITHM TO BEGIN UNDER THE ASSUMPTION
* THAT STORAGE HAS BEEN ALLOCATED.
*
* *** SUPPLYING R AND DR ***
*
* SRN2G USES IV AND V IN THE SAME WAY AS NL2SOL, WITH A SMALL
* NUMBER OF OBVIOUS CHANGES. ONE DIFFERENCE BETWEEN SRN2G AND
* NL2ITR IS THAT INITIAL FUNCTION AND GRADIENT INFORMATION NEED NOT
* BE SUPPLIED IN THE VERY FIRST CALL ON SRN2G, THE ONE WITH
* IV(1) = 0 OR 12. ANOTHER DIFFERENCE IS THAT SRN2G RETURNS WITH
* IV(1) = -2 WHEN IT WANTS ANOTHER LOOK AT THE OLD JACOBIAN MATRIX
* AND THE CURRENT RESIDUAL -- THE ONE CORRESPONDING TO X AND
* IV(NFGCAL). IT THEN RETURNS WITH IV(1) = -3 WHEN IT WANTS TO SEE
* BOTH THE NEW RESIDUAL AND THE NEW JACOBIAN MATRIX AT ONCE. NOTE
* THAT IV(NFGCAL) = IV(7) CONTAINS THE VALUE THAT IV(NFCALL) = IV(6)
* HAD WHEN THE CURRENT RESIDUAL WAS EVALUATED. ALSO NOTE THAT THE
* VALUE OF X CORRESPONDING TO THE OLD JACOBIAN MATRIX IS STORED IN
* V, STARTING AT V(IV(X0)) = V(IV(43)).
* ANOTHER NEW RETURN... SRN2G IV(1) = -1 WHEN IT WANTS BOTH THE
* RESIDUAL AND THE JACOBIAN TO BE EVALUATED AT X.
* A NEW RESIDUAL VECTOR MUST BE SUPPLIED WHEN SRN2G RETURNS WITH
* IV(1) = 1 OR -1. THIS TAKES THE FORM OF VALUES OF R(I,X) PASSED
* IN R(I-N1+1), I = N1(1)N2. YOU MAY PASS ALL THESE VALUES AT ONCE
* (I.E., N1 = 1 AND N2 = N) OR IN PIECES BY MAKING SEVERAL CALLS ON
* SRN2G. EACH TIME SRN2G RETURNS WITH IV(1) = 1, N1 WILL HAVE
* BEEN SET TO THE INDEX OF THE NEXT RESIDUAL THAT SRN2G EXPECTS TO
* SEE, AND N2 WILL BE SET TO THE INDEX OF THE HIGHEST RESIDUAL THAT
* COULD BE GIVEN ON THE NEXT CALL, I.E., N2 = N1 + ND - 1. (THUS
* WHEN SRN2G FIRST RETURNS WITH IV(1) = 1 FOR A NEW X, IT WILL
* HAVE SET N1 TO 1 AND N2 TO MIN(ND,N).) THE CALLER MAY PROVIDE
* FEWER THAN N2-N1+1 RESIDUALS ON THE NEXT CALL BY SETTING N2 TO
* A SMALLER VALUE. SRN2G ASSUMES IT HAS SEEN ALL THE RESIDUALS
* FOR THE CURRENT X WHEN IT IS CALLED WITH N2 .GE. N.
* EXAMPLE... SUPPOSE N = 80 AND THAT R IS TO BE PASSED IN 8
* BLOCKS OF SIZE 10. THE FOLLOWING CODE WOULD DO THE JOB.
*
* N = 80
* ND = 10
* ...
* DO 10 K = 1, 8
* *** COMPUTE R(I,X) FOR I = 10*K-9 TO 10*K ***
* *** AND STORE THEM IN R(1),...,R(10) ***
* CALL SRN2G(..., R, ...)
* 10 CONTINUE
*
* THE SITUATION IS SIMILAR WHEN GRADIENT INFORMATION IS
* REQUIRED, I.E., WHEN SRN2G RETURNS WITH IV(1) = 2, -1, OR -2.
* NOTE THAT SRN2G OVERWRITES R, BUT THAT IN THE SPECIAL CASE OF
* N1 = 1 AND N2 = N ON PREVIOUS CALLS, SRN2G NEVER RETURNS WITH
* IV(1) = -2. IT SHOULD BE CLEAR THAT THE PARTIAL DERIVATIVE OF
* R(I,X) WITH RESPECT TO X(L) IS TO BE STORED IN DR(I-N1+1,L),
* L = 1(1)P, I = N1(1)N2. IT IS ESSENTIAL THAT R(I) AND DR(I,L)
* ALL CORRESPOND TO THE SAME RESIDUALS WHEN IV(1) = -1 OR -2.
*
* *** COVARIANCE MATRIX ***
*
* IV(RDREQ) = IV(57) TELLS WHETHER TO COMPUTE A COVARIANCE */
/* MATRIX AND/OR REGRESSION DIAGNOSTICS... 0 MEANS NEITHER,
* 1 MEANS COVARIANCE MATRIX ONLY, 2 MEANS REG. DIAGNOSTICS ONLY,
* 3 MEANS BOTH. AS WITH NL2SOL, IV(COVREQ) = IV(15) TELLS WHAT
* HESSIAN APPROXIMATION TO USE IN THIS COMPUTING.
*
* *** REGRESSION DIAGNOSTICS ***
*
* SEE THE COMMENTS IN SUBROUTINE SN2G.
*
* *** GENERAL ***
*
* CODED BY DAVID M. GAY.
*
* ++++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++
*
* *** EXTERNAL FUNCTIONS AND SUBROUTINES ***
* */
/* ------------------------------------------------------------------
* SC7VFN... FINISHES COVARIANCE COMPUTATION.
* SIVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS.
* SD7TPR... COMPUTES INNER PRODUCT OF TWO VECTORS.
* SD7UPD... UPDATES SCALE VECTOR D.
* SG7LIT.... PERFORMS BASIC MINIMIZATION ALGORITHM.
* SITSUM.... PRINTS ITERATION SUMMARY, INFO ABOUT INITIAL AND FINAL X.
* SL7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX.
* SN2CVP... PRINTS COVARIANCE MATRIX.
* SN2LRD... COMPUTES REGRESSION DIAGNOSTICS.
* SQ7APL... APPLIES QR TRANSFORMATIONS STORED BY SQ7RAD.
* SQ7RAD.... ADDS A NEW BLOCK OF ROWS TO QR DECOMPOSITION.
* SV7CPY.... COPIES ONE VECTOR TO ANOTHER.
* SV7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR.
*
* *** LOCAL VARIABLES ***
* */
/* *** SUBSCRIPTS FOR IV AND V ***
* */
/* *** IV SUBSCRIPT VALUES ***
* */
/* *** V SUBSCRIPT VALUES *** */
/* ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++
* */
lh = p*(p + 1)/2;
if (Iv[1] == 0)
sivset( 1, iv, liv, lv, v );
iv1 = Iv[1];
if (iv1 > 2)
goto L_10;
nn = *n2 - *n1 + 1;
Iv[RESTOR] = 0;
i = iv1 + 4;
if (Iv[TOOBIG] == 0)
switch (i)
{
case 1: goto L_150;
case 2: goto L_130;
case 3: goto L_150;
case 4: goto L_120;
case 5: goto L_120;
case 6: goto L_150;
}
if (i != 5)
Iv[1] = 2;
goto L_40;
/* *** FRESH START OR RESTART -- CHECK INPUT INTEGERS ***
* */
L_10:
if (nd <= 0)
goto L_210;
if (p <= 0)
goto L_210;
if (n <= 0)
goto L_210;
if (iv1 == 14)
goto L_30;
if (iv1 > 16)
goto L_300;
if (iv1 < 12)
goto L_40;
if (iv1 == 12)
Iv[1] = 13;
if (Iv[1] != 13)
goto L_20;
Iv[IVNEED] += p;
Iv[VNEED] += p*(p + 13)/2;
L_20:
sg7lit( d, x, iv, liv, lv, p, p, v, x, x );
if (Iv[1] != 14)
goto L_999;
/* *** STORAGE ALLOCATION ***
* */
Iv[IPIVOT] = Iv[NEXTIV];
Iv[NEXTIV] = Iv[IPIVOT] + p;
Iv[Y] = Iv[NEXTV];
Iv[G] = Iv[Y] + p;
Iv[JCN] = Iv[G] + p;
Iv[RMAT] = Iv[JCN] + p;
Iv[QTR] = Iv[RMAT] + lh;
Iv[JTOL] = Iv[QTR] + p;
Iv[NEXTV] = Iv[JTOL] + 2*p;
if (iv1 == 13)
goto L_999;
L_30:
jtol1 = Iv[JTOL];
if (V[DINIT] >= ZERO)
sv7scp( p, d, V[DINIT] );
if (V[DTINIT] > ZERO)
sv7scp( p, &V[jtol1], V[DTINIT] );
i = jtol1 + p;
if (V[D0INIT] > ZERO)
sv7scp( p, &V[i], V[D0INIT] );
Iv[NF0] = 0;
Iv[NF1] = 0;
if (nd >= n)
goto L_40;
/* *** SPECIAL CASE HANDLING OF FIRST FUNCTION AND GRADIENT EVALUATION
* *** -- ASK FOR BOTH RESIDUAL AND JACOBIAN AT ONCE
* */
g1 = Iv[G];
y1 = Iv[Y];
sg7lit( d, &V[g1], iv, liv, lv, p, p, v, x, &V[y1] );
if (Iv[1] != 1)
goto L_220;
V[F] = ZERO;
sv7scp( p, &V[g1], ZERO );
Iv[1] = -1;
qtr1 = Iv[QTR];
sv7scp( p, &V[qtr1], ZERO );
Iv[REGD] = 0;
rmat1 = Iv[RMAT];
goto L_100;
L_40:
g1 = Iv[G];
y1 = Iv[Y];
sg7lit( d, &V[g1], iv, liv, lv, p, p, v, x, &V[y1] );
switch (IARITHIF(Iv[1] - 2))
{
case -1: goto L_50;
case 0: goto L_60;
case 1: goto L_220;
}
L_50:
V[F] = ZERO;
if (Iv[NF1] == 0)
goto L_260;
if (Iv[RESTOR] != 2)
goto L_260;
Iv[NF0] = Iv[NF1];
sv7cpy( n, rd, r );
Iv[REGD] = 0;
goto L_260;
L_60:
sv7scp( p, &V[g1], ZERO );
if (Iv[MODE] > 0)
goto L_230;
rmat1 = Iv[RMAT];
qtr1 = Iv[QTR];
sv7scp( p, &V[qtr1], ZERO );
Iv[REGD] = 0;
if (nd < n)
goto L_90;
if (*n1 != 1)
goto L_90;
if (Iv[MODE] < 0)
goto L_100;
if (Iv[NF1] == Iv[NFGCAL])
goto L_70;
if (Iv[NF0] != Iv[NFGCAL])
goto L_90;
sv7cpy( n, r, rd );
goto L_80;
L_70:
sv7cpy( n, rd, r );
L_80:
sq7apl( nd, n, p, dr, rd, 0 );
sl7vml( p, &V[y1], &V[rmat1], rd );
goto L_110;
L_90:
Iv[1] = -2;
if (Iv[MODE] < 0)
Iv[1] = -1;
L_100:
sv7scp( p, &V[y1], ZERO );
L_110:
sv7scp( lh, &V[rmat1], ZERO );
goto L_260;
/* *** COMPUTE F(X) ***
* */
L_120:
t = sv2nrm( nn, r );
if (t > V[RLIMIT])
goto L_200;
V[F] += HALF*SQ(t);
if (*n2 < n)
goto L_270;
if (*n1 == 1)
Iv[NF1] = Iv[NFCALL];
goto L_40;
/* *** COMPUTE Y ***
* */
L_130:
y1 = Iv[Y];
yi = y1;
for (l = 1; l <= p; l++)
{
V[yi] += sd7tpr( nn, &DR(l - 1,0), r );
yi += 1;
}
if (*n2 < n)
goto L_270;
Iv[1] = 2;
if (*n1 > 1)
Iv[1] = -3;
goto L_260;
/* *** COMPUTE GRADIENT INFORMATION ***
* */
L_150:
if (Iv[MODE] > p)
goto L_240;
g1 = Iv[G];
ivmode = Iv[MODE];
if (ivmode < 0)
goto L_170;
if (ivmode == 0)
goto L_180;
Iv[1] = 2;
/* *** COMPUTE GRADIENT ONLY (FOR USE IN COVARIANCE COMPUTATION) ***
* */
gi = g1;
for (l = 1; l <= p; l++)
{
V[gi] += sd7tpr( nn, r, &DR(l - 1,0) );
gi += 1;
}
goto L_190;
/* *** COMPUTE INITIAL FUNCTION VALUE WHEN ND .LT. N ***
* */
L_170:
if (n <= nd)
goto L_180;
t = sv2nrm( nn, r );
if (t > V[RLIMIT])
goto L_200;
V[F] += HALF*SQ(t);
/* *** UPDATE D IF DESIRED ***
* */
L_180:
if (Iv[DTYPE] > 0)
sd7upd( d, dr, iv, liv, lv, n, nd, nn, *n2, p, v );
/* *** COMPUTE RMAT AND QTR ***
* */
qtr1 = Iv[QTR];
rmat1 = Iv[RMAT];
sq7rad( nn, nd, p, &V[qtr1], TRUE, &V[rmat1], dr, r );
Iv[NF1] = 0;
L_190:
if (*n2 < n)
goto L_270;
if (ivmode > 0)
goto L_40;
Iv[NF00] = Iv[NFGCAL];
/* *** COMPUTE G FROM RMAT AND QTR ***
* */
sl7vml( p, &V[g1], &V[rmat1], &V[qtr1] );
Iv[1] = 2;
if (ivmode == 0)
goto L_40;
if (n <= nd)
goto L_40;
/* *** FINISH SPECIAL CASE HANDLING OF FIRST FUNCTION AND GRADIENT
* */
y1 = Iv[Y];
Iv[1] = 1;
sg7lit( d, &V[g1], iv, liv, lv, p, p, v, x, &V[y1] );
if (Iv[1] != 2)
goto L_220;
goto L_40;
/* *** MISC. DETAILS ***
*
* *** X IS OUT OF RANGE (OVERSIZE STEP) ***
* */
L_200:
Iv[TOOBIG] = 1;
goto L_40;
/* *** BAD N, ND, OR P ***
* */
L_210:
Iv[1] = 66;
goto L_300;
/* *** CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE ***
* */
L_220:
if (Iv[COVMAT] != 0)
goto L_290;
if (Iv[REGD] != 0)
goto L_290;
/* *** SEE IF CHOLESKY FACTOR OF HESSIAN IS AVAILABLE ***
* */
k = Iv[FDH];
if (k <= 0)
goto L_280;
if (Iv[RDREQ] <= 0)
goto L_290;
/* *** COMPUTE REGRESSION DIAGNOSTICS AND DEFAULT COVARIANCE IF
* DESIRED ***
* */
i = 0;
if ((Iv[RDREQ]%4) >= 2)
i = 1;
if ((Iv[RDREQ]%2) == 1 && labs( Iv[COVREQ] ) <= 1)
i += 2;
if (i == 0)
goto L_250;
Iv[MODE] = p + i;
Iv[NGCALL] += 1;
Iv[NGCOV] += 1;
Iv[CNVCOD] = Iv[1];
if (i < 2)
goto L_230;
l = labs( Iv[H] );
sv7scp( lh, &V[l], ZERO );
L_230:
Iv[NFCOV] += 1;
Iv[NFCALL] += 1;
Iv[NFGCAL] = Iv[NFCALL];
Iv[1] = -1;
goto L_260;
L_240:
l = Iv[LMAT];
sn2lrd( dr, iv, &V[l], lh, liv, lv, nd, nn, p, r, rd, v );
if (*n2 < n)
goto L_270;
if (*n1 > 1)
goto L_250;
/* *** ENSURE WE CAN RESTART -- AND MAKE RETURN STATE OF DR
* *** INDEPENDENT OF WHETHER REGRESSION DIAGNOSTICS ARE COMPUTED.
* *** USE STEP VECTOR (ALLOCATED BY SG7LIT) FOR SCRATCH.
* */
rmat1 = Iv[RMAT];
sv7scp( lh, &V[rmat1], ZERO );
sq7rad( nn, nd, p, r, FALSE, &V[rmat1], dr, r );
Iv[NF1] = 0;
/* *** FINISH COMPUTING COVARIANCE ***
* */
L_250:
l = Iv[LMAT];
sc7vfn( iv, &V[l], lh, liv, lv, n, p, v );
goto L_290;
/* *** RETURN FOR MORE FUNCTION OR GRADIENT INFORMATION ***
* */
L_260:
*n2 = 0;
L_270:
*n1 = *n2 + 1;
*n2 += nd;
if (*n2 > n)
*n2 = n;
goto L_999;
/* *** COME HERE FOR INDEFINITE FINITE-DIFFERENCE HESSIAN ***
* */
L_280:
Iv[COVMAT] = k;
Iv[REGD] = k;
/* *** PRINT SUMMARY OF FINAL ITERATION AND OTHER REQUESTED ITEMS ***
* */
L_290:
g1 = Iv[G];
L_300:
sitsum( d, &V[g1], iv, liv, lv, p, v, x );
if (Iv[1] <= 6 && Iv[RDREQ] > 0)
sn2cvp( iv, liv, lv, p, v );
L_999:
return;
/* *** LAST LINE OF SRN2G FOLLOWS *** */
#undef DR
} /* end of function */
/* ================================================================== */
/* PARAMETER translations */
#define COSMIN 47
#define DGNORM 1
#define DIG 37
#define DSTNRM 2
#define F0 13
#define FDIF 11
#define FUZZ 45
#define GTSTEP 4
#define HC 71
#define IERR 75
#define INCFAC 23
#define INITS 25
#define IRC 29
#define KAGQT 33
#define KALM 34
#define LMAX0 35
#define LMAXS 36
#define MODEL 5
#define MXFCAL 17
#define MXITER 18
#define NEGONE (-1.e0)
#define NITER 31
#define NVSAVE 9
#define ONE 1.e0
#define ONEP2 1.2e0
#define PHMXFC 21
#define PREDUC 7
#define RAD0 9
#define RADFAC 16
#define RADINC 8
#define RADIUS 8
#define RCOND 53
#define RELDX 17
#define S 62
#define SIZE 55
#define STEP 40
#define STGLIM 11
#define STLSTG 41
#define STPPAR 5
#define SUSED 64
#define SWITCH_ 12
#define TUNER4 29
#define TUNER5 30
#define VSAVE 60
#define W 65
#define WSCALE 56
#define X0 43
#define XIRC 13
/* end of PARAMETER translations */
void /*FUNCTION*/ sg7lit(
float d[],
float gg[],
long iv[],
long liv,
long lv,
long p,
long ps,
float v[],
float x[],
float yy[])
{
long int dig1, g01, h1, hc1, i, ipiv1, j, k, l, lmat1, lstgst,
pp1o2, qtr1, rmat1, rstrst, s1, step1, stpmod, temp1, temp2,
w1, x01;
float e, sttsst, t, t1, tp;
/* OFFSET Vectors w/subscript range: 1 to dimension */
float *const D = &d[0] - 1;
float *const Gg = &gg[0] - 1;
long *const Iv = &iv[0] - 1;
float *const V = &v[0] - 1;
float *const X = &x[0] - 1;
float *const Yy = &yy[0] - 1;
/* end of OFFSET VECTORS */
/*>> 1990-06-12 CLL @ JPL
*>> 1990-04-23 CLL (Recent revision by DMG)
*** from netlib, Mon Apr 23 20:37:24 EDT 1990 ***
*>> 1990-02-20 CLL @ JPL
*
* *** CARRY OUT NL2SOL-LIKE ITERATIONS FOR GENERALIZED LINEAR ***
* *** REGRESSION PROBLEMS (AND OTHERS OF SIMILAR STRUCTURE) ***
*
* *** PARAMETER DECLARATIONS ***
* */
/* ------------------------- PARAMETER USAGE --------------------------
*
* D.... SCALE VECTOR.
* IV... INTEGER VALUE ARRAY.
* LIV.. LENGTH OF IV. MUST BE AT LEAST 82.
* LH... LENGTH OF H = P*(P+1)/2.
* LV... LENGTH OF V. MUST BE AT LEAST P*(3*P + 19)/2 + 7.
* GG... GRADIENT AT X (WHEN IV(1) = 2).
* P.... NUMBER OF PARAMETERS (COMPONENTS IN X).
* PS... NUMBER OF NONZERO ROWS AND COLUMNS IN S.
* V.... FLOATING-POINT VALUE ARRAY.
* X.... PARAMETER VECTOR.
* YY... PART OF YIELD VECTOR (WHEN IV(1)= 2, SCRATCH OTHERWISE).
*
* *** DISCUSSION ***
*
* SG7LIT PERFORMS NL2SOL-LIKE ITERATIONS FOR A VARIETY OF
* REGRESSION PROBLEMS THAT ARE SIMILAR TO NONLINEAR LEAST-SQUARES
* IN THAT THE HESSIAN IS THE SUM OF TWO TERMS, A READILY-COMPUTED
* FIRST-ORDER TERM AND A SECOND-ORDER TERM. THE CALLER SUPPLIES
* THE FIRST-ORDER TERM OF THE HESSIAN IN HC (LOWER TRIANGLE, STORED
* COMPACTLY BY ROWS IN V, STARTING AT IV(HC)), AND SG7LIT BUILDS AN
* APPROXIMATION, S, TO THE SECOND-ORDER TERM. THE CALLER ALSO
* PROVIDES THE FUNCTION VALUE, GRADIENT, AND PART OF THE YIELD
* VECTOR USED IN UPDATING S. SG7LIT DECIDES DYNAMICALLY WHETHER OR
* NOT TO USE S WHEN CHOOSING THE NEXT STEP TO TRY... THE HESSIAN
* APPROXIMATION USED IS EITHER HC ALONE (GAUSS-NEWTON MODEL) OR
* HC + S (AUGMENTED MODEL).
*
* IF PS .LT. P, THEN ROWS AND COLUMNS PS+1...P OF S ARE KEPT
* CONSTANT. THEY WILL BE ZERO UNLESS THE CALLER SETS IV(INITS) TO
* 1 OR 2 AND SUPPLIES NONZERO VALUES FOR THEM, OR THE CALLER SETS
* IV(INITS) TO 3 OR 4 AND THE FINITE-DIFFERENCE INITIAL S THEN
* COMPUTED HAS NONZERO VALUES IN THESE ROWS.
*
* IF IV(INITS) IS 3 OR 4, THEN THE INITIAL S IS COMPUTED BY
* FINITE DIFFERENCES. 3 MEANS USE FUNCTION DIFFERENCES, 4 MEANS
* USE GRADIENT DIFFERENCES. FINITE DIFFERENCING IS DONE THE SAME
* WAY AS IN COMPUTING A COVARIANCE MATRIX (WITH IV(COVREQ) = -1, -2,
* 1, OR 2).
*
* FOR UPDATING S,SG7LIT ASSUMES THAT THE GRADIENT HAS THE FORM
* OF A SUM OVER I OF RHO(I,X)*GRAD(R(I,X)), WHERE GRAD DENOTES THE
* GRADIENT WITH RESPECT TO X. THE TRUE SECOND-ORDER TERM THEN IS
* THE SUM OVER I OF RHO(I,X)*HESSIAN(R(I,X)). IF X = X0 + STEP,
* THEN WE WISH TO UPDATE S SO THAT S*STEP IS THE SUM OVER I OF
* RHO(I,X)*(GRAD(R(I,X)) - GRAD(R(I,X0))). THE CALLER MUST SUPPLY
* PART OF THIS IN YY, NAMELY THE SUM OVER I OF
* RHO(I,X)*GRAD(R(I,X0)), WHEN CALLING SG7LIT WITH IV(1) = 2 AND
* IV(MODE) = 0 (WHERE MODE = 38). GG THEN CONTANS THE OTHER PART,
* SO THAT THE DESIRED YIELD VECTOR IS GG - YY. IF PS .LT. P, THEN
* THE ABOVE DISCUSSION APPLIES ONLY TO THE FIRST PS COMPONENTS OF
* GRAD(R(I,X)), STEP, AND YY.
*
* PARAMETERS IV, P, V, AND X ARE THE SAME AS THE CORRESPONDING
* ONES TO NL2SOL (WHICH SEE), EXCEPT THAT V CAN BE SHORTER
* (SINCE THE PART OF V THAT NL2SOL USES FOR STORING D, J, AND R IS
* NOT NEEDED). MOREOVER, COMPARED WITH NL2SOL, IV(1) MAY HAVE THE
* TWO ADDITIONAL OUTPUT VALUES 1 AND 2, WHICH ARE EXPLAINED BELOW,
* AS IS THE USE OF IV(TOOBIG) AND IV(NFGCAL). THE VALUES IV(D),
* IV(J), AND IV(R), WHICH ARE OUTPUT VALUES FROM NL2SOL (AND
* NL2SNO), ARE NOT REFERENCED BY SG7LIT OR THE SUBROUTINES IT CALLS.
*
* WHEN SG7LIT IS FIRST CALLED, I.E., WHEN SG7LIT IS CALLED WITH
* IV(1) = 0 OR 12, V(F), GG, AND HC NEED NOT BE INITIALIZED. TO
* OBTAIN THESE STARTING VALUES,SG7LIT RETURNS FIRST WITH IV(1) = 1,
* THEN WITH IV(1) = 2, WITH IV(MODE) = -1 IN BOTH CASES. ON
* SUBSEQUENT RETURNS WITH IV(1) = 2, IV(MODE) = 0 IMPLIES THAT
* YY MUST ALSO BE SUPPLIED. (NOTE THAT YY IS USED FOR SCRATCH --
* ITS INPUT CONTENTS ARE LOST. BY CONTRAST, HC IS NEVER CHANGED.)
* ONCE CONVERGENCE HAS BEEN OBTAINED, IV(RDREQ) AND IV(COVREQ) MAY
* IMPLY THAT A FINITE-DIFFERENCE HESSIAN SHOULD BE COMPUTED FOR USE
* IN COMPUTING A COVARIANCE MATRIX. IN THIS CASE SG7LIT WILL MAKE A
* NUMBER OF RETURNS WITH IV(1) = 1 OR 2 AND IV(MODE) POSITIVE.
* WHEN IV(MODE) IS POSITIVE, YY SHOULD NOT BE CHANGED.
*
* IV(1) = 1 MEANS THE CALLER SHOULD SET V(F) (I.E., V(10)) TO F(X), THE
* FUNCTION VALUE AT X, AND CALL SG7LIT AGAIN, HAVING CHANGED
* NONE OF THE OTHER PARAMETERS. AN EXCEPTION OCCURS IF F(X)
* CANNOT BE EVALUATED (E.G. IF OVERFLOW WOULD OCCUR), WHICH
* MAY HAPPEN BECAUSE OF AN OVERSIZED STEP. IN THIS CASE
* THE CALLER SHOULD SET IV(TOOBIG) = IV(2) TO 1, WHICH WILL
* CAUSE SG7LIT TO IGNORE V(F) AND TRY A SMALLER STEP. NOTE
* THAT THE CURRENT FUNCTION EVALUATION COUNT IS AVAILABLE
* IN IV(NFCALL) = IV(6). THIS MAY BE USED TO IDENTIFY
* WHICH COPY OF SAVED INFORMATION SHOULD BE USED IN COM-
* PUTING GG, HC, AND YY THE NEXT TIME SG7LIT RETURNS WITH
* IV(1) = 2. SEE MLPIT FOR AN EXAMPLE OF THIS.
* IV(1) = 2 MEANS THE CALLER SHOULD SET GG TO G(X), THE GRADIENT OF F AT
* X. THE CALLER SHOULD ALSO SET HC TO THE GAUSS-NEWTON
* HESSIAN AT X. IF IV(MODE) = 0, THEN THE CALLER SHOULD
* ALSO COMPUTE THE PART OF THE YIELD VECTOR DESCRIBED ABOVE.
* THE CALLER SHOULD THEN CALL SG7LIT AGAIN (WITH IV(1) = 2).
* THE CALLER MAY ALSO CHANGE D AT THIS TIME, BUT SHOULD NOT
* CHANGE X. NOTE THAT IV(NFGCAL) = IV(7) CONTAINS THE
* VALUE THAT IV(NFCALL) HAD DURING THE RETURN WITH
* IV(1) = 1 IN WHICH X HAD THE SAME VALUE AS IT NOW HAS.
* IV(NFGCAL) IS EITHER IV(NFCALL) OR IV(NFCALL) - 1. MLPIT
* IS AN EXAMPLE WHERE THIS INFORMATION IS USED. IF GG OR HC */
/* CANNOT BE EVALUATED AT X, THEN THE CALLER MAY SET
* IV(TOOBIG) TO 1, IN WHICH CASE SG7LIT WILL RETURN WITH
* IV(1) = 15.
*
* *** GENERAL ***
*
* CODED BY DAVID M. GAY.
* THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
* SUPPORTED IN PART BY D.O.E. GRANT EX-76-A-01-2295 TO MIT/CCREMS.
*
* (SEE NL2SOL FOR REFERENCES.)
*
* ------------------------------------------------------------------
* References to the function STOPX have been commented out of this
* subroutine. If one wishes to be able to terminate this package
* gracefully using a keybord "Break" key, one can provide a STOPX
* function that returns .true. if the Break key has been pressed
* since the last call to STOPX, and otherwise returns .false., and
* then uncomment the references to STOPX in this subr.
* -- CLL 6/12/90
* ++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++++
*
* *** LOCAL VARIABLES ***
*
* integer DUMMY */
/* *** CONSTANTS ***
* */
/* *** EXTERNAL FUNCTIONS AND SUBROUTINES ***
*
* external STOPX
* LOGICAL STOPX */
/* ------------------------------------------------------------------
* SA7SST.... ASSESSES CANDIDATE STEP.
* SD7TPR... RETURNS INNER PRODUCT OF TWO VECTORS.
* SF7HES.... COMPUTE FINITE-DIFFERENCE HESSIAN (FOR COVARIANCE).
* SG7QTS.... COMPUTES GOLDFELD-QUANDT-TROTTER STEP (AUGMENTED MODEL).
* SITSUM.... PRINTS ITERATION SUMMARY AND INFO ON INITIAL AND FINAL X.
* SL7MST... COMPUTES LEVENBERG-MARQUARDT STEP (GAUSS-NEWTON MODEL).
* SL7SRT.... COMPUTES CHOLESKY FACTOR OF (LOWER TRIANG. OF) SYM. MATRIX.
* SL7SQR... COMPUTES L * L**T FROM LOWER TRIANGULAR MATRIX L.
* SL7TVM... COMPUTES L**T * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX.
* SL7SVX... ESTIMATES LARGEST SING. VALUE OF LOWER TRIANG. MATRIX.
* SL7SVN... ESTIMATES SMALLEST SING. VALUE OF LOWER TRIANG. MATRIX.
* SL7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX.
* SPARCK.... CHECK VALIDITY OF IV AND V INPUT COMPONENTS.
* SRLDST... COMPUTES V(RELDX) = RELATIVE STEP SIZE.
* SR7MDC... RETURNS MACHINE-DEPENDENT CONSTANTS.
* SS7LUP... PERFORMS QUASI-NEWTON UPDATE ON COMPACTLY STORED LOWER TRI-
* ANGLE OF A SYMMETRIC MATRIX.
* STOPX.... RETURNS .TRUE. IF THE BREAK KEY HAS BEEN PRESSED.
* Call to STOPX commented out. -- CLL 6/12/90
* SV2AXY.... COMPUTES SCALAR TIMES ONE VECTOR PLUS ANOTHER.
* SV7CPY.... COPIES ONE VECTOR TO ANOTHER.
* SV7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR.
* SV2NRM... RETURNS THE 2-NORM OF A VECTOR.
*
* *** SUBSCRIPTS FOR IV AND V ***
* */
/* *** IV SUBSCRIPT VALUES ***
* */
/* *** V SUBSCRIPT VALUES ***
* */
/* ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++
* */
i = Iv[1];
if (i == 1)
goto L_40;
if (i == 2)
goto L_50;
if (i == 12 || i == 13)
Iv[VNEED] += p*(3*p + 19)/2 + 7;
sparck( 1, d, iv, liv, lv, p, v );
i = Iv[1] - 2;
if (i > 12)
goto L_999;
switch (i)
{
case 1: goto L_290;
case 2: goto L_290;
case 3: goto L_290;
case 4: goto L_290;
case 5: goto L_290;
case 6: goto L_290;
case 7: goto L_170;
case 8: goto L_120;
case 9: goto L_170;
case 10: goto L_10;
case 11: goto L_10;
case 12: goto L_20;
}
/* *** STORAGE ALLOCATION ***
* */
L_10:
pp1o2 = p*(p + 1)/2;
Iv[S] = Iv[LMAT] + pp1o2;
Iv[X0] = Iv[S] + pp1o2;
Iv[STEP] = Iv[X0] + p;
Iv[STLSTG] = Iv[STEP] + p;
Iv[DIG] = Iv[STLSTG] + p;
Iv[W] = Iv[DIG] + p;
Iv[H] = Iv[W] + 4*p + 7;
Iv[NEXTV] = Iv[H] + pp1o2;
if (Iv[1] != 13)
goto L_20;
Iv[1] = 14;
goto L_999;
/* *** INITIALIZATION ***
* */
L_20:
Iv[NITER] = 0;
Iv[NFCALL] = 1;
Iv[NGCALL] = 1;
Iv[NFGCAL] = 1;
Iv[MODE] = -1;
Iv[STGLIM] = 2;
Iv[TOOBIG] = 0;
Iv[CNVCOD] = 0;
Iv[COVMAT] = 0;
Iv[NFCOV] = 0;
Iv[NGCOV] = 0;
Iv[RADINC] = 0;
Iv[RESTOR] = 0;
Iv[FDH] = 0;
V[RAD0] = ZERO;
V[STPPAR] = ZERO;
V[RADIUS] = V[LMAX0]/(ONE + V[PHMXFC]);
/* *** SET INITIAL MODEL AND S MATRIX ***
* */
Iv[MODEL] = 1;
if (Iv[S] < 0)
goto L_999;
if (Iv[INITS] > 1)
Iv[MODEL] = 2;
s1 = Iv[S];
if (Iv[INITS] == 0 || Iv[INITS] > 2)
sv7scp( p*(p + 1)/2, &V[s1], ZERO );
Iv[1] = 1;
j = Iv[IPIVOT];
if (j <= 0)
goto L_999;
for (i = 1; i <= p; i++)
{
Iv[j] = i;
j += 1;
}
goto L_999;
/* *** NEW FUNCTION VALUE ***
* */
L_40:
if (Iv[MODE] == 0)
goto L_290;
if (Iv[MODE] > 0)
goto L_520;
Iv[1] = 2;
if (Iv[TOOBIG] == 0)
goto L_999;
Iv[1] = 63;
goto L_999;
/* *** NEW GRADIENT ***
* */
L_50:
Iv[KALM] = -1;
Iv[KAGQT] = -1;
Iv[FDH] = 0;
if (Iv[MODE] > 0)
goto L_520;
/* *** MAKE SURE GRADIENT COULD BE COMPUTED ***
* */
if (Iv[TOOBIG] == 0)
goto L_60;
Iv[1] = 65;
goto L_999;
L_60:
if (Iv[HC] <= 0 && Iv[RMAT] <= 0)
goto L_610;
/* *** COMPUTE D**-1 * GRADIENT ***
* */
dig1 = Iv[DIG];
k = dig1;
for (i = 1; i <= p; i++)
{
V[k] = Gg[i]/D[i];
k += 1;
}
V[DGNORM] = sv2nrm( p, &V[dig1] );
if (Iv[CNVCOD] != 0)
goto L_510;
if (Iv[MODE] == 0)
goto L_440;
Iv[MODE] = 0;
V[F0] = V[F];
if (Iv[INITS] <= 2)
goto L_100;
/* *** ARRANGE FOR FINITE-DIFFERENCE INITIAL S ***
* */
Iv[XIRC] = Iv[COVREQ];
Iv[COVREQ] = -1;
if (Iv[INITS] > 3)
Iv[COVREQ] = 1;
Iv[CNVCOD] = 70;
goto L_530;
/* *** COME TO NEXT STMT AFTER COMPUTING F.D. HESSIAN FOR INIT. S ***
* */
L_80:
Iv[CNVCOD] = 0;
Iv[MODE] = 0;
Iv[NFCOV] = 0;
Iv[NGCOV] = 0;
Iv[COVREQ] = Iv[XIRC];
s1 = Iv[S];
pp1o2 = ps*(ps + 1)/2;
hc1 = Iv[HC];
if (hc1 <= 0)
goto L_90;
sv2axy( pp1o2, &V[s1], NEGONE, &V[hc1], &V[h1] );
goto L_100;
L_90:
rmat1 = Iv[RMAT];
sl7sqr( ps, &V[s1], &V[rmat1] );
sv2axy( pp1o2, &V[s1], NEGONE, &V[s1], &V[h1] );
L_100:
Iv[1] = 2;
/* ---------------------------- MAIN LOOP -----------------------------
*
*
* *** PRINT ITERATION SUMMARY, CHECK ITERATION LIMIT ***
* */
L_110:
sitsum( d, gg, iv, liv, lv, p, v, x );
L_120:
k = Iv[NITER];
if (k < Iv[MXITER])
goto L_130;
Iv[1] = 10;
goto L_999;
L_130:
Iv[NITER] = k + 1;
/* *** UPDATE RADIUS ***
* */
if (k == 0)
goto L_150;
step1 = Iv[STEP];
for (i = 1; i <= p; i++)
{
V[step1] *= D[i];
step1 += 1;
}
step1 = Iv[STEP];
t = V[RADFAC]*sv2nrm( p, &V[step1] );
if (V[RADFAC] < ONE || t > V[RADIUS])
V[RADIUS] = t;
/* *** INITIALIZE FOR START OF NEXT ITERATION ***
* */
L_150:
x01 = Iv[X0];
V[F0] = V[F];
Iv[IRC] = 4;
Iv[H] = -labs( Iv[H] );
Iv[SUSED] = Iv[MODEL];
/* *** COPY X TO X0 ***
* */
sv7cpy( p, &V[x01], x );
/* *** CHECK STOPX AND FUNCTION EVALUATION LIMIT ***
* */
L_160:
;
/* if (STOPX(DUMMY)) then
* IV(1) = 11
* GO TO 190
* else */
goto L_180;
/* endif
*
* *** COME HERE WHEN RESTARTING AFTER FUNC. EVAL. LIMIT OR STOPX.
* */
L_170:
if (V[F] >= V[F0])
goto L_180;
V[RADFAC] = ONE;
k = Iv[NITER];
goto L_130;
L_180:
if (Iv[NFCALL] < Iv[MXFCAL] + Iv[NFCOV])
goto L_200;
Iv[1] = 9;
/* 190 continue */
if (V[F] >= V[F0])
goto L_999;
/* *** IN CASE OF STOPX OR FUNCTION EVALUATION LIMIT WITH
* *** IMPROVED V(F), EVALUATE THE GRADIENT AT X.
* */
Iv[CNVCOD] = Iv[1];
goto L_430;
/*. . . . . . . . . . . . . COMPUTE CANDIDATE STEP . . . . . . . . . .
* */
L_200:
step1 = Iv[STEP];
w1 = Iv[W];
h1 = Iv[H];
t1 = ONE;
if (Iv[MODEL] == 2)
goto L_210;
t1 = ZERO;
/* *** COMPUTE LEVENBERG-MARQUARDT STEP IF POSSIBLE...
* */
rmat1 = Iv[RMAT];
if (rmat1 <= 0)
goto L_210;
qtr1 = Iv[QTR];
if (qtr1 <= 0)
goto L_210;
ipiv1 = Iv[IPIVOT];
sl7mst( d, gg, Iv[IERR], &Iv[ipiv1], &Iv[KALM], p, &V[qtr1], &V[rmat1],
&V[step1], v, &V[w1] );
/* *** H IS STORED IN THE END OF W AND HAS JUST BEEN OVERWRITTEN,
* *** SO WE MARK IT INVALID... */
Iv[H] = -labs( h1 );
/* *** EVEN IF H WERE STORED ELSEWHERE, IT WOULD BE NECESSARY TO
* *** MARK INVALID THE INFORMATION SG7QTS MAY HAVE STORED IN V... */
Iv[KAGQT] = -1;
goto L_260;
L_210:
if (h1 > 0)
goto L_250;
/* *** SET H TO D**-1 * (HC + T1*S) * D**-1. ***
* */
h1 = -h1;
Iv[H] = h1;
Iv[FDH] = 0;
j = Iv[HC];
if (j > 0)
goto L_220;
j = h1;
rmat1 = Iv[RMAT];
sl7sqr( p, &V[h1], &V[rmat1] );
L_220:
s1 = Iv[S];
for (i = 1; i <= p; i++)
{
t = ONE/D[i];
for (k = 1; k <= i; k++)
{
V[h1] = t*(V[j] + t1*V[s1])/D[k];
j += 1;
h1 += 1;
s1 += 1;
}
}
h1 = Iv[H];
Iv[KAGQT] = -1;
/* *** COMPUTE ACTUAL GOLDFELD-QUANDT-TROTTER STEP ***
* */
L_250:
dig1 = Iv[DIG];
lmat1 = Iv[LMAT];
sg7qts( d, &V[dig1], &V[h1], &Iv[KAGQT], &V[lmat1], p, &V[step1],
v, &V[w1] );
if (Iv[KALM] > 0)
Iv[KALM] = 0;
L_260:
if (Iv[IRC] != 6)
goto L_270;
if (Iv[RESTOR] != 2)
goto L_290;
rstrst = 2;
goto L_300;
/* *** CHECK WHETHER EVALUATING F(X0 + STEP) LOOKS WORTHWHILE ***
* */
L_270:
Iv[TOOBIG] = 0;
if (V[DSTNRM] <= ZERO)
goto L_290;
if (Iv[IRC] != 5)
goto L_280;
if (V[RADFAC] <= ONE)
goto L_280;
if (V[PREDUC] > ONEP2*V[FDIF])
goto L_280;
if (Iv[RESTOR] != 2)
goto L_290;
rstrst = 0;
goto L_300;
/* *** COMPUTE F(X0 + STEP) ***
* */
L_280:
x01 = Iv[X0];
step1 = Iv[STEP];
sv2axy( p, x, ONE, &V[step1], &V[x01] );
Iv[NFCALL] += 1;
Iv[1] = 1;
goto L_999;
/*. . . . . . . . . . . . . ASSESS CANDIDATE STEP . . . . . . . . . . .
* */
L_290:
rstrst = 3;
L_300:
x01 = Iv[X0];
V[RELDX] = srldst( p, d, x, &V[x01] );
sa7sst( iv, liv, lv, v );
step1 = Iv[STEP];
lstgst = Iv[STLSTG];
i = Iv[RESTOR] + 1;
switch (i)
{
case 1: goto L_340;
case 2: goto L_310;
case 3: goto L_320;
case 4: goto L_330;
}
L_310:
sv7cpy( p, x, &V[x01] );
goto L_340;
L_320:
sv7cpy( p, &V[lstgst], &V[step1] );
goto L_340;
L_330:
sv7cpy( p, &V[step1], &V[lstgst] );
sv2axy( p, x, ONE, &V[step1], &V[x01] );
V[RELDX] = srldst( p, d, x, &V[x01] );
Iv[RESTOR] = rstrst;
/* *** IF NECESSARY, SWITCH MODELS ***
* */
L_340:
if (Iv[SWITCH_] == 0)
goto L_350;
Iv[H] = -labs( Iv[H] );
Iv[SUSED] += 2;
l = Iv[VSAVE];
sv7cpy( NVSAVE, v, &V[l] );
L_350:
l = Iv[IRC] - 4;
stpmod = Iv[MODEL];
if (l > 0)
switch (l)
{
case 1: goto L_370;
case 2: goto L_380;
case 3: goto L_390;
case 4: goto L_390;
case 5: goto L_390;
case 6: goto L_390;
case 7: goto L_390;
case 8: goto L_390;
case 9: goto L_500;
case 10: goto L_440;
}
/* *** DECIDE WHETHER TO CHANGE MODELS ***
* */
e = V[PREDUC] - V[FDIF];
s1 = Iv[S];
ss7lvm( ps, yy, &V[s1], &V[step1] );
sttsst = HALF*sd7tpr( ps, &V[step1], yy );
if (Iv[MODEL] == 1)
sttsst = -sttsst;
if (fabsf( e + sttsst )*V[FUZZ] >= fabsf( e ))
goto L_360;
/* *** SWITCH MODELS ***
* */
Iv[MODEL] = 3 - Iv[MODEL];
if (-2 < l)
goto L_400;
Iv[H] = -labs( Iv[H] );
Iv[SUSED] += 2;
l = Iv[VSAVE];
sv7cpy( NVSAVE, &V[l], v );
goto L_160;
L_360:
if (-3 < l)
goto L_400;
/* *** RECOMPUTE STEP WITH NEW RADIUS ***
* */
L_370:
V[RADIUS] = V[RADFAC]*V[DSTNRM];
goto L_160;
/* *** COMPUTE STEP OF LENGTH V(LMAXS) FOR SINGULAR CONVERGENCE TEST
* */
L_380:
V[RADIUS] = V[LMAXS];
goto L_200;
/* *** CONVERGENCE OR FALSE CONVERGENCE ***
* */
L_390:
Iv[CNVCOD] = l;
if (V[F] >= V[F0])
goto L_510;
if (Iv[XIRC] == 14)
goto L_510;
Iv[XIRC] = 14;
/*. . . . . . . . . . . . PROCESS ACCEPTABLE STEP . . . . . . . . . . .
* */
L_400:
Iv[COVMAT] = 0;
Iv[REGD] = 0;
/* *** SEE WHETHER TO SET V(RADFAC) BY GRADIENT TESTS ***
* */
if (Iv[IRC] != 3)
goto L_430;
step1 = Iv[STEP];
temp1 = Iv[STLSTG];
temp2 = Iv[W];
/* *** SET TEMP1 = HESSIAN * STEP FOR USE IN GRADIENT TESTS ***
* */
hc1 = Iv[HC];
if (hc1 <= 0)
goto L_410;
ss7lvm( p, &V[temp1], &V[hc1], &V[step1] );
goto L_420;
L_410:
rmat1 = Iv[RMAT];
sl7tvm( p, &V[temp1], &V[rmat1], &V[step1] );
sl7vml( p, &V[temp1], &V[rmat1], &V[temp1] );
L_420:
if (stpmod == 1)
goto L_430;
s1 = Iv[S];
ss7lvm( ps, &V[temp2], &V[s1], &V[step1] );
sv2axy( ps, &V[temp1], ONE, &V[temp2], &V[temp1] );
/* *** SAVE OLD GRADIENT AND COMPUTE NEW ONE ***
* */
L_430:
Iv[NGCALL] += 1;
g01 = Iv[W];
sv7cpy( p, &V[g01], gg );
Iv[1] = 2;
Iv[TOOBIG] = 0;
goto L_999;
/* *** INITIALIZATIONS -- G0 = GG - G0, ETC. ***
* */
L_440:
g01 = Iv[W];
sv2axy( p, &V[g01], NEGONE, &V[g01], gg );
step1 = Iv[STEP];
temp1 = Iv[STLSTG];
temp2 = Iv[W];
if (Iv[IRC] != 3)
goto L_470;
/* *** SET V(RADFAC) BY GRADIENT TESTS ***
*
* *** SET TEMP1 = D**-1 * (HESSIAN * STEP + (G(X0) - G(X))) ***
* */
k = temp1;
l = g01;
for (i = 1; i <= p; i++)
{
V[k] = (V[k] - V[l])/D[i];
k += 1;
l += 1;
}
/* *** DO GRADIENT TESTS ***
* */
if (sv2nrm( p, &V[temp1] ) <= V[DGNORM]*V[TUNER4])
goto L_460;
if (sd7tpr( p, gg, &V[step1] ) >= V[GTSTEP]*V[TUNER5])
goto L_470;
L_460:
V[RADFAC] = V[INCFAC];
/* *** COMPUTE YY VECTOR NEEDED FOR UPDATING S ***
* */
L_470:
sv2axy( ps, yy, NEGONE, yy, gg );
/* *** DETERMINE SIZING FACTOR V(SIZE) ***
*
* *** SET TEMP1 = S * STEP *** */
s1 = Iv[S];
ss7lvm( ps, &V[temp1], &V[s1], &V[step1] );
t1 = fabsf( sd7tpr( ps, &V[step1], &V[temp1] ) );
t = fabsf( sd7tpr( ps, &V[step1], yy ) );
V[SIZE] = ONE;
if (t < t1)
V[SIZE] = t/t1;
/* *** SET G0 TO WCHMTD CHOICE OF FLETCHER AND AL-BAALI ***
* */
hc1 = Iv[HC];
if (hc1 <= 0)
goto L_480;
ss7lvm( ps, &V[g01], &V[hc1], &V[step1] );
goto L_490;
L_480:
rmat1 = Iv[RMAT];
sl7tvm( ps, &V[g01], &V[rmat1], &V[step1] );
sl7vml( ps, &V[g01], &V[rmat1], &V[g01] );
L_490:
sv2axy( ps, &V[g01], ONE, yy, &V[g01] );
/* *** UPDATE S ***
* */
ss7lup( &V[s1], V[COSMIN], ps, V[SIZE], &V[step1], &V[temp1],
&V[temp2], &V[g01], &V[WSCALE], yy );
Iv[1] = 2;
goto L_110;
/*. . . . . . . . . . . . . . MISC. DETAILS . . . . . . . . . . . . . .
*
* *** BAD PARAMETERS TO ASSESS ***
* */
L_500:
Iv[1] = 64;
goto L_999;
/* *** CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE ***
* */
L_510:
if (Iv[RDREQ] == 0)
goto L_600;
if (Iv[FDH] != 0)
goto L_600;
if (Iv[CNVCOD] >= 7)
goto L_600;
if (Iv[REGD] > 0)
goto L_600;
if (Iv[COVMAT] > 0)
goto L_600;
if (labs( Iv[COVREQ] ) >= 3)
goto L_560;
if (Iv[RESTOR] == 0)
Iv[RESTOR] = 2;
goto L_530;
/* *** COMPUTE FINITE-DIFFERENCE HESSIAN FOR COMPUTING COVARIANCE ***
* */
L_520:
Iv[RESTOR] = 0;
L_530:
sf7hes( d, gg, &i, iv, liv, lv, p, v, x );
switch (i)
{
case 1: goto L_540;
case 2: goto L_550;
case 3: goto L_580;
}
L_540:
Iv[NFCOV] += 1;
Iv[NFCALL] += 1;
Iv[1] = 1;
goto L_999;
L_550:
Iv[NGCOV] += 1;
Iv[NGCALL] += 1;
Iv[NFGCAL] = Iv[NFCALL] + Iv[NGCOV];
Iv[1] = 2;
goto L_999;
L_560:
h1 = labs( Iv[H] );
Iv[H] = -h1;
pp1o2 = p*(p + 1)/2;
rmat1 = Iv[RMAT];
if (rmat1 <= 0)
goto L_570;
lmat1 = Iv[LMAT];
sv7cpy( pp1o2, &V[lmat1], &V[rmat1] );
V[RCOND] = ZERO;
goto L_590;
L_570:
hc1 = Iv[HC];
Iv[FDH] = h1;
sv7cpy( p*(p + 1)/2, &V[h1], &V[hc1] );
/* *** COMPUTE CHOLESKY FACTOR OF FINITE-DIFFERENCE HESSIAN
* *** FOR USE IN CALLER*S COVARIANCE CALCULATION...
* */
L_580:
lmat1 = Iv[LMAT];
h1 = Iv[FDH];
if (h1 <= 0)
goto L_600;
if (Iv[CNVCOD] == 70)
goto L_80;
sl7srt( 1, p, &V[lmat1], &V[h1], &i );
Iv[FDH] = -1;
V[RCOND] = ZERO;
if (i != 0)
goto L_600;
L_590:
Iv[FDH] = -1;
step1 = Iv[STEP];
t = sl7svn( p, &V[lmat1], &V[step1], &V[step1] );
if (t <= ZERO)
goto L_600;
tp = sl7svx( p, &V[lmat1], &V[step1], &V[step1] );
if (tp != ZERO)
{
t /= tp;
if (t > sr7mdc( 4 ))
Iv[FDH] = h1;
V[RCOND] = t;
}
L_600:
Iv[MODE] = 0;
Iv[1] = Iv[CNVCOD];
Iv[CNVCOD] = 0;
goto L_999;
/* *** SPECIAL RETURN FOR MISSING HESSIAN INFORMATION -- BOTH
* *** IV(HC) .LE. 0 AND IV(RMAT) .LE. 0
* */
L_610:
Iv[1] = 1400;
L_999:
return;
/* *** LAST LINE OF SG7LIT FOLLOWS *** */
} /* end of function */
/* ================================================================== */
void /*FUNCTION*/ sn2lrd(
float *dr,
long iv[],
float l[],
long lh,
long liv,
long lv,
long nd,
long nn,
long p,
float r[],
float rd[],
float v[])
{
#define DR(I_,J_) (*(dr+(I_)*(nd)+(J_)))
long int cov, i, j, m, step1, _i, _r;
float a, s, t;
static float onev[1];
static int _aini = 1;
/* OFFSET Vectors w/subscript range: 1 to dimension */
long *const Iv = &iv[0] - 1;
float *const L = &l[0] - 1;
float *const Onev = &onev[0] - 1;
float *const R = &r[0] - 1;
float *const Rd = &rd[0] - 1;
float *const V = &v[0] - 1;
/* end of OFFSET VECTORS */
if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
Onev[1] = 1.e0;
_aini = 0;
}
/* *** COMPUTE REGRESSION DIAGNOSTIC AND DEFAULT COVARIANCE MATRIX FOR
* SRN2G ***
*
* *** PARAMETERS ***
* */
/* *** CODED BY DAVID M. GAY (WINTER 1982, FALL 1983) ***
*
* *** EXTERNAL FUNCTIONS AND SUBROUTINES ***
* */
/* *** LOCAL VARIABLES ***
* */
/* *** CONSTANTS ***
* */
/* *** IV SUBSCRIPTS ***
* */
/* +++++++++++++++++++++++++++++++ BODY +++++++++++++++++++++++++++++++
* */
step1 = Iv[STEP];
i = Iv[RDREQ];
if (i <= 0)
goto L_999;
if ((i%4) < 2)
goto L_30;
sv7scp( nn, rd, NEGONE );
for (i = 1; i <= nn; i++)
{
a = SQ(R[i]);
m = step1;
for (j = 1; j <= p; j++)
{
V[m] = DR(j - 1,i - 1);
m += 1;
}
sl7ivm( p, &V[step1], l, &V[step1] );
s = sd7tpr( p, &V[step1], &V[step1] );
t = ONE - s;
if (t <= ZERO)
goto L_20;
a = a*s/t;
Rd[i] = sqrtf( a );
L_20:
;
}
L_30:
if (Iv[MODE] - p < 2)
goto L_999;
/* *** COMPUTE DEFAULT COVARIANCE MATRIX ***
* */
cov = labs( Iv[H] );
for (i = 1; i <= nn; i++)
{
m = step1;
for (j = 1; j <= p; j++)
{
V[m] = DR(j - 1,i - 1);
m += 1;
}
sl7ivm( p, &V[step1], l, &V[step1] );
sl7itv( p, &V[step1], l, &V[step1] );
so7prd( 1, lh, p, &V[cov], onev, &V[step1], &V[step1] );
}
L_999:
return;
/* *** LAST CARD OF SN2LRD FOLLOWS *** */
#undef DR
} /* end of function */
void /*FUNCTION*/ sc7vfn(
long iv[],
float l[],
long lh,
long liv,
long lv,
long n,
long p,
float v[])
{
long int cov, i;
static float half = 0.5e0;
/* OFFSET Vectors w/subscript range: 1 to dimension */
long *const Iv = &iv[0] - 1;
float *const L = &l[0] - 1;
float *const V = &v[0] - 1;
/* end of OFFSET VECTORS */
/* *** FINISH COVARIANCE COMPUTATION FOR SRN2G, SRNSG ***
* */
/* *** LOCAL VARIABLES ***
* */
/* *** SUBSCRIPTS FOR IV AND V ***
* */
/* *** BODY ***
* */
Iv[1] = Iv[CNVCOD];
i = Iv[MODE] - p;
Iv[MODE] = 0;
Iv[CNVCOD] = 0;
if (Iv[FDH] <= 0)
goto L_999;
if (ipow(i - 2,2) == 1)
Iv[REGD] = 1;
if ((Iv[RDREQ]%2) != 1)
goto L_999;
/* *** FINISH COMPUTING COVARIANCE MATRIX = INVERSE OF F.D. HESSIAN.
* */
cov = labs( Iv[H] );
Iv[FDH] = 0;
if (Iv[COVMAT] != 0)
goto L_999;
if (i >= 2)
goto L_10;
sl7nvr( p, &V[cov], l );
sl7tsq( p, &V[cov], &V[cov] );
L_10:
sv7scl( lh, &V[cov], V[F]/(half*(float)( max( 1, n - p ) )), &V[cov] );
Iv[COVMAT] = cov;
L_999:
return;
/* *** LAST LINE OF SC7VFN FOLLOWS *** */
} /* end of function */
/* PARAMETER translations */
#define DELTA 52
#define DELTA0 44
#define DLTFDC 42
#define FX 53
#define NEGPT5 (-0.5e0)
#define SAVEI 63
#define TWO 2.e0
#define XMSAVE 51
/* end of PARAMETER translations */
void /*FUNCTION*/ sf7hes(
float d[],
float gg[],
long *irt,
long iv[],
long liv,
long lv,
long p,
float v[],
float x[])
{
long int gsave1, hes, hmi, hpi, hpm, i, k, kind, l, m, mm1, mm1o2,
pp1o2, stp0, stpi, stpm;
float del;
/* OFFSET Vectors w/subscript range: 1 to dimension */
float *const D = &d[0] - 1;
float *const Gg = &gg[0] - 1;
long *const Iv = &iv[0] - 1;
float *const V = &v[0] - 1;
float *const X = &x[0] - 1;
/* end of OFFSET VECTORS */
/* *** COMPUTE FINITE-DIFFERENCE HESSIAN, STORE IT IN V STARTING
* *** AT V(IV(FDH)) = V(-IV(H)).
*
* *** IF IV(COVREQ) .GE. 0 THEN SF7HES USES GRADIENT DIFFERENCES,
* *** OTHERWISE FUNCTION DIFFERENCES. STORAGE IN V IS AS IN SG7LIT.
*
* IRT VALUES...
* 1 = COMPUTE FUNCTION VALUE, I.E., V(F).
* 2 = COMPUTE G.
* 3 = DONE.
*
*
* *** PARAMETER DECLARATIONS ***
* */
/* *** LOCAL VARIABLES ***
* */
/* *** EXTERNAL SUBROUTINES ***
* */
/* SV7CPY.... COPY ONE VECTOR TO ANOTHER.
*
* *** SUBSCRIPTS FOR IV AND V ***
* */
/* ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++
* */
*irt = 4;
kind = Iv[COVREQ];
m = Iv[MODE];
if (m > 0)
goto L_10;
Iv[H] = -labs( Iv[H] );
Iv[FDH] = 0;
Iv[KAGQT] = -1;
V[FX] = V[F];
L_10:
if (m > p)
goto L_999;
if (kind < 0)
goto L_110;
/* *** COMPUTE FINITE-DIFFERENCE HESSIAN USING BOTH FUNCTION AND
* *** GRADIENT VALUES.
* */
gsave1 = Iv[W] + p;
if (m > 0)
goto L_20;
/* *** FIRST CALL ON SF7HES. SET GSAVE = G, TAKE FIRST STEP *** */
sv7cpy( p, &V[gsave1], gg );
Iv[SWITCH_] = Iv[NFGCAL];
goto L_90;
L_20:
del = V[DELTA];
X[m] = V[XMSAVE];
if (Iv[TOOBIG] == 0)
goto L_40;
/* *** HANDLE OVERSIZE V(DELTA) ***
* */
if (del*X[m] > ZERO)
goto L_30;
/* *** WE ALREADY TRIED SHRINKING V(DELTA), SO QUIT *** */
Iv[FDH] = -2;
goto L_220;
/* *** TRY SHRINKING V(DELTA) *** */
L_30:
del *= NEGPT5;
goto L_100;
L_40:
hes = -Iv[H];
/* *** SET GG = (GG - GSAVE)/DEL ***
* */
for (i = 1; i <= p; i++)
{
Gg[i] = (Gg[i] - V[gsave1])/del;
gsave1 += 1;
}
/* *** ADD GG AS NEW COL. TO FINITE-DIFF. HESSIAN MATRIX ***
* */
k = hes + m*(m - 1)/2;
l = k + m - 2;
if (m == 1)
goto L_70;
/* *** SET H(I,M) = 0.5 * (H(I,M) + GG(I)) FOR I = 1 TO M-1 ***
* */
mm1 = m - 1;
for (i = 1; i <= mm1; i++)
{
V[k] = HALF*(V[k] + Gg[i]);
k += 1;
}
/* *** ADD H(I,M) = GG(I) FOR I = M TO P ***
* */
L_70:
l += 1;
for (i = m; i <= p; i++)
{
V[l] = Gg[i];
l += i;
}
L_90:
m += 1;
Iv[MODE] = m;
if (m > p)
goto L_210;
/* *** CHOOSE NEXT FINITE-DIFFERENCE STEP, RETURN TO GET GG THERE ***
* */
del = V[DELTA0]*fmaxf( ONE/D[m], fabsf( X[m] ) );
if (X[m] < ZERO)
del = -del;
V[XMSAVE] = X[m];
L_100:
X[m] += del;
V[DELTA] = del;
*irt = 2;
goto L_999;
/* *** COMPUTE FINITE-DIFFERENCE HESSIAN USING FUNCTION VALUES ONLY.
* */
L_110:
stp0 = Iv[W] + p - 1;
mm1 = m - 1;
mm1o2 = m*mm1/2;
if (m > 0)
goto L_120;
/* *** FIRST CALL ON SF7HES. *** */
Iv[SAVEI] = 0;
goto L_200;
L_120:
i = Iv[SAVEI];
hes = -Iv[H];
if (i > 0)
goto L_180;
if (Iv[TOOBIG] == 0)
goto L_140;
/* *** HANDLE OVERSIZE STEP ***
* */
stpm = stp0 + m;
del = V[stpm];
if (del*X[XMSAVE] > ZERO)
goto L_130;
/* *** WE ALREADY TRIED SHRINKING THE STEP, SO QUIT *** */
Iv[FDH] = -2;
goto L_220;
/* *** TRY SHRINKING THE STEP *** */
L_130:
del *= NEGPT5;
X[m] = X[XMSAVE] + del;
V[stpm] = del;
*irt = 1;
goto L_999;
/* *** SAVE F(X + STP(M)*E(M)) IN H(P,M) ***
* */
L_140:
pp1o2 = p*(p - 1)/2;
hpm = hes + pp1o2 + mm1;
V[hpm] = V[F];
/* *** START COMPUTING ROW M OF THE FINITE-DIFFERENCE HESSIAN H. ***
* */
hmi = hes + mm1o2;
if (mm1 == 0)
goto L_160;
hpi = hes + pp1o2;
for (i = 1; i <= mm1; i++)
{
V[hmi] = V[FX] - (V[F] + V[hpi]);
hmi += 1;
hpi += 1;
}
L_160:
V[hmi] = V[F] - TWO*V[FX];
/* *** COMPUTE FUNCTION VALUES NEEDED TO COMPLETE ROW M OF H. ***
* */
i = 1;
L_170:
Iv[SAVEI] = i;
stpi = stp0 + i;
V[DELTA] = X[i];
X[i] += V[stpi];
if (i == m)
X[i] = V[XMSAVE] - V[stpi];
*irt = 1;
goto L_999;
L_180:
X[i] = V[DELTA];
if (Iv[TOOBIG] == 0)
goto L_190;
/* *** PUNT IN THE EVENT OF AN OVERSIZE STEP *** */
Iv[FDH] = -2;
goto L_220;
/* *** FINISH COMPUTING H(M,I) ***
* */
L_190:
stpi = stp0 + i;
hmi = hes + mm1o2 + i - 1;
stpm = stp0 + m;
V[hmi] = (V[hmi] + V[F])/(V[stpi]*V[stpm]);
i += 1;
if (i <= m)
goto L_170;
Iv[SAVEI] = 0;
X[m] = V[XMSAVE];
L_200:
m += 1;
Iv[MODE] = m;
if (m > p)
goto L_210;
/* *** PREPARE TO COMPUTE ROW M OF THE FINITE-DIFFERENCE HESSIAN H.
* *** COMPUTE M-TH STEP SIZE STP(M), THEN RETURN TO OBTAIN
* *** F(X + STP(M)*E(M)), WHERE E(M) = M-TH STD. UNIT VECTOR.
* */
del = V[DLTFDC]*fmaxf( ONE/D[m], fabsf( X[m] ) );
if (X[m] < ZERO)
del = -del;
V[XMSAVE] = X[m];
X[m] += del;
stpm = stp0 + m;
V[stpm] = del;
*irt = 1;
goto L_999;
/* *** RESTORE V(F), ETC. ***
* */
L_210:
Iv[FDH] = hes;
L_220:
V[F] = V[FX];
*irt = 3;
if (kind < 0)
goto L_999;
Iv[NFGCAL] = Iv[SWITCH_];
gsave1 = Iv[W] + p;
sv7cpy( p, gg, &V[gsave1] );
goto L_999;
L_999:
return;
/* *** LAST CARD OF SF7HES FOLLOWS *** */
} /* end of function */
/* PARAMETER translations */
#define COVPRT 14
#define NEEDHD 36
#define PRUNIT 21
#define STATPR 23
/* end of PARAMETER translations */
void /*FUNCTION*/ sn2cvp(
long iv[],
long liv,
long lv,
long p,
float v[])
{
long int cov1, i, i1, ii, pu;
float t;
/* OFFSET Vectors w/subscript range: 1 to dimension */
long *const Iv = &iv[0] - 1;
float *const V = &v[0] - 1;
/* end of OFFSET VECTORS */
/* *** PRINT COVARIANCE MATRIX FOR SRN2G ***
* 6/27/90 CLL changed 'SCALE' to 'VARFAC' in output labels.
* ------------------------------------------------------------------ */
long int j, k;
/* INTEGER J */
/* *** LOCAL VARIABLES ***
* */
/* *** IV SUBSCRIPTS ***
* */
/* *** BODY ***
* */
/*++(~.C.) Default UNITNO='(PU,'
*++(.C.) Default UNITNO='(*,'
*++ Replace "(*, " = UNITNO
* */
if (Iv[1] > 8)
goto L_999;
pu = Iv[PRUNIT];
if (pu == 0)
goto L_999;
cov1 = Iv[COVMAT];
if (-2 == cov1)
{
printf("\n ++++++ OVERSIZE STEPS IN COMPUTING COVARIANCE +++++\n");
}
if (Iv[STATPR] == 0)
goto L_30;
if (Iv[NFCOV] > 0)
{
printf("\n %4ld EXTRA FUNC. EVALS FOR COVARIANCE AND DIAGNOSTICS.\n", Iv[NFCOV]);
}
if (Iv[NGCOV] > 0)
{
printf(" %4ld EXTRA GRAD. EVALS FOR COVARIANCE AND DIAGNOSTICS.\n", Iv[NGCOV]);
}
L_30:
if (Iv[COVPRT] <= 0)
goto L_999;
if (Iv[REGD] <= 0 && cov1 <= 0)
goto L_70;
Iv[NEEDHD] = 1;
t = SQ(V[RCOND]);
if (labs( Iv[COVREQ] ) > 2)
goto L_50;
printf("\n RECIPROCAL CONDITION OF F.D. HESSIAN = AT MOST%10.2g\n", t);
goto L_70;
L_50:
printf("\n RECIPROCAL CONDITION OF (J**T)*J = AT LEAST%10.2g\n", t);
L_70:
if ((Iv[COVPRT]%2) == 0)
goto L_999;
Iv[NEEDHD] = 1;
switch (IARITHIF(cov1))
{
case -1: goto L_80;
case 0: goto L_110;
case 1: goto L_130;
}
L_80:
if (-1 == cov1)
{
printf("\n ++++++ INDEFINITE COVARIANCE MATRIX ++++++\n");
}
goto L_999;
L_110:
printf("\n ++++++ COVARIANCE MATRIX NOT COMPUTED ++++++\n");
goto L_999;
L_130:
i = labs( Iv[COVREQ] );
if (i <= 1)
{
printf("\n COVARIANCE = VARFAC * H**-1 * (J**T * J) * H**-1\n WHERE H = F.D. HESSIAN\n \n");
}
if (i == 2)
{
printf(" \n COVARIANCE = VARFAC * H**-1, WHERE H = FINITE-DIFFERENCE HESSIAN\n \n");
}
if (i > 2)
{
printf("\n COVARIANCE = VARFAC * J**T * J\n \n");
}
ii = cov1 - 1;
for (i = 1; i <= p; i++)
{
i1 = ii + 1;
ii += i;
printf(" ROW%3ld ", i);
for (j = i1; j <= ii; j+=5){
for (k = j; k <= (j <= ii - 5 ? j+4 : ii); k++)
printf("%12.3g", V[k]);
printf("\n");
if (j <= ii - 5) printf(" ");}
}
/* WRITE(*, '('' ROW'',I3,2X,5g12.3/(9X,5g12.3))')I,(V(J),J=I1,II)
* */
L_999:
return;
/* *** LAST CARD OF SN2CVP FOLLOWS *** */
} /* end of function */
void /*FUNCTION*/ sn2rdp(
long iv[],
long liv,
long n,
float rd[])
{
long int pu;
/* OFFSET Vectors w/subscript range: 1 to dimension */
long *const Iv = &iv[0] - 1;
float *const Rd = &rd[0] - 1;
/* end of OFFSET VECTORS */
/* *** PRINT REGRESSION DIAGNOSTICS FOR MLPSL AND NL2S1 ***
*
* ------------------------------------------------------------------
*++ Code for .C. is active */
long int j,k;
/*++ End */
/* *** IV SUBSCRIPTS ***
* */
/* DATA COVPRT/14/, NEEDHD/36/, PRUNIT/21/, RDREQ/57/, REGD/67/ */
/* ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++
* */
pu = Iv[PRUNIT];
if (pu == 0)
goto L_999;
if (Iv[COVPRT] < 2)
goto L_999;
if (Iv[REGD] <= 0)
goto L_999;
Iv[NEEDHD] = 1;
printf(" REGRESSION DIAGNOSTIC = SQRT(G(I)**T * H(I)**-1 *G(I))...\n\n");
for(j=0; j < n; j+=6){
for (k = j; k < (j < n - 6 ? j+6 : n); k++)
printf("%12.3g", rd[k]);
printf("\n");}
L_999:
return;
/* WRITE(*, '(6g12.3)') RD
*
* *** LAST CARD OF SN2RDP FOLLOWS *** */
} /* end of function */
/* PARAMETER translations */
#undef ZERO
#define ZERO 0.0e0
/* end of PARAMETER translations */
void /*FUNCTION*/ so7prd(
long l,
long ls,
long pp,
float ss[],
float ww[],
float *yy,
float *z)
{
#define YY(I_,J_) (*(yy+(I_)*(pp)+(J_)))
#define Z(I_,J_) (*(z+(I_)*(pp)+(J_)))
long int i, j, k, m;
float wk, yi;
/* OFFSET Vectors w/subscript range: 1 to dimension */
float *const Ss = &ss[0] - 1;
float *const Ww = &ww[0] - 1;
/* end of OFFSET VECTORS */
/* *** FOR I = 1..L, SET SS = SS + WW(I)*YY(.,I)*(Z(.,I)**T), I.E.,
* *** ADD WW(I) TIMES THE OUTER PRODUCT OF Y(.,I) AND Z(.,I).
*
* ------------------------------------------------------------------ */
/* DIMENSION SS(PP*(PP+1)/2)
* */
for (k = 1; k <= l; k++)
{
wk = Ww[k];
if (wk == ZERO)
goto L_30;
m = 1;
for (i = 1; i <= pp; i++)
{
yi = wk*YY(k - 1,i - 1);
for (j = 1; j <= i; j++)
{
Ss[m] += yi*Z(k - 1,j - 1);
m += 1;
}
}
L_30:
;
}
return;
/* *** LAST CARD OF SO7PRD FOLLOWS *** */
#undef Z
#undef YY
} /* end of function */
/* PARAMETER translations */
#undef ZERO
#define ZERO 0.e0
/* end of PARAMETER translations */
void /*FUNCTION*/ sl7nvr(
long n,
float lin[],
float l[])
{
long int i, ii, im1, j0, j1, jj, k, k0, np1;
float t;
/* OFFSET Vectors w/subscript range: 1 to dimension */
float *const L = &l[0] - 1;
float *const Lin = &lin[0] - 1;
/* end of OFFSET VECTORS */
/* *** COMPUTE LIN = L**-1, BOTH N X N LOWER TRIANG. STORED ***
* *** COMPACTLY BY ROWS. LIN AND L MAY SHARE THE SAME STORAGE. ***
*
* *** PARAMETERS ***
*
* ------------------------------------------------------------------ */
/* DIMENSION L(N*(N+1)/2), LIN(N*(N+1)/2)
*
* *** LOCAL VARIABLES ***
* */
/* *** BODY ***
* */
np1 = n + 1;
j0 = n*(np1)/2;
for (ii = 1; ii <= n; ii++)
{
i = np1 - ii;
Lin[j0] = ONE/L[j0];
if (i <= 1)
goto L_999;
j1 = j0;
im1 = i - 1;
for (jj = 1; jj <= im1; jj++)
{
t = ZERO;
j0 = j1;
k0 = j1 - jj;
for (k = 1; k <= jj; k++)
{
t += -L[k0]*Lin[j0];
j0 -= 1;
k0 += k - i;
}
Lin[j0] = t/L[k0];
}
j0 -= 1;
}
L_999:
return;
/* *** LAST CARD OF SL7NVR FOLLOWS *** */
} /* end of function */
void /*FUNCTION*/ sl7tsq(
long n,
float a[],
float l[])
{
long int i, i1, ii, iim1, j, k, m;
float lii, lj;
/* OFFSET Vectors w/subscript range: 1 to dimension */
float *const A = &a[0] - 1;
float *const L = &l[0] - 1;
/* end of OFFSET VECTORS */
/* *** SET A TO LOWER TRIANGLE OF (L**T) * L ***
*
* *** L = N X N LOWER TRIANG. MATRIX STORED ROWWISE. ***
* *** A IS ALSO STORED ROWWISE AND MAY SHARE STORAGE WITH L. ***
*
* ------------------------------------------------------------------ */
/* DIMENSION A(N*(N+1)/2), L(N*(N+1)/2)
* */
ii = 0;
for (i = 1; i <= n; i++)
{
i1 = ii + 1;
ii += i;
m = 1;
if (i == 1)
goto L_30;
iim1 = ii - 1;
for (j = i1; j <= iim1; j++)
{
lj = L[j];
for (k = i1; k <= j; k++)
{
A[m] += lj*L[k];
m += 1;
}
}
L_30:
lii = L[ii];
for (j = i1; j <= ii; j++)
{
A[j] = lii*L[j];
}
}
return;
/* *** LAST CARD OF SL7TSQ FOLLOWS *** */
} /* end of function */