/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:22 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "srnsgb.h" #include /* PARAMETER translations */ #define AR 110 #define CSAVE 105 #define D 27 #define IERS 108 #define IPIVS 109 #define IV1SAV 104 #define IVNEED 3 #define J 70 #define MODE 35 #define NEXTIV 46 #define NEXTV 47 #define NFCALL 6 #define NFGCAL 7 #define PERM 58 #define R 61 #define REGD 67 #define REGD0 82 #define TOOBIG 2 #define VNEED 4 /* end of PARAMETER translations */ void /*FUNCTION*/ srnsgb( float *a, float alf[], float b[][2], float c[], float *da, long in[][2], long iv[], long l, long l1, long la, long liv, long lv, long n, long nda, long p, float v[], float y[]) { #define A(I_,J_) (*(a+(I_)*(la)+(J_))) #define DA(I_,J_) (*(da+(I_)*(la)+(J_))) long int ar1, csave1, d1, dr1, dr1l, i, i1, ier, ipiv1, iv1, j1, jlen, k, ll1o2, md, n1, nml, nran, r1, r1l, rd1; float singtl, t; static float machep = -1.e0; static float negone = -1.e0; static float sngfac = 1.e2; static float zero = 0.e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Alf = &alf[0] - 1; float *const C = &c[0] - 1; long *const Iv = &iv[0] - 1; float *const V = &v[0] - 1; float *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. *>> 2000-12-03 SRNSGB Krogh Removed references to IV(RESTOR). *>> 1998-10-29 SRNSGB Krogh Moved external statement up for mangle. *>> 1994-10-20 SRNSGB Krogh Changes to use M77CON *>> 1990-06-12 SRNSGB CLL @ JPL *>> 1990-02-14 CLL @ JPL *** from netlib, Fri Feb 9 13:10:09 EST 1990 *** * * *** ITERATION DRIVER FOR SEPARABLE NONLINEAR LEAST SQUARES, * *** WITH SIMPLE BOUNDS ON THE NONLINEAR VARIABLES. * * *** PARAMETER DECLARATIONS *** * */ /* *** PURPOSE *** * * GIVEN A SET OF N OBSERVATIONS Y(1)....Y(N) OF A DEPENDENT VARIABLE * T(1)...T(N), SRNSGB 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 * * SUBJECT TO THE SIMPLE BOUND CONSTRAINTS B(1,I) .LE. X(I) .LE. B(2,I), * I = 1(1)P. * * 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 SRNSGB 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 * SRNSGB SHOULD RETURN FOR THEM. * IV (I/O) INTEGER PARAMETER AND SCRATCH VECTOR. SRNSGB 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), SRNSGB 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 + 4*P. * LV (IN) LENGTH OF V. MUST BE AT LEAST * 105 + 2*N + L*(L+3)/2 + P*(2*P + 21 + N). * 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 *** * */ /*--S replaces "?": ?RNSGB,?ITSUM,?IVSET,?L7ITV,?L7SVN,?L7SVX,?Q7APL *--& ?Q7RFH,?R7MDC,?RN2GB,?S7CPR,?V2AXY,?V7CPY,?V7PRM *--& ?V7SCP * * SIVSET.... SUPPLIES DEFAULT PARAMETER VALUES. * SITSUM.... PRINTS ITERATION SUMMARY, INITIAL AND FINAL ALF. * SL7ITV... APPLIES INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX. * SL7SVX... ESTIMATES LARGEST SING. VALUE OF LOWER TRIANG. MATRIX. * SL7SVN... ESTIMATES SMALLEST SING. VALUE OF LOWER TRIANG. MATRIX. * SRN2GB... UNDERLYING NONLINEAR LEAST-SQUARES SOLVER. * SQ7APL... APPLIES HOUSEHOLDER TRANSFORMS STORED BY SQ7RFH. * SQ7RFH.... COMPUTES QR FACT. VIA HOUSEHOLDER TRANSFORMS WITH PIVOTING. * SR7MDC... RETURNS MACHINE-DEP. CONSTANTS. * SS7CPR... PRINTS LINEAR PARAMETERS AT SOLUTION. * SV2AXY.... ADDS MULTIPLE OF ONE VECTOR TO ANOTHER. * SV7CPY.... COPIES ONE VECTOR TO ANOTHER. * SV7PRM.... PERMUTES VECTOR. * DV7SCL... SCALES AND COPIES ONE VECTOR TO ANOTHER. * * *** LOCAL VARIABLES *** * */ /* *** SUBSCRIPTS FOR IV AND V *** * */ /* *** IV SUBSCRIPT VALUES *** * */ /*++++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++ * * */ if (Iv[1] == 0) sivset( 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_240; if (l < 0) goto L_240; if (n <= l) goto L_240; if (la < n) goto L_240; 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_240; ll1o2 = l*(l + 1)/2; jlen = n*p; i = l + p; 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; srn2gb( b, 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; L_30: srn2gb( b, &V[d1], &V[dr1l], iv, liv, lv, nml, n, &n1, &nml, p, &V[r1l], &V[rd1], v, alf ); /* IF (abs(IV(RESTOR)-2) .EQ. 1 .AND. L .GT. 0) * 1 CALL SV7CPY(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) .EQ. 2 .AND. L .GT. 0) CALL SV7CPY(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) sv7cpy( n, &V[r1], y ); if (l1 > l) sv2axy( 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]; sq7rfh( &ier, &Iv[ipiv1], n, la, 0, l, a, &V[ar1], ll1o2, c ); /* *** DETERMINE NUMERICAL RANK OF A *** * */ if (machep <= zero) machep = sr7mdc( 3 ); singtl = sngfac*(float)( max( l, n ) )*machep; k = l; if (ier != 0) k = ier - 1; L_70: if (k <= 0) goto L_90; t = sl7svx( k, &V[ar1], c, c ); if (t > zero) t = sl7svn( 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; sv7scp( l - k, &C[k + 1], zero ); L_100: Iv[IERS] = ier; if (k <= 0) goto L_110; /* *** APPLY HOUSEHOLDER TRANSFORMATONS TO RESIDUALS... * */ sq7apl( la, n, k, a, &V[r1], ier ); /* *** COMPUTING C NOW MAY SAVE A FUNCTION EVALUATION AT * *** THE LAST ITERATION. * */ sl7itv( k, c, &V[ar1], &V[r1] ); sv7prm( 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++) { sv2axy( 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: sv7scp( n*p, &V[dr1], zero ); goto L_999; /* *** COMPUTE NEW JACOBIAN *** * */ L_170: ; if (nda <= 0) goto L_240; 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]; sv2axy( n, &V[k], t, &DA(i - 1,0), &V[k] ); L_180: ; } if (iv1 == 2) goto L_190; Iv[1] = iv1; goto L_999; L_190: if (l <= 0) goto L_30; 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++) { sq7apl( la, n, nran, a, &V[k], ier ); k += n; } L_210: sv7cpy( 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 (Iv[REGD] == 1) Iv[REGD] = rd1; if (Iv[1] <= 11) ss7cpr( c, iv, l, liv ); goto L_999; L_240: Iv[1] = 66; sitsum( v, v, iv, liv, lv, p, v, alf ); L_999: return; /* *** LAST CARD OF SRNSGB FOLLOWS *** */ #undef DA #undef A } /* end of function */