/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:32:03 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "sblse.h" #include #include /* PARAMETER translations */ #define HALF 0.5e0 #define ONE 1.0e0 #define TWO 2.0e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ sblse( float *a, long lda, long mtrue, long n, long m1true, float b[], float bnd[][2], float unbnd, long kprint, float tol, long *ierr, float x[], float *rnorm, long *nsets, float w[], float size[], float tnorm[], float z[], float cc[], float ss[], long index[], long jstat[], float xt[], float rt[], float diff[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) LOGICAL32 cmode, find, free, free1, free2, freejs[2], hitbnd, movejf, msg1; long int bmode, i, i1, ibound, id, if0, if1, if2, ifnew, ig1, igaus1, igaus2, igiv1, igiv2, ii, iis, imax, iput, iput1, iput2, ir, irlast, is, isnew, itemp, iter, itmax, j, jcol, jd, jf, jfnew, jj, jjs, jpiv, jput, js, jsnew, kstep, l, lbound, len, m, m1, m1p1, mdb, minmn, nf, npr004, npr005, npr009, npr012, npr014, npr021, npr022, nsp1; float a1, absa1, alpha, amax, b1, bndl, bsave, bsize, btemp, fac, range, t, temp, tmax, up, wmax, xtry; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const B = &b[0] - 1; float *const Cc = &cc[0] - 1; float *const Diff = &diff[0] - 1; LOGICAL32 *const Freejs = &freejs[0] - 1; long *const Index = &index[0] - 1; long *const Jstat = &jstat[0] - 1; float *const Rt = &rt[0] - 1; float *const Size = &size[0] - 1; float *const Ss = &ss[0] - 1; float *const Tnorm = &tnorm[0] - 1; float *const W = &w[0] - 1; float *const X = &x[0] - 1; float *const Xt = &xt[0] - 1; float *const Z = &z[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 SBLSE Krogh Minor change for making .f90 version. *>> 2001-01-16 SBLSE Krogh Fixed problem with conversion to C. *>> 2000-12-01 SBLSE Krogh Removed unused parameter Factor. *>> 1999-12-22 SBLSE Krogh Declared some 'NPR' variables. *>> 1996-03-30 SBLSE Krogh Added external statement. *>> 1994-10-19 SBLSE Krogh Changes to use M77CON *>> 1992-10-09 SBLSE CLL Solve for correction rather than new X. *>> 1992-09-02 CLL Quit doing look-ahead updating of the B() array. *>> 1992-09-02 CLL Add random choice to combat cycling. *>> 1992-09-01 CLL Change elimination method in constraint matrix. *>> 1992-09-01 CLL Add more updating of SIZE(). *>> 1992-06-17 CLL Minor editing of comments. *>> 1990-07-04 CLL Minor editing of comments. *>> 1988-04-26 CLL *>> 1988-01-04 CLL * * *** Bounded variables Least Squares with Equality constraints *** * * In the following descriptions M denotes the input value MTRUE, and * M1 denotes the input value M1TRUE. * Given an M by N matrix, A, and an M-vector, b, compute an * N-vector, x, that satisfies the first M1 equations of the linear * system Ax = b as equality constraints, and satisfies the remaining * rows of the system in the least squares sense, subject to bounds * on the components of x. * * The bounds are specified in the form: * * BND(1,J) .le. X(J) .le. BND(2,J) * * where the special value UNBND is used with the meaning that * if BND(1,J) = UNBND, X(J) is unbounded below, and * if BND(2,J) = UNBND, X(J) is unbounded above. * * A problem with general linear inequality constraints can be put in * the form acceptable to this subroutine by the introduction of * slack variables. For example suppose a problem has inequality * constraints represented as Cx .ge. d where C is an M3 x N matrix. * Add M3 components to the end of the x vector, impose nonnegativity * bounds on these new components, and write the equality con- * straints, [C : -I] x = d, where I denotes the identity matrix of * order M3. * * In anticipation of this usage, this subroutine searches * for columns in the constraint equations that have only a single * nonzero element (i.e., singletons) and uses these, if the signs * permit, for efficient construction of the initial * triangularization of the constraint equations. * ------------------------------------------------------------------ * This code is a substantial generalization of the subroutine, * NNLS, that solves the least squares problem, Ax = b, subject to * all x(j) .ge. 0. The subroutine NNLS appeared in SOLVING * LEAST SQUARES PROBLEMS, by Lawson and Hanson, Prentice-Hall, 1974. * * This algorithm has two major phases, the initial solution of * the constraint equations, subject to the bounds on the variables, * and then the solution of the least-squares problem, subject to the * equality constraints and the bounds on the variables. * * The following Sftran procedures are used to get started and * to finish: * * procedure( INITIALIZE ) * procedure( TERMINATION ) * * The following Sftran procedures are used only in Phase 1: * * procedure( TRIANGULARIZE USING SINGLETONS ) * procedure( ALGORITHM BVEQ ) * procedure( BVEQ: SELECT ANOTHER COEF TO SOLVE FOR ) * procedure( BVEQ: MOVE JPNEW FROM SET F TO SET S ) * procedure( MSG FOR IERR = -1 ) * procedure( REARRANGE ROWS ) * * The following Sftran procedures are used only in Phase 2: * * procedure( ALGORITHM BVLS ) * procedure( BVLS: SELECT COEF TO SOLVE FOR ) * procedure( BVLS: MOVE JPNEW FROM SET F TO SET S ) * procedure( GIVENS ROTATIONS, USING IGIV1 AND IGIV2 ) * procedure( HOUSEHOLDER TRIANGULARIZATION ) * procedure( ELIMINATE IN LOWER ROWS ) * * The following Sftran procs are used in both Phases 1 & 2: * * procedure( CONTROL CONSTRAINT TESTING FOR SET S ) * procedure( SEE IF ALL CONSTRAINED COEFFS ARE FEASIBLE ) * procedure( PERTURB ANY OUTLIERS IN SET S TO THE BOUNDARY ) * procedure( MOVE COEF JZNEW FROM SET S TO SET F ) * procedure( GAUSSIAN ELIMINATION, USING IGAUS1 AND IGAUS2 ) * procedure( SOLVE TRIANGULAR SYSTEM, INPUT AND OUTPUT IN Z() ) * procedure( DEBUG ) * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * BVLS was developed by C. L. Lawson and R. J. Hanson at * JPL, 1973 June 12. This version was identified as BVLS2. * Changes made by Lawson, JPL, 1982 Sept. * Added NSETS into the argument list. */ /* 1982 Sept 16. Converted code to SFTRAN3. New version was called * BVLS5. * 1986 June 25. C. L. Lawson. Adapting code to MATH77 style. * 1987 Oct. CLL Changing name to SBLS. Using new MATH77 versions * of second-level subroutines. Developing SBLSE from SBLS to * solve the equality-constrained bounded variable problem. * 1988 April CLL & RJH. Added code to delete constraint rows that * are inconsistent or dependent. * 1992 Sept/Oct CLL Major changes. Solving for change to X rather * than directly for X. This changes formulas for updating the * right-side vector in B(). No longer need array SIZE(). * ------------------------------------------------------------------ * SUBROUTINE ARGUMENTS * * A() [In/Out] On entry A() contains the MTRUE by N matrix, A. * On return A() contains the product matrix, G*A, * where G is an MTRUE by MTRUE matrix generated * implicitly by this subroutine. * * LDA [In] LDA is the first dimensioning parameter for the * array, A(). Require LDA .ge. MTRUE. * * MTRUE [in] Total number of rows in the matrix, A. Require * MTRUE > 0. * * N [in] Number of columns in the matrix, A. Require N > 0. * * M1TRUE [in] Number of leading rows of [A,b] to be treated as * equality constraints. User must set M1TRUE = 0 if * there are no equality constraint rows. * Require 0 .le. M1TRUE .le. min(MTRUE, N). * * B() [In/Out] On entry B() contains the MTRUE-vector, B. * On return, B() contains G*(B - A[Z] * X[Z] ) where A[Z] * and X[Z] denote the parts indexed in Set F. * * BND(,) [In] If BND(1,J) .eq. UNBND then X(J) is unbounded below. * Otherwise BND(1,J) is the lower bound for X(J). * If BND(2,J) .eq. UNBND then X(J) is unbounded above. * Otherwise BND(2,J) is the upper bound for X(J). * If BND(1,J) .ne. UNBND and BND(2,J) .ne. UNBND, require * BND(1,J) .le. BND(2,J). * * UNBND [In] Special value to be used in BND(,) to indicate * unboundedness. * * KPRINT [in] Larger values give more printing. * 0: No printing. * .ge. 1: Print initial values of each component of X() and * a line for each move of an index between Sets * Z and P. * .ge. 2: Print details showing how index is selected to * move between Sets Z and P. * .ge. 3: Use extra input items XT(), RT(), DIFF() to * compute residuals at various stages to verify * transformed problems stay consistent with the * original problem. * .ge. 4: Print full data array [A : b] at various stages. * * TOL [in] Tolerance used in tests of variables against * their bounds, and in testing the right side of constraint * equations for zero. Used in a relative sense in testing * against nonzero bounds, and in an absolute sense testing * against zero. Suggest a value such as EPS**(0.75), however * a larger value has been found necessary to force * triangularization of the constraint equations in some cases. * Here EPS denotes the machine precision, i.e. the smallest * number EPS such that 1. + EPS is not 1. * * IERR [Out] Status indicator on return. * = -1 Failed to fully triangularize first M1TRUE rows. * Subroutine attempts to complete the computation, * omitting the constraint rows not triangularized. * User should check residuals of constraint rows. * = 0 No errors detected. * = 1 MTRUE .le. 0 or N .le. 0 * = 2 Inconsistent setting of bounds. * = 3 Too many iterations. Nominal ITMAX = 5*N. * * X() [Out] On entry X() need not be initialized. On return, * X() will contain the solution N-vector. * * RNORM [Out] The euclidean norm of the residual vector. * * NSETS [Out] Indicates the No. of components of the solution * vector, X(), that are in Set S. These values are obtained * by solving a system of equations, and thus will generally * not be at a bound. * * W() [Out] An N-array of working space. On return, W() will * contain the dual solution vector. This will satisfy * W(i) = 0 for all i in Set S, * W(i) .le. 0 for all i in Set F, such that X(i) is at its * lower bound, and * W(i) .ge. 0 for all i in Set F, such that X(i) is at its * upper bound. * * SIZE() [Scratch] !!! No longer needed 1992-10-12 !!! * Working space dimensioned at least M1TRUE. * SIZE(i) holds scaled roundoff error estimates for B(i). */ /* The roundoff error estimate is TOL times SIZE(i). * * TNORM() [Scratch] An N-array of working space. * * Z() [Scratch] Work array of size min(MTRUE, N). Used during * solution of triangular system to hold right-side vector * initially and the solution vector at conclusion. Also * used during interpolation to restore feasibility. * * CC(), SS() [Scratch] Work arrays of size min(MTRUE,N) each. * Used to hold cosines and sines of the Givens * transformations in order to support a column-oriented * implementation. * * INDEX() [Out] An INTEGER working array of length at least N. * On exit the contents of this array define the Sets * S, F, and H. * * JSTAT() [Scratch] Of length N. * * XT(1:N), RT(1:MTRUE), DIFF(1:MTRUE) For use in testing this subr. * Could be deleted in a production version. These arrays are * accessed only if KPRINT .ge. 3. * On entry XT() and RT() should contain the known * solution and residual vectors respectively. Diff() is scratch. * ------------------------------------------------------------------ * Description of Sets S, F, and H. * * The N components of X() are classified into three disjoint sets. * * INDEX(1) thru INDEX(NSETS) = Set S. * INDEX(IF1) thru INDEX(IF2) = Set F. * INDEX(IF2+1) thru INDEX(N) = Set H. * * NF = IF2 - IF1 + 1 * IF1 = NSETS + 1 = NSP1 * NSETS + NF .le. N * * The sizes of sets S and F are given by NSETS and NF. * Set H consists of the components that are constrained to a * unique value by the given constraints. * Sets S and F contain components that a priori are allowed a * possible range of values of nonzero length. * During the computation, set F contains components that are * temporarily being held constant. * Set S contains components whose values are determined by solving * the currently triangularized subsystem, or by adjusting such a * solution set to shift its values from infeasible to feasible. * Members of set F are available for consideration for moving to * set S. * ------------------------------------------------------------------ *--S replaces "?": ?BLSE, ?ERV1, ?HTCC, ?ROTG, ?ROT, ?DOT *--& ?NRM2, ?AXPY, ?COPY, ?SWAP, ?RANU * Other subprograms referenced directly: ERMOR, ERMSG, IERM1, IERV1 * ------------------------------------------------------------------ * Interesting internal variables * * MSG1 [logical] Used to cause the error message associated with * IERR = -1 to be issued only the first time the error is * detected. False means the message has not been issued. * True means it has. * ------------------------------------------------------------------ */ /* NPR Variables not declared by SFTRAN */ /* ------------------------------------------------------------------ */ goto L_30001; L_20002: if (*ierr != 0) return; if (!(m1 > 0)) goto L_20005; goto L_30002; /* Here all singletons with the right signs have been brought * into Set S. If NSETS .lt. M1 we must build the triangle * further. * */ L_20006: if (!(*nsets < m1)) goto L_20008; cmode = TRUE; goto L_30003; L_20009: ; L_20008: npr004 = 20010; goto L_30004; L_20010: ; L_20005: if ((*ierr > 0 || *nsets == minmn) || nf == 0) goto L_20003; for (jf = 1; jf <= n; jf++) { if (X[jf] != ZERO) { saxpy( m - m1, -X[jf], &A(jf - 1,m1p1 - 1), 1, &B[m1p1], 1 ); } } cmode = FALSE; if (!(m1 > 0)) goto L_20017; ig1 = 1; npr005 = 20018; goto L_30005; L_20018: ; L_20017: npr004 = 20019; goto L_30004; L_20019: ; goto L_30006; L_20020: npr004 = 20021; goto L_30004; L_20021: ; L_20003: goto L_30007; L_20022: return; /* ------------------------------------------------------------------ */ /* procedure( INITIALIZE ) * */ L_30001: m = mtrue; m1 = m1true; if (m <= 0 || n <= 0) { *ierr = 1; ierm1( "SBLSE", *ierr, 0, "Require M > 0 and N > 0.", "M", m, ',' ); ierv1( "N", n, '.' ); return; } *ierr = 0; iter = 0; itmax = 6*n; msg1 = FALSE; m1p1 = m1 + 1; minmn = min( m, n ); if2 = n; if1 = 1; *nsets = 0; nsp1 = 1; /* . Initialize INDEX() and set W() = zero. Setting W() = * . zero is not needed for the functioning of this * . subroutine, it is just done so the W() array will not * . have garbage in it if debugging printing is used. * */ for (j = 1; j <= n; j++) { Index[j] = j; W[j] = ZERO; } /* . BEGIN.. LOOP ON IF0 TO INITIALIZE X() */ if0 = if1; L_20028: if (if0 <= if2) { j = Index[if0]; if (bnd[j - 1][0] == unbnd) { if (bnd[j - 1][1] == unbnd) { X[j] = ZERO; } else { X[j] = fminf( ZERO, bnd[j - 1][1] ); } } else if (bnd[j - 1][1] == unbnd) { X[j] = fmaxf( ZERO, bnd[j - 1][0] ); } else { range = bnd[j - 1][1] - bnd[j - 1][0]; if (range == ZERO) { /* . Here X(J) is constrained to a single value. * */ Index[if0] = Index[if2]; Index[if2] = j; if0 -= 1; if2 -= 1; W[j] = ZERO; X[j] = bnd[j - 1][0]; } else if (range > ZERO) { /* . The following statement sets X(J) to 0 if the * . constraint interval includes 0, and otherwise sets * . X(J) to the endpoint of the constraint interval * . that is closest to 0. * */ X[j] = fmaxf( bnd[j - 1][0], fminf( bnd[j - 1][1], ZERO ) ); } else { /* ****** Error: IERR = 2 ****** */ *ierr = 2; ierm1( "SBLSE", *ierr, 0, "Inconsistent setting of bound for X(J)." , "J", j, ',' ); serv1( "UNBND", unbnd, ',' ); serv1( "BND(1,J)", bnd[j - 1][0], ',' ); serv1( "BND(2,J)", bnd[j - 1][1], ',' ); goto L_31001; } } if (kprint >= 1) { printf(" Initializing: J=%4ld, X(J)=%14.6g\n", j, X[j]); } if0 += 1; goto L_20028; } /* . END.. LOOP ON IF0 TO INITIALIZE X() */ nf = if2 - if1 + 1; L_31001: goto L_20002; /* ------------------------------------------------------------------ */ /* procedure( TRIANGULARIZE USING SINGLETONS ) * * Search for singletons. Use them to start the triangle if * their signs are right. * */ L_30002: itemp = if1; for (if0 = itemp; if0 <= if2; if0++) { if (*nsets == m1) goto L_20039; jf = Index[if0]; /* . The following loop sets II. * . If A(1..M, JF) contains exactly one nonzero, II will * . be set to the row index of that element. * . Otherwise II will be set to zero. * */ ii = 0; for (i = 1; i <= m; i++) { if (A(jf - 1,i - 1) != ZERO) { if (ii != 0) { ii = 0; goto L_20042; } ii = i; } } L_20042: ; if (ii > *nsets && ii <= m1) { btemp = B[ii] - sdot( n, &A(0,ii - 1), lda, x, 1 ); xtry = X[jf] + btemp/A(jf - 1,ii - 1); if ((bnd[jf - 1][0] == unbnd || xtry >= bnd[jf - 1][0]) && (bnd[jf - 1][1] == unbnd || xtry <= bnd[jf - 1][1])) { movejf = TRUE; } else { bsize = fabsf( B[ii] ); for (j = 1; j <= n; j++) { bsize += fabsf( A(j - 1,ii - 1)*X[j] ); } if (fabsf( btemp ) < tol*bsize) { movejf = TRUE; xtry = X[jf]; } else { movejf = FALSE; } } if (movejf) { /* . Move JF from Set F to Set S. * */ X[jf] = xtry; if (ii != nsp1) { sswap( n, &A(0,ii - 1), lda, &A(0,nsp1 - 1), lda ); B[ii] = B[nsp1]; } B[nsp1] = ZERO; if (if0 != if1) { Index[if0] = Index[if1]; Index[if1] = jf; } *nsets = nsp1; nsp1 = *nsets + 1; if1 = nsp1; nf -= 1; } } } L_20039: ; if (kprint >= 1) { printf(" Started using singletons. NSETS=%4ld\n", *nsets); } goto L_20006; /* ------------------------------------------------------------------ */ /* procedure( ALGORITHM BVEQ ) * * BVEQ == Bounded Variables solution of linear EQuations. * Solve this problem using rows 1 through M1 and the * variables indexed in Set S and Set F. * This proc will decrease the values of M1 and M if rows of the * equality constraint system are linearly dependent on earlier rows. * */ L_30003: if (kprint >= 1) { printf(" \n Entering Proc BVEQ.. NSETS=%5ld, NF=%5ld\n", *nsets, nf); } kstep = nsp1; L_20063: if (!(kstep <= m1)) goto L_20064; goto L_30008; L_20065: ; L_20066: ; /* . Quit on error flag, or * . if Set F is empty. * */ if (!(nf == 0)) goto L_20069; if (!(*nsets < m1)) goto L_20070; npr009 = 20070; goto L_30009; L_20070: goto L_20064; L_20069: if (!(*ierr > 0)) goto L_20071; goto L_20064; L_20071: ; goto L_30010; L_20072: if (!(find)) goto L_20074; goto L_30011; /* . Delta_x is now in Z(). * * . The following proc may move one variable from Set S * . to Set F. All variables remaining in Set S when it * . completes will be feasible. May set IERR positive. * */ L_20075: npr012 = 20076; goto L_30012; L_20076: if (!hitbnd) { kstep += 1; goto L_20067; } goto L_20073; /* . The following proc deletes a constraint row that appears * . to be dependant on preceeding constraint rows or * . inconsistent with preceeding constraint rows. It also * . decrements M1 and M. Will set IERR negative. * */ L_20074: npr009 = 20079; goto L_30009; L_20079: goto L_30013; L_20080: goto L_20067; L_20073: goto L_20066; L_20067: ; goto L_20063; L_20064: goto L_20009; /* ------------------------------------------------------------------ */ /* procedure( BVEQ: SELECT ANOTHER COEF TO SOLVE FOR ) * * . Search thru Set F for a new coef to solve for. * Pick column having elt of largest magnitude in row NSP1 * among those having acceptable sign. * On entry BSAVE contains the original input value of B(NSP1), and * B(NSP1) contains the residual for row NSP1 using the current X(). * Index JF = INDEX(IF0) is acceptable * if X(JF) is free to increase and B(NSP1)/A(NSP1,JF) .ge. zero, or * if X(JF) is free to decrease and B(NSP1)/A(NSP1,JF) .le. zero. * * . If no new coeff is found, then check abs(B(NSP1)) for being * small enough that it might be zero apart from roundoff error. * If not, set FIND = false and exit. * If so, set B(NSP1) = zero and try the search again. * * . We wish to force the maximum number of * components into the solution, and thus we accept components that * will enter with a zero change as well as those that will enter * with a change in the acceptable direction. * . If a coef is selected, set FIND = true. * The index of the selected coef is JSNEW = INDEX(ISNEW). * . If no coef is selected, set FIND = false. * */ L_30010: find = FALSE; amax = tol; isnew = 0; b1 = B[nsp1]; if (kprint >= 4) { printf(" \n Before selecting next JSNEW. [A:b] =\n \n"); for (i = 1; i <= (*nsets + 1); i++) { printf(" %3ld", i); printf(" "); for (j = 1; j <= n; j++) { printf("%13.6g", A(j - 1,i - 1)); } printf("%13.6g\n", B[i]); } } for (bmode = 1; bmode <= 2; bmode++) { if (kprint >= 2) { printf(" Scan Set F with B(NSP1) = %14.6g, BMODE = %2ld\n " " NSP1 = %4ld\n " " -- JF,A(NSP1,JF),FREE1,FREE2,FREE, JFbest --\n", b1, bmode, nsp1); } /* . Begin loop through Set F. */ for (if0 = if1; if0 <= if2; if0++) { jf = Index[if0]; a1 = A(jf - 1,nsp1 - 1); absa1 = fabsf( a1 ); if (absa1 <= amax) { if (kprint >= 2) { printf(" %4ld%14.6g\n", jf, a1); } } else { /* . Set FREE1 = true if X(JF) is not at the left end-point of * . its constraint region. * . Set FREE2 = true if X(JF) is not at the right end-point of * . its constraint region. * . Set FREE = true if X(JF) is not at either end-point of its * . constraint region. * */ free1 = bnd[jf - 1][0] == unbnd || X[jf] != bnd[jf - 1][0]; free2 = bnd[jf - 1][1] == unbnd || X[jf] != bnd[jf - 1][1]; free = free1 && free2; if (free) { amax = absa1; isnew = if0; } else if (free1) { if (b1/a1 <= ZERO) { amax = absa1; isnew = if0; } } else if (free2) { if (b1/a1 >= ZERO) { amax = absa1; isnew = if0; } } if (kprint >= 2) { if (isnew == 0) { printf(" %4ld%14.6g%2c%2c%2c\n", jf, a1, TorF(free1), TorF(free2), TorF(free)); } else { printf(" %4ld%14.6g%2c%2c%2c %4ld\n", jf, a1, TorF(free1), TorF(free2), TorF(free), Index[isnew]); } } } } if (amax != tol) goto L_20087; if (bmode == 1 && b1 != ZERO) { bsize = fabsf( bsave ); for (j = 1; j <= n; j++) { bsize += fabsf( A(j - 1,nsp1 - 1)*X[j] ); } if (fabsf( b1 ) <= tol*bsize) { if (kprint >= 2) { printf(" Setting B(NSP1) = 0., NSP1=%4ld\n Old B(NSP1) =%14.6g, TOL*BSIZE =%14.6g\n", nsp1, B[nsp1], tol*bsize); } B[nsp1] = ZERO; b1 = ZERO; goto L_20086; } } find = FALSE; goto L_31010; L_20086: ; } L_20087: ; /* . Selected index is JSNEW = INDEX(ISNEW) */ find = TRUE; jsnew = Index[isnew]; Freejs[1] = bnd[jsnew - 1][0] == unbnd || X[jsnew] != bnd[jsnew - 1][0]; Freejs[2] = bnd[jsnew - 1][1] == unbnd || X[jsnew] != bnd[jsnew - 1][1]; L_31010: goto L_20072; /* ------------------------------------------------------------------ */ /* procedure( BVEQ: MOVE JSNEW FROM SET F TO SET S ) * * The index JSNEW=INDEX(IFNEW) has been selected to be moved from * Set F to Set S. * B(1..M) has been adjusted as though X(JSNEW) = 0. * Update indices moving JSNEW from Set F to Set S. * Do one stage of Gaussian elimination in the upper matrix. * Copy B() to Z() and solve for Delta_x, leaving it in Z(). * Set W(JSNEW) = 0. * */ L_30011: if (kprint >= 1) { printf(" ==> Move index from F to S. JSNEW =%4ld\n", jsnew); } Index[isnew] = Index[if1]; Index[if1] = jsnew; if1 += 1; *nsets = nsp1; nsp1 += 1; nf -= 1; len = m1 - *nsets; W[jsnew] = ZERO; scopy( *nsets, b, 1, z, 1 ); npr014 = 20115; goto L_30014; L_20115: goto L_20075; /* ------------------------------------------------------------------ * procedure( MSG FOR IERR = -1 ) */ L_30009: *ierr = -1; if (!msg1) { msg1 = TRUE; ermsg( "SBLSE", *ierr, 0, "Failed to fully triangularize the equality constraint rows." , ',' ); ermor( "Could be due to dependent or inconsistent constraints." , ',' ); ermor( "The computation will be completed, ignoring the", ',' ); ermor( "un-triangularized constraint rows. User should check" , ',' ); ermor( "residuals of constraint equations.", '.' ); } switch (npr009) { case 20070: goto L_20070; case 20079: goto L_20079; } /* ------------------------------------------------------------------ * procedure( REARRANGE ROWS ) */ /* Here when Row NSP1 of A() appears to be either dependent or * inconsistent relative to preceeding constraint rows. * Delete Row NSP1. * Have NSP1 .le. M1. * If NSP1 < M1 copy row M1 to row NSP1. * Copy row M to row M1. * Decrease M1 and M each by 1. * * . If NSP1 < M1 copy row M1 to row NSP1. * */ L_30013: if (nsp1 < m1) { scopy( n, &A(0,m1 - 1), lda, &A(0,nsp1 - 1), lda ); B[nsp1] = B[m1]; Rt[nsp1] = Rt[m1]; } /* . Copy row M to row M1. * */ scopy( n, &A(0,m - 1), lda, &A(0,m1 - 1), lda ); B[m1] = B[m]; Rt[m1] = Rt[m]; /* . Decrement M1 and M. */ m -= 1; m1 -= 1; m1p1 = m1 + 1; minmn = min( m, n ); goto L_20080; /* ------------------------------------------------------------------ */ /* procedure( ALGORITHM BVLS ) * * . BVLS == Bounded Variables Least Squares * Solve this problem using rows 1 through M and the * variables indexed in Set S and Set F. * */ L_30006: if (kprint >= 1) { printf(" \n Entering Proc BVLS.. CMODE=%2c, M=%5ld,\n NSETS=%5ld" ", NF=%5ld\n", TorF(cmode), m, *nsets, nf); } L_20120: ; /* . Quit if Set F is empty, or * . if M cols have been triangularized, or * . on error flag. */ if (nf == 0 || *nsets == m) goto L_20121; goto L_30015; /* . If no index was found to be moved from set Z to set P, * . then go to termination. * */ L_20122: if (!find) goto L_20121; goto L_30016; L_20123: ; L_20124: ; /* When this loop is entered, NSETS will be > M1. * NSETS decreases by one on each repetition of this loop, but * will still be > M1 at the beginning of every * repetition of this loop. * */ scopy( *nsets, b, 1, z, 1 ); npr014 = 20126; goto L_30014; L_20126: npr012 = 20127; goto L_30012; /* . The above proc will set IERR nonzero if * . the iteration count is exceeded. * . Otherwise the proc will either set HITBND := true and * . decrement NSETS, or set HITBND := false and leave * . NSETS unchanged. NSETS will be .ge. M1. * */ L_20127: if (*ierr > 0) goto L_20121; if (!hitbnd || *nsets == m1) goto L_20125; goto L_20124; L_20125: ; /* . All coeffs in Set S are feasible. Loop back. */ goto L_20120; L_20121: ; goto L_20020; /* ------------------------------------------------------------------ */ /* procedure( BVLS: SELECT COEF TO SOLVE FOR ) * * . Search thru Set F for a new coef to solve for. * First compute dual coeff, W(J) = negative gradient. * Next classify the members of Set F by setting ISTAT(). * . ISTAT() Meaning * . -1 X(J) is rejected for moving to Set S. Either * . W(J) = 0, or X(J) is at a bound and the sign of W(J) * . indicates the objective function will increase if * . X(J) is moved away from the bound. * . 1 X(J) is accepted for further consideration * . as a candidate for moving to Set S. W(J) is * . nonzero and if X(J) is at a bound, the sign of W(J) * . indicates moving away from the bound will decrease * . the objective function. * * Among those with ISTAT(J) > 0, choose one having column norm * exceeding a locally determined tolerance, and among these the * largest abs(W(J)). * . Note that when the sign of W(J) is uncertain we reject * moving X(J) to Set S. This biases the method toward keeping * variables in Set F when in doubt. * . If a coef is selected, set FIND = true. * The index of the selected coef is JSNEW = INDEX(ISNEW). * . If no coef is selected, set FIND = false. * */ L_30015: find = FALSE; tmax = ZERO; len = m - *nsets; if (kprint >= 4) { printf(" \n Before computing W(): [A:b] =\n \n"); for (i = 1; i <= m; i++) { printf(" %3ld", i); printf(" "); for (j = 1; j <= n; j++) { printf("%13.6g", A(j - 1,i - 1)); } printf("%13.6g\n", B[i]); } } if (kprint >= 2) { printf(" Scan Set F:\n " " -- J,FREE1,FREE2,FREE, W(J), JSTAT(J), TNORM(J), W(+) --\n"); } /* . Begin loop through Set F. */ for (if0 = if1; if0 <= if2; if0++) { j = Index[if0]; /* . Set FREE1 = true if X(J) is not at the left end-point of * . its constraint region. * . Set FREE2 = true if X(J) is not at the right end-point of * . its constraint region. * . Set FREE = true if X(J) is not at either end-point of its * . constraint region. * */ free1 = bnd[j - 1][0] == unbnd || X[j] != bnd[j - 1][0]; free2 = bnd[j - 1][1] == unbnd || X[j] != bnd[j - 1][1]; free = free1 && free2; /* . Compute dual coefficient W(J). * */ W[j] = sdot( len, &A(j - 1,nsp1 - 1), 1, &B[nsp1], 1 ); /* . The following code sets JSTAT(J) and computes TNORM(J) * . and TMAX. * . The columns whose dual signs are ok to be moved into Set S * . will be identified by JSTAT(J) = 1. Those not ok * . have JSTAT(J) = -1. * . For the ok columns the norm is computed and stored in * . TNORM(J) and the max of these is accumulated in TMAX. * */ Jstat[j] = -1; if (W[j] != ZERO) { if (free) { Jstat[j] = 1; Tnorm[j] = snrm2( len, &A(j - 1,nsp1 - 1), 1 ); tmax = fmaxf( tmax, Tnorm[j] ); } else if (free1) { if (W[j] < ZERO) { Tnorm[j] = snrm2( len, &A(j - 1,nsp1 - 1), 1 ); tmax = fmaxf( tmax, Tnorm[j] ); Jstat[j] = 1; } } else if (free2) { if (W[j] > ZERO) { Tnorm[j] = snrm2( len, &A(j - 1,nsp1 - 1), 1 ); tmax = fmaxf( tmax, Tnorm[j] ); Jstat[j] = 1; } } } if (kprint >= 2) { if (Jstat[j] == 1) { printf(" %4ld%2c%2c%2c%14.6g%4ld%14.6g\n", j, TorF(free1), TorF(free2), TorF(free), W[j], Jstat[j], Tnorm[j]); } else { printf(" %4ld%2c%2c%2c%14.6g%4ld\n", j, TorF(free1), TorF(free2), TorF(free), W[j], Jstat[j]); } } } if (tmax == ZERO) { find = FALSE; goto L_31015; } find = TRUE; /* . Among the ok columns having norm .ge. half*TMAX * . find one having largest abs(W(J)). * */ wmax = ZERO; for (if0 = if1; if0 <= if2; if0++) { j = Index[if0]; if (Jstat[j] > 0) { if (Tnorm[j] >= HALF*tmax) { if (fabsf( W[j] ) >= wmax) { isnew = if0; wmax = fabsf( W[j] ); } } } } /* . Selected index is JSNEW = INDEX(ISNEW) */ jsnew = Index[isnew]; Freejs[1] = bnd[jsnew - 1][0] == unbnd || X[jsnew] != bnd[jsnew - 1][0]; Freejs[2] = bnd[jsnew - 1][1] == unbnd || X[jsnew] != bnd[jsnew - 1][1]; L_31015: goto L_20122; /* ------------------------------------------------------------------ */ /* procedure( BVLS: MOVE JSNEW FROM SET F TO SET S ) * * The index JSNEW=INDEX(IF0) has been selected to be moved from * Set F to Set S. * Compute new Householder transformation. Apply it to A and b. * Update indices defining Sets P and Z. * Zero subdiagonal elts in col JSNEW, Set W(JSNEW)=0. * */ L_30016: if (kprint >= 1) { printf(" ==> Move index from F to S. JSNEW =%4ld\n", jsnew); } /* . Compute new Householder transformation and apply it to B(). * */ shtcc( 1, nsp1, nsp1 + 1, m, &A(jsnew - 1,0), &up, b, m, 1 ); /* . Update pointers. */ Index[isnew] = Index[if1]; Index[if1] = jsnew; if1 += 1; *nsets = nsp1; nsp1 += 1; nf -= 1; for (jf = if1; jf <= if2; jf++) { jj = Index[jf]; shtcc( 2, *nsets, nsp1, m, &A(jsnew - 1,0), &up, &A(jj - 1,0), lda, 1 ); } if (kprint >= 3) shtcc( 2, *nsets, nsp1, m, &A(jsnew - 1,0), &up, &Rt[1], m, 1 ); for (l = nsp1; l <= m; l++) { A(jsnew - 1,l - 1) = ZERO; } W[jsnew] = ZERO; goto L_20123; /* ------------------------------------------------------------------ */ /* procedure( CONTROL CONSTRAINT TESTING FOR SET S ) * * On entry, the Delta_x obtained by solving the current Set S * is in the array Z(). * If X + Delta_x satisfies the bounds, this proc returns with * X := X + Delta_x, HITBND := false, * B(min(M1,NSETS):NSETS) := 0. * Otherwise this proc returns with * X := X + ALPHA * Delta_x, HITBND := true, * B(min(M1,NSETS):NSETS) := (1 - ALPHA)*B(min(M1,NSETS):NSETS), * assigns values to IFNEW and IBOUND, and decrements NSETS. * */ L_30012: if (kprint >= 2) { printf(" Solved Set S obtaining Delta_x:\n"); for (i = 1; i <= *nsets; i++) { printf(" J,X(J)=%4ld%14.6g, D_X=%14.6g\n", Index[i], X[Index[i]] + Z[i], Z[i]); } } iter += 1; if (iter > itmax) { /* . ****** Error: IERR = 3 ****** */ *ierr = 3; ierm1( "SBLSE", *ierr, 0, "Iteration count ITER exceeds ITMAX." , "ITER", iter, ',' ); ierv1( "ITMAX", itmax, '.' ); goto L_31012; } goto L_30017; /* The above proc call sets HITBND. If HITBND = true then it * also sets ALPHA, IFNEW, and IBOUND. * */ L_20174: if (!(hitbnd)) goto L_20176; jfnew = Index[ifnew]; if (kprint >= 2) { printf(" Hit a bound. IFNEW=%4ld, JFNEW=INDEX(IFNEW)=%4ld,\n " "ALPHA=%14.6g\n", ifnew, jfnew, alpha); } /* Here ALPHA will be between 0 and 1 for use in computing * X() := X() + ALPHA * Delta_x. * */ for (is = 1; is <= *nsets; is++) { l = Index[is]; X[l] += alpha*Z[is]; if (kprint >= 2) { printf(" Interpolate to feasibility: L,X(L)=%4ld%14.6g\n", l, X[l]); } } for (i = min( m1p1, *nsets ); i <= *nsets; i++) { B[i] *= ONE - alpha; } goto L_30018; L_20183: goto L_20175; L_20176: for (is = 1; is <= *nsets; is++) { l = Index[is]; X[l] += Z[is]; } for (i = min( m1p1, *nsets ); i <= *nsets; i++) { B[i] = ZERO; } L_20175: if (!(*nsets > 0)) goto L_20191; goto L_30019; L_20192: ; L_20191: ; L_31012: switch (npr012) { case 20076: goto L_20076; case 20127: goto L_20127; } /* ------------------------------------------------------------------ */ /* procedure( SEE IF ALL CONSTRAINED COEFFS ARE FEASIBLE ) * * The Delta_x vector is currently in the array Z(). * See if each x + Delta_x in Set S is in the closure of its * constraint region. * If so, set HITBND = false. * If not, set HITBND = true, and also set ALPHA, IFNEW, and IBOUND. * Then ALPHA will satisfy 0. .le. ALPHA .lt. 1. * */ L_30017: alpha = TWO; for (is = 1; is <= *nsets; is++) { l = Index[is]; xtry = X[l] + Z[is]; if (bnd[l - 1][0] != unbnd && xtry < bnd[l - 1][0]) { /* . XTRY crosses lower bound. */ lbound = 1; } else if (bnd[l - 1][1] != unbnd && xtry > bnd[l - 1][1]) { /* . XTRY crosses upper bound. */ lbound = 2; } else { lbound = 0; } if (lbound != 0) { bndl = bnd[l - 1][lbound - 1]; t = (bndl - X[l])/(xtry - X[l]); if (kprint >= 2) { printf(" L=%4ld, LBOUND=%4ld, BND(LBOUND,L)=%15.7g\n X(L)=%15.7g, XTRY=%15.7g" ", T=%15.7g\n", l, lbound, bndl, X[l], xtry, t); } if (t <= alpha) { if (((t > ONE - tol || (l == jsnew && !Freejs[lbound] )) || (bndl == ZERO && fabsf( xtry ) <= tol)) || (bndl != ZERO && fabsf( xtry - bndl ) <= tol*fabsf( bndl ))) { /* . If the test of T is true * . we assume XTRY differs from its bound only due * . to roundoff. To avoid having this cause a move * . from Set S to Set F we set XTRY to its bound and * . ignore this value of T. * . Similarly, if the tests of L and FREEJS() are * . passed, it means the * . coeff JSNEW, which is the coeff most recently * . moved into P, has crossed the bound it was on just * . before it was moved into P. We assume this is due * . to roundoff error and adjust the value back to this * . bound. * . Similarly if the other tests are passed we assume * . XTRY differs from its bound only due to roundoff. * */ if (kprint >= 1) { printf(" While testing against bounds, adjust XTRY for roundoff error.\n " "IS=%4ld, L=INDEX(IS)=%4ld, Old XTRY=%16.8g\n", is, l, xtry); } Z[is] = bnd[l - 1][lbound - 1] - X[l]; if (kprint >= 1) { printf(" , New XTRY=%16.8g\n", X[l] + Z[is]); } } else { if (alpha == ZERO) { /* . Here T must be zero, and we have already saved a * . previous value of IS (into IFNEW) associated with * . a zero value of T. * . To reduce further the unlikely chance of * . cycling, we make a random selection between the * . current IS and the earlier saved one. As there is * . no need to continue the IS loop, we exit here. * */ if (sranu() < HALF) { ifnew = is; ibound = lbound; } goto L_20194; } else { alpha = t; ifnew = is; ibound = lbound; } } } } } L_20194: ; hitbnd = alpha != TWO; goto L_20174; /* ------------------------------------------------------------------ */ /* procedure( PERTURB ANY OUTLIERS IN SET S TO THE BOUNDARY ) * * . See if the remaining coeffs in set P are feasible. * They should be because of the way alpha was determined. * If any are infeasible it is due to round-off error. * Any that are infeasible or on a boundary will be set to * the boundary value and kept in set P. * */ L_30019: for (is = 1; is <= *nsets; is++) { ibound = 0; js = Index[is]; if (bnd[js - 1][0] != unbnd && X[js] < bnd[js - 1][0]) { ibound = 1; } else if (bnd[js - 1][1] != unbnd && X[js] > bnd[js - 1][1]) { ibound = 2; } if (ibound != 0) { /* . The current coeff that appears to have crossed a * . bound should, at worst be on its bound. * . We assume this is due to roundoff error and adjust * . the value back to the bound. */ if (kprint >= 1) { printf(" While testing against bounds, adjust X(JS) for roundoff error.\n " "JS=%4ld, Old X(JS)=%16.8g\n", js, X[js]); printf(" , New X(JS)=%16.8g\n", bnd[js - 1][ibound - 1]); } X[js] = bnd[js - 1][ibound - 1]; } } goto L_20192; /* ------------------------------------------------------------------ */ /* procedure( MOVE COEF JFNEW FROM SET S TO SET F ) * * Here we have JFNEW = INDEX(IFNEW) and this coeff has just hit its * bound indexed by IBOUND. * This proc is ref'd by Proc CONTROL CONSTRAINT TESTING, which is * ref'd by Proc BVEQ and Proc BVLS. * When reached via BVEQ we have CMODE = true * When reached via BVLS we have CMODE = false * */ L_30018: X[jfnew] = bnd[jfnew - 1][ibound - 1]; if (kprint >= 1) { printf(" ==> Move index from S to F. JFNEW =%5ld\n", jfnew); if (kprint >= 2) { printf(" X(JFNEW) =%14.6g\n", X[jfnew]); } } nsp1 = *nsets; *nsets -= 1; nf += 1; if1 -= 1; /* . If IFNEW .eq. NSETS+1 there is nothing more to do. * */ if (!(ifnew != nsp1)) goto L_20224; /* . Shift indices in Set S. */ for (is = ifnew; is <= *nsets; is++) { Index[is] = Index[is + 1]; } Index[if1] = jfnew; /* . Here the columns INDEX(IFNEW .. NSETS) are of lengths * . IFNEW+1 .. NSETS+1, respectively. We must eliminate the * . last element of each of these columns. * */ if (!(ifnew < m1)) goto L_20229; igaus1 = ifnew; igaus2 = min( *nsets - 1, m1 - 1 ); goto L_30020; L_20230: if (!(*nsets <= m1)) goto L_20232; /* . Use diag elt (NSETS,NSETS) to eliminate elt below it. * */ js = Index[*nsets]; if (A(js - 1,nsp1 - 1) != ZERO) { fac = -A(js - 1,nsp1 - 1)/A(js - 1,*nsets - 1); saxpy( n, fac, &A(0,*nsets - 1), lda, &A(0,nsp1 - 1), lda ); A(js - 1,nsp1 - 1) = ZERO; B[nsp1] += fac*B[*nsets]; Rt[nsp1] += fac*Rt[*nsets]; } goto L_20231; L_20232: npr021 = 20235; goto L_30021; L_20235: ; L_20231: goto L_20228; L_20229: if (!(ifnew == m1)) goto L_20236; npr021 = 20237; goto L_30021; L_20237: goto L_20228; L_20236: igiv1 = ifnew; igiv2 = *nsets; npr022 = 20238; goto L_30022; L_20238: ; L_20228: ; L_20224: goto L_20183; /* ------------------------------------------------------------------ * procedure( COMPLETE RETRIANGULARIZATION FROM ROW M1 ON DOWN ) * * . Here the columns INDEX(M1 .. NSETS) are of lengths * . M1+1 .. NSETS+1, respectively. We must eliminate the * . last element of each of these columns. * * . Find element of largest magnitude in row M1 among column * . indices INDEX(M1 .. NSETS). * */ L_30021: amax = ZERO; for (is = m1; is <= *nsets; is++) { js = Index[is]; if (fabsf( A(js - 1,m1 - 1) ) > amax) { imax = is; amax = fabsf( A(js - 1,m1 - 1) ); } } /* . Interchange columns indexed by M1 and IMAX. * */ if (imax != m1) { jpiv = Index[imax]; Index[imax] = Index[m1]; Index[m1] = jpiv; } /* . Do one stage of Gaussian elimination. * */ ig1 = m1; npr005 = 20246; goto L_30005; /* . Use Householder transformations to * . retriangularization in columns M1+1 through IMAX-1 * */ L_20246: if (!(m1p1 <= imax - 1)) goto L_20248; iput1 = m1p1; iput2 = imax - 1; irlast = imax + 1; goto L_30023; /* . Use Givens transformations to * . retriangularization in columns max(M1+1,IMAX) through NSETS * */ L_20249: ; L_20248: if (!(max( m1p1, imax ) <= *nsets)) goto L_20251; igiv1 = max( m1p1, imax ); igiv2 = *nsets; npr022 = 20252; goto L_30022; L_20252: ; L_20251: switch (npr021) { case 20235: goto L_20235; case 20237: goto L_20237; } /* ------------------------------------------------------------------ * procedure( GIVENS ROTATIONS, USING IGIV1 AND IGIV2 ) * * This proc applies Givens rotations to pairs of rows, IS and IS+1, * for IS = IGIV1, ..., IGIV2. * This version applies these to all columns from 1 to N, even though * in general this causes arithmetic to be done on elements that * are known to be zero. * * do for IS = IGIV1, IGIV2 * JS=INDEX(IS) * call SROTG (A(IS,JS),A(IS+1,JS),CC,SS) * A(IS+1,JS)=ZERO * LEN = JS-1 * if(LEN .gt. 0) * * call SROT(LEN,A(IS,1),LDA,A(IS+1,1),LDA,CC,SS) * LEN = N - JS * if(LEN .gt. 0) * * call SROT(LEN,A(IS,JS+1),LDA,A(IS+1,JS+1),LDA,CC,SS) * call SROT (1, B(IS),1, B(IS+1),1, CC,SS) * if(KPRINT .ge. 3) * * call SROT (1, RT(IS),1, RT(IS+1),1, CC,SS) * end for ! IS * end proc !( GIVENS ROTATIONS, USING IGIV1 AND IGIV2 ) * ------------------------------------------------------------------ */ /* procedure( GIVENS ROTATIONS, USING IGIV1 AND IGIV2 ) * * This proc applies Givens rotations to pairs of rows, IS and IS+1, * for IS = IGIV1, ..., IGIV2. * This version applies these only to the columns that might contain * nonzeros in the affected rows. * */ L_30022: for (l = igiv1; l <= if2; l++) { jcol = Index[l]; /* . Apply previously determined Givens transformations. * */ for (ir = igiv1; ir <= min( l - 1, igiv2 ); ir++) { temp = Cc[ir]*A(jcol - 1,ir - 1) + Ss[ir]*A(jcol - 1,ir); A(jcol - 1,ir) = -Ss[ir]*A(jcol - 1,ir - 1) + Cc[ir]*A(jcol - 1,ir); A(jcol - 1,ir - 1) = temp; } /* . Compute a new Givens transformation. * */ if (l <= igiv2) { srotg( &A(jcol - 1,l - 1), &A(jcol - 1,l), &Cc[l], &Ss[l] ); A(jcol - 1,l) = ZERO; } } /* . Apply previous transformations to B() * */ for (ir = igiv1; ir <= igiv2; ir++) { temp = Cc[ir]*B[ir] + Ss[ir]*B[ir + 1]; B[ir + 1] = -Ss[ir]*B[ir] + Cc[ir]*B[ir + 1]; B[ir] = temp; } /* . Apply previous transformations to RT() for code testing. * */ if (kprint >= 3) { for (ir = igiv1; ir <= igiv2; ir++) { temp = Cc[ir]*Rt[ir] + Ss[ir]*Rt[ir + 1]; Rt[ir + 1] = -Ss[ir]*Rt[ir] + Cc[ir]*Rt[ir + 1]; Rt[ir] = temp; } } switch (npr022) { case 20238: goto L_20238; case 20252: goto L_20252; } /* ----------------------------------------------------------------- */ /* procedure( GAUSSIAN ELIMINATION, USING IGAUS1 AND IGAUS2 ) * * This proc applies Gaussian elimination to pairs of rows, * IS and IS+1, for IS = IGAUS1 .. IGAUS2. * */ L_30020: ; /*++ CODE for ~.C. is inactive * if(KPRINT .ge. 2) print '(1x,a,2i4)', * * 'Gauss elim in pairs of rows. IGAUS1, IGAUS2 =',IGAUS1,IGAUS2 *++ CODE for .C. is active */ if (kprint >= 2) printf( " Gauss elim in pairs of rows. IGAUS1, IGAUS2 =%4ld %4ld\n", igaus1, igaus2); for (is = igaus1; is <= igaus2; is++) { /*++ END */ js = Index[is]; if (fabsf( A(js - 1,is) ) > fabsf( A(js - 1,is - 1) )) { sswap( n, &A(0,is - 1), lda, &A(0,is), lda ); temp = B[is]; B[is] = B[is + 1]; B[is + 1] = temp; temp = Rt[is]; Rt[is] = Rt[is + 1]; Rt[is + 1] = temp; } fac = -A(js - 1,is)/A(js - 1,is - 1); saxpy( n, fac, &A(0,is - 1), lda, &A(0,is), lda ); A(js - 1,is) = ZERO; B[is + 1] += fac*B[is]; Rt[is + 1] += fac*Rt[is]; } goto L_20230; /* ------------------------------------------------------------------ */ /* procedure( HOUSEHOLDER TRIANGULARIZATION ) * * Build new columns of triangular matrix indexed by IPUT1 thru IPUT2 * Last row to be affected is IRLAST. * */ L_30023: for (iput = iput1; iput <= iput2; iput++) { jput = Index[iput]; shtcc( 1, iput, iput + 1, irlast, &A(jput - 1,0), &up, b, m, 1 ); if (kprint >= 4) { printf(" Hous Tri.. IPUT,IRLAST,B(IPUT:IRLAST)=%4ld\n %14.6ld", iput, irlast); for (i = iput; i <= irlast; i++) { printf("%14.6g", B[i]); } printf("\n"); } for (i = iput + 1; i <= if2; i++) { shtcc( 2, iput, iput + 1, irlast, &A(jput - 1,0), &up, &A(Index[i] - 1,0), lda, 1 ); } if (kprint >= 3) shtcc( 2, iput, iput + 1, irlast, &A(jput - 1,0), &up, &Rt[1], m, 1 ); for (i = iput + 1; i <= irlast; i++) { A(jput - 1,i - 1) = ZERO; } } if (!(kprint >= 3)) goto L_20283; npr004 = 20283; goto L_30004; L_20283: goto L_20249; /* ------------------------------------------------------------------ */ /* procedure( PREPARE ROW NSP1 TO BE OBJECTIVE FUNCTION ) * * Save current B(NSP1) into BSAVE. * Set B(NSP1) to be the residual in row NSP1 for the current * X vector. * Do Gaussian elimination in row NSP1 using preceeding rows as * pivot rows. Don't need to apply the Gaussian elim in the B() * vector because all of B(1:NSETS) is zero. * */ L_30008: bsave = B[nsp1]; B[nsp1] -= sdot( n, &A(0,nsp1 - 1), lda, x, 1 ); for (is = 1; is <= *nsets; is++) { js = Index[is]; if (is > 1) { A(js - 1,nsp1 - 1) -= sdot( is - 1, &Z[1], 1, &A(js - 1,0), 1 ); } Z[is] = A(js - 1,nsp1 - 1)/A(js - 1,is - 1); A(js - 1,nsp1 - 1) = ZERO; } for (if0 = if1; if0 <= if2; if0++) { jf = Index[if0]; A(jf - 1,nsp1 - 1) -= sdot( *nsets, &Z[1], 1, &A(jf - 1,0), 1 ); } goto L_20065; /* ------------------------------------------------------------------ */ /* procedure( ELIMINATE IN LOWER ROWS ) * * The equality constraint rows 1 thru M1 have been triangularized. * Elimination in columns 1 thru IG1-1 in the lower matrix has been * previously done. Now eliminate in columns IG1 thru M1 of the * lower matrix, i.e. in rows M1+1 through M. * */ L_30005: ; if (kprint >= 3) { printf(" \n Gaussian elimination. IG1,M1=%4ld%4ld\n", ig1, m1); } for (is = ig1; is <= m1; is++) { js = Index[is]; len = 0; for (i = m; i >= m1p1; i--) { if (A(js - 1,i - 1) != ZERO) { len = i - m1; goto L_20296; } } L_20296: ; if (len != 0) { fac = -ONE/A(js - 1,is - 1); for (i = is + 1; i <= if2; i++) { j = Index[i]; temp = fac*A(j - 1,is - 1); saxpy( len, temp, &A(js - 1,m1p1 - 1), 1, &A(j - 1,m1p1 - 1), 1 ); } temp = fac*B[is]; saxpy( len, temp, &A(js - 1,m1p1 - 1), 1, &B[m1p1], 1 ); if (kprint >= 3) saxpy( len, fac*Rt[is], &A(js - 1,m1p1 - 1), 1, &Rt[m1p1], 1 ); for (i = m1p1; i <= (m1 + len); i++) { A(js - 1,i - 1) = ZERO; } } } if (!(kprint >= 3)) goto L_20308; npr004 = 20308; goto L_30004; L_20308: switch (npr005) { case 20018: goto L_20018; case 20246: goto L_20246; } /* ------------------------------------------------------------------ * procedure( SOLVE TRIANGULAR SYSTEM, INPUT AND OUTPUT IN Z() ) * * . Solve the triangular system defined by the column vectors * . indexed in Set S. The right-side vector is given in * . Z(), and the solution vector will be returned in Z(). * . It is assured that all diagonal terms are nonzero. * */ L_30014: for (iis = *nsets; iis >= 1; iis--) { if (iis != *nsets) saxpy( iis, -Z[iis + 1], &A(jjs - 1,0), 1, z, 1 ); jjs = Index[iis]; Z[iis] /= A(jjs - 1,iis - 1); } switch (npr014) { case 20115: goto L_20115; case 20126: goto L_20126; } /* ------------------------------------------------------------------ */ /* procedure( TERMINATION ) * * Compute the norm of the final residual vector. * Clean up W(). */ L_30007: *rnorm = ZERO; if (*ierr <= 0 || *ierr >= 3) { i1 = max( m1p1, nsp1 ); len = m - i1 + 1; if (len > ZERO) { *rnorm = snrm2( len, &B[i1], 1 ); } else { for (j = 1; j <= n; j++) { W[j] = ZERO; } } } goto L_20022; /* ------------------------------------------------------------------ */ /* procedure( DEBUG ) * */ L_30004: if (kprint >= 3) { if (cmode) { mdb = m1; } else { mdb = m; } scopy( mdb, b, 1, diff, 1 ); for (jd = 1; jd <= n; jd++) { saxpy( mdb, X[jd] - Xt[jd], &A(jd - 1,0), 1, diff, 1 ); } saxpy( mdb, -ONE, rt, 1, diff, 1 ); printf(" \n Consistency test. DIFF()=\n"); printf(" "); for (id = 1; id <= mdb; id++) { printf("%14.6g", Diff[id]); } printf("\n"); } switch (npr004) { case 20010: goto L_20010; case 20019: goto L_20019; case 20021: goto L_20021; case 20283: goto L_20283; case 20308: goto L_20308; } /* ------------------------------------------------------------------ */ #undef A } /* end of function */