/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:18 */
/*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 "drnsg.h"
#include <stdlib.h>
/* PARAMETER translations */
#define AR 110
#define CNVCOD 55
#define COVMAT 26
#define COVREQ 15
#define CSAVE 105
#define CVRQSV 106
#define D 27
#define FDH 74
#define H 56
#define IERS 108
#define IPIVS 109
#define IV1SAV 104
#define IVNEED 3
#define J 70
#define LMAT 42
#define MODE 35
#define NEXTIV 46
#define NEXTV 47
#define NFCALL 6
#define NFCOV 52
#define NFGCAL 7
#define NGCALL 30
#define NGCOV 53
#define PERM 58
#define R 61
#define RCOND 53
#define RDREQ 57
#define RDRQSV 107
#define REGD 67
#define REGD0 82
#define RESTOR 9
#define TOOBIG 2
#define VNEED 4
/* end of PARAMETER translations */
void /*FUNCTION*/ drnsg(
double *a,
double alf[],
double c[],
double *da,
long in[][2],
long iv[],
long l,
long l1,
long la,
long liv,
long lv,
long n,
long nda,
long p,
double v[],
double y[])
{
#define A(I_,J_) (*(a+(I_)*(la)+(J_)))
#define DA(I_,J_) (*(da+(I_)*(la)+(J_)))
LOGICAL32 nocov;
long int ar1, csave1, d1, dr1, dr1l, dri, dri1, fdh0, hsave, i,
i1, ier, ipiv1, iv1, j1, jlen, k, lh, li, ll1o2, md, n1, nml,
nran, pp, pp1, r1, r1l, rd1, temp1;
double singtl, t;
static double machep = -1.e0;
static double negone = -1.e0;
static double sngfac = 1.e2;
static double zero = 0.e0;
/* OFFSET Vectors w/subscript range: 1 to dimension */
double *const Alf = &alf[0] - 1;
double *const C = &c[0] - 1;
long *const Iv = &iv[0] - 1;
double *const V = &v[0] - 1;
double *const Y = &y[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.
*>> 1998-10-29 DRNSG Krogh Moved external statement up for mangle.
*>> 1994-10-20 DRNSG Krogh Changes to use M77CON
*>> 1990-06-12 DRNSG 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
* (Calls DN2RDP with 4, rather than 6, arguments. -- CLL)
*
* *** ITERATION DRIVER FOR SEPARABLE NONLINEAR LEAST SQUARES.
*
* *** PARAMETER DECLARATIONS ***
* */
/* *** PURPOSE ***
*
* GIVEN A SET OF N OBSERVATIONS Y(1)....Y(N) OF A DEPENDENT VARIABLE
* T(1)...T(N), DRNSG ATTEMPTS TO COMPUTE A LEAST SQUARES FIT
* TO A FUNCTION ETA (THE MODEL) WHICH IS A LINEAR COMBINATION
*
* L
* ETA(C,ALF,T) = SUM C * PHI(ALF,T) +PHI (ALF,T)
* J=1 J J L+1
*
* OF NONLINEAR FUNCTIONS PHI(J) DEPENDENT ON T AND ALF(1),...,ALF(P)
* (.E.G. A SUM OF EXPONENTIALS OR GAUSSIANS). THAT IS, IT DETERMINES
* NONLINEAR PARAMETERS ALF WHICH MINIMIZE
*
* 2 N 2
* NORM(RESIDUAL) = SUM (Y - ETA(C,ALF,T )).
* I=1 I I
*
* THE (L+1)ST TERM IS OPTIONAL.
*
*
* *** PARAMETERS ***
*
* A (IN) MATRIX PHI(ALF,T) OF THE MODEL.
* ALF (I/O) NONLINEAR PARAMETERS.
* INPUT = INITIAL GUESS,
* OUTPUT = BEST ESTIMATE FOUND.
* C (OUT) LINEAR PARAMETERS (ESTIMATED).
* DA (IN) DERIVATIVES OF COLUMNS OF A WITH RESPECT TO COMPONENTS
* OF ALF, AS SPECIFIED BY THE IN ARRAY...
* IN (IN) WHEN DRNSG IS CALLED WITH IV(1) = 2 OR -2, THEN FOR
* I = 1(1)NDA, COLUMN I OF DA IS THE PARTIAL
* DERIVATIVE WITH RESPECT TO ALF(IN(1,I)) OF COLUMN
* IN(2,I) OF A, UNLESS IV(1,I) IS NOT POSITIVE (IN
* WHICH CASE COLUMN I OF DA IS IGNORED. IV(1) = -2
* MEANS THERE ARE MORE COLUMNS OF DA TO COME AND
* DRNSG SHOULD RETURN FOR THEM.
* IV (I/O) INTEGER PARAMETER AND SCRATCH VECTOR. DRNSG RETURNS
* WITH IV(1) = 1 WHEN IT WANTS A TO BE EVALUATED AT
* ALF AND WITH IV(1) = 2 WHEN IT WANTS DA TO BE
* EVALUATED AT ALF. WHEN CALLED WITH IV(1) = -2
* (AFTER A RETURN WITH IV(1) = 2), DRNSG RETURNS
* WITH IV(1) = -2 TO GET MORE COLUMNS OF DA.
* L (IN) NUMBER OF LINEAR PARAMETERS TO BE ESTIMATED.
* L1 (IN) L+1 IF PHI(L+1) IS IN THE MODEL, L IF NOT.
* LA (IN) LEAD DIMENSION OF A. MUST BE AT LEAST N.
* LIV (IN) LENGTH OF IV. MUST BE AT LEAST 110 + L + P.
* LV (IN) LENGTH OF V. MUST BE AT LEAST
* 105 + 2*N + JLEN + L*(L+3)/2 + P*(2*P + 17),
* WHERE JLEN = (L+P)*(N+L+P+1), UNLESS NEITHER A
* COVARIANCE MATRIX NOR REGRESSION DIAGNOSTICS ARE
* REQUESTED, IN WHICH CASE JLEN = N*P.
* N (IN) NUMBER OF OBSERVATIONS.
* NDA (IN) NUMBER OF COLUMNS IN DA AND IN.
* P (IN) NUMBER OF NONLINEAR PARAMETERS TO BE ESTIMATED.
* V (I/O) FLOATING-POINT PARAMETER AND SCRATCH VECTOR.
* Y (IN) RIGHT-HAND SIDE VECTOR.
*
* *** EXTERNAL SUBROUTINES ***
* */
/*--D replaces "?": ?RNSG ,?C7VFN,?D7TPR,?ITSUM,?IVSET,?L7ITV,?L7SRT
*--& ?L7SVN,?L7SVX,?N2CVP,?N2LRD,?N2RDP,?Q7APL,?Q7RAD
*--& ?Q7RFH,?R7MDC,?RN2G ,?S7CPR,?V2AXY,?V7CPY,?V7PRM
*--& ?V7SCL,?V7SCP
* DC7VFN... FINISHES COVARIANCE COMPUTATION.
* DIVSET.... SUPPLIES DEFAULT PARAMETER VALUES.
* DD7TPR... RETURNS INNER PRODUCT OF TWO VECTORS.
* DITSUM.... PRINTS ITERATION SUMMARY, INITIAL AND FINAL ALF.
* DL7ITV... APPLIES INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX.
* DL7SRT.... COMPUTES (PARTIAL) CHOLESKY FACTORIZATION.
* DL7SVX... ESTIMATES LARGEST SING. VALUE OF LOWER TRIANG. MATRIX.
* DL7SVN... ESTIMATES SMALLEST SING. VALUE OF LOWER TRIANG. MATRIX.
* DN2CVP... PRINTS COVARIANCE MATRIX.
* DN2LRD... COMPUTES COVARIANCE AND REGRESSION DIAGNOSTICS.
* DN2RDP... PRINTS REGRESSION DIAGNOSTICS.
* DRN2G... UNDERLYING NONLINEAR LEAST-SQUARES SOLVER.
* DQ7APL... APPLIES HOUSEHOLDER TRANSFORMS STORED BY DQ7RFH.
* DQ7RFH.... COMPUTES QR FACT. VIA HOUSEHOLDER TRANSFORMS WITH PIVOTING.
* DQ7RAD.... QR FACT., NO PIVOTING.
* DR7MDC... RETURNS MACHINE-DEP. CONSTANTS.
* DS7CPR... PRINTS LINEAR PARAMETERS AT SOLUTION.
* DV2AXY.... ADDS MULTIPLE OF ONE VECTOR TO ANOTHER.
* DV7CPY.... COPIES ONE VECTOR TO ANOTHER.
* DV7PRM.... PERMUTES A VECTOR.
* DV7SCL... SCALES AND COPIES ONE VECTOR TO ANOTHER.
* DV7SCP... SETS ALL COMPONENTS OF A VECTOR TO A SCALAR.
*
* *** LOCAL VARIABLES ***
* */
/* *** SUBSCRIPTS FOR IV AND V ***
* */
/* *** IV SUBSCRIPT VALUES ***
*
* /6
* DATA AR/110/, CNVCOD/55/, COVMAT/26/, COVREQ/15/, CSAVE/105/,
* 1 CVRQSV/106/, D/27/, FDH/74/, H/56/, IERS/108/, IPIVS/109/,
* 2 IV1SAV/104/, IVNEED/3/, J/70/, LMAT/42/, MODE/35/,
* 3 NEXTIV/46/, NEXTV/47/, NFCALL/6/, NFCOV/52/, NFGCAL/7/,
* 4 NGCALL/30/, NGCOV/53/, PERM/58/, R/61/, RCOND/53/, RDREQ/57/,
* 5 RDRQSV/107/, REGD/67/, REGD0/82/, RESTOR/9/, TOOBIG/2/,
* 6 VNEED/4/
* /7 */
/*/ */
/*++++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++
*
* */
if (Iv[1] == 0)
divset( 1, iv, liv, lv, v );
n1 = 1;
nml = n;
iv1 = Iv[1];
if (iv1 <= 2)
goto L_20;
/* *** CHECK INPUT INTEGERS ***
* */
if (p <= 0)
goto L_370;
if (l < 0)
goto L_370;
if (n <= l)
goto L_370;
if (la < n)
goto L_370;
if (iv1 < 12)
goto L_20;
if (iv1 == 14)
goto L_20;
if (iv1 == 12)
Iv[1] = 13;
/* *** FRESH START -- COMPUTE STORAGE REQUIREMENTS ***
* */
if (Iv[1] > 16)
goto L_370;
ll1o2 = l*(l + 1)/2;
jlen = n*p;
i = l + p;
if (Iv[RDREQ] > 0 && Iv[COVREQ] != 0)
jlen = i*(n + i + 1);
if (Iv[1] != 13)
goto L_10;
Iv[IVNEED] += l;
Iv[VNEED] += p + 2*n + jlen + ll1o2 + l;
L_10:
if (Iv[PERM] <= AR)
Iv[PERM] = AR + 1;
drn2g( v, v, iv, liv, lv, n, n, &n1, &nml, p, v, v, v, alf );
if (Iv[1] != 14)
goto L_999;
/* *** STORAGE ALLOCATION ***
* */
Iv[IPIVS] = Iv[NEXTIV];
Iv[NEXTIV] += l;
Iv[D] = Iv[NEXTV];
Iv[REGD0] = Iv[D] + p;
Iv[AR] = Iv[REGD0] + n;
Iv[CSAVE] = Iv[AR] + ll1o2;
Iv[J] = Iv[CSAVE] + l;
Iv[R] = Iv[J] + jlen;
Iv[NEXTV] = Iv[R] + n;
Iv[IERS] = 0;
if (iv1 == 13)
goto L_999;
/* *** SET POINTERS INTO IV AND V ***
* */
L_20:
ar1 = Iv[AR];
d1 = Iv[D];
dr1 = Iv[J];
dr1l = dr1 + l;
r1 = Iv[R];
r1l = r1 + l;
rd1 = Iv[REGD0];
csave1 = Iv[CSAVE];
nml = n - l;
if (iv1 <= 2)
goto L_50;
/* *** IF F.D. HESSIAN WILL BE NEEDED (FOR COVARIANCE OR REG.
* *** DIAGNOSTICS), HAVE DRN2G COMPUTE ONLY THE PART CORRESP.
* *** TO ALF WITH C FIXED...
* */
if (l <= 0)
goto L_30;
Iv[CVRQSV] = Iv[COVREQ];
if (labs( Iv[COVREQ] ) >= 3)
Iv[COVREQ] = 0;
Iv[RDRQSV] = Iv[RDREQ];
if (Iv[RDREQ] > 0)
Iv[RDREQ] = -1;
L_30:
drn2g( &V[d1], &V[dr1l], iv, liv, lv, nml, n, &n1, &nml, p, &V[r1l],
&V[rd1], v, alf );
if (labs( Iv[RESTOR] - 2 ) == 1 && l > 0)
dv7cpy( l, c, &V[csave1] );
iv1 = Iv[1];
switch (IARITHIF(iv1 - 2))
{
case -1: goto L_40;
case 0: goto L_150;
case 1: goto L_230;
}
/* *** NEW FUNCTION VALUE (RESIDUAL) NEEDED ***
* */
L_40:
Iv[IV1SAV] = Iv[1];
Iv[1] = labs( iv1 );
if (Iv[RESTOR] == 2 && l > 0)
dv7cpy( l, &V[csave1], c );
goto L_999;
/* *** COMPUTE NEW RESIDUAL OR GRADIENT ***
* */
L_50:
Iv[1] = Iv[IV1SAV];
md = Iv[MODE];
if (md <= 0)
goto L_60;
nml = n;
dr1l = dr1;
r1l = r1;
L_60:
if (Iv[TOOBIG] != 0)
goto L_30;
if (labs( iv1 ) == 2)
goto L_170;
/* *** COMPUTE NEW RESIDUAL ***
* */
if (l1 <= l)
dv7cpy( n, &V[r1], y );
if (l1 > l)
dv2axy( n, &V[r1], negone, &A(l1 - 1,0), y );
if (md > 0)
goto L_120;
ier = 0;
if (l <= 0)
goto L_110;
ll1o2 = l*(l + 1)/2;
ipiv1 = Iv[IPIVS];
dq7rfh( &ier, &Iv[ipiv1], n, la, 0, l, a, &V[ar1], ll1o2, c );
/* *** DETERMINE NUMERICAL RANK OF A ***
* */
if (machep <= zero)
machep = dr7mdc( 3 );
singtl = sngfac*(double)( max( l, n ) )*machep;
k = l;
if (ier != 0)
k = ier - 1;
L_70:
if (k <= 0)
goto L_90;
t = dl7svx( k, &V[ar1], c, c );
if (t > zero)
t = dl7svn( k, &V[ar1], c, c )/t;
if (t > singtl)
goto L_80;
k -= 1;
goto L_70;
/* *** RECORD RANK IN IV(IERS)... IV(IERS) = 0 MEANS FULL RANK,
* *** IV(IERS) .GT. 0 MEANS RANK IV(IERS) - 1.
* */
L_80:
if (k >= l)
goto L_100;
L_90:
ier = k + 1;
dv7scp( l - k, &C[k + 1], zero );
L_100:
Iv[IERS] = ier;
if (k <= 0)
goto L_110;
/* *** APPLY HOUSEHOLDER TRANSFORMATONS TO RESIDUALS...
* */
dq7apl( la, n, k, a, &V[r1], ier );
/* *** COMPUTING C NOW MAY SAVE A FUNCTION EVALUATION AT
* *** THE LAST ITERATION.
* */
dl7itv( k, c, &V[ar1], &V[r1] );
dv7prm( l, &Iv[ipiv1], c );
L_110:
if (Iv[1] < 2)
goto L_220;
goto L_999;
/* *** RESIDUAL COMPUTATION FOR F.D. HESSIAN ***
* */
L_120:
if (l <= 0)
goto L_140;
for (i = 1; i <= l; i++)
{
dv2axy( n, &V[r1], -C[i], &A(i - 1,0), &V[r1] );
}
L_140:
if (Iv[1] > 0)
goto L_30;
Iv[1] = 2;
goto L_160;
/* *** NEW GRADIENT (JACOBIAN) NEEDED ***
* */
L_150:
Iv[IV1SAV] = iv1;
if (Iv[NFGCAL] != Iv[NFCALL])
Iv[1] = 1;
L_160:
dv7scp( n*p, &V[dr1], zero );
goto L_999;
/* *** COMPUTE NEW JACOBIAN ***
* */
L_170:
nocov = md <= p || labs( Iv[COVREQ] ) >= 3;
fdh0 = dr1 + n*(p + l);
if (nda <= 0)
goto L_370;
for (i = 1; i <= nda; i++)
{
i1 = in[i - 1][0] - 1;
if (i1 < 0)
goto L_180;
j1 = in[i - 1][1];
k = dr1 + i1*n;
t = negone;
if (j1 <= l)
t = -C[j1];
dv2axy( n, &V[k], t, &DA(i - 1,0), &V[k] );
if (nocov)
goto L_180;
if (j1 > l)
goto L_180;
/* *** ADD IN (L,P) PORTION OF SECOND-ORDER PART OF HESSIAN
* *** FOR COVARIANCE OR REG. DIAG. COMPUTATIONS... */
j1 += p;
k = fdh0 + j1*(j1 - 1)/2 + i1;
V[k] -= dd7tpr( n, &V[r1], &DA(i - 1,0) );
L_180:
;
}
if (iv1 == 2)
goto L_190;
Iv[1] = iv1;
goto L_999;
L_190:
if (l <= 0)
goto L_30;
if (md > p)
goto L_240;
if (md > 0)
goto L_30;
k = dr1;
ier = Iv[IERS];
nran = l;
if (ier > 0)
nran = ier - 1;
if (nran <= 0)
goto L_210;
for (i = 1; i <= p; i++)
{
dq7apl( la, n, nran, a, &V[k], ier );
k += n;
}
L_210:
dv7cpy( l, &V[csave1], c );
L_220:
if (ier == 0)
goto L_30;
/* *** ADJUST SUBSCRIPTS DESCRIBING R AND DR...
* */
nran = ier - 1;
dr1l = dr1 + nran;
nml = n - nran;
r1l = r1 + nran;
goto L_30;
/* *** CONVERGENCE OR LIMIT REACHED ***
* */
L_230:
if (l <= 0)
goto L_350;
Iv[COVREQ] = Iv[CVRQSV];
Iv[RDREQ] = Iv[RDRQSV];
if (Iv[1] > 6)
goto L_360;
if ((Iv[RDREQ]%4) == 0)
goto L_360;
if (Iv[FDH] <= 0 && labs( Iv[COVREQ] ) < 3)
goto L_360;
if (Iv[REGD] > 0)
goto L_360;
if (Iv[COVMAT] > 0)
goto L_360;
/* *** PREPARE TO FINISH COMPUTING COVARIANCE MATRIX AND REG. DIAG. ***
* */
pp = l + p;
i = 0;
if ((Iv[RDREQ]%4) >= 2)
i = 1;
if ((Iv[RDREQ]%2) == 1 && labs( Iv[COVREQ] ) == 1)
i += 2;
Iv[MODE] = pp + i;
i = dr1 + n*pp;
k = p*(p + 1)/2;
i1 = Iv[LMAT];
dv7cpy( k, &V[i], &V[i1] );
i += k;
dv7scp( pp*(pp + 1)/2 - k, &V[i], zero );
Iv[NFCOV] += 1;
Iv[NFCALL] += 1;
Iv[NFGCAL] = Iv[NFCALL];
Iv[CNVCOD] = Iv[1];
Iv[IV1SAV] = -1;
Iv[1] = 1;
Iv[NGCALL] += 1;
Iv[NGCOV] += 1;
goto L_999;
/* *** FINISH COVARIANCE COMPUTATION ***
* */
L_240:
i = dr1 + n*p;
for (i1 = 1; i1 <= l; i1++)
{
dv7scl( n, &V[i], negone, &A(i1 - 1,0) );
i += n;
}
pp = l + p;
hsave = Iv[H];
k = dr1 + n*pp;
lh = pp*(pp + 1)/2;
if (labs( Iv[COVREQ] ) < 3)
goto L_270;
i = Iv[MODE] - 4;
if (i >= pp)
goto L_260;
dv7scp( lh, &V[k], zero );
dq7rad( n, n, pp, v, FALSE, &V[k], &V[dr1], v );
Iv[MODE] = i + 8;
Iv[1] = 2;
Iv[NGCALL] += 1;
Iv[NGCOV] += 1;
goto L_160;
L_260:
Iv[MODE] = i;
goto L_300;
L_270:
pp1 = p + 1;
dri = dr1 + n*p;
li = k + p*pp1/2;
for (i = pp1; i <= pp; i++)
{
dri1 = dr1;
for (i1 = 1; i1 <= i; i1++)
{
V[li] += dd7tpr( n, &V[dri], &V[dri1] );
li += 1;
dri1 += n;
}
dri += n;
}
dl7srt( pp1, pp, &V[k], &V[k], &i );
if (i != 0)
goto L_310;
L_300:
temp1 = k + lh;
t = dl7svn( pp, &V[k], &V[temp1], &V[temp1] );
if (t <= zero)
goto L_310;
t /= dl7svx( pp, &V[k], &V[temp1], &V[temp1] );
V[RCOND] = t;
if (t > dr7mdc( 4 ))
goto L_320;
L_310:
Iv[REGD] = -1;
Iv[COVMAT] = -1;
Iv[FDH] = -1;
goto L_340;
L_320:
Iv[H] = temp1;
Iv[FDH] = labs( hsave );
if (Iv[MODE] - pp < 2)
goto L_330;
i = Iv[H];
dv7scp( lh, &V[i], zero );
L_330:
dn2lrd( &V[dr1], iv, &V[k], lh, liv, lv, n, n, pp, &V[r1], &V[rd1],
v );
L_340:
dc7vfn( iv, &V[k], lh, liv, lv, n, pp, v );
Iv[H] = hsave;
L_350:
if (Iv[REGD] == 1)
Iv[REGD] = rd1;
L_360:
if (Iv[1] <= 11)
ds7cpr( c, iv, l, liv );
if (Iv[1] > 6)
goto L_999;
dn2cvp( iv, liv, lv, p + l, v );
dn2rdp( iv, liv, n, &V[rd1] );
goto L_999;
L_370:
Iv[1] = 66;
ditsum( v, v, iv, liv, lv, p, v, alf );
L_999:
return;
/* *** LAST CARD OF DRNSG FOLLOWS *** */
#undef DA
#undef A
} /* end of function */