/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:47 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "snqsol.h" #include #include #include /* PARAMETER translations */ #define FAC1 0.75e0 #define IWDIAG 4 #define IWTOLT 3 #define MAXF1 200 /* end of PARAMETER translations */ void /*FUNCTION*/ snqsol( void (*snqfj)(long,float[],float[],float*,long*), long n, float x[], float fvec[], float xtol, long iopt[], float w[], long idimw) { LOGICAL32 haved, havej, jpos, trace; long int _l0, dmode, iwa1, iwa2, iwa3, iwa4, iwfjac, iwgnst, iwqtf, iwr, j, jabs, k, maxfev, ml, mu, ni, nprint; float epsfcn, factor; static float epsmch = 0.0e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Fvec = &fvec[0] - 1; long *const Iopt = &iopt[0] - 1; float *const W = &w[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. *>> 2001-05-25 SNQSOL Krogh Minor change for making .f90 version. *>> 2000-12-01 SNQSOL Krogh Removed unused parameter P001. *>> 1996-05-16 SNQSOL Krogh Changes to use .C. and C%%. *>> 1996-03-30 SNQSOL Krogh Added external stmts. SIN => VSIN, etc. *>> 1994-11-02 SNQSOL Krogh Changes to use M77CON *>> 1992-04-27 SNQSOL CLL Deleted unreferenced stmt label. *>> 1992-04-07 CAO Extra comma in Print removed (error from VAX compile) *>> 1992-01-15 CLL *>> 1991-12-18 CLL & FTK Adding treatment of slow convergence to 0. *>> 1991-12-05 CLL & FTK Adding Option vector interface. *>> 1990-04-20 CLL@JPL Adapting code from Minpack for MATH77 * * Solves a system of N nonlinear equations in N unknowns. * SNQSOL is the the user-interface subroutine. It calls SNQSL1 which * contains the top-level logic of the solution algorithm. * SNQSOL & SNQSL1 also need: * Other subroutines that are in this file: * SNQFDJ, SNQDOG, SNQQFM, SNQQRF, SNQUPD. * Other subprograms from the MATH77 library: SNRM2, SERV1, * [D/R]1MACH (Fortan 77 only), IERV1, & IERM1. * A user-provided subroutine: SNQFJ. * * Most of these subprograms are derived from MINPACK-1. * MINPACK-1, 1980, was developed by Jorge J. More', * Burton S. Garbow, and Kenneth E. Hillstrom, Argonne Nat'l Lab. * The MINPACK-1 code was obtained as FILE05 from MINPACK/EX from * Netlib, downloaded to JPL on Tue Feb 6 12:17:45 EST 1990. * * Old Name New Name * -------- -------- * HYBRJ1, HYBRD1 SNQSOL (Completely redesigned.) * HYBRJ, HYBRD SNQSL1 (Algorithm and code changes.) * DOGLEG SNQDOG * ENORM SNRM2 in BLAS and MATH77 * FDJAC1 SNQFDJ * QFORM SNQQFM * QRFAC SNQQRF * R1MPYQ SNQAQ * R1UPDT SNQUPD * [D/S]PMPAR [D/R]1MACH in file amach.f (Fortran 77 only) * FCN SNQFJ * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Arguments for SNQSOL * * call SNQSOL(SNQFJ, N, X, FVEC, XTOL, IOPT, W, IDIMW) * * SNQFJ Name of user-supplied subroutine. * * N [in] Problem size * X(N) [inout] Initial and final x-vector. * FVEC(N) [out] Final F values. * XTOL [in] Rel. Conv. tolerance on weighted X * IOPT() [inout] First 3 elements contain output values. * IOPT(1) = INFO. Output status. * IOPT(2) = NFEV. No. of F evals used. * IOPT(3) = NJEV. No. of evals of Jacobian. * * Ramaining elements in IOPT() select options. * * Option No. of Affected variables Affected variables * Number arguments in SNQSOL. in SNQSL1. * 1 0 HAVEJ HAVEJ * 2 1 DMODE, HAVED, W(4:3+N) HAVED, DIAG(1:N) * 3 1 NPRINT NPRINT * 4 1 MAXFEV MAXFEV * 5 2 ML, MU ML, MU * 6 0 W(1) EPSFCN * 7 0 W(2) FACTOR * 8 0 TRACE TRACE * * Functionality of options, listed by option numbers in square brackets. * [1] If set, user is not computing a Jacobian. * This subr sets HAVEJ = .false. * [2] Arg: DMODE = 1, 2, or 3. * 1. This subr sets DIAG() to all ones and HAVED = .true. * 2. User has set DIAG(). This subr sets HAVED = .true. * 3. This subr sets HAVED = .false. so SNQSL1 will set * DIAG() dynamically. * [3] Arg: NPRINT Print control. * [4] Arg: MAXFEV Limit on no. of F evals. * [5] Args: ML & MU Band structure. * [6] If set means EPSFCN has been set in W(1). * [7] If set means FACTOR has been set in W(2). * [8] If set, this subr sets TRACE = .true., else this subr * sets TRACE = .false. When TRACE is .true., SNQSL1 prints * detailed intermediate results. * * W() [inout] W(1) and W(2) may be used to pass EPSFCN and FACTOR * to the subroutine. W(3) contains TOLTST on return. * W( 4 : 3+(15*N + 3*N**2)/2 ) is used as work space. * * EPSFCN W(1) Error in F evals. Used in computing * approx derivs. * FACTOR W(2) Algorithm parameter. * TOLTST W(3) Output. Final value of quantity compared * with XTOL for convergence test. * DIAG(N) W(4:N+3) Scaling values. May be input or */ /* computed. See option 2. * WA1(N) W() Work space of length N. * WA2(N) W() Work space of length N. * WA3(N) W() Work space of length N. * WA4(N) W() Work space of length N. * GNSTEP(N) W() Work space of length N. * QTF(N) W() Wrk space. At end has (Q**t)*F. * FJAC(N,N) W() Work space for Jacobian. At end has Q of * QR factorization. * R( (N + N**2)/2 ) W() Wrk spc. At end has Packed R of * QR factorization. * IDIMW [in] Dimension of W(). Require IDIMW .ge. 3+(15*N+3*N**2)/2 * ------------------------------------------------------------------ *--S replaces "?": ?NQSOL,?NQSL1,?ERV1,?NQFJ,?NQDOG,?NRM2,?NQFDJ *--& ?NQQFM,?NQQRF,?NQAQ,?NQUPD * Also uses IERM1, IERV1 * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ * */ if (epsmch == 0.0e0) epsmch = FLT_EPSILON; ni = n; Iopt[1] = 1; if (ni <= 0) { ierm1( "SNQSOL", Iopt[1], 0, "Require N > 0", "N", ni, '.' ); goto L_900; } if (idimw < 3 + (ni*(15 + 3*ni))/2) { ierm1( "SNQSOL", Iopt[1], 0, "Require IDIMW .ge. NEED", "IDIMW" , idimw, ',' ); ierv1( "NEED =", 3 + (ni*(15 + 3*ni))/2, '.' ); goto L_900; } /* Set default values. */ havej = TRUE; dmode = 1; nprint = 0; maxfev = MAXF1*(ni + 1); ml = ni - 1; mu = ml; epsfcn = epsmch; factor = FAC1; trace = FALSE; /* Loop on K beginning with K = 4 and * terminating when an option code, J, is zero. */ k = 4; L_20: ; j = Iopt[k]; jabs = labs( j ); jpos = j > 0; switch (jabs + 1) { case 1: goto L_40; case 2: goto L_31; case 3: goto L_32; case 4: goto L_33; case 5: goto L_34; case 6: goto L_35; case 7: goto L_36; case 8: goto L_37; case 9: goto L_38; } /* ANSI Standard Fortran 77 drops thru to here if JABS is * larger than 7. This is an error condition. * */ ierm1( "SNQSOL", Iopt[1], 0, "IOPT(K) must be in [-7..7]", "K" , k, ',' ); ierv1( "IOPT(K)", j, '.' ); goto L_900; L_31: havej = !jpos; k += 1; goto L_20; /* Option 2. Argument = 1, 2, or 3. Default = 1. * 1. This subr sets DIAG() to all ones. * 2. User has set DIAG(). * 3. Subr SNQSL1 sets DIAG() dynamically. */ L_32: if (jpos && Iopt[k + 1] == 2) { dmode = 2; } else if (jpos && Iopt[k + 1] == 3) { dmode = 3; } else if (!jpos || Iopt[k + 1] == 1) { dmode = 1; } else { /* Error. */ ierm1( "SNQSOL", Iopt[1], 0, "Bad argument for Option 2.", "Argument", Iopt[k + 1], '.' ); goto L_900; } k += 2; goto L_20; L_33: if (jpos) { nprint = Iopt[k + 1]; } else { nprint = 0; } k += 2; goto L_20; L_34: if (jpos) { maxfev = Iopt[k + 1]; } else { maxfev = MAXF1*(ni + 1); } k += 2; goto L_20; L_35: if (jpos) { ml = Iopt[k + 1]; mu = Iopt[k + 2]; } else { ml = ni + 1; mu = ml; } k += 3; goto L_20; L_36: if (jpos) { epsfcn = W[1]; } else { epsfcn = epsmch; } k += 1; goto L_20; L_37: if (jpos) { factor = W[2]; } else { factor = FAC1; } k += 1; goto L_20; L_38: if (jpos) { trace = TRUE; } else { trace = FALSE; } k += 1; goto L_20; /* End loop on K */ L_40: ; /* Option 2. DMODE = 1, 2, or 3. * 1. This subr sets DIAG() to all ones. * 2. User has set DIAG(). * 3. Subr SNQSL1 sets DIAG() dynamically. */ if (dmode == 1) { haved = TRUE; for (k = IWDIAG; k <= (IWDIAG + ni - 1); k++) { W[k] = 1.0e0; } } else { haved = dmode == 2; } iwa1 = IWDIAG + ni; iwa2 = iwa1 + ni; iwa3 = iwa2 + ni; iwa4 = iwa3 + ni; iwgnst = iwa4 + ni; iwqtf = iwgnst + ni; iwfjac = iwqtf + ni; iwr = iwfjac + ni*ni; /* IWNEXT = IWR + (N * (N+1)) / 2 Next available loc in W(). * */ snqsl1( snqfj, ni, x, fvec, xtol, &Iopt[1], &Iopt[2], &Iopt[3], nprint, havej, maxfev, haved, ml, mu, epsfcn, factor, trace, &W[IWTOLT], &W[IWDIAG], &W[iwa1], &W[iwa2], &W[iwa3], &W[iwa4], &W[iwgnst], &W[iwqtf], &W[iwfjac], &W[iwr] ); return; /* Error return */ L_900: ; Iopt[2] = 0; Iopt[3] = 0; W[3] = 0.0e0; return; } /* end of function */ /* ================================================================== */ /* PARAMETER translations */ #define COMPUT 1 #define KEEP 3 #define ONE 1.0e0 #define P0001 0.0001e0 #define P1 0.1e0 #define UPDATE 2 #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ snqsl1( void (*snqfj)(long,float[],float[],float*,long*), long n, float x[], float fvec[], float xtol, long *info, long *nfev, long *njev, long nprint, LOGICAL32 havej, long maxfev, LOGICAL32 haved, long ml, long mu, float epsfcn, float factor, LOGICAL32 trace, float *toltst, float diag[], float wa1[], float wa2[], float wa3[], float wa4[], float gnstep[], float qtf[], float *fjac, float r[]) { #define FJAC(I_,J_) (*(fjac+(I_)*(n)+(J_))) LOGICAL32 jeval, newtok, newx, sing, tryzer; long int _l0, i, iflag, iwa[1], j, jact, jact0, l, ldfjac, lr, msum, nbest, ncfail, ncsuc, nextpr, nloop, nslow1, nslow2, numnwt, nupdat; float actred, delta, fnorm, fnorm1, hlim0, hlim1, pnorm, prered, ratio, sum, temp, xnorm; static float epsmch = 0.0e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Diag = &diag[0] - 1; float *const Fvec = &fvec[0] - 1; float *const Gnstep = &gnstep[0] - 1; long *const Iwa = &iwa[0] - 1; float *const Qtf = &qtf[0] - 1; float *const R = &r[0] - 1; float *const Wa1 = &wa1[0] - 1; float *const Wa2 = &wa2[0] - 1; float *const Wa3 = &wa3[0] - 1; float *const Wa4 = &wa4[0] - 1; float *const X = &x[0] - 1; /* end of OFFSET VECTORS */ /*>> 1991-12-04 CLL *>> 1991-12-02 CLL *>> 1991-06-18 CLL@JPL Adapting code from Minpack for MATH77 */ /* 26 arguments. * Dimension of R() must be (N + N**2)/2. * Total space occupied by EPSFCN, FACTOR, and TOLTST through R is * 3 + (15*N + 3*N**2)/2 */ /* ********** * * SUBROUTINE SNQSL1 * * THE PURPOSE OF SNQSL1 IS TO FIND A ZERO OF A SYSTEM OF * N NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION * OF THE POWELL HYBRID METHOD. THE USER MUST PROVIDE A * SUBROUTINE WHICH CALCULATES THE FUNCTIONS and THE JACOBIAN. * * ------------------------------------------------------------------ * Arguments * * SNQFJ is THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH * CALCULATES THE FUNCTIONS and THE JACOBIAN. SNQFJ MUST * BE DECLARED IN AN EXTERNAL STATEMENT IN THE USER * CALLING PROGRAM. SNQFJ will not be called with IFLAG = 2 * if HAVEJ is .false. SNQFJ will not be called with IFLAG = 0 * if NPRINT is <= 0. * SNQFJ is specified as follows: * * subroutine SNQFJ(N, X, FVEC, FJAC, IFLAG) * integer N, IFLAG * real X(N), FVEC(N), FJAC(N,N) * ---------- * if IFLAG = 0, Print X() and FVEC() and return. * IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND * RETURN THIS VECTOR IN FVEC. DO NOT ALTER FJAC. * IF IFLAG = 2 CALCULATE THE JACOBIAN AT X AND * RETURN THIS MATRIX IN FJAC. DO NOT ALTER FVEC. * Set IFLAG to a negative value to force an immediate * termination of the solution procedure. Otherwise do not * alter IFLAG. * --------- * RETURN * END * * N is A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER * OF FUNCTIONS and VARIABLES. * * X is AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN * AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X * CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR. * * FVEC is AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS * THE FUNCTIONS EVALUATED AT THE OUTPUT X. * * XTOL is A NONNEGATIVE INPUT VARIABLE. TERMINATION * OCCURS WHEN THE RELATIVE ERROR BETWEEN TWO CONSECUTIVE * ITERATES is AT MOST XTOL. * * INFO [integer,out] If the user has terminated execution by setting * IFLAG negative in SNQFJ, INFO is set to IFLAG. * Otherwise, INFO is set as follows: * * INFO = 0 Successful termination. Radius of trust region has * been reduced to at most max(XTOL, machine precision). * * INFO = 1 IMPROPER INPUT PARAMETERS. * * INFO = 2 Number of calls to SNQFJ for function evaluations has * reached MAXFEV. * * INFO = 3 XTOL is TOO SMALL. NO FURTHER IMPROVEMENT IN * THE APPROXIMATE SOLUTION X is POSSIBLE. * * INFO = 4 Iteration is not making good progress, as * measured by the improvement through the last * five Jacobian evaluations. * * INFO = 5 Iteration is not making good progress, as * measured by the improvement through last * ten function evaluations. * * NFEV [out,integer] The number of calls to SNQFJ with IFLAG = 1. * * NJEV [out,integer] The number of evaluations of the Jacobian matrix. * If HAVEJ is .true. this will be the number of calls to SNQFJ with * IFLAG = 2. Otherwise it is the number of times the Jacobian has * been approximately computed by differencing. * * NPRINT [in, integer] Enables controlled printing of iterates if it * is positive. In this case, SNQFJ is called with IFLAG = 0 at the * beginning of the first iteration and every NPRINTth time a new X * vector is accepted as an improvement, and at termination. * On these calls the new best X and FVEC are made available for * printing. FVEC and FJAC should not be altered. * If NPRINT is not positive, no special calls to SNQFJ with * IFLAG = 0 will be made. * * HAVEJ [in, logical] True means the user subroutine SNQFJ contains * code for computing the Jacobian matrix, and false means it does * not. * * MAXFEV is A POSITIVE INTEGER INPUT VARIABLE. TERMINATION * OCCURS WHEN THE NUMBER OF CALLS TO SNQFJ WITH IFLAG = 1 * HAS REACHED MAXFEV. * * HAVED = true means initial values of DIAG() are given by the * calling program. False means this subroutine must compute * initial values of DIAG(). It will set DIAG(j) = the euclidean */ /* norm of column j, unless this is zero, in which case it will * set DIAG(j) = 0.0. * * ML and MU specify the band structure, if any, of the Jacobian * matrix. All nonzero elements of the Jacobian matrix lie * within the first ML subdiagonals, the main diagonal, and the * first MU superdiagonals. * ML and MU are only used when HAVEJ is .false. and are only useful * if ML+MU+1 < N. In this case they are used to * reduce the number of function evaluations in estimating * derivatives. If the Jacobian has no band structure set * ML = MU = N-1. * * EPSFCN is an input variable used in determining a suitable * step length for the forward-difference approximation. This * approximation assumes that the relative errors in the * functions are of the order of max(EPSFCN, Machine precision). * * FACTOR is a positive input variable used in determining the * initial step bound. This bound is set to the product of * FACTOR and the euclidean norm of DIAG*X if nonzero, or else * to FACTOR itself. In most cases FACTOR should lie in the * interval (0.1, 10.0). Default: FACTOR = 0.75. * * TRACE [in, logical] If true, this subr will print detailed * intermediate output. Otherwise it will not. * * TOLTST [out] Final value of quantity that is compared with * XTOL for convergence test. * * DIAG is an array of length N. If HAVED = false, * DIAG is internally set. If HAVED = true, DIAG() * MUST CONTAIN POSITIVE ENTRIES THAT SERVE AS * MULTIPLICATIVE SCALE FACTORS FOR THE VARIABLES. * * WA1, WA2, WA3, and WA4 are work arrays of length N. * * GNSTEP() [scratch] Work array of length N to save the * Gauss-Newton step vector computed in SNQDOG. * * QTF is AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS * THE VECTOR (Q TRANSPOSE)*FVEC. * * FJAC is AN OUTPUT N BY N ARRAY WHICH CONTAINS THE * ORTHOGONAL MATRIX Q PRODUCED BY THE QR FACTORIZATION * OF THE FINAL APPROXIMATE JACOBIAN. * * R is AN OUTPUT ARRAY OF LENGTH LR WHICH CONTAINS THE packed * UPPER TRIANGULAR MATRIX PRODUCED BY THE QR FACTORIZATION * OF THE FINAL APPROXIMATE JACOBIAN, STORED ROWWISE. * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * SUBPROGRAMS CALLED * * USER-SUPPLIED ...... SNQFJ * * MINPACK-SUPPLIED ... SNQDOG,R1MACH,SNRM2,SNQFDJ, * SNQQFM,SNQQRF,SNQAQ,SNQUPD * * FORTRAN-SUPPLIED ... abs,max,min,mod * * ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. * BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE' * Argonne Reports: ANL-80-68 and ANL-80-74, 1980. * * 1991-12-09 CLL at JPL. Replacing integer argument MODE that had * values 2 or 1 with logical argument HAVED related to MODE by * HAVED = MODE .eq. 2. Thus the user must set HAVED = .true. when * supplying the DIAG() values, and .false. otherwise. * ********** * ------------------------------------------------------------------ * Description of some of the local variables. * * DELTA [flpt] Diameter of trust region. * HLIM0 [flpt] Upper limit on DELTA when working with a computed * Jacobian. * HLIM1 [flpt] Upper limit on DELTA when working with an updated * Jacobian. * JACT [integer] Can have values of COMPUT, UPDATE, or KEEP. * Initially set to COMPUT. At the beginning of the main loop * we either compute a new Jacobian, update the Jacobian, or keep * the old Jacobian, depending on the setting of JACT. * JACT0 [integer] Saves the value of JACT at the beginning of the * main loop. As JACT gets changed in the loop, JACT0 is still * available as a record of what it was at the beginning of the loop. * JEVAL [logical] Set true whenever the Jacobian is computed, and set * false when it is updated. * NBEST [integer] Counter, incremented each time an x-vector is * accepted as being a better approximation to the solution. Used * in connection with NPRINT to trigger calles to SNQFJ for printing. * NCFAIL [integer] Counts consecutive "failed" steps since the last * computation of the Jacobian. NCFAIL is set to 0 when the Jacobian * is computed or when a step "succeeds" in the sense that * RATIO .ge. 0.1. It is incremented when RATIO .lt. 0.1. * NLOOP [integer] Counter for main iteration loop. * NUPDAT [integer] Counts number of consecutive times the Jacobian * matrix is updated. * TRYZER [logical] Initially set to true. While true, the algorithm * will monitor X's to see if they seem to be all approaching zero. * If so will try setting them all to zero. If this gives an exactly * zero function vector then we are finished. If not, we set TRYZER */ /* to false and restore X to its previous value (even if the function * value at X = 0 was an improvement) and omit any further testing * for X's approaching zero. (We tryed accepting the X reached by * this exceptional step if the function value was an improvement, * but in one test case this caused the algorithm to end at a local * nonzero minimum rather than finding a zero.) * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ * Set EPSMCH to the machine precision. * */ if (epsmch == 0.0e0) epsmch = FLT_EPSILON; /* Initialize values of output arguments. * */ *info = 1; *nfev = 0; *njev = 0; *toltst = 0.0e0; tryzer = TRUE; /* CHECK THE INPUT PARAMETERS FOR ERRORS. * We assume the condition N > 0 has already been checked in * the user-interface subroutine that called this one. */ if ((xtol < ZERO || maxfev <= 0) || factor <= ZERO) { ierm1( "SNQSL1", *info, 0, "Require MAXFEV > 0, XTOL .gt. 0.0, FACTOR > 0.0" , "MAXFEV", maxfev, ',' ); serv1( "XTOL", xtol, ',' ); serv1( "FACTOR", factor, '.' ); goto L_300; } if (!havej && (ml < 0 || mu < 0)) { ierm1( "SNQSL1", *info, 0, "With HAVEJ false, require ML .ge. 0 and MU .ge. 0" , "ML", ml, ',' ); ierv1( "MU", mu, ',' ); goto L_300; } /* HAVED = true means the user has set DIAG(). */ if (haved) { for (j = 1; j <= n; j++) { if (Diag[j] <= ZERO) { ierm1( "SNQSL1", *info, 0, "With HAVED = .true., require all DIAG(J) > 0.0" , "J", j, ',' ); serv1( "DIAG(J)", Diag[j], '.' ); goto L_300; } } } /* Initialize algorithm variables. */ *info = 0; jact = COMPUT; ldfjac = n; lr = (n*(n + 1))/2; msum = min( ml + mu + 1, n ); nbest = 1; ncsuc = 0; nextpr = 1; nloop = 0; nslow1 = 0; nslow2 = 0; numnwt = 0; /* Evaluate the function at the starting point. * Calculate and test its norm. * */ iflag = 1; (*snqfj)( n, x, fvec, fjac, &iflag ); *nfev = 1; /* CALL SNQFJ(N, X, FVEC, FJAC, IFLAG) */ if (iflag < 0) goto L_300; fnorm = snrm2( n, fvec, 1 ); if (trace) { printf(" %5ld Initial X:", nloop); printf("\n "); for (j = 1; j <= n; j++) { printf("%15.6g", X[j]); } printf("\n"); printf(" Initial FNORM:%15.6g\n", fnorm); } if (fnorm == 0.0e0) { goto L_300; } /* Beginning of main loop. * */ L_30: ; nloop += 1; jact0 = jact; /* Compute, Update, or Keep Jacobian, depending on JACT. * */ if (jact == COMPUT) { jeval = TRUE; nupdat = 0; ncfail = 0; /* CALCULATE THE JACOBIAN MATRIX. * */ if (trace) { printf(" %5ld Computing new Jacobian matrix.\n", nloop); } *njev += 1; if (havej) { iflag = 2; (*snqfj)( n, x, fvec, fjac, &iflag ); } else { /* CALL SNQFJ(N, X, FVEC, FJAC, IFLAG) */ snqfdj( snqfj, n, x, fvec, fjac, ldfjac, &iflag, ml, mu, epsfcn, wa1, wa2 ); *nfev += msum; } if (iflag < 0) goto L_300; /* COMPUTE THE QR FACTORIZATION OF THE JACOBIAN. * */ snqqrf( n, n, fjac, ldfjac, FALSE, iwa, 1, wa1, wa2, wa3 ); /* On the first iteration and if HAVED is .false., scale according * to the norms of the columns of the initial Jacobian. * Also on the first iteration calculate the norm of the scaled X * and initialize the trust region diameter, DELTA. * */ if (nloop == 1) { if (!haved) { for (j = 1; j <= n; j++) { Diag[j] = Wa2[j]; if (Wa2[j] == ZERO) Diag[j] = ONE; } } for (j = 1; j <= n; j++) { Wa3[j] = Diag[j]*X[j]; } xnorm = snrm2( n, wa3, 1 ); delta = factor*xnorm; if (delta == ZERO) delta = factor; } /* FORM (Q TRANSPOSE)*FVEC and STORE IN QTF. * */ for (i = 1; i <= n; i++) { Qtf[i] = Fvec[i]; } for (j = 1; j <= n; j++) { if (FJAC(j - 1,j - 1) == ZERO) goto L_110; sum = ZERO; for (i = j; i <= n; i++) { sum += FJAC(j - 1,i - 1)*Qtf[i]; } temp = -sum/FJAC(j - 1,j - 1); for (i = j; i <= n; i++) { Qtf[i] += FJAC(j - 1,i - 1)*temp; } L_110: ; } /* COPY THE TRIANGULAR FACTOR OF THE QR FACTORIZATION INTO R. * The diagonal elts come from WA1(). The strictly upper * triangular elts come from FJAC(,). The upper triangular matrix * will be stored, packed by rows, in R(). * */ sing = FALSE; for (j = 1; j <= n; j++) { l = j; for (i = 1; i <= (j - 1); i++) { R[l] = FJAC(j - 1,i - 1); l += n - i; } R[l] = Wa1[j]; if (Wa1[j] == ZERO) sing = TRUE; } /* ACCUMULATE THE ORTHOGONAL FACTOR IN FJAC. * */ snqqfm( n, n, fjac, ldfjac, wa1 ); /* RESCALE IF NECESSARY. * */ if (!haved) { for (j = 1; j <= n; j++) { Diag[j] = fmaxf( Diag[j], Wa2[j] ); } } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ } else if (jact == UPDATE) { /* CALCULATE THE RANK ONE MODIFICATION TO THE JACOBIAN * and UPDATE QTF IF NECESSARY. * */ if (trace) { printf(" %5ld Updating Jacobian matrix.\n", nloop); } nupdat += 1; jeval = FALSE; for (j = 1; j <= n; j++) { sum = ZERO; for (i = 1; i <= n; i++) { sum += FJAC(j - 1,i - 1)*Wa4[i]; } Wa2[j] = (sum - Wa3[j])/pnorm; Wa1[j] = Diag[j]*((Diag[j]*Wa1[j])/pnorm); if (ratio >= P0001) Qtf[j] = sum; } /* COMPUTE THE QR FACTORIZATION OF THE UPDATED JACOBIAN. * */ snqupd( n, n, r, lr, wa1, wa2, wa3, &sing ); snqaq( n, n, fjac, ldfjac, wa2, wa3 ); snqaq( 1, n, qtf, 1, wa2, wa3 ); } else { if (trace) { printf(" %5ld Keeping Jacobian matrix unchanged.\n", nloop); } } /* Now have a new or updated or retained Jacobian matrix. * * IF REQUESTED, CALL SNQFJ TO ENABLE PRINTING OF ITERATES. * */ if (nprint > 0) { if (nbest == nextpr) { iflag = 0; (*snqfj)( n, x, fvec, fjac, &iflag ); if (iflag < 0) goto L_300; /* CALL SNQFJ(N, X, FVEC, FJAC, IFLAG) */ nextpr += nprint; } } /* Determine the direction P, using a dogleg method, and * returning -P in WA1(). * */ snqdog( n, r, lr, diag, qtf, delta, wa1, &newtok, wa2, wa3, jact0 == KEEP, gnstep ); /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ if (tryzer) { /* NUMNWT counts number of consecutive * full Newton steps. */ if (newtok) { numnwt += 1; } else { numnwt = 0; } /* Test for convergence of some x components toward 0. * If this seems to be happening try setting such * components to 0. * */ if (numnwt >= 5 && ncsuc >= 4) { numnwt = 0; for (j = 1; j <= n; j++) { Wa2[j] = X[j] - Wa1[j]; if (fabsf( Wa2[j] ) <= 0.75e0*fabsf( X[j] )) { Wa2[j] = 0.0e0; } else { goto L_203; } } if (trace) { printf(" %5ld Trial setting of X() to zero.\n", nloop); } /* EVALUATE THE FUNCTION AT WA2() and CALCULATE ITS NORM. * */ iflag = 1; (*snqfj)( n, wa2, wa4, fjac, &iflag ); *nfev += 1; /* CALL SNQFJ(N, WA2, WA4, FJAC, IFLAG) */ if (iflag < 0) goto L_300; fnorm1 = snrm2( n, wa4, 1 ); if (trace) { printf(" %5ld FNORM1 = %15.6g\n", nloop, fnorm1); } if (fnorm1 == 0.0e0) { /* Accept new point as final solution. * Update X() and FVEC() and go to termination. * */ *info = 0; *toltst = 0.0e0; for (j = 1; j <= n; j++) { X[j] = Wa2[j]; Fvec[j] = Wa4[j]; } if (trace) { printf(" %5ld Accepting X = all zeros. \n", nloop); } goto L_300; } else { tryzer = FALSE; } } /* The following "endif" matches "if(TRYZER)then" */ } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * * STORE THE DIRECTION P and X + P. CALCULATE THE NORM OF P. * */ L_203: ; for (j = 1; j <= n; j++) { Wa1[j] = -Wa1[j]; Wa2[j] = X[j] + Wa1[j]; Wa3[j] = Diag[j]*Wa1[j]; } pnorm = snrm2( n, wa3, 1 ); if (trace) { printf(" %5ld DELTA PNORM\n %15.6g%15.6g\n", nloop, delta, pnorm); printf(" Trial X:"); printf("\n "); for (j = 1; j <= n; j++) { printf("%15.6g", Wa2[j]); } printf("\n"); } /* ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND. * */ if (nloop == 1) { delta = fminf( delta, pnorm ); hlim0 = delta; hlim1 = delta; } /* EVALUATE THE FUNCTION AT X + P and CALCULATE ITS NORM. * */ iflag = 1; (*snqfj)( n, wa2, wa4, fjac, &iflag ); *nfev += 1; /* CALL SNQFJ(N, WA2, WA4, FJAC, IFLAG) */ if (iflag < 0) goto L_300; fnorm1 = snrm2( n, wa4, 1 ); /* COMPUTE THE SCALED ACTUAL REDUCTION. * */ actred = -ONE; if (fnorm1 < fnorm) actred = ONE - powif(fnorm1/fnorm,2); /* COMPUTE THE SCALED PREDICTED REDUCTION. * */ l = 1; for (i = 1; i <= n; i++) { sum = ZERO; for (j = i; j <= n; j++) { sum += R[l]*Wa1[j]; l += 1; } Wa3[i] = Qtf[i] + sum; } temp = snrm2( n, wa3, 1 ); prered = ZERO; if (temp < fnorm) prered = ONE - powif(temp/fnorm,2); /* COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED * REDUCTION. * */ ratio = ZERO; if (prered > ZERO) ratio = actred/prered; if (trace) { printf(" %5ld FNORM1 ACTRED PRERED RATIO\n %15.6g%15.6g%15.6g%15.6g\n", nloop, fnorm1, actred, prered, ratio); } /* Analyze RATIO, NCSUC and JEVAL to decide on accepting or * rejecting the new X, and assigning new values to * NCSUC, JACT, and DELTA. * */ if (ratio < 0.0000e0) { ncsuc = 0; ncfail += 1; newx = FALSE; if (jeval) hlim0 = fminf( hlim0, 0.707107e0*pnorm ); hlim1 = fminf( hlim1, 0.707107e0*pnorm ); if (jeval || (ncfail <= 1 && nupdat <= 2)) { jact = KEEP; delta = 0.5e0*pnorm; } else { jact = COMPUT; delta = hlim0; } } else if (ratio < 0.1e0) { ncsuc = 0; ncfail += 1; newx = TRUE; if (ncfail <= 1 && nupdat <= 2) { jact = UPDATE; delta = 0.5e0*pnorm; } else { jact = COMPUT; delta = hlim0; } } else { /* Here we have RATIO .ge. 0.1 */ ncsuc += 1; ncfail = 0; newx = TRUE; jact = UPDATE; if (ratio < 0.5e0) { if (ncsuc >= 5) hlim1 = fmaxf( hlim1, 1.414214e0*pnorm ); if (ncsuc >= 2) delta = fminf( hlim1, fmaxf( delta, 1.414214e0*pnorm ) ); } else if (ratio < 0.9e0) { if (jact0 == COMPUT) hlim0 = fmaxf( hlim0, 1.414214e0*pnorm ); if (ncsuc >= 4) hlim1 = fmaxf( hlim1, 1.414214e0*pnorm ); if (ncsuc >= 2) delta = fminf( hlim1, fmaxf( delta, 1.414214e0*pnorm ) ); } else if (ratio < 1.1e0) { if (jact0 == COMPUT) hlim0 = fmaxf( hlim0, 2.0e0*pnorm ); if (ncsuc == 1) { delta = 1.414214e0*pnorm; } else { delta = 2.0e0*pnorm; } hlim1 = fmaxf( hlim1, delta ); } } hlim0 = fmaxf( hlim0, hlim1 ); if (trace) { printf(" %5ld NCSUC NCFAIL NUPDAT DELTA HLIM0 HLIM1\n %8ld%8ld%8ld%13.4g%13.4g%13.4g\n", nloop, ncsuc, ncfail, nupdat, delta, hlim0, hlim1); } if (newx) { /* Accept new X, FVEC, and their norms. */ for (j = 1; j <= n; j++) { X[j] = Wa2[j]; Wa2[j] = Diag[j]*X[j]; Fvec[j] = Wa4[j]; } xnorm = snrm2( n, wa2, 1 ); fnorm = fnorm1; nbest += 1; if (trace) { printf(" %5ld Accepting new X with XNORM = %15.6g\n", nloop, xnorm); } } /* DETERMINE THE PROGRESS OF THE ITERATION. * */ if (actred >= 0.001e0) { nslow1 = 0; } else { nslow1 += 1; } if (actred >= 0.1e0) { nslow2 = 0; } else if (jact0 == COMPUT) { nslow2 += 1; } if (trace) { printf(" %5ld NSLOW1 NSLOW2\n %11ld %11ld \n", nloop, nslow1, nslow2); } /* TEST FOR CONVERGENCE. * */ if (delta <= xtol*xnorm || fnorm == ZERO) { *info = 0; if (trace) { printf(" %5ld INFO XNORM\n %14ld%15.6g\n", nloop, *info, xnorm); } goto L_295; } /* TESTS FOR TERMINATION and STRINGENT TOLERANCES. * */ if (*nfev >= maxfev) *info = 2; if (P1*fmaxf( P1*delta, pnorm ) <= epsmch*xnorm) *info = 3; if (nslow2 == 5) *info = 4; if (nslow1 == 10) *info = 5; if (*info != 0) { if (trace) { printf(" %5ld INFO XNORM\n %14ld%15.6g\n", nloop, *info, xnorm); } ierm1( "SNQSL1", *info, 0, "Unsuccessful termination.", "INFO" , *info, '.' ); goto L_295; } goto L_30; /* End of main loop. * * Come to following stmt when INFO has been set to * 2, 3, 4, or 5, or to 0 due to successful XTOL test. */ L_295: ; if (xnorm != 0.0e0) { *toltst = delta/xnorm; } else { *toltst = delta; } /* Jump to following statement with IFLAG negative * or INFO = 1 or INFO = 0 due to FNORM being zero. * Here we have TOLTST = 0.0. */ L_300: ; /* TERMINATION, EITHER NORMAL OR USER IMPOSED. * */ if (iflag < 0) *info = iflag; if (trace) { printf(" %5ld Quitting with INFO = %3ld\n", nloop, *info); } iflag = 0; if (nprint > 0) (*snqfj)( n, x, fvec, fjac, &iflag ); if (*info < 0) { /* IF (NPRINT .gt. 0) CALL SNQFJ(N,X,FVEC,FJAC, IFLAG) */ ierm1( "SNQSL1", *info, 0, "Quitting because user code set IFLAG negative." , "IFLAG", *info, '.' ); } return; /* Last line of subroutine SNQSL1. * */ #undef FJAC } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ snqfdj( void (*snqfj)(long,float[],float[],float*,long*), long n, float x[], float fvec[], float *fjac, long ldfjac, long *iflag, long ml, long mu, float epsfcn, float wa1[], float wa2[]) { #define FJAC(I_,J_) (*(fjac+(I_)*(ldfjac)+(J_))) long int _d_l, _d_m, _do0, _do1, _do2, _do3, _l0, i, j, k, msum; float eps, epsmch, h, temp; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Fvec = &fvec[0] - 1; float *const Wa1 = &wa1[0] - 1; float *const Wa2 = &wa2[0] - 1; float *const X = &x[0] - 1; /* end of OFFSET VECTORS */ /*>> 1991-12-04 CLL Changed arg list of user supplied subroutine. *>> 1991-06-18 CLL@JPL Adapting code from Minpack for MATH77 */ /* ********** * * SUBROUTINE SNQFDJ * * THIS SUBROUTINE COMPUTES A FORWARD-DIFFERENCE APPROXIMATION * TO THE N BY N JACOBIAN MATRIX ASSOCIATED WITH A SPECIFIED * PROBLEM OF N FUNCTIONS IN N VARIABLES. IF THE JACOBIAN HAS * A BANDED FORM, THEN FUNCTION EVALUATIONS ARE SAVED BY ONLY * APPROXIMATING THE NONZERO TERMS. * * THE SUBROUTINE STATEMENT IS * * subroutine SNQFDJ(SNQFJ,N,X,FVEC,FJAC,LDFJAC,IFLAG,ML,MU,EPSFCN, * WA1,WA2) * * WHERE * * SNQFJ IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH * CALCULATES THE FUNCTIONS. SNQFJ MUST BE DECLARED * IN AN EXTERNAL STATEMENT IN THE USER CALLING * PROGRAM, and SHOULD BE WRITTEN AS FOLLOWS. * * subroutine SNQFJ(N,X,FVEC,IFLAG) * integer N,IFLAG * real X(N),FVEC(N) * ---------- * CALCULATE THE FUNCTIONS AT X AND * RETURN THIS VECTOR IN FVEC. * ---------- * RETURN * END * * THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY SNQFJ UNLESS * THE USER WANTS TO TERMINATE EXECUTION OF SNQFDJ. * IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER. * * N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER * OF FUNCTIONS and VARIABLES. * * X IS AN INPUT ARRAY OF LENGTH N. * * FVEC IS AN INPUT ARRAY OF LENGTH N WHICH MUST CONTAIN THE * FUNCTIONS EVALUATED AT X. * * FJAC IS AN OUTPUT N BY N ARRAY WHICH CONTAINS THE * APPROXIMATION TO THE JACOBIAN MATRIX EVALUATED AT X. * * LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N * WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC. * * IFLAG IS AN INTEGER VARIABLE WHICH CAN BE USED TO TERMINATE * THE EXECUTION OF SNQFDJ. SEE DESCRIPTION OF SNQFJ. * * ML IS A NONNEGATIVE INTEGER INPUT VARIABLE WHICH SPECIFIES * THE NUMBER OF SUBDIAGONALS WITHIN THE BAND OF THE * JACOBIAN MATRIX. IF THE JACOBIAN IS NOT BANDED, SET * ML TO AT LEAST N - 1. * * MU IS A NONNEGATIVE INTEGER INPUT VARIABLE WHICH SPECIFIES * THE NUMBER OF SUPERDIAGONALS WITHIN THE BAND OF THE * JACOBIAN MATRIX. IF THE JACOBIAN IS NOT BANDED, SET * MU TO AT LEAST N - 1. * * EPSFCN IS AN INPUT VARIABLE USED IN DETERMINING A SUITABLE * STEP LENGTH FOR THE FORWARD-DIFFERENCE APPROXIMATION. THIS * APPROXIMATION ASSUMES THAT THE RELATIVE ERRORS IN THE * FUNCTIONS ARE OF THE ORDER OF EPSFCN. IF EPSFCN IS LESS * THAN THE MACHINE PRECISION, IT IS ASSUMED THAT THE RELATIVE * ERRORS IN THE FUNCTIONS ARE OF THE ORDER OF THE MACHINE * PRECISION. * * WA1 and WA2 ARE WORK ARRAYS OF LENGTH N. IF ML + MU + 1 IS AT * LEAST N, THEN THE JACOBIAN IS CONSIDERED DENSE, and WA2 IS * NOT REFERENCED. * * SUBPROGRAMS CALLED * * MINPACK-SUPPLIED ... R1MACH * * FORTRAN-SUPPLIED ... abs,max,sqrt * * ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. * BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE * * ********** * ------------------------------------------------------------------ */ /*++ CODE for ~.C. is inactive * real DUMMY(1,1) *++ CODE for .C. & (.N. == 'S') is active */ float *dummy; /*++ CODE for .C. & (.N. == 'D') is inactive *%% double *dummy; *++ End */ /* EPSMCH IS THE MACHINE PRECISION. * */ epsmch = FLT_EPSILON; eps = sqrtf( fmaxf( epsfcn, epsmch ) ); *iflag = 1; msum = ml + mu + 1; if (msum < n) goto L_40; /* COMPUTATION OF DENSE APPROXIMATE JACOBIAN. * */ for (j = 1; j <= n; j++) { temp = X[j]; h = eps*fabsf( temp ); if (h == ZERO) h = eps; X[j] = temp + h; (*snqfj)( n, x, wa1, dummy, iflag ); if (*iflag < 0) goto L_30; /* CALL SNQFJ(N, X, WA1, DUMMY, IFLAG) */ X[j] = temp; for (i = 1; i <= n; i++) { FJAC(j - 1,i - 1) = (Wa1[i] - Fvec[i])/h; } } L_30: ; goto L_110; L_40: ; /* COMPUTATION OF BANDED APPROXIMATE JACOBIAN. * */ for (k = 1; k <= msum; k++) { for (j = k, _do0=DOCNT(j,n,_do1 = msum); _do0 > 0; j += _do1, _do0--) { Wa2[j] = X[j]; h = eps*fabsf( Wa2[j] ); if (h == ZERO) h = eps; X[j] = Wa2[j] + h; } (*snqfj)( n, x, wa1, dummy, iflag ); if (*iflag < 0) goto L_100; /* CALL SNQFJ(N, X, WA1, DUMMY, IFLAG) */ for (j = k, _do2=DOCNT(j,n,_do3 = msum); _do2 > 0; j += _do3, _do2--) { X[j] = Wa2[j]; h = eps*fabsf( Wa2[j] ); if (h == ZERO) h = eps; for (i = 1; i <= n; i++) { FJAC(j - 1,i - 1) = ZERO; if (i >= j - mu && i <= j + ml) FJAC(j - 1,i - 1) = (Wa1[i] - Fvec[i])/h; } } } L_100: ; L_110: ; return; /* Last line of subroutine SNQFDJ. * */ #undef FJAC } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ snqaq( long m, long n, float *a, long lda, float v[], float w[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) long int i, j, nm1, nmj; float temp, vcos, vsin; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const V = &v[0] - 1; float *const W = &w[0] - 1; /* end of OFFSET VECTORS */ /* ********** * * SUBROUTINE SNQAQ * * GIVEN AN M BY N MATRIX A, THIS SUBROUTINE COMPUTES A*Q WHERE * Q IS THE PRODUCT OF 2*(N - 1) TRANSFORMATIONS * * GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1) * * and GV(I), GW(I) ARE GIVENS ROTATIONS IN THE (I,N) PLANE WHICH * ELIMINATE ELEMENTS IN THE I-TH and N-TH PLANES, RESPECTIVELY. * Q ITSELF IS NOT GIVEN, RATHER THE INFORMATION TO RECOVER THE * GV, GW ROTATIONS IS SUPPLIED. * * THE SUBROUTINE STATEMENT IS * * subroutine SNQAQ(M,N,A,LDA,V,W) * * WHERE * * M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER * OF ROWS OF A. * * N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER * OF COLUMNS OF A. * * A IS AN M BY N ARRAY. ON INPUT A MUST CONTAIN THE MATRIX * TO BE POSTMULTIPLIED BY THE ORTHOGONAL MATRIX Q * DESCRIBED ABOVE. ON OUTPUT A*Q HAS REPLACED A. * * LDA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M * WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY A. * * V IS AN INPUT ARRAY OF LENGTH N. V(I) MUST CONTAIN THE * INFORMATION NECESSARY TO RECOVER THE GIVENS ROTATION GV(I) * DESCRIBED ABOVE. * * W IS AN INPUT ARRAY OF LENGTH N. W(I) MUST CONTAIN THE * INFORMATION NECESSARY TO RECOVER THE GIVENS ROTATION GW(I) * DESCRIBED ABOVE. * * SUBROUTINES CALLED * * FORTRAN-SUPPLIED ... abs,sqrt * * ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. * BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE * * ********** * ------------------------------------------------------------------ */ /* APPLY THE FIRST SET OF GIVENS ROTATIONS TO A. * */ nm1 = n - 1; if (nm1 < 1) goto L_50; for (nmj = 1; nmj <= nm1; nmj++) { j = n - nmj; if (fabsf( V[j] ) > ONE) vcos = ONE/V[j]; if (fabsf( V[j] ) > ONE) vsin = sqrtf( ONE - SQ(vcos) ); if (fabsf( V[j] ) <= ONE) vsin = V[j]; if (fabsf( V[j] ) <= ONE) vcos = sqrtf( ONE - SQ(vsin) ); for (i = 1; i <= m; i++) { temp = vcos*A(j - 1,i - 1) - vsin*A(n - 1,i - 1); A(n - 1,i - 1) = vsin*A(j - 1,i - 1) + vcos*A(n - 1,i - 1); A(j - 1,i - 1) = temp; } } /* APPLY THE SECOND SET OF GIVENS ROTATIONS TO A. * */ for (j = 1; j <= nm1; j++) { if (fabsf( W[j] ) > ONE) vcos = ONE/W[j]; if (fabsf( W[j] ) > ONE) vsin = sqrtf( ONE - SQ(vcos) ); if (fabsf( W[j] ) <= ONE) vsin = W[j]; if (fabsf( W[j] ) <= ONE) vcos = sqrtf( ONE - SQ(vsin) ); for (i = 1; i <= m; i++) { temp = vcos*A(j - 1,i - 1) + vsin*A(n - 1,i - 1); A(n - 1,i - 1) = -vsin*A(j - 1,i - 1) + vcos*A(n - 1,i - 1); A(j - 1,i - 1) = temp; } } L_50: ; return; /* Last line of subroutine SNQAQ. * */ #undef A } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ snqdog( long n, float r[], long lr, float diag[], float qtb[], float delta, float x[], LOGICAL32 *newtok, float wa1[], float wa2[], LOGICAL32 samej, float gnstep[]) { long int _l0, i, j, jj, jp1, k, l, np1; float alpha, bnorm, gnorm, qnorm, sgnorm, sum, temp; static float epsmch = 0.0e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Diag = &diag[0] - 1; float *const Gnstep = &gnstep[0] - 1; float *const Qtb = &qtb[0] - 1; float *const R = &r[0] - 1; float *const Wa1 = &wa1[0] - 1; float *const Wa2 = &wa2[0] - 1; float *const X = &x[0] - 1; /* end of OFFSET VECTORS */ /*>> 1992-01-03 CLL */ /* ********** * * subroutine SNQDOG * * Given an M by N matrix A, an N by N nonsingular diagonal * matrix D, an M-vector B, and a positive number DELTA, the * problem is to determine the convex combination X of the * gauss-newton and scaled gradient directions that minimizes * (A*X - B) in the least squares sense, subject to the * restriction that the euclidean norm of D*X be at most DELTA. * * This subroutine completes the solution of the problem * if it is provided with the necessary information from the * QR factorization of A. that is, if A = Q*R, where Q has * orthogonal columns and R is an upper triangular matrix, * then SNQDOG needs the full upper triangle of R and * the first N components of (Q transpose)*B. * * The subroutine statement is * * subroutine SNQDOG(N,R,LR,DIAG,QTB,DELTA,X,WA1,WA2) * * where * * N is A POSITIVE INTEGER INPUT VARIABLE SET TO THE ORDER OF R. * * R() [in] An ARRAY OF LENGTH LR WHICH MUST CONTAIN THE UPPER * TRIANGULAR MATRIX R STORED BY ROWS. * * LR is A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN * (N*(N+1))/2. * * DIAG() [in] An ARRAY OF LENGTH N WHICH MUST CONTAIN THE * DIAGONAL ELEMENTS OF THE MATRIX D. * * QTB() [in] An ARRAY OF LENGTH N WHICH MUST CONTAIN THE FIRST * N ELEMENTS OF THE VECTOR (Q TRANSPOSE)*B. * * DELTA is a POSITIVE INPUT VARIABLE WHICH SPECIFIES AN UPPER * BOUND ON THE EUCLIDEAN NORM OF D*X. * * X() [out] An ARRAY OF LENGTH N WHICH CONTAINS THE DESIRED * CONVEX COMBINATION OF THE GAUSS-NEWTON DIRECTION and THE * SCALED GRADIENT DIRECTION. * * NEWTOK [logical, out] True means the full Newton step was * used. False means a modified, shorter, step was used. * * WA1() and WA2() are work arrays of length N. * * SAMEJ [logical, in] True means we have the same Jacobian matrix as * on the previous call to this subr. The Gauss-Newton vector in * GNSTEP() can be reused. * * GNSTEP() [inout] On return holds the Gauss-Newton vector. On entry * with SAMEJ = .true., contains the GN vector from the previous call. * ------------------------------------------------------------------ * SUBPROGRAMS CALLED * * MINPACK-SUPPLIED ... R1MACH,SNRM2 * * FORTRAN-SUPPLIED ... abs,max,min,sqrt * * ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. * BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE * * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ * Set EPSMCH to the machine precision. * */ if (epsmch == 0.0e0) epsmch = FLT_EPSILON; if (!samej) { /* FIRST, CALCULATE THE GAUSS-NEWTON DIRECTION. * */ np1 = n + 1; jj = (n*(n + 1))/2 + 1; for (k = 1; k <= n; k++) { j = np1 - k; jp1 = j + 1; jj -= k; l = jj + 1; sum = ZERO; for (i = jp1; i <= n; i++) { sum += R[l]*Gnstep[i]; l += 1; } temp = R[jj]; if (temp == ZERO) { l = j; for (i = 1; i <= (j - 1); i++) { temp = fmaxf( temp, fabsf( R[l] ) ); l += n - i; } temp *= epsmch; } if (temp == ZERO) { Gnstep[j] = 0.0e0; } else { Gnstep[j] = (Qtb[j] - sum)/temp; } } } /* TEST WHETHER THE GAUSS-NEWTON DIRECTION is ACCEPTABLE. * */ for (j = 1; j <= n; j++) { Wa1[j] = ZERO; Wa2[j] = Diag[j]*Gnstep[j]; } qnorm = snrm2( n, wa2, 1 ); *newtok = qnorm <= delta; if (*newtok) { for (j = 1; j <= n; j++) { X[j] = Gnstep[j]; } goto L_140; } /* THE GAUSS-NEWTON DIRECTION is NOT ACCEPTABLE. * NEXT, CALCULATE THE SCALED GRADIENT DIRECTION. * */ l = 1; for (j = 1; j <= n; j++) { temp = Qtb[j]; for (i = j; i <= n; i++) { Wa1[i] += R[l]*temp; l += 1; } Wa1[j] /= Diag[j]; } /* CALCULATE THE NORM OF THE SCALED GRADIENT and TEST FOR * THE SPECIAL CASE IN WHICH THE SCALED GRADIENT is ZERO. * */ gnorm = snrm2( n, wa1, 1 ); if (gnorm == ZERO) { alpha = delta/qnorm; for (j = 1; j <= n; j++) { X[j] = alpha*Gnstep[j]; } goto L_140; } /* CALCULATE THE POINT ALONG THE SCALED GRADIENT * AT WHICH THE QUADRATIC is MINIMIZED. * */ for (j = 1; j <= n; j++) { Wa1[j] = (Wa1[j]/gnorm)/Diag[j]; } l = 1; for (j = 1; j <= n; j++) { sum = ZERO; for (i = j; i <= n; i++) { sum += R[l]*Wa1[i]; l += 1; } Wa2[j] = sum; } temp = snrm2( n, wa2, 1 ); sgnorm = (gnorm/temp)/temp; /* TEST WHETHER THE SCALED GRADIENT DIRECTION is ACCEPTABLE. * */ alpha = ZERO; if (sgnorm < delta) { /* THE SCALED GRADIENT DIRECTION is NOT ACCEPTABLE. * FINALLY, CALCULATE THE POINT ALONG THE dogleg * AT WHICH THE QUADRATIC is MINIMIZED. * */ bnorm = snrm2( n, qtb, 1 ); temp = (bnorm/gnorm)*(bnorm/qnorm)*(sgnorm/delta); temp += -(delta/qnorm)*powif(sgnorm/delta,2) + sqrtf( powif(temp - (delta/qnorm),2) + (ONE - powif(delta/qnorm,2))*(ONE - powif(sgnorm/ delta,2)) ); alpha = ((delta/qnorm)*(ONE - powif(sgnorm/delta,2)))/temp; } /* FORM APPROPRIATE CONVEX COMBINATION OF THE GAUSS-NEWTON * DIRECTION and THE SCALED GRADIENT DIRECTION. * */ temp = (ONE - alpha)*fminf( sgnorm, delta ); for (j = 1; j <= n; j++) { X[j] = temp*Wa1[j] + alpha*Gnstep[j]; } L_140: ; return; /* Last line of subroutine SNQDOG. * */ } /* end of function */ void /*FUNCTION*/ snqqfm( long m, long n, float *q, long ldq, float wa[]) { #define Q(I_,J_) (*(q+(I_)*(ldq)+(J_))) long int i, j, jm1, k, l, minmn, np1; float sum, temp; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Wa = &wa[0] - 1; /* end of OFFSET VECTORS */ /* ********** * * SUBROUTINE SNQQFM * * THIS SUBROUTINE PROCEEDS FROM THE COMPUTED QR FACTORIZATION OF * AN M BY N MATRIX A TO ACCUMULATE THE M BY M ORTHOGONAL MATRIX * Q FROM ITS FACTORED FORM. * * THE SUBROUTINE STATEMENT IS * * subroutine SNQQFM(M,N,Q,LDQ,WA) * * WHERE * * M is A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER * OF ROWS OF A and THE ORDER OF Q. * * N is A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER * OF COLUMNS OF A. * * Q is AN M BY M ARRAY. ON INPUT THE FULL LOWER TRAPEZOID IN * THE FIRST MIN(M,N) COLUMNS OF Q CONTAINS THE FACTORED FORM. * ON OUTPUT Q HAS BEEN ACCUMULATED INTO A SQUARE MATRIX. * * LDQ is A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M * WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY Q. * * WA is A WORK ARRAY OF LENGTH M. * * SUBPROGRAMS CALLED * * FORTRAN-SUPPLIED ... min * * ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. * BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE * * ------------------------------------------------------------------ */ /* ZERO OUT UPPER TRIANGLE OF Q IN THE FIRST MIN(M,N) COLUMNS. * */ minmn = min( m, n ); if (minmn < 2) goto L_30; for (j = 2; j <= minmn; j++) { jm1 = j - 1; for (i = 1; i <= jm1; i++) { Q(j - 1,i - 1) = ZERO; } } L_30: ; /* INITIALIZE REMAINING COLUMNS TO THOSE OF THE IDENTITY MATRIX. * */ np1 = n + 1; if (m < np1) goto L_60; for (j = np1; j <= m; j++) { for (i = 1; i <= m; i++) { Q(j - 1,i - 1) = ZERO; } Q(j - 1,j - 1) = ONE; } L_60: ; /* ACCUMULATE Q FROM ITS FACTORED FORM. * */ for (l = 1; l <= minmn; l++) { k = minmn - l + 1; for (i = k; i <= m; i++) { Wa[i] = Q(k - 1,i - 1); Q(k - 1,i - 1) = ZERO; } Q(k - 1,k - 1) = ONE; if (Wa[k] == ZERO) goto L_110; for (j = k; j <= m; j++) { sum = ZERO; for (i = k; i <= m; i++) { sum += Q(j - 1,i - 1)*Wa[i]; } temp = sum/Wa[k]; for (i = k; i <= m; i++) { Q(j - 1,i - 1) += -temp*Wa[i]; } } L_110: ; } return; /* Last line of subroutine SNQQFM. * */ #undef Q } /* end of function */ /* PARAMETER translations */ #define P05 0.05e0 /* end of PARAMETER translations */ void /*FUNCTION*/ snqqrf( long m, long n, float *a, long lda, LOGICAL32 pivot, long ipvt[], long lipvt, float rdiag[], float acnorm[], float wa[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) long int _l0, i, j, jp1, k, kmax, minmn; float ajnorm, epsmch, sum, temp; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Acnorm = &acnorm[0] - 1; long *const Ipvt = &ipvt[0] - 1; float *const Rdiag = &rdiag[0] - 1; float *const Wa = &wa[0] - 1; /* end of OFFSET VECTORS */ /* ********** * * SUBROUTINE SNQQRF * * THIS SUBROUTINE USES HOUSEHOLDER TRANSFORMATIONS WITH COLUMN * PIVOTING (OPTIONAL) TO COMPUTE A QR FACTORIZATION OF THE * M BY N MATRIX A. THAT IS, SNQQRF DETERMINES AN ORTHOGONAL * MATRIX Q, A PERMUTATION MATRIX P, and AN UPPER TRAPEZOIDAL * MATRIX R WITH DIAGONAL ELEMENTS OF NONINCREASING MAGNITUDE, * SUCH THAT A*P = Q*R. THE HOUSEHOLDER TRANSFORMATION FOR * COLUMN K, K = 1,2,...,MIN(M,N), is OF THE FORM * * T * I - (1/U(K))*U*U * * WHERE U HAS ZEROS IN THE FIRST K-1 POSITIONS. THE FORM OF * THIS TRANSFORMATION and THE METHOD OF PIVOTING FIRST * APPEARED IN THE CORRESPONDING LINPACK SUBROUTINE. * * THE SUBROUTINE STATEMENT IS * * subroutine SNQQRF(M,N,A,LDA,PIVOT,IPVT,LIPVT,RDIAG,ACNORM,WA) * * WHERE * * M is A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER * OF ROWS OF A. * * N is A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER * OF COLUMNS OF A. * * A is AN M BY N ARRAY. ON INPUT A CONTAINS THE MATRIX FOR * WHICH THE QR FACTORIZATION is TO BE COMPUTED. ON OUTPUT * THE STRICT UPPER TRAPEZOIDAL PART OF A CONTAINS THE STRICT * UPPER TRAPEZOIDAL PART OF R, and THE LOWER TRAPEZOIDAL * PART OF A CONTAINS A FACTORED FORM OF Q (THE NON-TRIVIAL * ELEMENTS OF THE U VECTORS DESCRIBED ABOVE). * * LDA is A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M * WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY A. * * PIVOT is A LOGICAL INPUT VARIABLE. IF PIVOT is SET TRUE, * THEN COLUMN PIVOTING is ENFORCED. IF PIVOT is SET FALSE, * THEN NO COLUMN PIVOTING is DONE. * * IPVT is AN INTEGER OUTPUT ARRAY OF LENGTH LIPVT. IPVT * DEFINES THE PERMUTATION MATRIX P SUCH THAT A*P = Q*R. * COLUMN J OF P is COLUMN IPVT(J) OF THE IDENTITY MATRIX. * IF PIVOT is FALSE, IPVT is NOT REFERENCED. * * LIPVT is A POSITIVE INTEGER INPUT VARIABLE. IF PIVOT is FALSE, * THEN LIPVT MAY BE AS SMALL AS 1. IF PIVOT is TRUE, THEN * LIPVT MUST BE AT LEAST N. * * RDIAG is AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE * DIAGONAL ELEMENTS OF R. * * ACNORM is AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE * NORMS OF THE CORRESPONDING COLUMNS OF THE INPUT MATRIX A. * IF THIS INFORMATION is NOT NEEDED, THEN ACNORM CAN COINCIDE * WITH RDIAG. * * WA is A WORK ARRAY OF LENGTH N. IF PIVOT is FALSE, THEN WA * CAN COINCIDE WITH RDIAG. * * SUBPROGRAMS CALLED * * MINPACK-SUPPLIED ... R1MACH,SNRM2 * * FORTRAN-SUPPLIED ... max,sqrt,min * * ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. * BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE * * ********** * ------------------------------------------------------------------ */ /* EPSMCH is THE MACHINE PRECISION. * */ epsmch = FLT_EPSILON; /* COMPUTE THE INITIAL COLUMN NORMS and INITIALIZE SEVERAL ARRAYS. * */ for (j = 1; j <= n; j++) { Acnorm[j] = snrm2( m, &A(j - 1,0), 1 ); Rdiag[j] = Acnorm[j]; Wa[j] = Rdiag[j]; if (pivot) Ipvt[j] = j; } /* REDUCE A TO R WITH HOUSEHOLDER TRANSFORMATIONS. * */ minmn = min( m, n ); for (j = 1; j <= minmn; j++) { if (!pivot) goto L_40; /* BRING THE COLUMN OF LARGEST NORM INTO THE PIVOT POSITION. * */ kmax = j; for (k = j; k <= n; k++) { if (Rdiag[k] > Rdiag[kmax]) kmax = k; } if (kmax == j) goto L_40; for (i = 1; i <= m; i++) { temp = A(j - 1,i - 1); A(j - 1,i - 1) = A(kmax - 1,i - 1); A(kmax - 1,i - 1) = temp; } Rdiag[kmax] = Rdiag[j]; Wa[kmax] = Wa[j]; k = Ipvt[j]; Ipvt[j] = Ipvt[kmax]; Ipvt[kmax] = k; L_40: ; /* COMPUTE THE HOUSEHOLDER TRANSFORMATION TO REDUCE THE * J-TH COLUMN OF A TO A MULTIPLE OF THE J-TH UNIT VECTOR. * */ ajnorm = snrm2( m - j + 1, &A(j - 1,j - 1), 1 ); if (ajnorm == ZERO) goto L_100; if (A(j - 1,j - 1) < ZERO) ajnorm = -ajnorm; for (i = j; i <= m; i++) { A(j - 1,i - 1) /= ajnorm; } A(j - 1,j - 1) += ONE; /* APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS * and UPDATE THE NORMS. * */ jp1 = j + 1; if (n < jp1) goto L_100; for (k = jp1; k <= n; k++) { sum = ZERO; for (i = j; i <= m; i++) { sum += A(j - 1,i - 1)*A(k - 1,i - 1); } temp = sum/A(j - 1,j - 1); for (i = j; i <= m; i++) { A(k - 1,i - 1) += -temp*A(j - 1,i - 1); } if (!pivot || Rdiag[k] == ZERO) goto L_80; temp = A(k - 1,j - 1)/Rdiag[k]; Rdiag[k] *= sqrtf( fmaxf( ZERO, ONE - SQ(temp) ) ); if (P05*powif(Rdiag[k]/Wa[k],2) > epsmch) goto L_80; Rdiag[k] = snrm2( m - j, &A(k - 1,jp1 - 1), 1 ); Wa[k] = Rdiag[k]; L_80: ; } L_100: ; Rdiag[j] = -ajnorm; } return; /* Last line of subroutine SNQQRF. * */ #undef A } /* end of function */ /* PARAMETER translations */ #define P25 0.25e0 #define P5 0.5e0 /* end of PARAMETER translations */ void /*FUNCTION*/ snqupd( long m, long n, float s[], long ls, float u[], float v[], float w[], LOGICAL32 *sing) { long int _l0, i, j, jj, l, nm1, nmj; float cotan, tau, temp, vcos, vsin, vtan; static float giant = 0.0e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const S = &s[0] - 1; float *const U = &u[0] - 1; float *const V = &v[0] - 1; float *const W = &w[0] - 1; /* end of OFFSET VECTORS */ /* ********** * * SUBROUTINE SNQUPD * * GIVEN AN M BY N LOWER TRAPEZOIDAL MATRIX S, AN M-VECTOR U, * and AN N-VECTOR V, THE PROBLEM is TO DETERMINE AN * ORTHOGONAL MATRIX Q SUCH THAT * * T * (S + U*V )*Q * * is AGAIN LOWER TRAPEZOIDAL. * * THIS SUBROUTINE DETERMINES Q AS THE PRODUCT OF 2*(N - 1) * TRANSFORMATIONS * * GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1) * * WHERE GV(I), GW(I) ARE GIVENS ROTATIONS IN THE (I,N) PLANE * WHICH ELIMINATE ELEMENTS IN THE I-TH and N-TH PLANES, * RESPECTIVELY. Q ITSELF is NOT ACCUMULATED, RATHER THE * INFORMATION TO RECOVER THE GV, GW ROTATIONS is RETURNED. * * THE SUBROUTINE STATEMENT IS * * subroutine SNQUPD(M,N,S,LS,U,V,W,SING) * * WHERE * * M is A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER * OF ROWS OF S. * * N is A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER * OF COLUMNS OF S. N MUST NOT EXCEED M. * * S is AN ARRAY OF LENGTH LS. ON INPUT S MUST CONTAIN THE LOWER * TRAPEZOIDAL MATRIX S STORED BY COLUMNS. ON OUTPUT S CONTAINS * THE LOWER TRAPEZOIDAL MATRIX PRODUCED AS DESCRIBED ABOVE. * * LS is A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN * (N*(2*M-N+1))/2. * * U is AN INPUT ARRAY OF LENGTH M WHICH MUST CONTAIN THE * VECTOR U. * * V is AN ARRAY OF LENGTH N. ON INPUT V MUST CONTAIN THE VECTOR * V. ON OUTPUT V(I) CONTAINS THE INFORMATION NECESSARY TO * RECOVER THE GIVENS ROTATION GV(I) DESCRIBED ABOVE. * * W is AN OUTPUT ARRAY OF LENGTH M. W(I) CONTAINS INFORMATION * NECESSARY TO RECOVER THE GIVENS ROTATION GW(I) DESCRIBED * ABOVE. * * SING is A LOGICAL OUTPUT VARIABLE. SING is SET TRUE IF ANY * OF THE DIAGONAL ELEMENTS OF THE OUTPUT S ARE ZERO. OTHERWISE * SING is SET FALSE. * * SUBPROGRAMS CALLED * * MINPACK-SUPPLIED ... R1MACH * * FORTRAN-SUPPLIED ... abs,sqrt * * ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. * BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE, * JOHN L. NAZARETH * * ********** * ------------------------------------------------------------------ */ /* GIANT is THE LARGEST MAGNITUDE. * */ if (giant == 0.0e0) giant = FLT_MAX; /* INITIALIZE THE DIAGONAL ELEMENT POINTER. * */ jj = (n*(2*m - n + 1))/2 - (m - n); /* MOVE THE NONTRIVIAL PART OF THE LAST COLUMN OF S INTO W. * */ l = jj; for (i = n; i <= m; i++) { W[i] = S[l]; l += 1; } /* ROTATE THE VECTOR V INTO A MULTIPLE OF THE N-TH UNIT VECTOR * IN SUCH A WAY THAT A SPIKE is INTRODUCED INTO W. * */ nm1 = n - 1; if (nm1 < 1) goto L_70; for (nmj = 1; nmj <= nm1; nmj++) { j = n - nmj; jj -= m - j + 1; W[j] = ZERO; if (V[j] == ZERO) goto L_50; /* DETERMINE A GIVENS ROTATION WHICH ELIMINATES THE * J-TH ELEMENT OF V. * */ if (fabsf( V[n] ) >= fabsf( V[j] )) goto L_20; cotan = V[n]/V[j]; vsin = P5/sqrtf( P25 + P25*SQ(cotan) ); vcos = vsin*cotan; tau = ONE; if (fabsf( vcos )*giant > ONE) tau = ONE/vcos; goto L_30; L_20: ; vtan = V[j]/V[n]; vcos = P5/sqrtf( P25 + P25*SQ(vtan) ); vsin = vcos*vtan; tau = vsin; L_30: ; /* APPLY THE TRANSFORMATION TO V and STORE THE INFORMATION * NECESSARY TO RECOVER THE GIVENS ROTATION. * */ V[n] = vsin*V[j] + vcos*V[n]; V[j] = tau; /* APPLY THE TRANSFORMATION TO S and EXTEND THE SPIKE IN W. * */ l = jj; for (i = j; i <= m; i++) { temp = vcos*S[l] - vsin*W[i]; W[i] = vsin*S[l] + vcos*W[i]; S[l] = temp; l += 1; } L_50: ; } L_70: ; /* ADD THE SPIKE FROM THE RANK 1 UPDATE TO W. * */ for (i = 1; i <= m; i++) { W[i] += V[n]*U[i]; } /* ELIMINATE THE SPIKE. * */ *sing = FALSE; if (nm1 < 1) goto L_140; for (j = 1; j <= nm1; j++) { if (W[j] == ZERO) goto L_120; /* DETERMINE A GIVENS ROTATION WHICH ELIMINATES THE * J-TH ELEMENT OF THE SPIKE. * */ if (fabsf( S[jj] ) >= fabsf( W[j] )) goto L_90; cotan = S[jj]/W[j]; vsin = P5/sqrtf( P25 + P25*SQ(cotan) ); vcos = vsin*cotan; tau = ONE; if (fabsf( vcos )*giant > ONE) tau = ONE/vcos; goto L_100; L_90: ; vtan = W[j]/S[jj]; vcos = P5/sqrtf( P25 + P25*SQ(vtan) ); vsin = vcos*vtan; tau = vsin; L_100: ; /* APPLY THE TRANSFORMATION TO S and REDUCE THE SPIKE IN W. * */ l = jj; for (i = j; i <= m; i++) { temp = vcos*S[l] + vsin*W[i]; W[i] = -vsin*S[l] + vcos*W[i]; S[l] = temp; l += 1; } /* STORE THE INFORMATION NECESSARY TO RECOVER THE * GIVENS ROTATION. * */ W[j] = tau; L_120: ; /* TEST FOR ZERO DIAGONAL ELEMENTS IN THE OUTPUT S. * */ if (S[jj] == ZERO) *sing = TRUE; jj += m - j + 1; } L_140: ; /* MOVE W BACK INTO THE LAST COLUMN OF THE OUTPUT S. * */ l = jj; for (i = n; i <= m; i++) { S[l] = W[i]; l += 1; } if (S[jj] == ZERO) *sing = TRUE; return; /* Last line of subroutine SNQUPD. * */ } /* end of function */