diff --git a/include/cantera/equil/vcs_internal.h b/include/cantera/equil/vcs_internal.h index a2bb96a3b..da6ccbd08 100644 --- a/include/cantera/equil/vcs_internal.h +++ b/include/cantera/equil/vcs_internal.h @@ -17,10 +17,6 @@ #include "cantera/base/global.h" -#ifndef ALTLINPROG -#define ALTLINPROG -#endif - namespace VCSnonideal { using Cantera::npos; diff --git a/include/cantera/equil/vcs_solve.h b/include/cantera/equil/vcs_solve.h index 1bbc1befd..494707b2f 100644 --- a/include/cantera/equil/vcs_solve.h +++ b/include/cantera/equil/vcs_solve.h @@ -1058,7 +1058,6 @@ public: */ int vcs_inest_TP(); -#ifdef ALTLINPROG //! Estimate the initial mole numbers by constrained linear programming /*! * This is done by running @@ -1072,7 +1071,6 @@ public: * non-negativity constraints. */ int vcs_setMolesLinProg(); -#endif double vcs_Total_Gibbs(double* w, double* fe, double* tPhMoles); @@ -2072,11 +2070,6 @@ public: }; -#ifdef ALTLINPROG -#else -int linprogmax(double*, double*, double*, double*, size_t, size_t, size_t); -#endif - } #endif diff --git a/src/equil/Makefile.am b/src/equil/Makefile.am index b05aba20c..b56713c8e 100644 --- a/src/equil/Makefile.am +++ b/src/equil/Makefile.am @@ -26,8 +26,7 @@ cc_sources = BasisOptimize.cpp vcs_Gibbs.cpp vcs_root1d.cpp \ vcs_elem_rearrange.cpp vcs_prob.cpp vcs_inest.cpp \ vcs_equilibrate.cpp vcs_rearrange.cpp vcs_util.cpp \ vcs_Exception.cpp vcs_report.cpp vcs_VolPhase.cpp \ - vcs_dbocls.c vcs_dbolsm.c vcs_dbols.c vcs_dvout.c \ - vcs_ivout.c vcs_dmout.c vcs_xerror.c vcs_solve_phaseStability.cpp + vcs_solve_phaseStability.cpp AM_CPPFLAGS = -I$(top_builddir)/include -I$(top_builddir)/ext/f2c_libs/ AM_CXXFLAGS = $(AM_CPPFLAGS) diff --git a/src/equil/vcs_dbocls.c b/src/equil/vcs_dbocls.c deleted file mode 100644 index 33fdf9127..000000000 --- a/src/equil/vcs_dbocls.c +++ /dev/null @@ -1,1338 +0,0 @@ -/* dbocls.f -- translated by f2c (version 20031025). - You must link the resulting object file with libf2c: - on Microsoft Windows system, link with libf2c.lib; - on Linux or Unix systems, link with .../path/to/libf2c.a -lm - or, if you install libf2c.a in a standard place, with -lf2c -lm - -- in that order, at the end of the command line, as in - cc *.o -lf2c -lm - Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., - - http://www.netlib.org/f2c/libf2c.zip -*/ - -#include "f2c.h" - -/* Table of constant values */ - -static integer c__1 = 1; -static integer c__0 = 0; -static integer c__2 = 2; -static integer c__4 = 4; -static long int lc__4 = 4; - -/* DECK DBOCLS */ -/* Subroutine */ -int dbocls_(doublereal* w, integer* mdw, integer* mcon, - integer* mrows, integer* ncols, doublereal* bl, doublereal* bu, - integer* ind, integer* iopt, doublereal* x, doublereal* rnormc, - doublereal* rnorm, integer* mode, doublereal* rw, integer* iw) -{ - /* Initialized data */ - - static integer igo = 0; - - /* System generated locals */ - integer w_dim1, w_offset, i__1, i__2; - doublereal d__1, d__2; - - /* Local variables */ - static integer i__, j; - static doublereal t, t1, t2; - static integer ip, jp, lp; - static doublereal wt; - static integer llb; - static doublereal one; - static integer lds, iiw, liw, llx, irw, lrw; - extern doublereal ddot_(integer*, doublereal*, integer*, doublereal*, - integer*); - static integer idum, lbou, lmdw, lndw, mdwl, nerr; - static real rdum; - static integer lenx, lliw, mnew, jopt[5], lopt; - static doublereal zero; - static integer mopt, llrw, mout; - extern doublereal dnrm2_(integer*, doublereal*, integer*); - static real rdum2; - static integer icase; - extern /* Subroutine */ int dscal_(integer*, doublereal*, doublereal*, - integer*); - static integer modec; - static logical accum; - static integer nchar; - extern /* Subroutine */ int dbols_(doublereal*, integer*, integer*, - integer*, doublereal*, doublereal*, integer*, integer*, - doublereal*, doublereal*, integer*, doublereal*, integer*); - static integer level; - extern doublereal dasum_(integer*, doublereal*, integer*); - static doublereal anorm, cnorm; - static integer lboum; - extern /* Subroutine */ int dcopy_(integer*, doublereal*, integer*, - doublereal*, integer*); - static integer liopt; - extern doublereal d1mach_(long int *); - static integer locacc; - static logical checkl; - static integer iscale, locdim; - static logical filter; - static doublereal drelpr; - static logical pretri; - static integer inrows; - extern /* Subroutine */ int xerrwv_(char*, integer*, integer*, integer - *, integer*, integer*, integer*, integer*, real*, real*, - ftnlen); - - /* ***BEGIN PROLOGUE DBOCLS */ - /* ***DATE WRITTEN 821220 (YYMMDD) */ - /* ***REVISION DATE 870803 (YYMMDD) */ - /* ***CATEGORY NO. K1A2A,G2E,G2H1,G2H2 */ - /* ***KEYWORDS LIBRARY=SLATEC,TYPE=DOUBLE PRECISION(SBOCLS-S DBOCLS-D), */ - /* BOUNDS,CONSTRAINTS,INEQUALITY,LEAST SQUARES,LINEAR */ - /* ***AUTHOR HANSON, R. J., SNLA */ - /* ***PURPOSE Solve the bounded and constrained least squares */ - /* problem consisting of solving the equation */ - /* E*X = F (in the least squares sense) */ - /* subject to the linear constraints */ - /* C*X = Y. */ - /* ***DESCRIPTION */ - - /* **** Double Precision Version of SBOCLS **** */ - /* **** All INPUT and OUTPUT real variables are DOUBLE PRECISION **** */ - - /* This subprogram solves the bounded and constrained least squares */ - /* problem. The problem statement is: */ - - /* Solve E*X = F (least squares sense), subject to constraints */ - /* C*X=Y. */ - - /* In this formulation both X and Y are unknowns, and both may */ - /* have bounds on any of their components. This formulation */ - /* of the problem allows the user to have equality and inequality */ - /* constraints as well as simple bounds on the solution components. */ - - /* This constrained linear least squares subprogram solves E*X=F */ - /* subject to C*X=Y, where E is MROWS by NCOLS, C is MCON by NCOLS. */ - - /* The user must have dimension statements of the form */ - - /* DIMENSION W(MDW,NCOLS+MCON+1), BL(NCOLS+MCON), BU(NCOLS+MCON), */ - /* * X(2*(NCOLS+MCON)+2+NX), RW(6*NCOLS+5*MCON) */ - /* INTEGER IND(NCOLS+MCON), IOPT(17+NI), IW(2*(NCOLS+MCON)) */ - - /* (here NX=number of extra locations required for the options; NX=0 */ - /* if no options are in use. Also NI=number of extra locations */ - /* for options 1-9.) */ - - /* INPUT */ - /* ----- */ - - /* ------------------------- */ - /* W(MDW,*),MCON,MROWS,NCOLS */ - /* ------------------------- */ - /* The array W contains the (possibly null) matrix [C:*] followed by */ - /* [E:F]. This must be placed in W as follows: */ - /* [C : *] */ - /* W = [ ] */ - /* [E : F] */ - /* The (*) after C indicates that this data can be undefined. The */ - /* matrix [E:F] has MROWS rows and NCOLS+1 columns. The matrix C is */ - /* placed in the first MCON rows of W(*,*) while [E:F] */ - /* follows in rows MCON+1 through MCON+MROWS of W(*,*). The vector F */ - /* is placed in rows MCON+1 through MCON+MROWS, column NCOLS+1. The */ - /* values of MDW and NCOLS must be positive; the value of MCON must */ - /* be nonnegative. An exception to this occurs when using option 1 */ - /* for accumulation of blocks of equations. In that case MROWS is an */ - /* OUTPUT variable only, and the matrix data for [E:F] is placed in */ - /* W(*,*), one block of rows at a time. See IOPT(*) contents, option */ - /* number 1, for further details. The row dimension, MDW, of the */ - /* array W(*,*) must satisfy the inequality: */ - - /* If using option 1, */ - /* MDW .ge. MCON + max(max. number of */ - /* rows accumulated, NCOLS) + 1. */ - /* If using option 8, */ - /* MDW .ge. MCON + MROWS. */ - /* Else */ - /* MDW .ge. MCON + max(MROWS, NCOLS). */ - - /* Other values are errors, but this is checked only when using */ - /* option=2. The value of MROWS is an output parameter when */ - /* using option number 1 for accumulating large blocks of least */ - /* squares equations before solving the problem. */ - /* See IOPT(*) contents for details about option 1. */ - - /* ------------------ */ - /* BL(*),BU(*),IND(*) */ - /* ------------------ */ - /* These arrays contain the information about the bounds that the */ - /* solution values are to satisfy. The value of IND(J) tells the */ - /* type of bound and BL(J) and BU(J) give the explicit values for */ - /* the respective upper and lower bounds on the unknowns X and Y. */ - /* The first NVARS entries of IND(*), BL(*) and BU(*) specify */ - /* bounds on X; the next MCON entries specify bounds on Y. */ - - /* 1. For IND(J)=1, require X(J) .ge. BL(J); */ - /* IF J.gt.NCOLS, Y(J-NCOLS) .ge. BL(J). */ - /* (the value of BU(J) is not used.) */ - /* 2. For IND(J)=2, require X(J) .le. BU(J); */ - /* IF J.gt.NCOLS, Y(J-NCOLS) .le. BU(J). */ - /* (the value of BL(J) is not used.) */ - /* 3. For IND(J)=3, require X(J) .ge. BL(J) and */ - /* X(J) .le. BU(J); */ - /* IF J.gt.NCOLS, Y(J-NCOLS) .ge. BL(J) and */ - /* Y(J-NCOLS) .le. BU(J). */ - /* (to impose equality constraints have BL(J)=BU(J)= */ - /* constraining value.) */ - /* 4. For IND(J)=4, no bounds on X(J) or Y(J-NCOLS) are required. */ - /* (the values of BL(J) and BU(J) are not used.) */ - - /* Values other than 1,2,3 or 4 for IND(J) are errors. In the case */ - /* IND(J)=3 (upper and lower bounds) the condition BL(J) .gt. BU(J) */ - /* is an error. The values BL(J), BU(J), J .gt. NCOLS, will be */ - /* changed. Significant changes mean that the constraints are */ - /* infeasible. (Users must make this decision themselves.) */ - /* The new values for BL(J), BU(J), J .gt. NCOLS, define a */ - /* region such that the perturbed problem is feasible. If users */ - /* know that their problem is feasible, this step can be skipped */ - /* by using option number 8 described below. */ - /* See IOPT(*) description. */ - - - /* ------- */ - /* IOPT(*) */ - /* ------- */ - /* This is the array where the user can specify nonstandard options */ - /* for DBOCLS( ). Most of the time this feature can be ignored by */ - /* setting the input value IOPT(1)=99. Occasionally users may have */ - /* needs that require use of the following subprogram options. For */ - /* details about how to use the options see below: IOPT(*) CONTENTS. */ - - /* Option Number Brief Statement of Purpose */ - /* ------ ------ ----- --------- -- ------- */ - /* 1 Return to user for accumulation of blocks */ - /* of least squares equations. The values */ - /* of IOPT(*) are changed with this option. */ - /* The changes are updates to pointers for */ - /* placing the rows of equations into position */ - /* for processing. */ - /* 2 Check lengths of all arrays used in the */ - /* subprogram. */ - /* 3 Column scaling of the data matrix, [C]. */ - /* [E] */ - /* 4 User provides column scaling for matrix [C]. */ - /* [E] */ - /* 5 Provide option array to the low-level */ - /* subprogram SBOLS( ). */ - /* 6 Provide option array to the low-level */ - /* subprogram SBOLSM( ). */ - /* 7 Move the IOPT(*) processing pointer. */ - /* 8 Do not preprocess the constraints to */ - /* resolve infeasibilities. */ - /* 9 Do not pretriangularize the least squares matrix. */ - /* 99 No more options to change. */ - - /* ---- */ - /* X(*) */ - /* ---- */ - /* This array is used to pass data associated with options 4,5 and */ - /* 6. Ignore this parameter (on input) if no options are used. */ - /* Otherwise see below: IOPT(*) CONTENTS. */ - - - /* OUTPUT */ - /* ------ */ - - /* ----------------- */ - /* X(*),RNORMC,RNORM */ - /* ----------------- */ - /* The array X(*) contains a solution (if MODE .ge.0 or .eq.-22) for */ - /* the constrained least squares problem. The value RNORMC is the */ - /* minimum residual vector length for the constraints C*X - Y = 0. */ - /* The value RNORM is the minimum residual vector length for the */ - /* least squares equations. Normally RNORMC=0, but in the case of */ - /* inconsistent constraints this value will be nonzero. */ - /* The values of X are returned in the first NVARS entries of X(*). */ - /* The values of Y are returned in the last MCON entries of X(*). */ - - /* ---- */ - /* MODE */ - /* ---- */ - /* The sign of MODE determines whether the subprogram has completed */ - /* normally, or encountered an error condition or abnormal status. A */ - /* value of MODE .ge. 0 signifies that the subprogram has completed */ - /* normally. The value of mode (.ge. 0) is the number of variables */ - /* in an active status: not at a bound nor at the value zero, for */ - /* the case of free variables. A negative value of MODE will be one */ - /* of the cases (-57)-(-41), (-37)-(-22), (-19)-(-2). Values .lt. -1 */ - /* correspond to an abnormal completion of the subprogram. These */ - /* error messages are in groups for the subprograms DBOCLS(), */ - /* SBOLSM(), and SBOLS(). An approximate solution will be returned */ - /* to the user only when max. iterations is reached, MODE=-22. */ - - /* ----------- */ - /* RW(*),IW(*) */ - /* ----------- */ - /* These are working arrays. (normally the user can ignore the */ - /* contents of these arrays.) */ - - /* IOPT(*) CONTENTS */ - /* ------- -------- */ - /* The option array allows a user to modify some internal variables */ - /* in the subprogram without recompiling the source code. A central */ - /* goal of the initial software design was to do a good job for most */ - /* people. Thus the use of options will be restricted to a select */ - /* group of users. The processing of the option array proceeds as */ - /* follows: a pointer, here called LP, is initially set to the value */ - /* 1. At the pointer position the option number is extracted and */ - /* used for locating other information that allows for options to be */ - /* changed. The portion of the array IOPT(*) that is used for each */ - /* option is fixed; the user and the subprogram both know how many */ - /* locations are needed for each option. The value of LP is updated */ - /* for each option based on the amount of storage in IOPT(*) that is */ - /* required. A great deal of error checking is done by the */ - /* subprogram on the contents of the option array. Nevertheless it */ - /* is still possible to give the subprogram optional input that is */ - /* meaningless. For example option 4 uses the locations */ - /* X(NCOLS+IOFF),...,X(NCOLS+IOFF+NCOLS-1) for passing scaling data. */ - /* The user must manage the allocation of these locations. */ - - /* 1 */ - /* - */ - /* This option allows the user to solve problems with a large number */ - /* of rows compared to the number of variables. The idea is that the */ - /* subprogram returns to the user (perhaps many times) and receives */ - /* new least squares equations from the calling program unit. */ - /* Eventually the user signals "that's all" and a solution is then */ - /* computed. The value of MROWS is an output variable when this */ - /* option is used. Its value is always in the range 0 .le. MROWS */ - /* .le. NCOLS+1. It is the number of rows after the */ - /* triangularization of the entire set of equations. If LP is the */ - /* processing pointer for IOPT(*), the usage for the sequential */ - /* processing of blocks of equations is */ - - - /* IOPT(LP)=1 */ - /* Move block of equations to W(*,*) starting at */ - /* the first row of W(*,*). */ - /* IOPT(LP+3)=# of rows in the block; user defined */ - - /* The user now calls DBOCLS( ) in a loop. The value of IOPT(LP+1) */ - /* directs the user's action. The value of IOPT(LP+2) points to */ - /* where the subsequent rows are to be placed in W(*,*). Both of */ - /* these values are first defined in the subprogram. The user */ - /* changes the value of IOPT(LP+1) (to 2) as a signal that all of */ - /* the rows have been processed. */ - - - /* . 4) { - nerr = 56; - nchar = 46; - xerrwv_("DBOCLS(). FOR J=(I1), IND(J)=(I2) MUST BE 1-4.", & - nchar, &nerr, &level, &c__2, &j, &ind[j], &c__0, & - rdum, &rdum, (ftnlen)46); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L260; - } - /* L10: */ - } - - /* SEE THAT BOUNDS ARE CONSISTENT. */ - i__1 = *ncols + *mcon; - for (j = 1; j <= i__1; ++j) { - if (ind[j] == 3) { - if (bl[j] > bu[j]) { - nerr = 57; - nchar = 58; - rdum2 = (real) bl[j]; - rdum = (real) bu[j]; - xerrwv_("DBOCLS(). FOR J=(I1), BOUND BL(J)=(R1) IS .GT. " - "BU(J)=(R2).", &nchar, &nerr, &level, &c__1, &j, & - idum, &c__2, &rdum2, &rdum, (ftnlen)58); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L260; - } - } - /* L20: */ - } - /* END PROCEDURE */ - /* DO(PROCESS OPTION ARRAY) */ - /* PROCEDURE(PROCESS OPTION ARRAY) */ - zero = 0.; - one = 1.; - drelpr = d1mach_(&lc__4); - checkl = FALSE_; - filter = TRUE_; - /*lenx = (*ncols + *mcon << 1) + 2;*/ - lenx = ((*ncols + *mcon) << 1) + 2; - iscale = 1; - igo = 1; - accum = FALSE_; - pretri = TRUE_; - lopt = 0; - mopt = 0; - lp = 0; - lds = 0; - /* DO FOREVER */ -L30: - lp += lds; - ip = iopt[lp + 1]; - jp = abs(ip); - - /* TEST FOR NO MORE OPTIIONS TO CHANGE. */ - if (ip == 99) { - if (lopt == 0) { - lopt = -(lp + 2); - } - if (mopt == 0) { - mopt = -(abs(lopt) + 7); - } - if (lopt < 0) { - lbou = abs(lopt); - } else { - lbou = lopt - 15; - } - - /* SEND COL. SCALING TO DBOLS(). */ - iopt[lbou] = 4; - iopt[lbou + 1] = 1; - - /* PASS AN OPTION ARRAY FOR DBOLSM(). */ - iopt[lbou + 2] = 5; - - /* LOC. OF OPTION ARRAY FOR DBOLSM( ). */ - iopt[lbou + 3] = 8; - - /* SKIP TO START OF USER-GIVEN OPTION ARRAY FOR DBOLS(). */ - iopt[lbou + 4] = 6; - iopt[lbou + 6] = 99; - if (lopt > 0) { - iopt[lbou + 5] = lopt - lbou + 1; - } else { - iopt[lbou + 4] = -iopt[lbou + 4]; - } - if (mopt < 0) { - lboum = abs(mopt); - } else { - lboum = mopt - 8; - } - - /* CHANGE PRETRIANGULARIZATION FACTOR IN DBOLSM(). */ - iopt[lboum] = 5; - iopt[lboum + 1] = *ncols + *mcon + 1; - - /* PASS WEIGHT TO DBOLSM() FOR RANK TEST. */ - iopt[lboum + 2] = 6; - iopt[lboum + 3] = *ncols + *mcon + 2; - iopt[lboum + 4] = *mcon; - - /* SKIP TO USER-GIVEN OPTION ARRAY FOR DBOLSM( ). */ - iopt[lboum + 5] = 1; - iopt[lboum + 7] = 99; - if (mopt > 0) { - iopt[lboum + 6] = mopt - lboum + 1; - } else { - iopt[lboum + 5] = -iopt[lboum + 5]; - } - /* EXIT FOREVER */ - goto L50; - } else if (jp == 99) { - lds = 1; - /* CYCLE FOREVER */ - goto L50; - } else if (jp == 1) { - if (ip > 0) { - - /* SET UP DIRECTION FLAG LOCATION, ROW STACKING POINTER */ - /* LOCATION, AND LOCATION FOR NUMBER OF NEW ROWS. */ - locacc = lp + 2; - - /* IOPT(LOCACC-1)=OPTION NUMBER FOR SEQ. ACCUMULATION. */ - /* CONTENTS.. IOPT(LOCACC )=USER DIRECTION FLAG, 1 OR 2. */ - /* IOPT(LOCACC+1)=ROW STACKING POINTER. */ - /* IOPT(LOCACC+2)=NUMBER OF NEW ROWS TO PROCESS. */ - /* USER ACTION WITH THIS OPTION.. */ - /* (SET UP OPTION DATA FOR SEQ. ACCUMULATION IN IOPT(*).) */ - /* (MOVE BLOCK OF EQUATIONS INTO W(*,*) STARTING AT FIRST */ - /* ROW OF W(*,*) BELOW THE ROWS FOR THE CONSTRAINT MATRIX C. */ - /* SET IOPT(LOCACC+2)=NO. OF LEAST SQUARES EQUATIONS IN BLOCK. */ - /* LOOP */ - /* CALL DBOCLS() */ - - /* IF(IOPT(LOCACC) .EQ. 1) THEN */ - /* STACK EQUAS. INTO W(*,*), STARTING AT */ - /* ROW IOPT(LOCACC+1). */ - /* INTO W(*,*). */ - /* SET IOPT(LOCACC+2)=NO. OF EQUAS. */ - /* IF LAST BLOCK OF EQUAS., SET IOPT(LOCACC)=2. */ - /* ELSE IF IOPT(LOCACC) .EQ. 2) THEN */ - /* (PROCESS IS OVER. EXIT LOOP.) */ - /* ELSE */ - /* (ERROR CONDITION. SHOULD NOT HAPPEN.) */ - /* END IF */ - /* END LOOP */ - iopt[locacc + 1] = *mcon + 1; - accum = TRUE_; - iopt[locacc] = igo; - } - lds = 4; - /* CYCLE FOREVER */ - goto L30; - } else if (jp == 2) { - if (ip > 0) { - - /* GET ACTUAL LENGTHS OF ARRAYS FOR CHECKING AGAINST NEEDS. */ - locdim = lp + 2; - - /* LMDW.GE.MCON+MAX(MOUT,NCOLS), IF MCON.GT.0 .AND FILTER */ - /* LMDW.GE.MCON+MOUT, OTHERWISE */ - - /* LNDW.GE.NCOLS+MCON+1 */ - /* LLB .GE.NCOLS+MCON */ - /* LLX .GE.2*(NCOLS+MCON)+2+EXTRA REQD. IN OPTIONS. */ - /* LLRW.GE.6*NCOLS+5*MCON */ - /* LLIW.GE.2*(NCOLS+MCON) */ - /* LIOP.GE. AMOUNT REQD. FOR OPTION ARRAY. */ - lmdw = iopt[locdim]; - lndw = iopt[locdim + 1]; - llb = iopt[locdim + 2]; - llx = iopt[locdim + 3]; - llrw = iopt[locdim + 4]; - lliw = iopt[locdim + 5]; - liopt = iopt[locdim + 6]; - checkl = TRUE_; - } - lds = 8; - /* CYCLE FOREVER */ - goto L30; - - /* OPTION TO MODIFY THE COLUMN SCALING. */ - } else if (jp == 3) { - if (ip > 0) { - iscale = iopt[lp + 2]; - - /* SEE THAT ISCALE IS 1 THRU 3. */ - if (iscale < 1 || iscale > 3) { - nerr = 48; - nchar = 41; - xerrwv_("DBOCLS(). ISCALE OPTION=(I1) MUST BE 1-3.", & - nchar, &nerr, &level, &c__1, &iscale, &idum, & - c__0, &rdum, &rdum, (ftnlen)41); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L260; - } - } - lds = 2; - /* CYCLE FOREVER */ - goto L30; - - /* IN THIS OPTION THE USER HAS PROVIDED SCALING. THE */ - /* SCALE FACTORS FOR THE COLUMNS BEGIN IN X(NCOLS+IOPT(LP+2)). */ - } else if (jp == 4) { - if (ip > 0) { - iscale = 4; - if (iopt[lp + 2] <= 0) { - nerr = 49; - nchar = 86; - xerrwv_("DBOCLS(). OFFSET PAST X(NCOLS) (I1) FOR USER-PR" - "OVIDED COLUMN SCALING MUST BE POSITIVE.", &nchar, - &nerr, &level, &c__1, &iopt[lp + 2], &idum, &c__0, - &rdum, &rdum, (ftnlen)86); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L260; - } - dcopy_(ncols, &x[*ncols + iopt[lp + 2]], &c__1, &rw[1], &c__1) - ; - lenx += *ncols; - i__1 = *ncols; - for (j = 1; j <= i__1; ++j) { - if (rw[j] <= zero) { - nerr = 50; - nchar = 84; - rdum2 = (real) rw[j]; - xerrwv_("DBOCLS(). EACH PROVIDED COL. SCALE FACTOR M" - "UST BE POSITIVE. COMP. (I1) NOW = (R1).", & - nchar, &nerr, &level, &c__1, &j, &idum, &c__1, - &rdum2, &rdum, (ftnlen)84); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L260; - } - /* L40: */ - } - } - lds = 2; - /* CYCLE FOREVER */ - goto L30; - - /* IN THIS OPTION AN OPTION ARRAY IS PROVIDED TO DBOLS(). */ - } else if (jp == 5) { - if (ip > 0) { - lopt = iopt[lp + 2]; - } - lds = 2; - /* CYCLE FOREVER */ - goto L30; - - /* IN THIS OPTION AN OPTION ARRAY IS PROVIDED TO DBOLSM(). */ - } else if (jp == 6) { - if (ip > 0) { - mopt = iopt[lp + 2]; - } - lds = 2; - /* CYCLE FOREVER */ - goto L30; - - /* THIS OPTION USES THE NEXT LOC OF IOPT(*) AS A */ - /* POINTER VALUE TO SKIP TO NEXT. */ - } else if (jp == 7) { - if (ip > 0) { - lp = iopt[lp + 2] - 1; - lds = 0; - } else { - lds = 2; - } - /* CYCLE FOREVER */ - goto L30; - - /* THIS OPTION AVOIDS THE CONSTRAINT RESOLVING PHASE FOR */ - /* THE LINEAR CONSTRAINTS C*X=Y. */ - } else if (jp == 8) { - filter = !(ip > 0); - lds = 1; - /* CYCLE FOREVER */ - goto L30; - - /* THIS OPTION SUPPRESSES PRETIRANGULARIZATION OF THE LEAST */ - /* SQUARES EQATIONS. */ - } else if (jp == 9) { - pretri = !(ip > 0); - lds = 1; - /* CYCLE FOREVER */ - goto L30; - - /* NO VALID OPTION NUMBER WAS NOTED. THIS IS AN ERROR CONDITION. */ - } else { - nerr = 51; - nchar = 48; - rdum = (real) idum; - xerrwv_("DBOCLS(). THE OPTION NUMBER=(I1) IS NOT DEFINED.", & - nchar, &nerr, &level, &c__1, &jp, &idum, &c__0, &rdum, & - rdum, (ftnlen)48); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L260; - } - /* END FOREVER */ - /* END PROCEDURE */ -L50: - if (checkl) { - /* DO(CHECK LENGTHS OF ARRAYS) */ - /* PROCEDURE(CHECK LENGTHS OF ARRAYS) */ - - /* THIS FEATURE ALLOWS THE USER TO MAKE SURE THAT THE */ - /* ARRAYS ARE LONG ENOUGH FOR THE INTENDED PROBLEM SIZE AND USE. */ - if (filter && ! accum) { - mdwl = *mcon + max(*mrows,*ncols); - } else { - mdwl = *mcon + *ncols + 1; - } - if (lmdw < mdwl) { - nerr = 41; - nchar = 88; - xerrwv_("DBOCLS(). THE ROW DIMENSION OF W(,)=(I1) MUST BE .G" - "E. THE NUMBER OF EFFECTIVE ROWS=(I2).", &nchar, &nerr, - &level, &c__2, &lmdw, &mdwl, &c__0, &rdum, &rdum, ( - ftnlen)88); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L260; - } - if (lndw < *ncols + *mcon + 1) { - nerr = 42; - nchar = 75; - i__1 = *ncols + *mcon + 1; - xerrwv_("DBOCLS(). THE COLUMN DIMENSION OF W(,)=(I1) MUST BE" - " .GE. NCOLS+MCON+1=(I2).", &nchar, &nerr, &level, & - c__2, &lndw, &i__1, &c__0, &rdum, &rdum, (ftnlen)75); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L260; - } - if (llb < *ncols + *mcon) { - nerr = 43; - nchar = 94; - i__1 = *ncols + *mcon; - xerrwv_("DBOCLS(). THE DIMENSIONS OF THE ARRAYS BL(),BU(), A" - "ND IND()=(I1) MUST BE .GE. NCOLS+MCON=(I2).", &nchar, - &nerr, &level, &c__2, &llb, &i__1, &c__0, &rdum, & - rdum, (ftnlen)94); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L260; - } - if (llx < lenx) { - nerr = 44; - nchar = 71; - xerrwv_("DBOCLS(). THE DIMENSION OF X()=(I1) MUST BE .GE. TH" - "E REQD. LENGTH=(I2).", &nchar, &nerr, &level, &c__2, & - llx, &lenx, &c__0, &rdum, &rdum, (ftnlen)71); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L260; - } - if (llrw < *ncols * 6 + *mcon * 5) { - nerr = 45; - nchar = 70; - i__1 = *ncols * 6 + *mcon * 5; - xerrwv_("DBOCLS(). THE DIMENSION OF RW()=(I1) MUST BE .GE. 6" - "*NCOLS+5*MCON=(I2).", &nchar, &nerr, &level, &c__2, & - llrw, &i__1, &c__0, &rdum, &rdum, (ftnlen)70); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L260; - } - if (lliw < (*ncols << 1) + (*mcon << 1)) { - nerr = 46; - nchar = 69; - i__1 = (*ncols << 1) + (*mcon << 1); - xerrwv_("DBOCLS() THE DIMENSION OF IW()=(I1) MUST BE .GE. 2*" - "NCOLS+2*MCON=(I2).", &nchar, &nerr, &level, &c__2, & - lliw, &i__1, &c__0, &rdum, &rdum, (ftnlen)69); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L260; - } - if (liopt < lp + 17) { - nerr = 47; - nchar = 72; - i__1 = lp + 17; - xerrwv_("DBOCLS(). THE DIMENSION OF IOPT()=(I1) MUST BE .GE." - " THE REQD. LEN.=(I2).", &nchar, &nerr, &level, &c__2, - &liopt, &i__1, &c__0, &rdum, &rdum, (ftnlen)72); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L260; - } - /* END PROCEDURE */ - } - } - - /* OPTIONALLY GO BACK TO THE USER FOR ACCUMULATION OF LEAST SQUARES */ - /* EQUATIONS AND DIRECTIONS FOR PROCESSING THESE EQUATIONS. */ - /* DO(ACCUMULATE LEAST SQUARES EQUATIONS) */ - /* PROCEDURE(ACCUMULATE LEAST SQUARES EQUATIONS) */ - if (accum) { - *mrows = iopt[locacc + 1] - 1 - *mcon; - inrows = iopt[locacc + 2]; - mnew = *mrows + inrows; - if (mnew < 0 || mnew + *mcon > *mdw) { - nerr = 52; - nchar = 66; - i__1 = *mdw - *mcon; - xerrwv_("DBOCLS(). NO. OF ROWS=(I1) MUST BE .GE. 0 .AND. .LE.MDW" - "-MCON=(I2)", &nchar, &nerr, &level, &c__2, &mnew, &i__1, & - c__0, &rdum, &rdum, (ftnlen)65); - /* (RETURN TO USER PROGRAM UNIT) */ - goto L260; - } - } - - /* USE THE SOFTWARE OF DBOLS( ) FOR THE TRIANGULARIZATION OF THE */ - /* LEAST SQUARES MATRIX. THIS MAY INVOLVE A SYSTALTIC INTERCHANGE */ - /* OF PROCESSING POINTERS BETWEEN THE CALLING AND CALLED (DBOLS()) */ - /* PROGRAM UNITS. */ - jopt[0] = 1; - jopt[1] = 2; - jopt[3] = *mrows; - jopt[4] = 99; - irw = *ncols + 1; - iiw = 1; - if (accum || pretri) { - dbols_(&w[*mcon + 1 + w_dim1], mdw, &mout, ncols, &bl[1], &bu[1], & - ind[1], jopt, &x[1], rnorm, mode, &rw[irw], &iw[iiw]); - } else { - mout = *mrows; - } - if (accum) { - accum = iopt[locacc] == 1; - iopt[locacc + 1] = jopt[2] + *mcon; - /* Computing MIN */ - i__1 = *ncols + 1; - *mrows = min(i__1,mnew); - } - /* END PROCEDURE */ - if (accum) { - return 0; - } - /* DO(SOLVE CONSTRAINED AND BOUNDED LEAST SQUARES PROBLEM) */ - /* PROCEDURE(SOLVE CONSTRAINED AND BOUNDED LEAST SQUARES PROBLEM) */ - - /* MOVE RIGHT HAND SIDE OF LEAST SQUARES EQUATIONS. */ - dcopy_(&mout, &w[*mcon + 1 + (*ncols + 1) * w_dim1], &c__1, &w[*mcon + 1 - + (*ncols + *mcon + 1) * w_dim1], &c__1); - if (*mcon > 0 && filter) { - - /* PROJECT THE LINEAR CONSTRAINTS INTO A REACHABLE SET. */ - i__1 = *mcon; - for (i__ = 1; i__ <= i__1; ++i__) { - dcopy_(ncols, &w[i__ + w_dim1], mdw, &w[*mcon + 1 + (*ncols + i__) - * w_dim1], &c__1); - /* L60: */ - } - - /* PLACE (-)IDENTITY MATRIX AFTER CONSTRAINT DATA. */ - i__1 = *ncols + *mcon + 1; - for (j = *ncols + 1; j <= i__1; ++j) { - w[j * w_dim1 + 1] = zero; - dcopy_(mcon, &w[j * w_dim1 + 1], &c__0, &w[j * w_dim1 + 1], &c__1) - ; - /* L70: */ - } - w[(*ncols + 1) * w_dim1 + 1] = -one; - i__1 = *mdw + 1; - dcopy_(mcon, &w[(*ncols + 1) * w_dim1 + 1], &c__0, &w[(*ncols + 1) * - w_dim1 + 1], &i__1); - - /* OBTAIN A 'FEASIBLE POINT' FOR THE LINEAR CONSTRAINTS. */ - jopt[0] = 99; - irw = *ncols + 1; - iiw = 1; - i__1 = *ncols + *mcon; - dbols_(&w[w_offset], mdw, mcon, &i__1, &bl[1], &bu[1], &ind[1], jopt, - &x[1], rnormc, &modec, &rw[irw], &iw[iiw]); - - /* ENLARGE THE BOUNDS SET, IF REQUIRED, TO INCLUDE POINTS THAT */ - /* CAN BE REACHED. */ - i__1 = *ncols + *mcon; - for (j = *ncols + 1; j <= i__1; ++j) { - icase = ind[j]; - if (icase < 4) { - t = ddot_(ncols, &w[*mcon + 1 + j * w_dim1], &c__1, &x[1], & - c__1); - } - switch (icase) { - case 1: - goto L80; - case 2: - goto L90; - case 3: - goto L100; - case 4: - goto L110; - } - goto L120; - /* CASE 1 */ -L80: - /* Computing MIN */ - d__1 = t, d__2 = bl[j]; - bl[j] = min(d__1,d__2); - goto L120; - /* CASE 2 */ -L90: - /* Computing MAX */ - d__1 = t, d__2 = bu[j]; - bu[j] = max(d__1,d__2); - goto L120; - /* CASE 3 */ -L100: - /* Computing MIN */ - d__1 = t, d__2 = bl[j]; - bl[j] = min(d__1,d__2); - /* Computing MAX */ - d__1 = t, d__2 = bu[j]; - bu[j] = max(d__1,d__2); - goto L120; - /* CASE 4 */ -L110: -L120: - /* L130: */ - ; - } - - /* MOVE CONSTRAINT DATA BACK TO THE ORIGINAL AREA. */ - i__1 = *ncols + *mcon; - for (j = *ncols + 1; j <= i__1; ++j) { - dcopy_(ncols, &w[*mcon + 1 + j * w_dim1], &c__1, &w[j - *ncols + - w_dim1], mdw); - /* L140: */ - } - } - if (*mcon > 0) { - i__1 = *ncols + *mcon; - for (j = *ncols + 1; j <= i__1; ++j) { - w[*mcon + 1 + j * w_dim1] = zero; - dcopy_(&mout, &w[*mcon + 1 + j * w_dim1], &c__0, &w[*mcon + 1 + j - * w_dim1], &c__1); - /* L150: */ - } - - /* PUT IN (-)IDENTITY MATRIX (POSSIBLY) ONCE AGAIN. */ - i__1 = *ncols + *mcon + 1; - for (j = *ncols + 1; j <= i__1; ++j) { - w[j * w_dim1 + 1] = zero; - dcopy_(mcon, &w[j * w_dim1 + 1], &c__0, &w[j * w_dim1 + 1], &c__1) - ; - /* L160: */ - } - w[(*ncols + 1) * w_dim1 + 1] = -one; - i__1 = *mdw + 1; - dcopy_(mcon, &w[(*ncols + 1) * w_dim1 + 1], &c__0, &w[(*ncols + 1) * - w_dim1 + 1], &i__1); - } - - /* COMPUTE NOMINAL COLUMN SCALING FOR THE UNWEIGHTED MATRIX. */ - cnorm = zero; - anorm = zero; - i__1 = *ncols; - for (j = 1; j <= i__1; ++j) { - t1 = dasum_(mcon, &w[j * w_dim1 + 1], &c__1); - t2 = dasum_(&mout, &w[*mcon + 1 + w_dim1], &c__1); - t = t1 + t2; - if (t == zero) { - t = one; - } - cnorm = max(cnorm,t1); - anorm = max(anorm,t2); - x[*ncols + *mcon + j] = one / t; - /* L170: */ - } - switch (iscale) { - case 1: - goto L180; - case 2: - goto L190; - case 3: - goto L210; - case 4: - goto L220; - } - goto L230; - /* CASE 1 */ -L180: - goto L230; - /* CASE 2 */ - - /* SCALE COLS. (BEFORE WEIGHTING) TO HAVE LENGTH ONE. */ -L190: - i__1 = *ncols; - for (j = 1; j <= i__1; ++j) { - i__2 = *mcon + mout; - t = dnrm2_(&i__2, &w[j * w_dim1 + 1], &c__1); - if (t == zero) { - t = one; - } - x[*ncols + *mcon + j] = one / t; - /* L200: */ - } - goto L230; - /* CASE 3 */ - - /* SUPPRESS SCALING (USE UNIT MATRIX). */ -L210: - x[*ncols + *mcon + 1] = one; - dcopy_(ncols, &x[*ncols + *mcon + 1], &c__0, &x[*ncols + *mcon + 1], & - c__1); - goto L230; - /* CASE 4 */ - - /* THE USER HAS PROVIDED SCALING. */ -L220: - dcopy_(ncols, &rw[1], &c__1, &x[*ncols + *mcon + 1], &c__1); -L230: - i__1 = *ncols + *mcon; - for (j = *ncols + 1; j <= i__1; ++j) { - x[*ncols + *mcon + j] = one; - /* L240: */ - } - - /* WEIGHT THE LEAST SQUARES EQUATIONS. */ - wt = drelpr; - if (anorm > zero) { - wt /= anorm; - } - if (cnorm > zero) { - wt *= cnorm; - } - i__1 = mout; - for (i__ = 1; i__ <= i__1; ++i__) { - dscal_(ncols, &wt, &w[i__ + *mcon + w_dim1], mdw); - /* L250: */ - } - dscal_(&mout, &wt, &w[*mcon + 1 + (*mcon + *ncols + 1) * w_dim1], &c__1); - lrw = 1; - liw = 1; - - /* SET THE NEW TRIANGULARIZATION FACTOR. */ - /*x[(*ncols + *mcon << 1) + 1] = zero;*/ - x[((*ncols + *mcon) << 1) + 1] = zero; - - /* SET THE WEIGHT TO USE IN COMPONENTS .GT. MCON, */ - /* WHEN MAKING LINEAR INDEPENDENCE TEST. */ - /*x[(*ncols + *mcon << 1) + 2] = one / wt; */ - x[((*ncols + *mcon) << 1) + 2] = one / wt; - i__1 = mout + *mcon; - i__2 = *ncols + *mcon; - dbols_(&w[w_offset], mdw, &i__1, &i__2, &bl[1], &bu[1], &ind[1], &iopt[ - lbou], &x[1], rnorm, mode, &rw[lrw], &iw[liw]); - *rnorm /= wt; - /* END PROCEDURE */ - /* PROCEDURE(RETURN TO USER PROGRAM UNIT) */ -L260: - if (*mode >= 0) { - *mode = -nerr; - } - igo = 0; - return 0; - /* END PROGRAM */ -} /* dbocls_ */ - diff --git a/src/equil/vcs_dbols.c b/src/equil/vcs_dbols.c deleted file mode 100644 index a33f342b3..000000000 --- a/src/equil/vcs_dbols.c +++ /dev/null @@ -1,964 +0,0 @@ -/* dbols.f -- translated by f2c (version 20031025). - You must link the resulting object file with libf2c: - on Microsoft Windows system, link with libf2c.lib; - on Linux or Unix systems, link with .../path/to/libf2c.a -lm - or, if you install libf2c.a in a standard place, with -lf2c -lm - -- in that order, at the end of the command line, as in - cc *.o -lf2c -lm - Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., - - http://www.netlib.org/f2c/libf2c.zip -*/ - -#include "f2c.h" - -/* Table of constant values */ - -static integer c__1 = 1; -static integer c__0 = 0; -static integer c__2 = 2; - -/* DECK DBOLS */ -/* Subroutine */ -int dbols_(doublereal* w, integer* mdw, integer* mrows, - integer* ncols, doublereal* bl, doublereal* bu, integer* ind, integer - *iopt, doublereal* x, doublereal* rnorm, integer* mode, doublereal * - rw, integer* iw) -{ - /* Initialized data */ - - static integer igo = 0; - - /* System generated locals */ - integer w_dim1, w_offset, i__1, i__2, i__3; - doublereal d__1; - - /* Local variables */ - static integer i__, j; - static doublereal sc; - static integer ip, jp, lp; - static doublereal ss; - static integer llb; - static doublereal one; - static integer lds, llx, ibig, idum, lmdw, lndw, nerr; - static real rdum; - static integer lenx, lliw, mnew; - extern /* Subroutine */ int drot_(integer*, doublereal*, integer*, - doublereal*, integer*, doublereal*, doublereal*); - static integer lopt; - static doublereal zero; - static integer llrw; - extern doublereal dnrm2_(integer*, doublereal*, integer*); - static real rdum2; - static integer nchar, level; - extern /* Subroutine */ int dcopy_(integer*, doublereal*, integer*, - doublereal*, integer*), drotg_(doublereal*, doublereal*, - doublereal*, doublereal*); - static integer liopt, locacc; - static logical checkl; - static integer iscale; - extern integer idamax_(integer*, doublereal*, integer*); - static integer locdim; - extern /* Subroutine */ int dbolsm_(doublereal*, integer*, integer*, - integer*, doublereal*, doublereal*, integer*, integer*, - doublereal*, doublereal*, integer*, doublereal*, doublereal*, - doublereal*, integer*, integer*); - static integer inrows; - extern /* Subroutine */ int xerrwv_(char*, integer*, integer*, integer - *, integer*, integer*, integer*, integer*, real*, real*, - ftnlen); - - /* ***BEGIN PROLOGUE DBOLS */ - /* ***DATE WRITTEN 821220 (YYMMDD) */ - /* ***REVISION DATE 861211 (YYMMDD) */ - /* ***CATEGORY NO. K1A2A,G2E,G2H1,G2H2 */ - /* ***KEYWORDS LIBRARY=SLATEC,TYPE=DOUBLE PRECISION(SBOLS-S DBOLS-D), */ - /* BOUNDS,CONSTRAINTS,INEQUALITY,LEAST SQUARES,LINEAR */ - /* ***AUTHOR HANSON, R. J., SNLA */ - /* ***PURPOSE Solve the problem */ - /* E*X = F (in the least squares sense) */ - /* with bounds on selected X values. */ - /* ***DESCRIPTION */ - - /* **** Double Precision Version of SBOLS **** */ - /* **** All INPUT and OUTPUT real variables are DOUBLE PRECISION **** */ - - /* The user must have dimension statements of the form: */ - - /* DIMENSION W(MDW,NCOLS+1), BL(NCOLS), BU(NCOLS), */ - /* * X(NCOLS+NX), RW(5*NCOLS) */ - /* INTEGER IND(NCOLS), IOPT(1+NI), IW(2*NCOLS) */ - - /* (here NX=number of extra locations required for option 4; NX=0 */ - /* for no options; NX=NCOLS if this option is in use. Here NI=number */ - /* of extra locations required for options 1-6; NI=0 for no */ - /* options.) */ - - /* INPUT */ - /* ----- */ - - /* -------------------- */ - /* W(MDW,*),MROWS,NCOLS */ - /* -------------------- */ - /* The array W(*,*) contains the matrix [E:F] on entry. The matrix */ - /* [E:F] has MROWS rows and NCOLS+1 columns. This data is placed in */ - /* the array W(*,*) with E occupying the first NCOLS columns and the */ - /* right side vector F in column NCOLS+1. The row dimension, MDW, of */ - /* the array W(*,*) must satisfy the inequality MDW .ge. MROWS. */ - /* Other values of MDW are errrors. The values of MROWS and NCOLS */ - /* must be positive. Other values are errors. There is an exception */ - /* to this when using option 1 for accumulation of blocks of */ - /* equations. In that case MROWS is an OUTPUT variable ONLY, and the */ - /* matrix data for [E:F] is placed in W(*,*), one block of rows at a */ - /* time. MROWS contains the number of rows in the matrix after */ - /* triangularizing several blocks of equations. This is an OUTPUT */ - /* parameter ONLY when option 1 is used. See IOPT(*) CONTENTS */ - /* for details about option 1. */ - - /* ------------------ */ - /* BL(*),BU(*),IND(*) */ - /* ------------------ */ - /* These arrays contain the information about the bounds that the */ - /* solution values are to satisfy. The value of IND(J) tells the */ - /* type of bound and BL(J) and BU(J) give the explicit values for */ - /* the respective upper and lower bounds. */ - - /* 1. For IND(J)=1, require X(J) .ge. BL(J). */ - /* (the value of BU(J) is not used.) */ - /* 2. For IND(J)=2, require X(J) .le. BU(J). */ - /* (the value of BL(J) is not used.) */ - /* 3. For IND(J)=3, require X(J) .ge. BL(J) and */ - /* X(J) .le. BU(J). */ - /* 4. For IND(J)=4, no bounds on X(J) are required. */ - /* (the values of BL(J) and BU(J) are not used.) */ - - /* Values other than 1,2,3 or 4 for IND(J) are errors. In the case */ - /* IND(J)=3 (upper and lower bounds) the condition BL(J) .gt. BU(J) */ - /* is an error. */ - - /* ------- */ - /* IOPT(*) */ - /* ------- */ - /* This is the array where the user can specify nonstandard options */ - /* for DBOLSM( ). Most of the time this feature can be ignored by */ - /* setting the input value IOPT(1)=99. Occasionally users may have */ - /* needs that require use of the following subprogram options. For */ - /* details about how to use the options see below: IOPT(*) CONTENTS. */ - - /* Option Number Brief Statement of Purpose */ - /* ------ ------ ----- --------- -- ------- */ - /* 1 Return to user for accumulation of blocks */ - /* of least squares equations. */ - /* 2 Check lengths of all arrays used in the */ - /* subprogram. */ - /* 3 Standard scaling of the data matrix, E. */ - /* 4 User provides column scaling for matrix E. */ - /* 5 Provide option array to the low-level */ - /* subprogram DBOLSM( ). */ - /* 6 Move the IOPT(*) processing pointer. */ - /* 99 No more options to change. */ - - /* ---- */ - /* X(*) */ - /* ---- */ - /* This array is used to pass data associated with option 4. Ignore */ - /* this parameter if this option is not used. Otherwise see below: */ - /* IOPT(*) CONTENTS. */ - - /* OUTPUT */ - /* ------ */ - - /* ---------- */ - /* X(*),RNORM */ - /* ---------- */ - /* The array X(*) contains a solution (if MODE .ge.0 or .eq.-22) for */ - /* the constrained least squares problem. The value RNORM is the */ - /* minimum residual vector length. */ - - /* ---- */ - /* MODE */ - /* ---- */ - /* The sign of MODE determines whether the subprogram has completed */ - /* normally, or encountered an error condition or abnormal status. A */ - /* value of MODE .ge. 0 signifies that the subprogram has completed */ - /* normally. The value of MODE (.GE. 0) is the number of variables */ - /* in an active status: not at a bound nor at the value ZERO, for */ - /* the case of free variables. A negative value of MODE will be one */ - /* of the cases -37,-36,...,-22, or -17,...,-2. Values .lt. -1 */ - /* correspond to an abnormal completion of the subprogram. To */ - /* understand the abnormal completion codes see below: ERROR */ - /* MESSAGES for DBOLS( ). AN approximate solution will be returned */ - /* to the user only when max. iterations is reached, MODE=-22. */ - /* Values for MODE=-37,...,-22 come from the low-level subprogram */ - /* DBOLSM(). See the section ERROR MESSAGES for DBOLSM() in the */ - /* documentation for DBOLSM(). */ - - /* ----------- */ - /* RW(*),IW(*) */ - /* ----------- */ - /* These are working arrays with 5*NCOLS and 2*NCOLS entries. */ - /* (normally the user can ignore the contents of these arrays, */ - /* but they must be dimensioned properly.) */ - - /* IOPT(*) CONTENTS */ - /* ------- -------- */ - /* The option array allows a user to modify internal variables in */ - /* the subprogram without recompiling the source code. A central */ - /* goal of the initial software design was to do a good job for most */ - /* people. Thus the use of options will be restricted to a select */ - /* group of users. The processing of the option array proceeds as */ - /* follows: a pointer, here called LP, is initially set to the value */ - /* 1. This value is updated as each option is processed. At the */ - /* pointer position the option number is extracted and used for */ - /* locating other information that allows for options to be changed. */ - /* The portion of the array IOPT(*) that is used for each option is */ - /* fixed; the user and the subprogram both know how many locations */ - /* are needed for each option. A great deal of error checking is */ - /* done by the subprogram on the contents of the option array. */ - /* Nevertheless it is still possible to give the subprogram optional */ - /* input that is meaningless. For example option 4 uses the */ - /* locations X(NCOLS+IOFF),...,X(NCOLS+IOFF+NCOLS-1) for passing */ - /* scaling data. The user must manage the allocation of these */ - /* locations. */ - - /* 1 */ - /* - */ - /* This option allows the user to solve problems with a large number */ - /* of rows compared to the number of variables. The idea is that the */ - /* subprogram returns to the user (perhaps many times) and receives */ - /* new least squares equations from the calling program unit. */ - /* Eventually the user signals "that's all" and then computes the */ - /* solution with one final call to subprogram DBOLS( ). The value of */ - /* MROWS is an OUTPUT variable when this option is used. Its value */ - /* is always in the range 0 .le. MROWS .le. NCOLS+1. It is equal to */ - /* the number of rows after the triangularization of the entire set */ - /* of equations. If LP is the processing pointer for IOPT(*), the */ - /* usage for the sequential processing of blocks of equations is */ - - /* IOPT(LP)=1 */ - /* Move block of equations to W(*,*) starting at */ - /* the first row of W(*,*). */ - /* IOPT(LP+3)=# of rows in the block; user defined */ - - /* The user now calls DBOLS( ) in a loop. The value of IOPT(LP+1) */ - /* directs the user's action. The value of IOPT(LP+2) points to */ - /* where the subsequent rows are to be placed in W(*,*). */ - - /* . 4) { - nerr = 4; - nchar = 45; - xerrwv_("DBOLS(). FOR J=(I1), IND(J)=(I2) MUST BE 1-4.", & - nchar, &nerr, &level, &c__2, &j, &ind[j], &c__0, & - rdum, &rdum, (ftnlen)45); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L190; - } - /* L10: */ - } - - /* SEE THAT BOUNDS ARE CONSISTENT. */ - i__1 = *ncols; - for (j = 1; j <= i__1; ++j) { - if (ind[j] == 3) { - if (bl[j] > bu[j]) { - nerr = 5; - nchar = 57; - rdum2 = (real) bl[j]; - rdum = (real) bu[j]; - xerrwv_("DBOLS(). FOR J=(I1), BOUND BL(J)=(R1) IS .GT. B" - "U(J)=(R2).", &nchar, &nerr, &level, &c__1, &j, & - idum, &c__2, &rdum2, &rdum, (ftnlen)57); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L190; - } - } - /* L20: */ - } - /* END PROCEDURE */ - /* DO(PROCESS OPTION ARRAY) */ - /* PROCEDURE(PROCESS OPTION ARRAY) */ - zero = 0.; - one = 1.; - checkl = FALSE_; - lenx = *ncols; - iscale = 1; - igo = 2; - lopt = 0; - lp = 0; - lds = 0; -L30: - lp += lds; - ip = iopt[lp + 1]; - jp = abs(ip); - - /* TEST FOR NO MORE OPTIONS. */ - if (ip == 99) { - if (lopt == 0) { - lopt = lp + 1; - } - goto L50; - } else if (jp == 99) { - lds = 1; - goto L30; - } else if (jp == 1) { - if (ip > 0) { - - /* SET UP DIRECTION FLAG, ROW STACKING POINTER */ - /* LOCATION, AND LOCATION FOR NUMBER OF NEW ROWS. */ - locacc = lp + 2; - - /* IOPT(LOCACC-1)=OPTION NUMBER FOR SEQ. ACCUMULATION. */ - /* CONTENTS.. IOPT(LOCACC )=USER DIRECTION FLAG, 1 OR 2. */ - /* IOPT(LOCACC+1)=ROW STACKING POINTER. */ - /* IOPT(LOCACC+2)=NUMBER OF NEW ROWS TO PROCESS. */ - /* USER ACTION WITH THIS OPTION.. */ - /* (SET UP OPTION DATA FOR SEQ. ACCUMULATION IN IOPT(*). */ - /* MUST ALSO START PROCESS WITH IOPT(LOCACC)=1.) */ - /* (MOVE BLOCK OF EQUATIONS INTO W(*,*) STARTING AT FIRST */ - /* ROW OF W(*,*). SET IOPT(LOCACC+2)=NO. OF ROWS IN BLOCK.) */ - /* LOOP */ - /* CALL DBOLS() */ - - /* IF(IOPT(LOCACC) .EQ. 1) THEN */ - /* STACK EQUAS., STARTING AT ROW IOPT(LOCACC+1), */ - /* INTO W(*,*). */ - /* SET IOPT(LOCACC+2)=NO. OF EQUAS. */ - /* IF LAST BLOCK OF EQUAS., SET IOPT(LOCACC)=2. */ - /* ELSE IF IOPT(LOCACC) .EQ. 2) THEN */ - /* (PROCESS IS OVER. EXIT LOOP.) */ - /* ELSE */ - /* (ERROR CONDITION. SHOULD NOT HAPPEN.) */ - /* END IF */ - /* END LOOP */ - /* SET IOPT(LOCACC-1)=-OPTION NUMBER FOR SEQ. ACCUMULATION. */ - /* CALL DBOLS( ) */ - iopt[locacc + 1] = 1; - igo = 1; - } - lds = 4; - goto L30; - } else if (jp == 2) { - if (ip > 0) { - - /* GET ACTUAL LENGTHS OF ARRAYS FOR CHECKING AGAINST NEEDS. */ - locdim = lp + 2; - - /* LMDW.GE.MROWS */ - /* LNDW.GE.NCOLS+1 */ - /* LLB .GE.NCOLS */ - /* LLX .GE.NCOLS+EXTRA REQD. IN OPTIONS. */ - /* LLRW.GE.5*NCOLS */ - /* LLIW.GE.2*NCOLS */ - /* LIOP.GE. AMOUNT REQD. FOR IOPTION ARRAY. */ - lmdw = iopt[locdim]; - lndw = iopt[locdim + 1]; - llb = iopt[locdim + 2]; - llx = iopt[locdim + 3]; - llrw = iopt[locdim + 4]; - lliw = iopt[locdim + 5]; - liopt = iopt[locdim + 6]; - checkl = TRUE_; - } - lds = 8; - goto L30; - - /* OPTION TO MODIFY THE COLUMN SCALING. */ - } else if (jp == 3) { - if (ip > 0) { - iscale = iopt[lp + 2]; - - /* SEE THAT ISCALE IS 1 THRU 3. */ - if (iscale < 1 || iscale > 3) { - nerr = 7; - nchar = 40; - xerrwv_("DBOLS(). ISCALE OPTION=(I1) MUST BE 1-3.", & - nchar, &nerr, &level, &c__1, &iscale, &idum, & - c__0, &rdum, &rdum, (ftnlen)40); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L190; - } - } - lds = 2; - /* CYCLE FOREVER */ - goto L30; - - /* IN THIS OPTION THE USER HAS PROVIDED SCALING. THE */ - /* SCALE FACTORS FOR THE COLUMNS BEGIN IN X(NCOLS+IOPT(LP+2)). */ - } else if (jp == 4) { - if (ip > 0) { - iscale = 4; - if (iopt[lp + 2] <= 0) { - nerr = 8; - nchar = 85; - xerrwv_("DBOLS(). OFFSET PAST X(NCOLS) (I1) FOR USER-PRO" - "VIDED COLUMN SCALING MUST BE POSITIVE.", &nchar, & - nerr, &level, &c__1, &iopt[lp + 2], &idum, &c__0, - &rdum, &rdum, (ftnlen)85); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L190; - } - dcopy_(ncols, &x[*ncols + iopt[lp + 2]], &c__1, &rw[1], &c__1) - ; - lenx += *ncols; - i__1 = *ncols; - for (j = 1; j <= i__1; ++j) { - if (rw[j] <= zero) { - nerr = 9; - nchar = 85; - rdum2 = (real) rw[j]; - xerrwv_("DBOLS(). EACH PROVIDED COL. SCALE FACTOR MU" - "ST BE POSITIVE. COMPONENT (I1) NOW = (R1).", & - nchar, &nerr, &level, &c__1, &j, &idum, &c__1, - &rdum2, &rdum, (ftnlen)85); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L190; - } - /* L40: */ - } - } - lds = 2; - /* CYCLE FOREVER */ - goto L30; - - /* IN THIS OPTION AN OPTION ARRAY IS PROVIDED TO DBOLSM(). */ - } else if (jp == 5) { - if (ip > 0) { - lopt = iopt[lp + 2]; - } - lds = 2; - /* CYCLE FOREVER */ - goto L30; - - /* THIS OPTION USES THE NEXT LOC OF IOPT(*) AS AN */ - /* INCREMENT TO SKIP. */ - } else if (jp == 6) { - if (ip > 0) { - lp = iopt[lp + 2] - 1; - lds = 0; - } else { - lds = 2; - } - /* CYCLE FOREVER */ - goto L30; - - /* NO VALID OPTION NUMBER WAS NOTED. THIS IS AN ERROR CONDITION. */ - } else { - nerr = 6; - nchar = 47; - rdum2 = (real) idum; - xerrwv_("DBOLS(). THE OPTION NUMBER=(I1) IS NOT DEFINED.", &nchar, - &nerr, &level, &c__1, &jp, &idum, &c__0, &rdum2, &rdum2, - (ftnlen)47); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L190; - } -L50: - /* END PROCEDURE */ - if (checkl) { - /* DO(CHECK LENGTHS OF ARRAYS) */ - /* PROCEDURE(CHECK LENGTHS OF ARRAYS) */ - - /* THIS FEATURE ALLOWS THE USER TO MAKE SURE THAT THE */ - /* ARRAYS ARE LONG ENOUGH FOR THE INTENDED PROBLEM SIZE AND USE. */ - if (lmdw < *mrows) { - nerr = 11; - nchar = 76; - xerrwv_("DBOLS(). THE ROW DIMENSION OF W(,)=(I1) MUST BE .GE" - ".THE NUMBER OF ROWS=(I2).", &nchar, &nerr, &level, & - c__2, &lmdw, mrows, &c__0, &rdum, &rdum, (ftnlen)76); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L190; - } - if (lndw < *ncols + 1) { - nerr = 12; - nchar = 69; - i__1 = *ncols + 1; - xerrwv_("DBOLS(). THE COLUMN DIMENSION OF W(,)=(I1) MUST BE " - ".GE. NCOLS+1=(I2).", &nchar, &nerr, &level, &c__2, & - lndw, &i__1, &c__0, &rdum, &rdum, (ftnlen)69); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L190; - } - if (llb < *ncols) { - nerr = 13; - nchar = 88; - xerrwv_("DBOLS(). THE DIMENSIONS OF THE ARRAYS BL(),BU(), AN" - "D IND()=(I1) MUST BE .GE. NCOLS=(I2).", &nchar, &nerr, - &level, &c__2, &llb, ncols, &c__0, &rdum, &rdum, ( - ftnlen)88); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L190; - } - if (llx < lenx) { - nerr = 14; - nchar = 70; - xerrwv_("DBOLS(). THE DIMENSION OF X()=(I1) MUST BE .GE. THE" - " REQD. LENGTH=(I2).", &nchar, &nerr, &level, &c__2, & - llx, &lenx, &c__0, &rdum, &rdum, (ftnlen)70); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L190; - } - if (llrw < *ncols * 5) { - nerr = 15; - nchar = 62; - i__1 = *ncols * 5; - xerrwv_("DBOLS(). THE DIMENSION OF RW()=(I1) MUST BE .GE. 5*" - "NCOLS=(I2).", &nchar, &nerr, &level, &c__2, &llrw, & - i__1, &c__0, &rdum, &rdum, (ftnlen)62); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L190; - } - if (lliw < *ncols << 1) { - nerr = 16; - nchar = 61; - i__1 = *ncols << 1; - xerrwv_("DBOLS() THE DIMENSION OF IW()=(I1) MUST BE .GE. 2*N" - "COLS=(I2).", &nchar, &nerr, &level, &c__2, &lliw, & - i__1, &c__0, &rdum, &rdum, (ftnlen)61); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L190; - } - if (liopt < lp + 1) { - nerr = 17; - nchar = 71; - i__1 = lp + 1; - xerrwv_("DBOLS(). THE DIMENSION OF IOPT()=(I1) MUST BE .GE. " - "THE REQD. LEN.=(I2).", &nchar, &nerr, &level, &c__2, & - liopt, &i__1, &c__0, &rdum, &rdum, (ftnlen)71); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L190; - } - /* END PROCEDURE */ - } - } - switch (igo) { - case 1: - goto L60; - case 2: - goto L90; - } - goto L180; - - /* GO BACK TO THE USER FOR ACCUMULATION OF LEAST SQUARES */ - /* EQUATIONS AND DIRECTIONS TO QUIT PROCESSING. */ - /* CASE 1 */ -L60: - /* DO(ACCUMULATE LEAST SQUARES EQUATIONS) */ - /* PROCEDURE(ACCUMULATE LEAST SQUARES EQUATIONS) */ - *mrows = iopt[locacc + 1] - 1; - inrows = iopt[locacc + 2]; - mnew = *mrows + inrows; - if (mnew < 0 || mnew > *mdw) { - nerr = 10; - nchar = 61; - xerrwv_("DBOLS(). NO. OF ROWS=(I1) MUST BE .GE. 0 .AND. .LE. MDW=(I2" - ").", &nchar, &nerr, &level, &c__2, &mnew, mdw, &c__0, &rdum, & - rdum, (ftnlen)61); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L190; - } - /* Computing MIN */ - i__2 = *ncols + 1; - i__1 = min(i__2,mnew); - for (j = 1; j <= i__1; ++j) { - i__2 = max(*mrows,j) + 1; - for (i__ = mnew; i__ >= i__2; --i__) { - i__3 = i__ - j; - ibig = idamax_(&i__3, &w[j + j * w_dim1], &c__1) + j - 1; - - /* PIVOT FOR INCREASED STABILITY. */ - drotg_(&w[ibig + j * w_dim1], &w[i__ + j * w_dim1], &sc, &ss); - i__3 = *ncols + 1 - j; - drot_(&i__3, &w[ibig + (j + 1) * w_dim1], mdw, &w[i__ + (j + 1) * - w_dim1], mdw, &sc, &ss); - w[i__ + j * w_dim1] = zero; - /* L70: */ - } - /* L80: */ - } - /* Computing MIN */ - i__1 = *ncols + 1; - *mrows = min(i__1,mnew); - iopt[locacc + 1] = *mrows + 1; - igo = iopt[locacc]; - /* END PROCEDURE */ - if (igo == 2) { - igo = 0; - } - goto L180; - /* CASE 2 */ -L90: - /* DO(INITIALIZE VARIABLES AND DATA VALUES) */ - /* PROCEDURE(INITIALIZE VARIABLES AND DATA VALUES) */ - i__1 = *ncols; - for (j = 1; j <= i__1; ++j) { - switch (iscale) { - case 1: - goto L100; - case 2: - goto L110; - case 3: - goto L120; - case 4: - goto L130; - } - goto L140; -L100: - /* CASE 1 */ - - /* THIS IS THE NOMINAL SCALING. EACH NONZERO */ - /* COL. HAS MAX. NORM EQUAL TO ONE. */ - ibig = idamax_(mrows, &w[j * w_dim1 + 1], &c__1); - rw[j] = (d__1 = w[ibig + j * w_dim1], abs(d__1)); - if (rw[j] == zero) { - rw[j] = one; - } else { - rw[j] = one / rw[j]; - } - goto L140; -L110: - /* CASE 2 */ - - /* THIS CHOICE OF SCALING MAKES EACH NONZERO COLUMN */ - /* HAVE EUCLIDEAN LENGTH EQUAL TO ONE. */ - rw[j] = dnrm2_(mrows, &w[j * w_dim1 + 1], &c__1); - if (rw[j] == zero) { - rw[j] = one; - } else { - rw[j] = one / rw[j]; - } - goto L140; -L120: - /* CASE 3 */ - - /* THIS CASE EFFECTIVELY SUPPRESSES SCALING BY SETTING */ - /* THE SCALING MATRIX TO THE IDENTITY MATRIX. */ - rw[1] = one; - dcopy_(ncols, &rw[1], &c__0, &rw[1], &c__1); - goto L160; -L130: - /* CASE 4 */ - goto L160; -L140: - /* L150: */ - ; - } -L160: - /* END PROCEDURE */ - /* DO(SOLVE BOUNDED LEAST SQUARES PROBLEM) */ - /* PROCEDURE(SOLVE BOUNDED LEAST SQUARES PROBLEM) */ - - /* INITIALIZE IBASIS(*), J=1,NCOLS, AND IBB(*), J=1,NCOLS, */ - /* TO =J,AND =1, FOR USE IN DBOLSM( ). */ - i__1 = *ncols; - for (j = 1; j <= i__1; ++j) { - iw[j] = j; - iw[j + *ncols] = 1; - rw[*ncols * 3 + j] = bl[j]; - rw[(*ncols << 2) + j] = bu[j]; - /* L170: */ - } - dbolsm_(&w[w_offset], mdw, mrows, ncols, &rw[*ncols * 3 + 1], &rw[(*ncols - << 2) + 1], &ind[1], &iopt[lopt], &x[1], rnorm, mode, &rw[*ncols - + 1], &rw[(*ncols << 1) + 1], &rw[1], &iw[1], &iw[*ncols + 1]); - /* END PROCEDURE */ - igo = 0; -L180: - return 0; - /* PROCEDURE(RETURN TO USER PROGRAM UNIT) */ -L190: - if (*mode >= 0) { - *mode = -nerr; - } - igo = 0; - return 0; - /* END PROCEDURE */ -} /* dbols_ */ - diff --git a/src/equil/vcs_dbolsm.c b/src/equil/vcs_dbolsm.c deleted file mode 100644 index 5e8f03033..000000000 --- a/src/equil/vcs_dbolsm.c +++ /dev/null @@ -1,1586 +0,0 @@ -/* dbolsm.f -- translated by f2c (version 20031025). - You must link the resulting object file with libf2c: - on Microsoft Windows system, link with libf2c.lib; - on Linux or Unix systems, link with .../path/to/libf2c.a -lm - or, if you install libf2c.a in a standard place, with -lf2c -lm - -- in that order, at the end of the command line, as in - cc *.o -lf2c -lm - Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., - - http://www.netlib.org/f2c/libf2c.zip -*/ - -#include "f2c.h" - -/* Table of constant values */ - -static integer c__1 = 1; -static integer c__0 = 0; -static integer c__2 = 2; -static long int lc__2 = 2; -static integer c_n4 = -4; -static integer c__4 = 4; -static long int lc__4 = 4; - -/* DECK DBOLSM */ - -/* 4/16/91 Harry K. Moffat Div. 1126 844-6912 */ -/* Fixed an error! */ -/* DC and DS were not declared double precision */ - -/* Subroutine */ -int dbolsm_(doublereal* w, integer* mdw, integer* minput, - integer* ncols, doublereal* bl, doublereal* bu, integer* ind, integer - *iopt, doublereal* x, doublereal* rnorm, integer* mode, doublereal * - rw, doublereal* ww, doublereal* scl, integer* ibasis, integer* ibb) -{ - /* System generated locals */ - integer w_dim1, w_offset, i__1, i__2, i__3, i__4; - doublereal d__1, d__2; - - /* Builtin functions */ - double sqrt(doublereal); - - /* Local variables */ - static integer i__, j; - static doublereal t, t1, t2, dc, sc, ds; - static integer ip, jp, lp; - static doublereal ss, wt, cl1, cl2, cl3, fac, big, one; - static integer lds; - static doublereal bou, two, beta; - static integer jbig, jmag, ioff, jcol; - static doublereal wbig; - extern doublereal ddot_(integer*, doublereal*, integer*, doublereal*, - integer*); - static doublereal wmag; - static integer idum, mval, iter, nerr; - static real rdum; - extern /* Subroutine */ int drot_(integer*, doublereal*, integer*, - doublereal*, integer*, doublereal*, doublereal*); - static doublereal zero, xnew; - extern doublereal dnrm2_(integer*, doublereal*, integer*); - static real rdum2; - static integer icase; - static doublereal alpha; - static integer nchar, level; - static logical found; - static integer nsetb; - extern /* Subroutine */ int dcopy_(integer*, doublereal*, integer*, - doublereal*, integer*); - static integer igopr, itemp, itmax; - extern /* Subroutine */ int drotg_(doublereal*, doublereal*, doublereal - *, doublereal*), dswap_(integer*, doublereal*, integer*, - doublereal*, integer*), daxpy_(integer*, doublereal*, - doublereal*, integer*, doublereal*, integer*); - static integer lgopr; - extern /* Subroutine */ int dmout_(integer*, integer*, integer*, - doublereal*, char*, integer*, ftnlen); - static integer jdrop; - extern doublereal d1mach_(long int*); - extern /* Subroutine */ int dvout_(integer*, doublereal*, char*, - integer*, ftnlen), ivout_(integer*, integer*, char*, integer * - , ftnlen); - static integer mrows, jdrop1, jdrop2, jlarge; - static doublereal colabv, colblo, wlarge; - static integer nlevel; - static doublereal tolind; - static integer iprint; - static logical constr; - static doublereal tolsze; - extern /* Subroutine */ int xerrwv_(char*, integer*, integer*, integer - *, integer*, integer*, integer*, integer*, real*, real*, - ftnlen); - - /* ***BEGIN PROLOGUE DBOLSM */ - /* ***REFER TO SBOCLS,SBOLS */ - /* ***ROUTINES CALLED D1MACH,DAXPY,DCOPY,DDOT,DMOUT,DNRM2,DROT,DROTG, */ - /* DSWAP,DVOUT,IVOUT,XERRWV */ - /* ***DESCRIPTION */ - - /* **** Double Precision Version of SBOLSM **** */ - /* **** All INPUT and OUTPUT real variables are DOUBLE PRECISION **** */ - - /* Solve E*X = F (least squares sense) with bounds on */ - /* selected X values. */ - /* The user must have dimension statements of the form: */ - - /* DIMENSION W(MDW,NCOLS+1), BL(NCOLS), BU(NCOLS), */ - /* * X(NCOLS+NX), RW(NCOLS), WW(NCOLS), SCL(NCOLS) */ - /* INTEGER IND(NCOLS), IOPT(1+NI), IBASIS(NCOLS), IBB(NCOLS) */ - - /* (here NX=number of extra locations required for options 1,...,7; */ - /* NX=0 for no options; here NI=number of extra locations possibly */ - /* required for options 1-7; NI=0 for no options; NI=14 if all the */ - /* options are simultaneously in use.) */ - - /* INPUT */ - /* ----- */ - - /* -------------------- */ - /* W(MDW,*),MROWS,NCOLS */ - /* -------------------- */ - /* The array w(*,*) contains the matrix [E:F] on entry. The matrix */ - /* [E:F] has MROWS rows and NCOLS+1 columns. This data is placed in */ - /* the array W(*,*) with E occupying the first NCOLS columns and the */ - /* right side vector F in column NCOLS+1. The row dimension, MDW, of */ - /* the array W(*,*) must satisfy the inequality MDW .ge. MROWS. */ - /* Other values of MDW are errors. The values of MROWS and NCOLS */ - /* must be positive. Other values are errors. */ - - /* ------------------ */ - /* BL(*),BU(*),IND(*) */ - /* ------------------ */ - /* These arrays contain the information about the bounds that the */ - /* solution values are to satisfy. The value of IND(J) tells the */ - /* type of bound and BL(J) and BU(J) give the explicit values for */ - /* the respective upper and lower bounds. */ - - /* 1. For IND(J)=1, require X(J) .ge. BL(J). */ - /* 2. For IND(J)=2, require X(J) .le. BU(J). */ - /* 3. For IND(J)=3, require X(J) .ge. BL(J) and */ - /* X(J) .le. BU(J). */ - /* 4. For IND(J)=4, no bounds on X(J) are required. */ - /* The values of BL(*),BL(*) are modified by the subprogram. Values */ - /* other than 1,2,3 or 4 for IND(J) are errors. In the case IND(J)=3 */ - /* (upper and lower bounds) the condition BL(J) .gt. BU(J) is an */ - /* error. */ - - /* ------- */ - /* IOPT(*) */ - /* ------- */ - /* This is the array where the user can specify nonstandard options */ - /* for DBOLSM( ). Most of the time this feature can be ignored by */ - /* setting the input value IOPT(1)=99. Occasionally users may have */ - /* needs that require use of the following subprogram options. For */ - /* details about how to use the options see below: IOPT(*) CONTENTS. */ - - /* Option Number Brief Statement of Purpose */ - /* ----- ------ ----- --------- -- ------- */ - /* 1 Move the IOPT(*) processing pointer. */ - /* 2 Change rank determination tolerance. */ - /* 3 Change blow-up factor that determines the */ - /* size of variables being dropped from active */ - /* status. */ - /* 4 Reset the maximum number of iterations to use */ - /* in solving the problem. */ - /* 5 The data matrix is triangularized before the */ - /* problem is solved whenever (NCOLS/MROWS) .lt. */ - /* FAC. Change the value of FAC. */ - /* 6 Redefine the weighting matrix used for */ - /* linear independence checking. */ - /* 7 Debug output is desired. */ - /* 99 No more options to change. */ - - /* ---- */ - /* X(*) */ - /* ---- */ - /* This array is used to pass data associated with options 1,2,3 and */ - /* 5. Ignore this input parameter if none of these options are used. */ - /* Otherwise see below: IOPT(*) CONTENTS. */ - - /* ---------------- */ - /* IBASIS(*),IBB(*) */ - /* ---------------- */ - /* These arrays must be initialized by the user. The values */ - /* IBASIS(J)=J, J=1,...,NCOLS */ - /* IBB(J) =1, J=1,...,NCOLS */ - /* are appropriate except when using nonstandard features. */ - - /* ------ */ - /* SCL(*) */ - /* ------ */ - /* This is the array of scaling factors to use on the columns of the */ - /* matrix E. These values must be defined by the user. To suppress */ - /* any column scaling set SCL(J)=1.0, J=1,...,NCOLS. */ - - /* OUTPUT */ - /* ------ */ - - /* ---------- */ - /* X(*),RNORM */ - /* ---------- */ - /* The array X(*) contains a solution (if MODE .ge.0 or .eq.-22) for */ - /* the constrained least squares problem. The value RNORM is the */ - /* minimum residual vector length. */ - - /* ---- */ - /* MODE */ - /* ---- */ - /* The sign of mode determines whether the subprogram has completed */ - /* normally, or encountered an error condition or abnormal status. */ - /* A value of MODE .ge. 0 signifies that the subprogram has completed */ - /* normally. The value of MODE (.ge. 0) is the number of variables */ - /* in an active status: not at a bound nor at the value ZERO, for */ - /* the case of free variables. A negative value of MODE will be one */ - /* of the 20 cases -40,-39,...,-22, or -1. Values .lt. -1 correspond */ - /* to an abnormal completion of the subprogram. To understand the */ - /* abnormal completion codes see below: ERROR MESSAGES for DBOLSM( ) */ - /* An approximate solution will be returned to the user only when */ - /* max. iterations is reached, MODE=-22. */ - - /* ----------- */ - /* RW(*),WW(*) */ - /* ----------- */ - /* These are working arrays each with NCOLS entries. The array RW(*) */ - /* contains the working (scaled, nonactive) solution values. The */ - /* array WW(*) contains the working (scaled, active) gradient vector */ - /* values. */ - - /* ---------------- */ - /* IBASIS(*),IBB(*) */ - /* ---------------- */ - /* These arrays contain information about the status of the solution */ - /* when MODE .ge. 0. The indices IBASIS(K), K=1,...,MODE, show the */ - /* nonactive variables; indices IBASIS(K), K=MODE+1,..., NCOLS are */ - /* the active variables. The value (IBB(J)-1) is the number of times */ - /* variable J was reflected from its upper bound. (normally the user */ - /* can ignore these parameters.) */ - - /* IOPT(*) CONTENTS */ - /* ------- -------- */ - /* The option array allows a user to modify internal variables in */ - /* the subprogram without recompiling the source code. A central */ - /* goal of the initial software design was to do a good job for most */ - /* people. Thus the use of options will be restricted to a select */ - /* group of users. The processing of the option array proceeds as */ - /* follows: a pointer, here called LP, is initially set to the value */ - /* 1. The value is updated as the options are processed. At the */ - /* pointer position the option number is extracted and used for */ - /* locating other information that allows for options to be changed. */ - /* The portion of the array IOPT(*) that is used for each option is */ - /* fixed; the user and the subprogram both know how many locations */ - /* are needed for each option. A great deal of error checking is */ - /* done by the subprogram on the contents of the option array. */ - /* Nevertheless it is still possible to give the subprogram optional */ - /* input that is meaningless. For example some of the options use */ - /* the location X(NCOLS+IOFF) for passing data. The user must manage */ - /* the allocation of these locations when more than one piece of */ - /* option data is being passed to the subprogram. */ - - /* 1 */ - /* - */ - /* Move the processing pointer (either forward or backward) to the */ - /* location IOPT(LP+1). The processing pointer is moved to location */ - /* LP+2 of IOPT(*) in case IOPT(LP)=-1. For example to skip over */ - /* locations 3,...,NCOLS+2 of IOPT(*), */ - - /* IOPT(1)=1 */ - /* IOPT(2)=NCOLS+3 */ - /* (IOPT(I), I=3,...,NCOLS+2 are not defined here.) */ - /* IOPT(NCOLS+3)=99 */ - /* CALL DBOLSM( ) */ - - /* CAUTION: Misuse of this option can yield some very hard */ - /* -to-find bugs. Use it with care. */ - - /* 2 */ - /* - */ - /* The algorithm that solves the bounded least squares problem */ - /* iteratively drops columns from the active set. This has the */ - /* effect of joining a new column vector to the QR factorization of */ - /* the rectangular matrix consisting of the partially triangularized */ - /* nonactive columns. After triangularizing this matrix a test is */ - /* made on the size of the pivot element. The column vector is */ - /* rejected as dependent if the magnitude of the pivot element is */ - /* .le. TOL* magnitude of the column in components strictly above */ - /* the pivot element. Nominally the value of this (rank) tolerance */ - /* is TOL = SQRT(SRELPR), where SRELPR is relative machine */ - /* precision. To change only the value of TOL, for example, */ - - /* X(NCOLS+1)=TOL */ - /* IOPT(1)=2 */ - /* IOPT(2)=1 */ - /* IOPT(3)=99 */ - /* CALL DBOLSM() */ - - /* Generally, if LP is the processing pointer for IOPT(*), */ - - /* X(NCOLS+IOFF)=TOL */ - /* IOPT(LP)=2 */ - /* IOPT(LP+1)=IOFF */ - /* . */ - /* CALL DBOLSM() */ - - /* The required length of IOPT(*) is increased by 2 if option 2 is */ - /* used; The required length of X(*) is increased by 1. A value of */ - /* IOFF .le. 0 is an error. A value of TOL .le. SRELPR gives a */ - /* warning message; it is not considered an error. */ - /* Here SRELPR is the relative machine precision. */ - - /* 3 */ - /* - */ - /* A solution component is left active (not used) if, roughly */ - /* speaking, it seems too large. Mathematically the new component is */ - /* left active if the magnitude is .ge.((vector norm of F)/(matrix */ - /* norm of E))/BLOWUP. Nominally the factor BLOWUP = SQRT(SRELPR) */ - /* where SRELPR is the relative machine precision. To change only */ - /* the value of BLOWUP, for example, */ - - /* X(NCOLS+2)=BLOWUP */ - /* IOPT(1)=3 */ - /* IOPT(2)=2 */ - /* IOPT(3)=99 */ - /* CALL DBOLSM() */ - - /* Generally, if LP is the processing pointer for IOPT(*), */ - - /* X(NCOLS+IOFF)=BLOWUP */ - /* IOPT(LP)=3 */ - /* IOPT(LP+1)=IOFF */ - /* . */ - /* CALL DBOLSM() */ - - /* The required length of IOPT(*) is increased by 2 if option 3 is */ - /* used; the required length of X(*) is increased by 1. A value of */ - /* IOFF .le. 0 is an error. A value of BLOWUP .le. 0.0 is an error. */ - - /* 4 */ - /* - */ - /* Normally the algorithm for solving the bounded least squares */ - /* problem requires between NCOLS/3 and NCOLS drop-add steps to */ - /* converge. (this remark is based on examining a small number of */ - /* test cases.) The amount of arithmetic for such problems is */ - /* typically about twice that required for linear least squares if */ - /* there are no bounds and if plane rotations are used in the */ - /* solution method. Convergence of the algorithm, while */ - /* mathematically certain, can be much slower than indicated. To */ - /* avoid this potential but unlikely event ITMAX drop-add steps are */ - /* permitted. Nominally ITMAX=5*(MAX(MROWS,NCOLS)). To change the */ - /* value of ITMAX, for example, */ - - /* IOPT(1)=4 */ - /* IOPT(2)=ITMAX */ - /* IOPT(3)=99 */ - /* CALL DBOLSM() */ - - /* Generally, if LP is the processing pointer for IOPT(*), */ - - /* IOPT(LP)=4 */ - /* IOPT(LP+1)=ITMAX */ - /* . */ - /* CALL DBOLSM() */ - - /* The value of ITMAX must be .gt. 0. Other values are errors. Use */ - /* of this option increases the required length of IOPT(*) by 2. */ - - /* 5 */ - /* - */ - /* For purposes of increased efficiency the MROWS by NCOLS+1 data */ - /* matrix [E:F] is triangularized as a first step whenever MROWS */ - /* satisfies FAC*MROWS .gt. NCOLS. Nominally FAC=0.75. To change the */ - /* value of FAC, */ - - /* X(NCOLS+3)=FAC */ - /* IOPT(1)=5 */ - /* IOPT(2)=3 */ - /* IOPT(3)=99 */ - /* CALL DBOLSM() */ - - /* Generally, if LP is the processing pointer for IOPT(*), */ - - /* X(NCOLS+IOFF)=FAC */ - /* IOPT(LP)=5 */ - /* IOPT(LP+1)=IOFF */ - /* . */ - /* CALL DBOLSM() */ - - /* The value of FAC must be nonnegative. Other values are errors. */ - /* Resetting FAC=0.0 suppresses the initial triangularization step. */ - /* Use of this option increases the required length of IOPT(*) by 2; */ - /* The required length of of X(*) is increased by 1. */ - - /* 6 */ - /* - */ - /* The norm used in testing the magnitudes of the pivot element */ - /* compared to the mass of the column above the pivot line can be */ - /* changed. The type of change that this option allows is to weight */ - /* the components with an index larger than MVAL by the parameter */ - /* WT. Normally MVAL=0 and WT=1. To change both the values MVAL and */ - /* WT, where LP is the processing pointer for IOPT(*), */ - - /* X(NCOLS+IOFF)=WT */ - /* IOPT(LP)=6 */ - /* IOPT(LP+1)=IOFF */ - /* IOPT(LP+2)=MVAL */ - - /* Use of this option increases the required length of IOPT(*) by 3. */ - /* The length of X(*) is increased by 1. Values of MVAL must be */ - /* nonnegative and not greater than MROWS. Other values are errors. */ - /* The value of WT must be positive. Any other value is an error. If */ - /* either error condition is present a message will be printed. */ - - /* 7 */ - /* - */ - /* Debug output, showing the detailed add-drop steps for the */ - /* constrained least squares problem, is desired. This option is */ - /* intended to be used to locate suspected bugs. */ - - /* 99 */ - /* -- */ - /* There are no more options to change. */ - - /* The values for options are 1,...,7,99, and are the only ones */ - /* permitted. Other values are errors. Options -99,-1,...,-7 mean */ - /* that the repective options 99,1,...,7 are left at their default */ - /* values. An example is the option to modify the (rank) tolerance: */ - - /* X(NCOLS+1)=TOL */ - /* IOPT(1)=-2 */ - /* IOPT(2)=1 */ - /* IOPT(3)=99 */ - - /* Error Messages for DBOLSM( ) */ - /* ----- -------- --- --------- */ - /* WARNING IN... */ - /* DBOLSM(). MORE THAN (I1)=ITMAX ITERATIONS SOLVING BOUNDED LEAST */ - /* SQUARES PROBLEM. */ - /* IN ABOVE MESSAGE, I1= 3 */ - /* ERROR NUMBER = 22 */ - /* (NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.) */ - - /* WARNING IN... */ - /* DBOLSM. THE OPTION NUMBER=(I1) IS NOT DEFINED. */ - /* IN ABOVE MESSAGE, I1= 0 */ - /* ERROR NUMBER = 23 */ - /* (NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.) */ - - /* WARNING IN... */ - /* DBOLSM(). THE OFFSET=(I1) BEYOND POSTION NCOLS=(I2) */ - /* MUST BE POSITIVE FOR OPTION NUMBER 2. */ - /* IN ABOVE MESSAGE, I1= 0 */ - /* IN ABOVE MESSAGE, I2= 1 */ - /* ERROR NUMBER = 24 */ - /* (NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.) */ - - /* WARNING IN... */ - /* DBOLSM(). THE TOLERANCE FOR RANK DETERMINATION=(R1) */ - /* IS LESS THAN MACHINE PRECISION=(R2). */ - /* IN ABOVE MESSAGE, R1= 0. */ - /* IN ABOVE MESSAGE, R2= .7105427358E-14 */ - /* ERROR NUMBER = 25 */ - - /* WARNING IN... */ - /* DBOLSM(). THE OFFSET=(I1) BEYOND POSITION NCOLS=(I2) MUST */ - /* BE POSTIVE FOR OPTION NUMBER 3. */ - /* IN ABOVE MESSAGE, I1= 0 */ - /* IN ABOVE MESSAGE, I2= 1 */ - /* ERROR NUMBER = 26 */ - /* (NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.) */ - - /* WARNING IN... */ - /* DBOLSM(). THE RECIPROCAL OF THE BLOW-UP FACTOR FOR REJECTING */ - /* VARIABLES MUST BE POSITIVE. NOW=(R1). */ - /* IN ABOVE MESSAGE, R1= 0. */ - /* ERROR NUMBER = 27 */ - /* (NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.) */ - - /* WARNING IN... */ - /* DBOLSM(). THE MAXIMUM NUMBER OF ITERATIONS=(I1) MUST BE POSITIVE. */ - /* IN ABOVE MESSAGE, I1= 0 */ - /* ERROR NUMBER = 28 */ - /* (NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.) */ - - /* WARNING IN... */ - /* DBOLSM(). THE OFFSET=(I1) BEYOND POSITION NCOLS=(I2) MUST BE */ - /* POSTIVE FOR OPTION NUMBER 5. */ - /* IN ABOVE MESSAGE, I1= 0 */ - /* IN ABOVE MESSAGE, I2= 1 */ - /* ERROR NUMBER = 29 */ - /* (NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.) */ - - /* WARNING IN... */ - /* DBOLSM(). THE FACTOR (NCOLS/MROWS) WHERE PRETRIANGULARIZING IS */ - /* PERFORMED MUST BE NONNEGATIVE. NOW=(R1). */ - /* IN ABOVE MESSAGE, R1= -.2500000000E+00 */ - /* ERROR NUMBER = 30 */ - /* (NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.) */ - - /* WARNING IN... */ - /* DBOLSM(). THE NUMBER OF ROWS=(I1) MUST BE POSITIVE. */ - /* IN ABOVE MESSAGE, I1= 0 */ - /* ERROR NUMBER = 31 */ - /* (NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.) */ - - /* WARNING IN... */ - /* DBOLSM(). THE NUMBER OF COLS.=(I1) MUST BE POSTIVE. */ - /* IN ABOVE MESSAGE, I1= 0 */ - /* ERROR NUMBER = 32 */ - /* (NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.) */ - - /* WARNING IN... */ - /* DBOLSM(). THE ROW DIMENSION OF W(,)=(I1) MUST BE .GE. THE */ - /* NUMBER OF ROWS =(I2). */ - /* IN ABOVE MESSAGE, I1= 0 */ - /* IN ABOVE MESSAGE, I2= 1 */ - /* ERROR NUMBER = 33 */ - /* (NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.) */ - - /* WARNING IN... */ - /* DBOLSM(). FOR J=(I1) THE CONSTRAINT INDICATOR MUST BE 1-4. */ - /* IN ABOVE MESSAGE, I1= 1 */ - /* IN ABOVE MESSAGE, I2= 0 */ - /* ERROR NUMBER = 34 */ - /* (NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.) */ - - /* WARNING IN... */ - /* DBOLSM(). FOR J=(I1) THE LOWER BOUND=(R1) IS .GT. THE UPPER */ - /* BOUND=(R2). */ - /* IN ABOVE MESSAGE, I1= 1 */ - /* IN ABOVE MESSAGE, R1= 0. */ - /* IN ABOVE MESSAGE, R2= -.1000000000E+01 */ - /* ERROR NUMBER = 35 */ - /* (NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.) */ - - /* WARNING IN... */ - /* DBOLSM(). THE INPUT ORDER OF COLUMNS=(I1) IS NOT BETWEEN 1 */ - /* AND NCOLS=(I2). */ - /* IN ABOVE MESSAGE, I1= 0 */ - /* IN ABOVE MESSAGE, I2= 1 */ - /* ERROR NUMBER = 36 */ - /* (NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.) */ - - /* WARNING IN... */ - /* DBOLSM(). THE BOUND POLARITY FLAG IN COMPONENT J=(I1) MUST */ - /* BE POSITIVE. NOW=(I2). */ - /* IN ABOVE MESSAGE, I1= 1 */ - /* IN ABOVE MESSAGE, I2= 0 */ - /* ERROR NUMBER = 37 */ - /* (NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.) */ - - /* WARNING IN... */ - /* DBOLSM(). THE ROW SEPARATOR TO APPLY WEIGHTING (I1) MUST LIE */ - /* BETWEEN 0 AND MROWS (I2). WEIGHT (R1) MUST BE POSITIVE. */ - /* IN ABOVE MESSAGE, I1= -1 */ - /* IN ABOVE MESSAGE, I2= 2 */ - /* IN ABOVE MESSAGE, R1= 0. */ - /* ERROR NUMBER = 38 */ - /* (NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.) */ - - /* WARNING IN... */ - /* DBOLSM(). THE OFFSET (I1) BEYOND POSITION NCOLS=(I2) MUST BE */ - /* POSITIVE FOR OPTION NUMBER 7. */ - /* IN ABOVE MESSAGE, I1= -1 */ - /* IN ABOVE MESSAGE, I2= 2 */ - /* ERROR NUMBER = 39 */ - /* (NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.) */ - - /* WARNING IN... */ - /* DBOLSM(). THE COLUMN PIVOTING THRESHOLD FACTOR MUST BE */ - /* POSITIVE. NOW=(R1). */ - /* IN ABOVE MESSAGE, R1= 0. */ - /* ERROR NUMBER = 40 */ - /* (NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.) */ - /* ***END PROLOGUE DBOLSM */ - - /* PURPOSE */ - /* ------- */ - /* THIS IS THE MAIN SUBPROGRAM THAT SOLVES THE BOUNDED */ - /* LEAST SQUARES PROBLEM. THE PROBLEM SOLVED HERE IS: */ - - /* SOLVE E*X = F (LEAST SQUARES SENSE) */ - /* WITH BOUNDS ON SELECTED X VALUES. */ - - /* REVISED 850422-0800 */ - /* REVISED YYMMDD-HHMM */ - - /* TO CHANGE THIS SUBPROGRAM FROM SINGLE TO DOUBLE PRECISION BEGIN */ - /* EDITING AT THE CARD 'C++'. */ - /* CHANGE THE SUBPROGRAM NAME TO DBOLSM AND THE STRINGS */ - /* /SAXPY/ TO /DAXPY/, /SCOPY/ TO /DCOPY/, */ - /* /SDOT/ TO /DDOT/, /SNRM2/ TO /DNRM2/, */ - /* /SROTG/ TO /DROTG/, /SROT/ TO /DROT/, /R1MACH/ TO /D1MACH/, */ - /* /SVOUT/ TO /DVOUT/, /SMOUT/ TO /DMOUT/, */ - /* /SSWAP/ TO /DSWAP/, /E0/ TO /D0/, */ - /* /REAL / TO /DOUBLE PRECISION/. */ - /* ++ */ - - /* THIS VARIABLE REMAINS TYPE REAL. */ - /* ***FIRST EXECUTABLE STATEMENT DBOLSM */ - /* Parameter adjustments */ - w_dim1 = *mdw; - w_offset = 1 + w_dim1; - w -= w_offset; - --bl; - --bu; - --ind; - --iopt; - --x; - --rw; - --ww; - --scl; - --ibasis; - --ibb; - - /* Function Body */ - level = 1; - /* -> HKM 5/21/04 Put this line below in to avoid a valgrind incident */ - mrows = 0; - - /* VERIFY THAT THE PROBLEM DIMENSIONS ARE DEFINED PROPERLY. */ - if (*minput <= 0) { - nerr = 31; - nchar = 51; - xerrwv_("DBOLSM(). THE NUMBER OF ROWS=(I1) MUST BE POSITIVE.", &nchar, - &nerr, &level, &c__1, minput, &idum, &c__0, &rdum, &rdum, ( - ftnlen)51); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L610; - } - if (*ncols <= 0) { - nerr = 32; - nchar = 51; - xerrwv_("DBOLSM(). THE NUMBER OF COLS.=(I1) MUST BE POSTIVE.", &nchar, - &nerr, &level, &c__1, ncols, &idum, &c__0, &rdum, &rdum, ( - ftnlen)51); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L610; - } - if (*mdw < *minput) { - nerr = 33; - nchar = 78; - xerrwv_("DBOLSM(). THE ROW DIMENSION OF W(,)=(I1) MUST BE .GE. THE N" - "UMBER OF ROWS=(I2).", &nchar, &nerr, &level, &c__2, mdw, & - mrows, &c__0, &rdum, &rdum, (ftnlen)78); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L610; - } - - /* VERIFY THAT BOUND INFORMATION IS CORRECT. */ - i__1 = *ncols; - for (j = 1; j <= i__1; ++j) { - if (ind[j] < 1 || ind[j] > 4) { - nerr = 34; - nchar = 58; - xerrwv_("DBOLSM(). FOR J=(I1) THE CONSTRAINT INDICATOR MUST BE 1" - "-4.", &nchar, &nerr, &level, &c__2, &j, &ind[j], &c__0, & - rdum, &rdum, (ftnlen)58); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L610; - } - /* L10: */ - } - i__1 = *ncols; - for (j = 1; j <= i__1; ++j) { - if (ind[j] == 3) { - if (bu[j] < bl[j]) { - nerr = 35; - nchar = 71; - rdum2 = (real) bl[j]; - rdum = (real) bu[j]; - xerrwv_("DBOLSM(). FOR J=(I1) THE LOWER BOUND=(R1) IS .GT. T" - "HE UPPER BOUND=(R2).", &nchar, &nerr, &level, &c__1, & - j, &idum, &c__2, &rdum2, &rdum, (ftnlen)71); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L610; - } - } - /* L20: */ - } - - /* CHECK THAT PERMUTATION AND POLARITY ARRAYS HAVE BEEN SET. */ - i__1 = *ncols; - for (j = 1; j <= i__1; ++j) { - if (ibasis[j] < 1 || ibasis[j] > *ncols) { - nerr = 36; - nchar = 74; - xerrwv_("DBOLSM(). THE INPUT ORDER OF COLUMNS=(I1) IS NOT BETWEE" - "N 1 AND NCOLS=(I2).", &nchar, &nerr, &level, &c__2, & - ibasis[j], ncols, &c__0, &rdum, &rdum, (ftnlen)74); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L610; - } - if (ibb[j] <= 0) { - nerr = 37; - nchar = 81; - xerrwv_("DBOLSM(). THE BOUND POLARITY FLAG IN COMPONENT J=(I1) M" - "UST BE POSITIVE. NOW=(I2).", &nchar, &nerr, &level, &c__2, - &j, &ibb[j], &c__0, &rdum, &rdum, (ftnlen)81); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L610; - } - /* L30: */ - } - /* DO(PROCESS OPTION ARRAY) */ - goto L580; -L40: - /* DO(INITIALIZE VARIABLES AND DATA VALUES) */ - goto L470; -L50: - if (iprint > 0) { - i__1 = *ncols + 1; - dmout_(&mrows, &i__1, mdw, &w[w_offset], "(' PRETRI. INPUT MATRIX')", - &c_n4, (ftnlen)25); - dvout_(ncols, &bl[1], "(' LOWER BOUNDS')", &c_n4, (ftnlen)17); - dvout_(ncols, &bu[1], "(' UPPER BOUNDS')", &c_n4, (ftnlen)17); - } -L60: - ++iter; - if (iter <= itmax) { - goto L80; - } - nerr = 22; - nchar = 80; - xerrwv_("DBOLSM(). MORE THAN (I1)=ITMAX ITERATIONS SOLVING BOUNDED LEAST" - " SQUARES PROBLEM.", &nchar, &nerr, &level, &c__1, &itmax, &idum, & - c__0, &rdum, &rdum, (ftnlen)80); - /* DO(RESCALE AND TRANSLATE VARIABLES) */ - igopr = 1; - goto L130; -L70: - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L610; -L80: - /* DO(FIND A VARIABLE TO BECOME NON-ACTIVE) */ - goto L190; -L90: - if (found) { - goto L110; - } - /* DO(RESCALE AND TRANSLATE VARIABLES) */ - igopr = 2; - goto L130; -L100: - *mode = nsetb; - return 0; -L110: - /* DO(MAKE MOVE AND UPDATE FACTORIZATION) */ - goto L290; -L120: - goto L60; - /* PROCEDURE(RESCALE AND TRANSLATE VARIABLES) */ -L130: - dcopy_(&nsetb, &x[1], &c__1, &rw[1], &c__1); - x[1] = zero; - dcopy_(ncols, &x[1], &c__0, &x[1], &c__1); - i__1 = nsetb; - for (j = 1; j <= i__1; ++j) { - jcol = (i__2 = ibasis[j], abs(i__2)); - x[jcol] = rw[j] * (d__1 = scl[jcol], abs(d__1)); - /* L140: */ - } - i__1 = *ncols; - for (j = 1; j <= i__1; ++j) { - if (ibb[j] % 2 == 0) { - x[j] = bu[j] - x[j]; - } - /* L150: */ - } - i__1 = *ncols; - for (j = 1; j <= i__1; ++j) { - jcol = ibasis[j]; - if (jcol < 0) { - x[-jcol] = bl[-jcol] + x[-jcol]; - } - /* L160: */ - } - i__1 = *ncols; - for (j = 1; j <= i__1; ++j) { - if (scl[j] < zero) { - x[j] = -x[j]; - } - /* L170: */ - } - i__ = max(nsetb,mval); - i__1 = mrows - i__; - /* Computing MIN */ - i__2 = i__ + 1; - *rnorm = dnrm2_(&i__1, &w[min(i__2,mrows) + (*ncols + 1) * w_dim1], &c__1) - ; - /* END PROCEDURE */ - switch (igopr) { - case 1: - goto L70; - case 2: - goto L100; - } - /* PROCEDURE(FIND A VARIABLE TO BECOME NON-ACTIVE) */ -L190: - - /* COMPUTE (NEGATIVE) OF GRADIENT VECTOR, W= */ - /* (TRANSPOSE OF E)*(F-E*X). */ - ww[1] = zero; - dcopy_(ncols, &ww[1], &c__0, &ww[1], &c__1); - i__1 = *ncols; - for (j = nsetb + 1; j <= i__1; ++j) { - jcol = (i__2 = ibasis[j], abs(i__2)); - i__2 = mrows - nsetb; - /* Computing MIN */ - i__3 = nsetb + 1; - /* Computing MIN */ - i__4 = nsetb + 1; - ww[j] = ddot_(&i__2, &w[min(i__3,mrows) + j * w_dim1], &c__1, &w[min( - i__4,mrows) + (*ncols + 1) * w_dim1], &c__1) * (d__1 = scl[ - jcol], abs(d__1)); - /* L200: */ - } - if (iprint > 0) { - dvout_(ncols, &ww[1], "(' GRADIENT VALUES')", &c_n4, (ftnlen)20); - ivout_(ncols, &ibasis[1], "(' INTERNAL VARIABLE ORDER')", &c_n4, ( - ftnlen)28); - ivout_(ncols, &ibb[1], "(' BOUND POLARITY')", &c_n4, (ftnlen)19); - } -L210: - - /* IF ACTIVE SET = NUMBER OF TOTAL ROWS, QUIT. */ - if (nsetb == mrows) { - found = FALSE_; - /* EXIT PROCEDURE */ - goto L90; - } - - /* CHOOSE AN EXTREMAL COMPONENT OF GRADIENT VECTOR */ - /* FOR A CANDIDATE TO BECOME NON-ACTIVE. */ - wlarge = -big; - wmag = -big; - i__1 = *ncols; - for (j = nsetb + 1; j <= i__1; ++j) { - t = ww[j]; - if (t == big) { - goto L220; - } - itemp = ibasis[j]; - jcol = abs(itemp); - i__2 = mval - nsetb; - /* Computing MIN */ - i__3 = nsetb + 1; - t1 = dnrm2_(&i__2, &w[min(i__3,mrows) + j * w_dim1], &c__1); - if (itemp < 0) { - if (ibb[jcol] % 2 == 0) { - t = -t; - } - if (t < zero) { - goto L220; - } - if (mval > nsetb) { - t = t1; - } - if (t > wlarge) { - wlarge = t; - jlarge = j; - } - } else { - if (mval > nsetb) { - t = t1; - } - if (abs(t) > wmag) { - wmag = abs(t); - jmag = j; - } - } -L220: - ; - } - - /* CHOOSE MAG. OF LARGEST COMPONENT OF GRADIENT FOR CANDIDATE. */ - jbig = 0; - wbig = zero; - if (wlarge > zero) { - jbig = jlarge; - wbig = wlarge; - } - if (wmag >= wbig) { - jbig = jmag; - wbig = wmag; - } - if (jbig == 0) { - found = FALSE_; - if (iprint > 0) { - ivout_(&c__0, &i__, "(' FOUND NO VARIABLE TO ENTER')", &c_n4, ( - ftnlen)31); - } - /* EXIT PROCEDURE */ - goto L90; - } - - /* SEE IF THE INCOMING COL. IS SUFFICIENTLY INDEPENDENT. */ - /* THIS TEST IS MADE BEFORE AN ELIMINATION IS PERFORMED. */ - if (iprint > 0) { - ivout_(&c__1, &jbig, "(' TRY TO BRING IN THIS COL.')", &c_n4, (ftnlen) - 30); - } - if (mval <= nsetb) { - cl1 = dnrm2_(&mval, &w[jbig * w_dim1 + 1], &c__1); - i__1 = nsetb - mval; - /* Computing MIN */ - i__2 = mval + 1; - cl2 = abs(wt) * dnrm2_(&i__1, &w[min(i__2,mrows) + jbig * w_dim1], & - c__1); - i__1 = mrows - nsetb; - /* Computing MIN */ - i__2 = nsetb + 1; - cl3 = abs(wt) * dnrm2_(&i__1, &w[min(i__2,mrows) + jbig * w_dim1], & - c__1); - drotg_(&cl1, &cl2, &dc, &ds); - colabv = abs(cl1); - colblo = cl3; - } else { - cl1 = dnrm2_(&nsetb, &w[jbig * w_dim1 + 1], &c__1); - i__1 = mval - nsetb; - /* Computing MIN */ - i__2 = nsetb + 1; - cl2 = dnrm2_(&i__1, &w[min(i__2,mrows) + jbig * w_dim1], &c__1); - i__1 = mrows - mval; - /* Computing MIN */ - i__2 = mval + 1; - cl3 = abs(wt) * dnrm2_(&i__1, &w[min(i__2,mrows) + jbig * w_dim1], & - c__1); - colabv = cl1; - drotg_(&cl2, &cl3, &dc, &ds); - colblo = abs(cl2); - } - if (colblo <= tolind * colabv) { - ww[jbig] = big; - if (iprint > 0) { - ivout_(&c__0, &i__, "(' VARIABLE IS DEPENDENT, NOT USED.')", & - c_n4, (ftnlen)37); - } - goto L210; - } - - /* SWAP MATRIX COLS. NSETB+1 AND JBIG, PLUS POINTER INFO., AND */ - /* GRADIENT VALUES. */ - ++nsetb; - if (nsetb != jbig) { - dswap_(&mrows, &w[nsetb * w_dim1 + 1], &c__1, &w[jbig * w_dim1 + 1], & - c__1); - dswap_(&c__1, &ww[nsetb], &c__1, &ww[jbig], &c__1); - itemp = ibasis[nsetb]; - ibasis[nsetb] = ibasis[jbig]; - ibasis[jbig] = itemp; - } - - /* ELIMINATE ENTRIES BELOW THE PIVOT LINE IN COL. NSETB. */ - if (mrows > nsetb) { - i__1 = nsetb + 1; - for (i__ = mrows; i__ >= i__1; --i__) { - if (i__ == mval + 1) { - goto L230; - } - drotg_(&w[i__ - 1 + nsetb * w_dim1], &w[i__ + nsetb * w_dim1], & - sc, &ss); - w[i__ + nsetb * w_dim1] = zero; - i__2 = *ncols - nsetb + 1; - drot_(&i__2, &w[i__ - 1 + (nsetb + 1) * w_dim1], mdw, &w[i__ + ( - nsetb + 1) * w_dim1], mdw, &sc, &ss); -L230: - ; - } - if (mval >= nsetb && mval < mrows) { - drotg_(&w[nsetb + nsetb * w_dim1], &w[mval + 1 + nsetb * w_dim1], - &sc, &ss); - w[mval + 1 + nsetb * w_dim1] = zero; - i__1 = *ncols - nsetb + 1; - drot_(&i__1, &w[nsetb + (nsetb + 1) * w_dim1], mdw, &w[mval + 1 + - (nsetb + 1) * w_dim1], mdw, &sc, &ss); - } - } - if (w[nsetb + nsetb * w_dim1] == zero) { - ww[nsetb] = big; - --nsetb; - if (iprint > 0) { - ivout_(&c__0, &i__, "(' PIVOT IS ZERO, NOT USED.')", &c_n4, ( - ftnlen)29); - } - goto L210; - } - - /* CHECK THAT NEW VARIABLE IS MOVING IN THE RIGHT DIRECTION. */ - itemp = ibasis[nsetb]; - jcol = abs(itemp); - xnew = w[nsetb + (*ncols + 1) * w_dim1] / w[nsetb + nsetb * w_dim1] / ( - d__1 = scl[jcol], abs(d__1)); - /* CONT: DO BLOCK */ - /* QUIT: DO BLOCK */ - if (itemp < 0) { - /* IF(WW(NSETB).GE.ZERO.AND.XNEW.LE.ZERO) EXIT(QUIT) */ - /* IF(WW(NSETB).LE.ZERO.AND.XNEW.GE.ZERO) EXIT(QUIT) */ - if (ww[nsetb] >= zero && xnew <= zero) { - goto L240; - } - if (ww[nsetb] <= zero && xnew >= zero) { - goto L240; - } - } - /* EXIT(CONT) */ - goto L250; - /* END BLOCK */ -L240: - ww[nsetb] = big; - --nsetb; - if (iprint > 0) { - ivout_(&c__0, &i__, "(' VARIABLE HAS BAD DIRECTION, NOT USED.')", & - c_n4, (ftnlen)42); - } - goto L210; - /* END BLOCK */ -L250: - found = TRUE_; - /* EXIT PROCEDURE */ - goto L260; -L260: - /* END PROCEDURE */ - goto L90; - /* PROCEDURE(SOLVE THE TRIANGULAR SYSTEM) */ -L270: - dcopy_(&nsetb, &w[(*ncols + 1) * w_dim1 + 1], &c__1, &rw[1], &c__1); - for (j = nsetb; j >= 1; --j) { - rw[j] /= w[j + j * w_dim1]; - jcol = (i__1 = ibasis[j], abs(i__1)); - t = rw[j]; - if (ibb[jcol] % 2 == 0) { - rw[j] = -rw[j]; - } - i__1 = j - 1; - d__1 = -t; - daxpy_(&i__1, &d__1, &w[j * w_dim1 + 1], &c__1, &rw[1], &c__1); - rw[j] /= (d__1 = scl[jcol], abs(d__1)); - /* L280: */ - } - if (iprint > 0) { - dvout_(&nsetb, &rw[1], "(' SOLN. VALUES')", &c_n4, (ftnlen)17); - ivout_(&nsetb, &ibasis[1], "(' COLS. USED')", &c_n4, (ftnlen)15); - } - /* END PROCEDURE */ - switch (lgopr) { - case 1: - goto L300; - case 2: - goto L440; - } - /* PROCEDURE(MAKE MOVE AND UPDATE FACTORIZATION) */ -L290: - /* DO(SOLVE THE TRIANGULAR SYSTEM) */ - lgopr = 1; - goto L270; -L300: - - /* SEE IF THE UNCONSTRAINED SOL. (OBTAINED BY SOLVING THE */ - /* TRIANGULAR SYSTEM) SATISFIES THE PROBLEM BOUNDS. */ - alpha = two; - beta = two; - x[nsetb] = zero; - i__1 = nsetb; - for (j = 1; j <= i__1; ++j) { - itemp = ibasis[j]; - jcol = abs(itemp); - t1 = two; - t2 = two; - if (itemp < 0) { - bou = zero; - } else { - bou = bl[jcol]; - } - if (-bou != big) { - bou /= (d__1 = scl[jcol], abs(d__1)); - } - if (rw[j] <= bou) { - t1 = (x[j] - bou) / (x[j] - rw[j]); - } - bou = bu[jcol]; - if (bou != big) { - bou /= (d__1 = scl[jcol], abs(d__1)); - } - if (rw[j] >= bou) { - t2 = (bou - x[j]) / (rw[j] - x[j]); - } - - /* IF NOT, THEN COMPUTE A STEP LENGTH SO THAT THE */ - /* VARIABLES REMAIN FEASIBLE. */ - if (t1 < alpha) { - alpha = t1; - jdrop1 = j; - } - if (t2 < beta) { - beta = t2; - jdrop2 = j; - } - /* L310: */ - } - constr = alpha < two || beta < two; - if (constr) { - goto L320; - } - - /* ACCEPT THE CANDIDATE BECAUSE IT SATISFIES THE STATED BOUNDS */ - /* ON THE VARIABLES. */ - dcopy_(&nsetb, &rw[1], &c__1, &x[1], &c__1); - goto L120; -L320: - - /* TAKE A STEP THAT IS AS LARGE AS POSSIBLE WITH ALL */ - /* VARIABLES REMAINING FEASIBLE. */ - i__1 = nsetb; - for (j = 1; j <= i__1; ++j) { - x[j] += min(alpha,beta) * (rw[j] - x[j]); - /* L330: */ - } - if (alpha <= beta) { - jdrop2 = 0; - } else { - jdrop1 = 0; - } -L340: - if (jdrop1 + jdrop2 > 0 && nsetb > 0) { - goto L350; - } - goto L460; -L350: - jdrop = jdrop1 + jdrop2; - itemp = ibasis[jdrop]; - jcol = abs(itemp); - if (jdrop2 > 0) { - - /* VARIABLE IS AT AN UPPER BOUND. SUBTRACT MULTIPLE OF THIS COL. */ - /* FROM RIGHT HAND SIDE. */ - t = bu[jcol]; - if (itemp > 0) { - bu[jcol] = t - bl[jcol]; - bl[jcol] = -t; - itemp = -itemp; - scl[jcol] = -scl[jcol]; - i__1 = jdrop; - for (i__ = 1; i__ <= i__1; ++i__) { - w[i__ + jdrop * w_dim1] = -w[i__ + jdrop * w_dim1]; - /* L360: */ - } - } else { - ++ibb[jcol]; - if (ibb[jcol] % 2 == 0) { - t = -t; - } - } - /* VARIABLE IS AT A LOWER BOUND. */ - } else { - if ((doublereal) itemp < zero) { - t = zero; - } else { - t = -bl[jcol]; - bu[jcol] += t; - itemp = -itemp; - } - } - daxpy_(&jdrop, &t, &w[jdrop * w_dim1 + 1], &c__1, &w[(*ncols + 1) * - w_dim1 + 1], &c__1); - - /* MOVE CERTAIN COLS. LEFT TO ACHIEVE UPPER HESSENBERG FORM. */ - dcopy_(&jdrop, &w[jdrop * w_dim1 + 1], &c__1, &rw[1], &c__1); - i__1 = nsetb; - for (j = jdrop + 1; j <= i__1; ++j) { - ibasis[j - 1] = ibasis[j]; - x[j - 1] = x[j]; - dcopy_(&j, &w[j * w_dim1 + 1], &c__1, &w[(j - 1) * w_dim1 + 1], &c__1) - ; - /* L370: */ - } - ibasis[nsetb] = itemp; - w[nsetb * w_dim1 + 1] = zero; - i__1 = mrows - jdrop; - dcopy_(&i__1, &w[nsetb * w_dim1 + 1], &c__0, &w[jdrop + 1 + nsetb * - w_dim1], &c__1); - dcopy_(&jdrop, &rw[1], &c__1, &w[nsetb * w_dim1 + 1], &c__1); - - /* TRANSFORM THE MATRIX FROM UPPER HESSENBERG FORM TO */ - /* UPPER TRIANGULAR FORM. */ - --nsetb; - /* SMLL: */ - /* *DO BLOCK */ - /* NRML: */ - /* *DO BLOCK */ - i__1 = nsetb; - for (i__ = jdrop; i__ <= i__1; ++i__) { - - /* LOOK FOR SMALL PIVOTS AND AVOID MIXING WEIGHTED AND */ - /* NONWEIGHTED ROWS. */ - if (i__ == mval) { - t = zero; - i__2 = nsetb; - for (j = i__; j <= i__2; ++j) { - jcol = (i__3 = ibasis[j], abs(i__3)); - t1 = (d__1 = w[i__ + j * w_dim1] * scl[jcol], abs(d__1)); - if (t1 > t) { - jbig = j; - t = t1; - } - /* L380: */ - } - /* EXIT(NRML) */ - goto L400; - } - drotg_(&w[i__ + i__ * w_dim1], &w[i__ + 1 + i__ * w_dim1], &sc, &ss); - w[i__ + 1 + i__ * w_dim1] = zero; - i__2 = *ncols - i__ + 1; - drot_(&i__2, &w[i__ + (i__ + 1) * w_dim1], mdw, &w[i__ + 1 + (i__ + 1) - * w_dim1], mdw, &sc, &ss); - /* L390: */ - } - /* EXIT (SMLL) */ - goto L430; - /* END BLOCK */ -L400: - - /* THE TRIANGULARIZATION IS COMPLETED BY GIVING UP */ - /* THE HESSENBERG FORM AND TRIANGULARIZING A RECTANGULAR MATRIX. */ - dswap_(&mrows, &w[i__ * w_dim1 + 1], &c__1, &w[jbig * w_dim1 + 1], &c__1); - dswap_(&c__1, &ww[i__], &c__1, &ww[jbig], &c__1); - dswap_(&c__1, &x[i__], &c__1, &x[jbig], &c__1); - itemp = ibasis[i__]; - ibasis[i__] = ibasis[jbig]; - ibasis[jbig] = itemp; - jbig = i__; - i__1 = nsetb; - for (j = jbig; j <= i__1; ++j) { - i__2 = mrows; - for (i__ = j + 1; i__ <= i__2; ++i__) { - drotg_(&w[j + j * w_dim1], &w[i__ + j * w_dim1], &sc, &ss); - w[i__ + j * w_dim1] = zero; - i__3 = *ncols - j + 1; - drot_(&i__3, &w[j + (j + 1) * w_dim1], mdw, &w[i__ + (j + 1) * - w_dim1], mdw, &sc, &ss); - /* L410: */ - } - /* L420: */ - } - /* END BLOCK */ -L430: - - /* SEE IF THE REMAINING COEFFICIENTS ARE FEASIBLE. THEY SHOULD */ - /* BE BECAUSE OF THE WAY MIN(ALPHA,BETA) WAS CHOSEN. ANY THAT ARE */ - /* NOT FEASIBLE WILL BE SET TO THEIR BOUNDS AND */ - /* APPROPRIATELY TRANSLATED. */ - jdrop1 = 0; - jdrop2 = 0; - /* DO(SOLVE THE TRIANGULAR SYSTEM) */ - lgopr = 2; - goto L270; -L440: - dcopy_(&nsetb, &rw[1], &c__1, &x[1], &c__1); - i__1 = nsetb; - for (j = 1; j <= i__1; ++j) { - itemp = ibasis[j]; - jcol = abs(itemp); - if (itemp < 0) { - bou = zero; - } else { - bou = bl[jcol]; - } - if (-bou != big) { - bou /= (d__1 = scl[jcol], abs(d__1)); - } - if (x[j] <= bou) { - jdrop1 = j; - goto L340; - } - bou = bu[jcol]; - if (bou != big) { - bou /= (d__1 = scl[jcol], abs(d__1)); - } - if (x[j] >= bou) { - jdrop2 = j; - goto L340; - } - /* L450: */ - } - goto L340; -L460: - /* END PROCEDURE */ - goto L120; - /* PROCEDURE(INITIALIZE VARIABLES AND DATA VALUES) */ -L470: - - /* PRETRIANGULARIZE RECTANGULAR ARRAYS OF CERTAIN SIZES */ - /* FOR INCREASED EFFICIENCY. */ - if (fac * *minput > (doublereal)(*ncols)) { - i__1 = *ncols + 1; - for (j = 1; j <= i__1; ++j) { - i__2 = j + mval + 1; - for (i__ = *minput; i__ >= i__2; --i__) { - drotg_(&w[i__ - 1 + j * w_dim1], &w[i__ + j * w_dim1], &sc, & - ss); - w[i__ + j * w_dim1] = zero; - i__3 = *ncols - j + 1; - drot_(&i__3, &w[i__ - 1 + (j + 1) * w_dim1], mdw, &w[i__ + (j - + 1) * w_dim1], mdw, &sc, &ss); - /* L480: */ - } - /* L490: */ - } - mrows = *ncols + mval + 1; - } else { - mrows = *minput; - } - - /* SET THE X(*) ARRAY TO ZERO SO ALL COMPONENTS ARE DEFINED. */ - x[1] = zero; - dcopy_(ncols, &x[1], &c__0, &x[1], &c__1); - - /* THE ARRAYS IBASIS(*), IBB(*) ARE INITIALIZED BY THE CALLING */ - /* PROGRAM UNIT. */ - /* THE COL. SCALING IS DEFINED IN THE CALLING PROGRAM UNIT. */ - /* 'BIG' IS PLUS INFINITY ON THIS MACHINE. */ - big = d1mach_(&lc__2); - i__1 = *ncols; - for (j = 1; j <= i__1; ++j) { - icase = ind[j]; - /* DO CASE(ICASE,4) */ - switch (icase) { - case 1: - goto L500; - case 2: - goto L510; - case 3: - goto L520; - case 4: - goto L530; - } - goto L540; - /* CASE 1 */ -L500: - bu[j] = big; - goto L540; - /* CASE 2 */ -L510: - bl[j] = -big; - goto L540; - /* CASE 3 */ -L520: - goto L540; - /* CASE 4 */ -L530: - bl[j] = -big; - bu[j] = big; - /* END CASE */ -L540: - /* L550: */ - ; - } - i__1 = *ncols; - for (j = 1; j <= i__1; ++j) { - if ((bl[j] <= zero && zero <= bu[j] && - (d__1 = bu[j], abs(d__1)) < (d__2 = bl[j], abs(d__2))) || - bu[j] < zero) { - t = bu[j]; - bu[j] = -bl[j]; - bl[j] = -t; - scl[j] = -scl[j]; - i__2 = mrows; - for (i__ = 1; i__ <= i__2; ++i__) { - w[i__ + j * w_dim1] = -w[i__ + j * w_dim1]; - /* L560: */ - } - } - - /* INDICES IN SET T(=TIGHT) ARE DENOTED BY NEGATIVE VALUES */ - /* OF IBASIS(*). */ - if (bl[j] >= zero) { - ibasis[j] = -ibasis[j]; - t = -bl[j]; - bu[j] += t; - daxpy_(&mrows, &t, &w[j * w_dim1 + 1], &c__1, &w[(*ncols + 1) * - w_dim1 + 1], &c__1); - } - /* L570: */ - } - nsetb = 0; - iter = 0; - /* END PROCEDURE */ - goto L50; - /* PROCEDURE(PROCESS OPTION ARRAY) */ -L580: - zero = 0.; - fac = .75; - one = 1.; - two = 2.; - tolind = sqrt(d1mach_(&lc__4)); - tolsze = sqrt(d1mach_(&lc__4)); - itmax = max(mrows,*ncols) * 5; - wt = one; - mval = 0; - iprint = 0; - - /* CHANGES TO SOME PARAMETERS CAN OCCUR THROUGH THE OPTION */ - /* ARRAY, IOPT(*). PROCESS THIS ARRAY LOOKING CAREFULLY */ - /* FOR INPUT DATA ERRORS. */ - lp = 0; - lds = 0; -L590: - lp += lds; - - /* TEST FOR NO MORE OPTIONS. */ - ip = iopt[lp + 1]; - jp = abs(ip); - if (ip == 99) { - goto L600; - } else if (jp == 99) { - lds = 1; - goto L590; - } else if (jp == 1) { - - /* MOVE THE IOPT(*) PROCESSING POINTER. */ - if (ip > 0) { - lp = iopt[lp + 2] - 1; - lds = 0; - } else { - lds = 2; - } - goto L590; - } else if (jp == 2) { - - /* CHANGE TOLERANCE FOR RANK DETERMINATION. */ - if (ip > 0) { - ioff = iopt[lp + 2]; - if (ioff <= 0) { - nerr = 24; - nchar = 89; - xerrwv_("DBOLSM(). THE OFFSET=(I1) BEYOND POSTION NCOLS=(I2)" - " MUST BE POSITIVE FOR OPTION NUMBER 2.", &nchar, & - nerr, &level, &c__2, &ioff, ncols, &c__0, &rdum, & - rdum, (ftnlen)89); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L610; - } - tolind = x[*ncols + ioff]; - if (tolind < d1mach_(&lc__4)) { - nerr = 25; - nlevel = 0; - nchar = 88; - rdum2 = (real) tolind; - rdum = (real) d1mach_(&lc__4); - xerrwv_("DBOLSM(). THE TOLERANCE FOR RANK DETERMINATION=(R1)" - " IS LESS THAN MACHINE PRECISION=(R2).", &nchar, &nerr, - &nlevel, &c__0, &idum, &idum, &c__2, &rdum2, &rdum, ( - ftnlen)88); - } - } - lds = 2; - goto L590; - } else if (jp == 3) { - - /* CHANGE BLOWUP FACTOR FOR ALLOWING VARIABLES TO BECOME */ - /* INACTIVE. */ - if (ip > 0) { - ioff = iopt[lp + 2]; - if (ioff <= 0) { - nerr = 26; - nchar = 89; - xerrwv_("DBOLSM(). THE OFFSET=(I1) BEYOND POSITION NCOLS=(I2" - ") MUST BE POSTIVE FOR OPTION NUMBER 3.", &nchar, & - nerr, &level, &c__2, &ioff, ncols, &c__0, &rdum, & - rdum, (ftnlen)89); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L610; - } - tolsze = x[*ncols + ioff]; - if (tolsze <= zero) { - nerr = 27; - rdum2 = (real) tolsze; - xerrwv_("DBOLSM(). THE RECIPROCAL OF THE BLOW-UP FACTOR FOR " - "REJECTING VARIABLES MUST BE POSITIVE. NOW=(R1).", & - nchar, &nerr, &level, &c__0, &idum, &idum, &c__1, & - rdum2, &rdum, (ftnlen)98); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L610; - } - } - lds = 2; - goto L590; - } else if (jp == 4) { - - /* CHANGE THE MAX. NO. OF ITERATIONS ALLOWED. */ - if (ip > 0) { - itmax = iopt[lp + 2]; - if (itmax <= 0) { - nerr = 28; - nchar = 65; - xerrwv_("DBOLSM(). THE MAXIMUM NUMBER OF ITERATIONS=(I1) MUS" - "T BE POSITIVE.", &nchar, &nerr, &level, &c__1, &itmax, - &idum, &c__0, &rdum, &rdum, (ftnlen)65); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L610; - } - } - lds = 2; - goto L590; - } else if (jp == 5) { - - /* CHANGE THE FACTOR FOR PRETRIANGULARIZING THE DATA MATRIX. */ - if (ip > 0) { - ioff = iopt[lp + 2]; - if (ioff <= 0) { - nerr = 29; - nchar = 89; - xerrwv_("DBOLSM(). THE OFFSET=(I1) BEYOND POSITION NCOLS=(I2" - ") MUST BE POSTIVE FOR OPTION NUMBER 5.", &nchar, & - nerr, &level, &c__2, &ioff, ncols, &c__0, &rdum, & - rdum, (ftnlen)89); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L610; - } - fac = x[*ncols + ioff]; - if (fac < zero) { - nerr = 30; - nlevel = 0; - nchar = 104; - rdum2 = (real) fac; - xerrwv_("DBOLSM(). THE FACTOR (NCOLS/MROWS) WHERE PRE-TRIANG" - "ULARIZING IS PERFORMED MUST BE NONNEGATIVE. NOW=(R1)." - , &nchar, &nerr, &nlevel, &c__0, &idum, &idum, &c__1, - &rdum2, &rdum, (ftnlen)104); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L610; - } - } - lds = 2; - goto L590; - } else if (jp == 6) { - - /* CHANGE THE WEIGHTING FACTOR (FROM ONE) TO APPLY TO COMPONENTS */ - /* NUMBERED .GT. MVAL (INITIALLY SET TO 1.) THIS TRICK IS NEEDED */ - /* FOR APPLICATIONS OF THIS SUBPROGRAM TO THE HEAVILY WEIGHTED */ - /* LEAST SQUARES PROBLEM THAT COME FROM EQUALITY CONSTRAINTS. */ - if (ip > 0) { - ioff = iopt[lp + 2]; - mval = iopt[lp + 3]; - wt = x[*ncols + ioff]; - } - if (mval < 0 || mval > *minput || wt <= zero) { - nerr = 38; - nlevel = 0; - nchar = 116; - rdum2 = (real) wt; - xerrwv_("DBOLSM(). THE ROW SEPARATOR TO APPLY WEIGHTING (I1) MUS" - "T LIE BETWEEN 0 AND MROWS (I2). WEIGHT (R1) MUST BE POSI" - "TIVE.", &nchar, &nerr, &nlevel, &c__2, &mval, minput, & - c__1, &rdum2, &rdum, (ftnlen)116); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L610; - } - lds = 3; - goto L590; - - /* TURN ON DEBUG OUTPUT. */ - } else if (jp == 7) { - if (ip > 0) { - iprint = 1; - } - lds = 2; - goto L590; - } else { - nerr = 23; - nchar = 46; - xerrwv_("DBOLSM. THE OPTION NUMBER=(I1) IS NOT DEFINED.", &nchar, & - nerr, &level, &c__1, &ip, &idum, &c__0, &rdum, &rdum, (ftnlen) - 46); - /* DO(RETURN TO USER PROGRAM UNIT) */ - goto L610; - } -L600: - /* END PROCEDURE */ - goto L40; - /* PROCEDURE(RETURN TO USER PROGRAM UNIT) */ -L610: - *mode = -nerr; - return 0; - /* END PROCEDURE */ - /* END PROGRAM */ -} /* dbolsm_ */ - diff --git a/src/equil/vcs_dmout.c b/src/equil/vcs_dmout.c deleted file mode 100644 index 27b173fea..000000000 --- a/src/equil/vcs_dmout.c +++ /dev/null @@ -1,477 +0,0 @@ -/* dmout.f -- translated by f2c (version 20031025). - You must link the resulting object file with libf2c: - on Microsoft Windows system, link with libf2c.lib; - on Linux or Unix systems, link with .../path/to/libf2c.a -lm - or, if you install libf2c.a in a standard place, with -lf2c -lm - -- in that order, at the end of the command line, as in - cc *.o -lf2c -lm - Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., - - http://www.netlib.org/f2c/libf2c.zip -*/ - -#include "f2c.h" - -/* Table of constant values */ - -static integer c__2 = 2; -static integer c__1 = 1; - -/* DECK DMOUT */ -/* Subroutine */ -int dmout_(integer* m, integer* n, integer* lda, doublereal * - a, char* ifmt, integer* idigit, ftnlen ifmt_len) -{ - /* Initialized data */ - - static char icol[3] = "COL"; - - /* Format strings */ - static char fmt_1010[] = "(10x,10(4x,a,i4,1x))"; - static char fmt_1009[] = "(1x,\002ROW\002,i4,2x,1p10d12.3)"; - static char fmt_1000[] = "(10x,8(5x,a,i4,2x))"; - static char fmt_1004[] = "(1x,\002ROW\002,i4,2x,1p8d14.5)"; - static char fmt_1001[] = "(10x,5(9x,a,i4,6x))"; - static char fmt_1005[] = "(1x,\002ROW\002,i4,2x,1p5d22.13)"; - static char fmt_1002[] = "(10x,4(12x,a,i4,9x))"; - static char fmt_1006[] = "(1x,\002ROW\002,i4,2x,1p4d28.19)"; - static char fmt_1003[] = "(10x,3(16x,a,i4,13x))"; - static char fmt_1007[] = "(1x,\002ROW\002,i4,2x,1p3d36.27)"; - - /* System generated locals */ - integer a_dim1, a_offset, i__1, i__2, i__3; - cilist ci__1; - - /* Builtin functions */ - integer s_wsfe(cilist*), e_wsfe(void), do_fio(integer*, char*, ftnlen); - - /* Local variables */ - static integer i__, j, k1, k2, lout; - extern integer i1mach_(integer*); - static integer ndigit; - - /* Fortran I/O blocks */ - static cilist io___6 = { 0, 0, 0, fmt_1010, 0 }; - static cilist io___8 = { 0, 0, 0, fmt_1009, 0 }; - static cilist io___10 = { 0, 0, 0, fmt_1000, 0 }; - static cilist io___11 = { 0, 0, 0, fmt_1004, 0 }; - static cilist io___12 = { 0, 0, 0, fmt_1001, 0 }; - static cilist io___13 = { 0, 0, 0, fmt_1005, 0 }; - static cilist io___14 = { 0, 0, 0, fmt_1002, 0 }; - static cilist io___15 = { 0, 0, 0, fmt_1006, 0 }; - static cilist io___16 = { 0, 0, 0, fmt_1003, 0 }; - static cilist io___17 = { 0, 0, 0, fmt_1007, 0 }; - static cilist io___18 = { 0, 0, 0, fmt_1000, 0 }; - static cilist io___19 = { 0, 0, 0, fmt_1009, 0 }; - static cilist io___20 = { 0, 0, 0, fmt_1000, 0 }; - static cilist io___21 = { 0, 0, 0, fmt_1004, 0 }; - static cilist io___22 = { 0, 0, 0, fmt_1001, 0 }; - static cilist io___23 = { 0, 0, 0, fmt_1005, 0 }; - static cilist io___24 = { 0, 0, 0, fmt_1002, 0 }; - static cilist io___25 = { 0, 0, 0, fmt_1006, 0 }; - static cilist io___26 = { 0, 0, 0, fmt_1003, 0 }; - static cilist io___27 = { 0, 0, 0, fmt_1007, 0 }; - - - /* ***BEGIN PROLOGUE DMOUT */ - /* ***REFER TO DBOCLS,DFC */ - /* ***ROUTINES CALLED I1MACH */ - /* ***DESCRIPTION */ - - /* DOUBLE PRECISION MATRIX OUTPUT ROUTINE. */ - - /* INPUT.. */ - - /* M,N,LDA,A(*,*) PRINT THE DOUBLE PRECISION ARRAY A(I,J),I = 1,...,M, */ - /* J=1,...,N, ON OUTPUT UNIT LOUT=6. LDA IS THE DECLARED */ - /* FIRST DIMENSION OF A(*,*) AS SPECIFIED IN THE CALLING */ - /* PROGRAM. THE HEADING IN THE FORTRAN FORMAT STATEMENT */ - /* IFMT(*), DESCRIBED BELOW, IS PRINTED AS A FIRST STEP. */ - /* THE COMPONENTS A(I,J) ARE INDEXED, ON OUTPUT, IN A */ - /* PLEASANT FORMAT. */ - /* IFMT(*) A FORTRAN FORMAT STATEMENT. THIS IS PRINTED ON */ - /* OUTPUT UNIT LOUT=6 WITH THE VARIABLE FORMAT FORTRAN */ - /* STATEMENT */ - /* WRITE(LOUT,IFMT). */ - /* IDIGIT PRINT AT LEAST IABS(IDIGIT) DECIMAL DIGITS PER NUMBER. */ - /* THE SUBPROGRAM WILL CHOOSE THAT INTEGER 4,6,14,20 OR */ - /* 28 WHICH WILL PRINT AT LEAST IABS(IDIGIT) NUMBER OF */ - /* PLACES. IF IDIGIT.LT.0, 72 PRINTING COLUMNS ARE */ - /* UTILIZED TO WRITE EACH LINE OF OUTPUT OF THE ARRAY */ - /* A(*,*). (THIS CAN BE USED ON MOST TIME-SHARING */ - /* TERMINALS). IF IDIGIT.GE.0, 133 PRINTING COLUMNS ARE */ - /* UTILIZED. (THIS CAN BE USED ON MOST LINE PRINTERS). */ - - /* EXAMPLE.. */ - - /* PRINT AN ARRAY CALLED (SIMPLEX TABLEAU ) OF SIZE 10 BY 20 SHOWING */ - /* 6 DECIMAL DIGITS PER NUMBER. THE USER IS RUNNING ON A TIME-SHARING */ - /* SYSTEM WITH A 72 COLUMN OUTPUT DEVICE. */ - - /* DOUBLE PRECISION TABLEU(20,20) */ - /* M = 10 */ - /* N = 20 */ - /* LDTABL = 20 */ - /* IDIGIT = -6 */ - /* CALL DMOUT(M,N,LDTABL,TABLEU,21H(16H1SIMPLEX TABLEAU),IDIGIT) */ - - - - /* AUTHORS JOHN A. WISNIEWSKI SANDIA LABS ALBUQUERQUE. */ - /* RICHARD J. HANSON SANDIA LABS ALBUQUERQUE. */ - /* DATE JULY 30,1978. */ - /* ***END PROLOGUE DMOUT */ - /* Parameter adjustments */ - a_dim1 = *lda; - a_offset = 1 + a_dim1; - a -= a_offset; - - /* Function Body */ - /* ***FIRST EXECUTABLE STATEMENT DMOUT */ - lout = i1mach_(&c__2); - ci__1.cierr = 0; - ci__1.ciunit = lout; - ci__1.cifmt = ifmt; - s_wsfe(&ci__1); - e_wsfe(); - if (*m <= 0 || *n <= 0 || *lda <= 0) { - return 0; - } - ndigit = *idigit; - if (*idigit == 0) { - ndigit = 4; - } - if (*idigit >= 0) { - goto L80; - } - - ndigit = -(*idigit); - if (ndigit > 4) { - goto L9; - } - - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 5) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 4; - k2 = min(i__2,i__3); - io___6.ciunit = lout; - s_wsfe(&io___6); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, icol, (ftnlen)3); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - } - e_wsfe(); - i__2 = *m; - for (i__ = 1; i__ <= i__2; ++i__) { - io___8.ciunit = lout; - s_wsfe(&io___8); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - i__3 = k2; - for (j = k1; j <= i__3; ++j) { - do_fio(&c__1, (char*)&a[i__ + j * a_dim1], (ftnlen)sizeof( - doublereal)); - } - e_wsfe(); - /* L5: */ - } - } - return 0; - -L9: - if (ndigit > 6) { - goto L20; - } - - i__2 = *n; - for (k1 = 1; k1 <= i__2; k1 += 4) { - /* Computing MIN */ - i__1 = *n, i__3 = k1 + 3; - k2 = min(i__1,i__3); - io___10.ciunit = lout; - s_wsfe(&io___10); - i__1 = k2; - for (i__ = k1; i__ <= i__1; ++i__) { - do_fio(&c__1, icol, (ftnlen)3); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - } - e_wsfe(); - i__1 = *m; - for (i__ = 1; i__ <= i__1; ++i__) { - io___11.ciunit = lout; - s_wsfe(&io___11); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - i__3 = k2; - for (j = k1; j <= i__3; ++j) { - do_fio(&c__1, (char*)&a[i__ + j * a_dim1], (ftnlen)sizeof( - doublereal)); - } - e_wsfe(); - /* L10: */ - } - } - return 0; - -L20: - if (ndigit > 14) { - goto L40; - } - - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 2) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 1; - k2 = min(i__2,i__3); - io___12.ciunit = lout; - s_wsfe(&io___12); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, icol, (ftnlen)3); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - } - e_wsfe(); - i__2 = *m; - for (i__ = 1; i__ <= i__2; ++i__) { - io___13.ciunit = lout; - s_wsfe(&io___13); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - i__3 = k2; - for (j = k1; j <= i__3; ++j) { - do_fio(&c__1, (char*)&a[i__ + j * a_dim1], (ftnlen)sizeof( - doublereal)); - } - e_wsfe(); - /* L30: */ - } - } - return 0; - -L40: - if (ndigit > 20) { - goto L60; - } - - i__2 = *n; - for (k1 = 1; k1 <= i__2; k1 += 2) { - /* Computing MIN */ - i__1 = *n, i__3 = k1 + 1; - k2 = min(i__1,i__3); - io___14.ciunit = lout; - s_wsfe(&io___14); - i__1 = k2; - for (i__ = k1; i__ <= i__1; ++i__) { - do_fio(&c__1, icol, (ftnlen)3); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - } - e_wsfe(); - i__1 = *m; - for (i__ = 1; i__ <= i__1; ++i__) { - io___15.ciunit = lout; - s_wsfe(&io___15); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - i__3 = k2; - for (j = k1; j <= i__3; ++j) { - do_fio(&c__1, (char*)&a[i__ + j * a_dim1], (ftnlen)sizeof( - doublereal)); - } - e_wsfe(); - /* L50: */ - } - } - return 0; - -L60: - i__1 = *n; - for (k1 = 1; k1 <= i__1; ++k1) { - k2 = k1; - io___16.ciunit = lout; - s_wsfe(&io___16); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, icol, (ftnlen)3); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - } - e_wsfe(); - i__2 = *m; - for (i__ = 1; i__ <= i__2; ++i__) { - io___17.ciunit = lout; - s_wsfe(&io___17); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - i__3 = k2; - for (j = k1; j <= i__3; ++j) { - do_fio(&c__1, (char*)&a[i__ + j * a_dim1], (ftnlen)sizeof( - doublereal)); - } - e_wsfe(); - /* L70: */ - } - } - return 0; - -L80: - if (ndigit > 4) { - goto L86; - } - - i__2 = *n; - for (k1 = 1; k1 <= i__2; k1 += 10) { - /* Computing MIN */ - i__1 = *n, i__3 = k1 + 9; - k2 = min(i__1,i__3); - io___18.ciunit = lout; - s_wsfe(&io___18); - i__1 = k2; - for (i__ = k1; i__ <= i__1; ++i__) { - do_fio(&c__1, icol, (ftnlen)3); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - } - e_wsfe(); - i__1 = *m; - for (i__ = 1; i__ <= i__1; ++i__) { - io___19.ciunit = lout; - s_wsfe(&io___19); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - i__3 = k2; - for (j = k1; j <= i__3; ++j) { - do_fio(&c__1, (char*)&a[i__ + j * a_dim1], (ftnlen)sizeof( - doublereal)); - } - e_wsfe(); - /* L85: */ - } - } - -L86: - if (ndigit > 6) { - goto L100; - } - - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 8) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 7; - k2 = min(i__2,i__3); - io___20.ciunit = lout; - s_wsfe(&io___20); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, icol, (ftnlen)3); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - } - e_wsfe(); - i__2 = *m; - for (i__ = 1; i__ <= i__2; ++i__) { - io___21.ciunit = lout; - s_wsfe(&io___21); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - i__3 = k2; - for (j = k1; j <= i__3; ++j) { - do_fio(&c__1, (char*)&a[i__ + j * a_dim1], (ftnlen)sizeof( - doublereal)); - } - e_wsfe(); - /* L90: */ - } - } - return 0; - -L100: - if (ndigit > 14) { - goto L120; - } - - i__2 = *n; - for (k1 = 1; k1 <= i__2; k1 += 5) { - /* Computing MIN */ - i__1 = *n, i__3 = k1 + 4; - k2 = min(i__1,i__3); - io___22.ciunit = lout; - s_wsfe(&io___22); - i__1 = k2; - for (i__ = k1; i__ <= i__1; ++i__) { - do_fio(&c__1, icol, (ftnlen)3); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - } - e_wsfe(); - i__1 = *m; - for (i__ = 1; i__ <= i__1; ++i__) { - io___23.ciunit = lout; - s_wsfe(&io___23); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - i__3 = k2; - for (j = k1; j <= i__3; ++j) { - do_fio(&c__1, (char*)&a[i__ + j * a_dim1], (ftnlen)sizeof( - doublereal)); - } - e_wsfe(); - /* L110: */ - } - } - return 0; - -L120: - if (ndigit > 20) { - goto L140; - } - - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 4) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 3; - k2 = min(i__2,i__3); - io___24.ciunit = lout; - s_wsfe(&io___24); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, icol, (ftnlen)3); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - } - e_wsfe(); - i__2 = *m; - for (i__ = 1; i__ <= i__2; ++i__) { - io___25.ciunit = lout; - s_wsfe(&io___25); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - i__3 = k2; - for (j = k1; j <= i__3; ++j) { - do_fio(&c__1, (char*)&a[i__ + j * a_dim1], (ftnlen)sizeof( - doublereal)); - } - e_wsfe(); - /* L130: */ - } - } - return 0; - -L140: - i__2 = *n; - for (k1 = 1; k1 <= i__2; k1 += 3) { - /* Computing MIN */ - i__1 = *n, i__3 = k1 + 2; - k2 = min(i__1,i__3); - io___26.ciunit = lout; - s_wsfe(&io___26); - i__1 = k2; - for (i__ = k1; i__ <= i__1; ++i__) { - do_fio(&c__1, icol, (ftnlen)3); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - } - e_wsfe(); - i__1 = *m; - for (i__ = 1; i__ <= i__1; ++i__) { - io___27.ciunit = lout; - s_wsfe(&io___27); - do_fio(&c__1, (char*)&i__, (ftnlen)sizeof(integer)); - i__3 = k2; - for (j = k1; j <= i__3; ++j) { - do_fio(&c__1, (char*)&a[i__ + j * a_dim1], (ftnlen)sizeof( - doublereal)); - } - e_wsfe(); - /* L150: */ - } - } - return 0; -} /* dmout_ */ - diff --git a/src/equil/vcs_dvout.c b/src/equil/vcs_dvout.c deleted file mode 100644 index 1e3d68e90..000000000 --- a/src/equil/vcs_dvout.c +++ /dev/null @@ -1,295 +0,0 @@ -/* dvout.f -- translated by f2c (version 20031025). - You must link the resulting object file with libf2c: - on Microsoft Windows system, link with libf2c.lib; - on Linux or Unix systems, link with .../path/to/libf2c.a -lm - or, if you install libf2c.a in a standard place, with -lf2c -lm - -- in that order, at the end of the command line, as in - cc *.o -lf2c -lm - Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., - - http://www.netlib.org/f2c/libf2c.zip -*/ - -#include "f2c.h" - -/* Table of constant values */ - -static integer c__2 = 2; -static integer c__1 = 1; - -/* DECK DVOUT */ -/* Subroutine */ -int dvout_(integer* n, doublereal* dx, char* ifmt, integer * - idigit, ftnlen ifmt_len) -{ - /* Format strings */ - static char fmt_1000[] = "(1x,i4,\002 - \002,i4,1x,1p8d14.5)"; - static char fmt_1001[] = "(1x,i4,\002 - \002,i4,1x,1p5d22.13)"; - static char fmt_1002[] = "(1x,i4,\002 - \002,i4,1x,1p4d28.19)"; - static char fmt_1003[] = "(1x,i4,\002 - \002,i4,1x,1p3d36.27)"; - - /* System generated locals */ - integer i__1, i__2, i__3; - cilist ci__1; - - /* Builtin functions */ - integer s_wsfe(cilist*), e_wsfe(void), do_fio(integer*, char*, ftnlen); - - /* Local variables */ - static integer i__, k1, k2, lout; - extern integer i1mach_(integer*); - static integer ndigit; - - /* Fortran I/O blocks */ - static cilist io___5 = { 0, 0, 0, fmt_1000, 0 }; - static cilist io___7 = { 0, 0, 0, fmt_1001, 0 }; - static cilist io___8 = { 0, 0, 0, fmt_1002, 0 }; - static cilist io___9 = { 0, 0, 0, fmt_1003, 0 }; - static cilist io___10 = { 0, 0, 0, fmt_1000, 0 }; - static cilist io___11 = { 0, 0, 0, fmt_1001, 0 }; - static cilist io___12 = { 0, 0, 0, fmt_1002, 0 }; - static cilist io___13 = { 0, 0, 0, fmt_1003, 0 }; - - - /* ***BEGIN PROLOGUE DVOUT */ - /* ***REFER TO DSPLP */ - /* ***ROUTINES CALLED I1MACH */ - /* ***DESCRIPTION */ - - /* REVISED FEB. 27, 1981. */ - - /* DOUBLE PRECISION VECTOR OUTPUT ROUTINE. */ - - /* INPUT.. */ - - /* N,DX(*) PRINT THE DOUBLE PRECISION ARRAY DX(I),I=1,...,N, ON */ - /* OUTPUT UNIT LOUT. THE HEADING IN THE FORTRAN FORMAT */ - /* STATEMENT IFMT(*), DESCRIBED BELOW, IS PRINTED AS A FIRST */ - /* STEP. THE COMPONENTS DX(I) ARE INDEXED, ON OUTPUT, */ - /* IN A PLEASANT FORMAT. */ - /* IFMT(*) A FORTRAN FORMAT STATEMENT. THIS IS PRINTED ON OUTPUT */ - /* UNIT LOUT WITH THE VARIABLE FORMAT FORTRAN STATEMENT */ - /* WRITE(LOUT,IFMT) */ - /* IDIGIT PRINT AT LEAST IABS(IDIGIT) DECIMAL DIGITS PER NUMBER. */ - /* THE SUBPROGRAM WILL CHOOSE THAT INTEGER 4,6,10 OR 14 */ - /* WHICH WILL PRINT AT LEAST IABS(IDIGIT) NUMBER OF */ - /* PLACES. IF IDIGIT.LT.0, 72 PRINTING COLUMNS ARE UTILIZED */ - /* TO WRITE EACH LINE OF OUTPUT OF THE ARRAY DX(*). (THIS */ - /* CAN BE USED ON MOST TIME-SHARING TERMINALS). IF */ - /* IDIGIT.GE.0, 133 PRINTING COLUMNS ARE UTILIZED. (THIS CAN */ - /* BE USED ON MOST LINE PRINTERS). */ - - /* EXAMPLE.. */ - - /* PRINT AN ARRAY CALLED (COSTS OF PURCHASES) OF LENGTH 100 SHOWING */ - /* 6 DECIMAL DIGITS PER NUMBER. THE USER IS RUNNING ON A TIME-SHARING */ - /* SYSTEM WITH A 72 COLUMN OUTPUT DEVICE. */ - - /* DOUBLE PRECISION COSTS(100) */ - /* N = 100 */ - /* IDIGIT = -6 */ - /* CALL DVOUT(N,COSTS,'(''1COSTS OF PURCHASES'')',IDIGIT) */ - - - - /* AUTHORS JOHN A. WISNIEWSKI SANDIA LABS ALBUQUERQUE. */ - /* RICHARD J. HANSON SANDIA LABS ALBUQUERQUE. */ - /* DATE JULY 27,1978. */ - /* ***END PROLOGUE DVOUT */ - /* ***FIRST EXECUTABLE STATEMENT DVOUT */ - /* Parameter adjustments */ - --dx; - - /* Function Body */ - lout = i1mach_(&c__2); - ci__1.cierr = 0; - ci__1.ciunit = lout; - ci__1.cifmt = ifmt; - s_wsfe(&ci__1); - e_wsfe(); - if (*n <= 0) { - return 0; - } - ndigit = *idigit; - if (*idigit == 0) { - ndigit = 6; - } - if (*idigit >= 0) { - goto L80; - } - - ndigit = -(*idigit); - if (ndigit > 6) { - goto L20; - } - - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 4) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 3; - k2 = min(i__2,i__3); - io___5.ciunit = lout; - s_wsfe(&io___5); - do_fio(&c__1, (char*)&k1, (ftnlen)sizeof(integer)); - do_fio(&c__1, (char*)&k2, (ftnlen)sizeof(integer)); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, (char*)&dx[i__], (ftnlen)sizeof(doublereal)); - } - e_wsfe(); - /* L10: */ - } - return 0; - -L20: - if (ndigit > 14) { - goto L40; - } - - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 2) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 1; - k2 = min(i__2,i__3); - io___7.ciunit = lout; - s_wsfe(&io___7); - do_fio(&c__1, (char*)&k1, (ftnlen)sizeof(integer)); - do_fio(&c__1, (char*)&k2, (ftnlen)sizeof(integer)); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, (char*)&dx[i__], (ftnlen)sizeof(doublereal)); - } - e_wsfe(); - /* L30: */ - } - return 0; - -L40: - if (ndigit > 20) { - goto L60; - } - - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 2) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 1; - k2 = min(i__2,i__3); - io___8.ciunit = lout; - s_wsfe(&io___8); - do_fio(&c__1, (char*)&k1, (ftnlen)sizeof(integer)); - do_fio(&c__1, (char*)&k2, (ftnlen)sizeof(integer)); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, (char*)&dx[i__], (ftnlen)sizeof(doublereal)); - } - e_wsfe(); - /* L50: */ - } - return 0; - -L60: - i__1 = *n; - for (k1 = 1; k1 <= i__1; ++k1) { - k2 = k1; - io___9.ciunit = lout; - s_wsfe(&io___9); - do_fio(&c__1, (char*)&k1, (ftnlen)sizeof(integer)); - do_fio(&c__1, (char*)&k2, (ftnlen)sizeof(integer)); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, (char*)&dx[i__], (ftnlen)sizeof(doublereal)); - } - e_wsfe(); - /* L70: */ - } - return 0; - -L80: - if (ndigit > 6) { - goto L100; - } - - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 8) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 7; - k2 = min(i__2,i__3); - io___10.ciunit = lout; - s_wsfe(&io___10); - do_fio(&c__1, (char*)&k1, (ftnlen)sizeof(integer)); - do_fio(&c__1, (char*)&k2, (ftnlen)sizeof(integer)); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, (char*)&dx[i__], (ftnlen)sizeof(doublereal)); - } - e_wsfe(); - /* L90: */ - } - return 0; - -L100: - if (ndigit > 14) { - goto L120; - } - - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 5) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 4; - k2 = min(i__2,i__3); - io___11.ciunit = lout; - s_wsfe(&io___11); - do_fio(&c__1, (char*)&k1, (ftnlen)sizeof(integer)); - do_fio(&c__1, (char*)&k2, (ftnlen)sizeof(integer)); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, (char*)&dx[i__], (ftnlen)sizeof(doublereal)); - } - e_wsfe(); - /* L110: */ - } - return 0; - -L120: - if (ndigit > 20) { - goto L140; - } - - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 4) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 3; - k2 = min(i__2,i__3); - io___12.ciunit = lout; - s_wsfe(&io___12); - do_fio(&c__1, (char*)&k1, (ftnlen)sizeof(integer)); - do_fio(&c__1, (char*)&k2, (ftnlen)sizeof(integer)); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, (char*)&dx[i__], (ftnlen)sizeof(doublereal)); - } - e_wsfe(); - /* L130: */ - } - return 0; - -L140: - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 3) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 2; - k2 = min(i__2,i__3); - io___13.ciunit = lout; - s_wsfe(&io___13); - do_fio(&c__1, (char*)&k1, (ftnlen)sizeof(integer)); - do_fio(&c__1, (char*)&k2, (ftnlen)sizeof(integer)); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, (char*)&dx[i__], (ftnlen)sizeof(doublereal)); - } - e_wsfe(); - /* L150: */ - } - return 0; -} /* dvout_ */ - diff --git a/src/equil/vcs_inest.cpp b/src/equil/vcs_inest.cpp index 3076a22f2..8f5980274 100644 --- a/src/equil/vcs_inest.cpp +++ b/src/equil/vcs_inest.cpp @@ -62,35 +62,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, * * Note, both of these programs do this. */ -#ifdef ALTLINPROG vcs_setMolesLinProg(); -#else - std::vector ax(m_numElemConstraints*nspecies, 0.0); - std::vector bb(m_numElemConstraints, 0.0); - std::vector cc(nspecies, 0.0); - - int neActive = 0; - size_t jj = 0; - for (size_t j = 0; j < m_numElemConstraints; j++) { - if (m_elementActive[j]) { - neActive++; - bb[jj] = m_elemAbundancesGoal[j]; - jj++; - } - } - for (kspec = 0; kspec < nspecies; ++kspec) { - cc[kspec] = -m_SSfeSpecies[kspec]; - jj = 0; - for (size_t j = 0; j < m_numElemConstraints; ++j) { - if (m_elementActive[j]) { - ax[jj + kspec * neActive] = m_formulaMatrix[j][kspec]; - jj++; - } - } - } - linprogmax(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(cc), VCS_DATA_PTR(ax), - VCS_DATA_PTR(bb), neActive, nspecies, neActive); -#endif #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { diff --git a/src/equil/vcs_ivout.c b/src/equil/vcs_ivout.c deleted file mode 100644 index 3dec630ad..000000000 --- a/src/equil/vcs_ivout.c +++ /dev/null @@ -1,299 +0,0 @@ -/* ivout.f -- translated by f2c (version 20031025). - You must link the resulting object file with libf2c: - on Microsoft Windows system, link with libf2c.lib; - on Linux or Unix systems, link with .../path/to/libf2c.a -lm - or, if you install libf2c.a in a standard place, with -lf2c -lm - -- in that order, at the end of the command line, as in - cc *.o -lf2c -lm - Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., - - http://www.netlib.org/f2c/libf2c.zip -*/ - -#include "f2c.h" - -/* Table of constant values */ - -static integer c__1 = 1; - -/* DECK IVOUT */ -/* Subroutine */ -int ivout_(integer* n, integer* ix, char* ifmt, integer * - idigit, ftnlen ifmt_len) -{ - /* Format strings */ - static char fmt_1000[] = "(1x,i4,\002 - \002,i4,20(1x,i5))"; - static char fmt_1001[] = "(1x,i4,\002 - \002,i4,15(1x,i7))"; - static char fmt_1002[] = "(1x,i4,\002 - \002,i4,10(1x,i11))"; - static char fmt_1003[] = "(1x,i4,\002 - \002,i4,7(1x,i15))"; - - /* System generated locals */ - integer i__1, i__2, i__3; - cilist ci__1; - - /* Builtin functions */ - integer s_wsfe(cilist*), e_wsfe(void), do_fio(integer*, char*, ftnlen); - - /* Local variables */ - static integer i__, j, k1, k2, lout; - extern integer i1mach_(integer*); - static integer ndigit; - - /* Fortran I/O blocks */ - static cilist io___6 = { 0, 0, 0, fmt_1000, 0 }; - static cilist io___8 = { 0, 0, 0, fmt_1001, 0 }; - static cilist io___9 = { 0, 0, 0, fmt_1002, 0 }; - static cilist io___10 = { 0, 0, 0, fmt_1003, 0 }; - static cilist io___11 = { 0, 0, 0, fmt_1000, 0 }; - static cilist io___12 = { 0, 0, 0, fmt_1001, 0 }; - static cilist io___13 = { 0, 0, 0, fmt_1002, 0 }; - static cilist io___14 = { 0, 0, 0, fmt_1003, 0 }; - - - /* ***BEGIN PROLOGUE IVOUT */ - /* ***REFER TO SPLP */ - /* ***ROUTINES CALLED I1MACH */ - /* ***DESCRIPTION */ - - /* REVISED FEB. 27, 1981. */ - - /* INTEGER VECTOR OUTPUT ROUTINE. */ - - /* INPUT.. */ - - /* N,IX(*) PRINT THE INTEGER ARRAY IX(I),I=1,...,N, ON OUTPUT */ - /* UNIT LOUT. THE HEADING IN THE FORTRAN FORMAT */ - /* STATEMENT IFMT(*), DESCRIBED BELOW, IS PRINTED AS A FIRST */ - /* STEP. THE COMPONENTS IX(I) ARE INDEXED, ON OUTPUT, */ - /* IN A PLEASANT FORMAT. */ - /* IFMT(*) A FORTRAN FORMAT STATEMENT. THIS IS PRINTED ON OUTPUT */ - /* UNIT LOUT WITH THE VARIABLE FORMAT FORTRAN STATEMENT */ - /* WRITE(LOUT,IFMT) */ - /* IDIGIT PRINT UP TO IABS(IDIGIT) DECIMAL DIGITS PER NUMBER. */ - /* THE SUBPROGRAM WILL CHOOSE THAT INTEGER 4,6,10 OR 14 */ - /* WHICH WILL PRINT AT LEAST IABS(IDIGIT) NUMBER OF */ - /* PLACES. IF IDIGIT.LT.0, 72 PRINTING COLUMNS ARE UTILIZED */ - /* TO WRITE EACH LINE OF OUTPUT OF THE ARRAY IX(*). (THIS */ - /* CAN BE USED ON MOST TIME-SHARING TERMINALS). IF */ - /* IDIGIT.GE.0, 133 PRINTING COLUMNS ARE UTILIZED. (THIS CAN */ - /* BE USED ON MOST LINE PRINTERS). */ - - /* EXAMPLE.. */ - - /* PRINT AN ARRAY CALLED (COSTS OF PURCHASES) OF LENGTH 100 SHOWING */ - /* 6 DECIMAL DIGITS PER NUMBER. THE USER IS RUNNING ON A TIME-SHARING */ - /* SYSTEM WITH A 72 COLUMN OUTPUT DEVICE. */ - - /* DIMENSION ICOSTS(100) */ - /* N = 100 */ - /* IDIGIT = -6 */ - /* CALL IVOUT(N,ICOSTS,'(''1COSTS OF PURCHASES'')',IDIGIT) */ - - - - /* AUTHORS JOHN A. WISNIEWSKI SANDIA LABS ALBUQUERQUE. */ - /* RICHARD J. HANSON SANDIA LABS ALBUQUERQUE. */ - /* DATE JULY 27,1978. */ - /* ***END PROLOGUE IVOUT */ - - /* GET THE UNIT NUMBER WHERE OUTPUTWILL BE WRITTEN. */ - /* ***FIRST EXECUTABLE STATEMENT IVOUT */ - /* Parameter adjustments */ - --ix; - - /* Function Body */ - j = 2; - lout = i1mach_(&j); - ci__1.cierr = 0; - ci__1.ciunit = lout; - ci__1.cifmt = ifmt; - s_wsfe(&ci__1); - e_wsfe(); - if (*n <= 0) { - return 0; - } - ndigit = *idigit; - if (*idigit == 0) { - ndigit = 4; - } - if (*idigit >= 0) { - goto L80; - } - - ndigit = -(*idigit); - if (ndigit > 4) { - goto L20; - } - - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 10) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 9; - k2 = min(i__2,i__3); - io___6.ciunit = lout; - s_wsfe(&io___6); - do_fio(&c__1, (char*)&k1, (ftnlen)sizeof(integer)); - do_fio(&c__1, (char*)&k2, (ftnlen)sizeof(integer)); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, (char*)&ix[i__], (ftnlen)sizeof(integer)); - } - e_wsfe(); - /* L10: */ - } - return 0; - -L20: - if (ndigit > 6) { - goto L40; - } - - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 7) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 6; - k2 = min(i__2,i__3); - io___8.ciunit = lout; - s_wsfe(&io___8); - do_fio(&c__1, (char*)&k1, (ftnlen)sizeof(integer)); - do_fio(&c__1, (char*)&k2, (ftnlen)sizeof(integer)); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, (char*)&ix[i__], (ftnlen)sizeof(integer)); - } - e_wsfe(); - /* L30: */ - } - return 0; - -L40: - if (ndigit > 10) { - goto L60; - } - - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 5) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 4; - k2 = min(i__2,i__3); - io___9.ciunit = lout; - s_wsfe(&io___9); - do_fio(&c__1, (char*)&k1, (ftnlen)sizeof(integer)); - do_fio(&c__1, (char*)&k2, (ftnlen)sizeof(integer)); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, (char*)&ix[i__], (ftnlen)sizeof(integer)); - } - e_wsfe(); - /* L50: */ - } - return 0; - -L60: - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 3) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 2; - k2 = min(i__2,i__3); - io___10.ciunit = lout; - s_wsfe(&io___10); - do_fio(&c__1, (char*)&k1, (ftnlen)sizeof(integer)); - do_fio(&c__1, (char*)&k2, (ftnlen)sizeof(integer)); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, (char*)&ix[i__], (ftnlen)sizeof(integer)); - } - e_wsfe(); - /* L70: */ - } - return 0; - -L80: - if (ndigit > 4) { - goto L100; - } - - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 20) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 19; - k2 = min(i__2,i__3); - io___11.ciunit = lout; - s_wsfe(&io___11); - do_fio(&c__1, (char*)&k1, (ftnlen)sizeof(integer)); - do_fio(&c__1, (char*)&k2, (ftnlen)sizeof(integer)); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, (char*)&ix[i__], (ftnlen)sizeof(integer)); - } - e_wsfe(); - /* L90: */ - } - return 0; - -L100: - if (ndigit > 6) { - goto L120; - } - - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 15) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 14; - k2 = min(i__2,i__3); - io___12.ciunit = lout; - s_wsfe(&io___12); - do_fio(&c__1, (char*)&k1, (ftnlen)sizeof(integer)); - do_fio(&c__1, (char*)&k2, (ftnlen)sizeof(integer)); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, (char*)&ix[i__], (ftnlen)sizeof(integer)); - } - e_wsfe(); - /* L110: */ - } - return 0; - -L120: - if (ndigit > 10) { - goto L140; - } - - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 10) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 9; - k2 = min(i__2,i__3); - io___13.ciunit = lout; - s_wsfe(&io___13); - do_fio(&c__1, (char*)&k1, (ftnlen)sizeof(integer)); - do_fio(&c__1, (char*)&k2, (ftnlen)sizeof(integer)); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, (char*)&ix[i__], (ftnlen)sizeof(integer)); - } - e_wsfe(); - /* L130: */ - } - return 0; - -L140: - i__1 = *n; - for (k1 = 1; k1 <= i__1; k1 += 7) { - /* Computing MIN */ - i__2 = *n, i__3 = k1 + 6; - k2 = min(i__2,i__3); - io___14.ciunit = lout; - s_wsfe(&io___14); - do_fio(&c__1, (char*)&k1, (ftnlen)sizeof(integer)); - do_fio(&c__1, (char*)&k2, (ftnlen)sizeof(integer)); - i__2 = k2; - for (i__ = k1; i__ <= i__2; ++i__) { - do_fio(&c__1, (char*)&ix[i__], (ftnlen)sizeof(integer)); - } - e_wsfe(); - /* L150: */ - } - return 0; -} /* ivout_ */ - diff --git a/src/equil/vcs_setMolesLinProg.cpp b/src/equil/vcs_setMolesLinProg.cpp index 10e14e190..033d74d87 100644 --- a/src/equil/vcs_setMolesLinProg.cpp +++ b/src/equil/vcs_setMolesLinProg.cpp @@ -18,23 +18,10 @@ #include #include -#ifdef hpux -#define dbocls_ dbocls -#endif -#ifdef DEBUG_MODE -//extern int vcs_debug_print_lvl; -#endif #ifndef MAX #define MAX(x,y) (( (x) > (y) ) ? (x) : (y)) #endif - -extern "C" void dbocls_(double* W, int* MDW, int* MCON, int* MROWS, - int* NCOLS, - double* BL, double* BU, int* IND, int* IOPT, - double* X, double* RNORMC, double* RNORM, - int* MODE, double* RW, int* IW); - using namespace std; namespace VCSnonideal @@ -58,7 +45,6 @@ static void printProgress(const vector &spName, } #endif -#ifdef ALTLINPROG //! Estimate the initial mole numbers. /*! * This is done by running @@ -249,133 +235,4 @@ int VCS_SOLVE::vcs_setMolesLinProg() return retn; } -#else // ALTLINPROG - -int linprogmax(double* XMOLES, double* CC, double* AX, double* BB, - size_t NE, size_t M, size_t NE0) - -/*----------------------------------------------------------------------- -* Find XMOLES(I), i = 1, M such that -* Maximize CC dot W, subject to the NE constraints: -* -* [AX] [XMOLES] = [BB] -* and XMOLES(i) > 0 -* -* Input -* --------- -* AX(NE, M) - matrix of constraints AX(I,J) = ax(i + j*ne0) -* BB(NE) - contraint values -* CC(M) - Vector of "Good Values" to maximize -* -* Output -* --------- -* XMOLES(M) - optimal value of XMOLES() -*----------------------------------------------------------------------*/ -{ - int MROWS, MCON, NCOLS, NX, NI, MDW, i, j, MODE; - double sum, F[1], RNORMC, RNORM, *W, *BL, *BU, *RW, *X; - int* IND, *IW, *IOPT; - - MROWS = 1; - MCON = (int) NE; - NCOLS = (int) M; - MDW = MCON + NCOLS; - NX = 0; - NI = 0; - - sum = 0.0; - for (i = 0; i < NCOLS; i++) { - sum += fabs(CC[i]); - } - F[0] = sum * 1000.; - if (F[0] <= 0.0) { - F[0] = 1000.; - } - - BL = (double*) malloc(2*(NCOLS+MCON) * sizeof(double)); - BU = BL + (NCOLS+MCON); - IND = (int*) malloc((NCOLS+MCON) * sizeof(int)); - RW = (double*) malloc((6*NCOLS + 5*MCON) * sizeof(double)); - IW = (int*) malloc((2*NCOLS + 2*MCON) * sizeof(int)); - IOPT = (int*) malloc((17 + NI) * sizeof(int)); - X = (double*) malloc((2*(NCOLS+MCON) + 2 + NX) * sizeof(double)); - W = (double*) malloc((MDW*(NCOLS+MCON+1)) * sizeof(double)); - if (W == NULL) { - plogf("linproxmax ERROR: can not malloc memory of size %d bytes\n", - (int)((MDW*(NCOLS+MCON+1)) * sizeof(double))); - if (BL != NULL) { - free((void*) BL); - } - if (IND != NULL) { - free((void*) IND); - } - if (RW != NULL) { - free((void*) RW); - } - if (IW != NULL) { - free((void*) IW); - } - if (IOPT != NULL) { - free((void*) IOPT); - } - if (W != NULL) { - free((void*) W); - } - return -1; - } - for (j = 0; j < MCON; j++) { - for (i = 0; i < NCOLS; i++) { - W[j + i*MDW] = AX[j + i*NE0]; - } - } - for (i = 0; i < NCOLS; i++) { - W[MCON + i*MDW] = CC[i]; - } - W[MCON + (NCOLS)*MDW] = F[0]; - IOPT[0] = 99; - - for (j = 0; j < NCOLS; j++) { - IND[j] = 1; - BL[j] = 0.0; - BU[j] = 1.0e200; - } - for (j = 0; j < MCON; j++) { - IND[j + NCOLS] = 3; - BL[j + NCOLS] = BB[j]; - BU[j + NCOLS] = BL[j + NCOLS]; - } - - - dbocls_(W, &MDW, &MCON, &MROWS, &NCOLS, BL, BU, IND, IOPT, - X, &RNORMC, &RNORM, &MODE, RW, IW); - if (MODE != 0) { - plogf("Return from DBOCLS was not normal, MODE = %d\n", MODE); - plogf(" refer to subroutine DBOCLS for resolution\n"); - plogf(" RNORMC = %g\n", RNORMC); - } - - for (j = 0; j < NCOLS; j++) { - XMOLES[j] = X[j]; - } -#ifdef DEBUG_MODE - //sum = 0.0; - //for (j = 0; j < NCOLS; j++) { - // sum += XMOLES[j] * CC[j]; - //} - //if (vcs_debug_print_lvl >= 2) { - // plogf(" -- linmaxc: Final Maximized Value = %g\n", sum); - //} -#endif - - free((void*)W); - free((void*)BL); - free((void*)IND); - free((void*)RW); - free((void*)IW); - free((void*)IOPT); - free((void*)X); - - return 0; -} -#endif // ALTLINPROG } diff --git a/src/equil/vcs_xerror.c b/src/equil/vcs_xerror.c deleted file mode 100644 index 312b79885..000000000 --- a/src/equil/vcs_xerror.c +++ /dev/null @@ -1,255 +0,0 @@ -/* xerror.f -- translated by f2c (version 20031025). - You must link the resulting object file with libf2c: - on Microsoft Windows system, link with libf2c.lib; - on Linux or Unix systems, link with .../path/to/libf2c.a -lm - or, if you install libf2c.a in a standard place, with -lf2c -lm - -- in that order, at the end of the command line, as in - cc *.o -lf2c -lm - Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., - - http://www.netlib.org/f2c/libf2c.zip -*/ - -#include "f2c.h" - -/* Table of constant values */ - -static integer c__2 = 2; -static integer c__1 = 1; - -/* Subroutine */ -int s88fmt_(integer* n, integer* ivalue, integer* ifmt) -{ - /* Initialized data */ - - static struct { - char e_1[40]; - integer e_2; - } equiv_5 = { "0 1 2 3 4 5 6 7 8 9 ", 0 }; - -#define idigit ((integer *)&equiv_5) - - - static integer it, nt, index; - - /* ----------------------------------------------------------------------- */ - /* S88FMT replaces ifmt(1), ... ,ifmt(n) with the */ - /* characters corresponding to the n least significant */ - /* digits of ivalue. */ - /* Taken from the Bell laboratories port library error handler */ - /* latest revision --- 7 June 1978. */ - - /* Jones R.E., *SLATEC common mathematical library error handling */ - /* package*, SAND78-1189, Sandia Laboratories, 1978. */ - /* ----------------------------------------------------------------------- */ - /* Parameter adjustments */ - --ifmt; - - /* Function Body */ - nt = *n; - it = *ivalue; -L10: - if (nt == 0) { - return 0; - } - index = it % 10; - ifmt[nt] = idigit[index]; - it /= 10; - --nt; - goto L10; -} /* s88fmt_ */ - -#undef idigit - - - -/* ----------------------------------------------------------------------- */ - -/* Subroutine */ int xerrwv_(char* msg, integer* nmes, integer* nerr, integer - *iert, integer* ni, integer* i1, integer* i2, integer* nr, real* r1, - real* r2, ftnlen msg_len) -{ - /* Initialized data */ - - static struct { - char e_1[16]; - integer fill_2[2]; - char e_3[4]; - integer fill_4[2]; - char e_5[8]; - integer fill_6[1]; - char e_7[4]; - } equiv_20 = { "( 1 X , ", {0}, "A ", {0}, ", A ", {0}, - ") " - }; - -#define nform ((integer *)&equiv_20) - - static integer ncpw = 4; - static integer mesflg = 1; - static integer lunit = 6; - - /* Format strings */ - static char fmt_20[] = "(6x,\002IN ABOVE MESSAGE, I1 =\002,i10)"; - static char fmt_30[] = "(6x,\002IN ABOVE MESSAGE, I1 =\002,i10,3x,\002I" - "2 =\002,i10)"; - static char fmt_40[] = "(6x,\002IN ABOVE MESSAGE, R1 =\002,d21.13)"; - static char fmt_50[] = "(6x,\002IN ABOVE, R1 =\002,d21.13,3x,\002R2 " - "=\002,d21.13)"; - - /* System generated locals */ - integer i__1; - - /* Builtin functions */ - integer s_wsfe(cilist*), do_fio(integer*, char*, ftnlen), e_wsfe(void); - /* Subroutine */ - int s_stop(char*, ftnlen); - - /* Local variables */ - static integer i__, nch, lun, nrem, nwds; - extern /* Subroutine */ int s88fmt_(integer*, integer*, integer*); - - /* Fortran I/O blocks */ - static cilist io___14 = { 0, 0, 0, (char*)nform, 0 }; - static cilist io___16 = { 0, 0, 0, fmt_20, 0 }; - static cilist io___17 = { 0, 0, 0, fmt_30, 0 }; - static cilist io___18 = { 0, 0, 0, fmt_40, 0 }; - static cilist io___19 = { 0, 0, 0, fmt_50, 0 }; - - - /* ----------------------------------------------------------------------- */ - /* Subroutine xerrwv, as given here, constitutes */ - /* a simplified version of the slatec error handling package */ - /* written by A. C. Hindmarsh at LLL. version of January 23, 1980, */ - /* modified by L. R. Petzold, April 1982. */ - - /* All arguments are input arguments. */ - - /* msg = the message (Hollerith literal or integer array). */ - /* nmes = the length of msg (number of characters). */ - /* nerr = the error number (not used). */ - /* iert = the error type.. */ - /* 1 means recoverable (control returns to caller). */ - /* 2 means fatal (run is aborted--see note below). */ - /* ni = number of integers (0, 1, or 2) to be printed with message. */ - /* i1,i2 = integers to be printed, depending on ni. */ - /* nr = number of reals (0, 1, or 2) to be printed with message. */ - /* r1,r2 = reals to be printed, depending on ni. */ - - /* note.. this routine is machine-dependent and specialized for use */ - /* in limited context, in the following ways.. */ - /* 1. the number of hollerith characters stored per word, denoted */ - /* by ncpw below, is set in a data statement below. */ - /* 2. the value of nmes is assumed to be at most 60. */ - /* (multi-line messages are generated by repeated calls.) */ - /* 3. if iert = 2, control passes to the statement stop */ - /* to abort the run. this statement may be machine-dependent. */ - /* 4. r1 and r2 are assumed to be in real and are printed */ - /* in d21.13 format. */ - /* 5. the data statement below contains default values of */ - /* mesflg = print control flag.. */ - /* 1 means print all messages (the default). */ - /* 0 means no printing. */ - /* lunit = logical unit number for messages. */ - /* the default is 6 (machine-dependent). */ - /* to change lunit, change the data statement */ - /* below. */ - /* ----------------------------------------------------------------------- */ - /* The following are instructions for installing this routine */ - /* in different machine environments. */ - - /* To change the default output unit, change the data statement */ - /* below. */ - - /* For a different number of characters per word, change the */ - /* data statement setting ncpw below. */ - /* Alternatives for various computers are shown in comment */ - /* cards. */ - - /* For a different run-abort command, change the statement following */ - /* statement 100 at the end. */ - /* ----------------------------------------------------------------------- */ - /* The following value of ncpw is valid for the cdc-6600 and */ - /* cdc-7600 computers. */ - /* data ncpw/10/ */ - /* The following is valid for the cray-1 computer. */ - /* data ncpw/8/ */ - /* The following is valid for the burroughs 6700 and 7800 computers. */ - /* data ncpw/6/ */ - /* The following is valid for the pdp-10 computer. */ - /* data ncpw/5/ */ - /* The following is valid for the vax computer with 4 bytes per integer, */ - /* and for the ibm-360, ibm-303x, and ibm-43xx computers. */ - /* data ncpw/4/ */ - /* The following is valid for the pdp-11, or vax with 2-byte integers. */ - /* data ncpw/2/ */ - /* ---------------------------------------------------------------------- */ - /* Parameter adjustments */ - msg -= msg_len; - - /* Function Body */ - - /* This may have to be modified according to the information */ - /* provided above */ - - - if (mesflg == 0) { - goto L100; - } - lun = lunit; - nch = min(*nmes,60); - nwds = nch / ncpw; - s88fmt_(&c__2, &nwds, &nform[4]); - s88fmt_(&c__2, &ncpw, &nform[7]); - nrem = nch - nwds * ncpw; - if (nrem > 0) { - ++nwds; - } - if (nrem < 1) { - nrem = 1; - } - s88fmt_(&c__1, &nrem, &nform[11]); - io___14.ciunit = lun; - s_wsfe(&io___14); - i__1 = nwds; - for (i__ = 1; i__ <= i__1; ++i__) { - do_fio(&c__1, msg + i__ * msg_len, msg_len); - } - e_wsfe(); - if (*ni == 1) { - io___16.ciunit = lun; - s_wsfe(&io___16); - do_fio(&c__1, (char*)&(*i1), (ftnlen)sizeof(integer)); - e_wsfe(); - } - if (*ni == 2) { - io___17.ciunit = lun; - s_wsfe(&io___17); - do_fio(&c__1, (char*)&(*i1), (ftnlen)sizeof(integer)); - do_fio(&c__1, (char*)&(*i2), (ftnlen)sizeof(integer)); - e_wsfe(); - } - if (*nr == 1) { - io___18.ciunit = lun; - s_wsfe(&io___18); - do_fio(&c__1, (char*)&(*r1), (ftnlen)sizeof(real)); - e_wsfe(); - } - if (*nr == 2) { - io___19.ciunit = lun; - s_wsfe(&io___19); - do_fio(&c__1, (char*)&(*r1), (ftnlen)sizeof(real)); - do_fio(&c__1, (char*)&(*r2), (ftnlen)sizeof(real)); - e_wsfe(); - } -L100: - if (*iert != 2) { - return 0; - } - s_stop("", (ftnlen)0); - return 0; -} /* xerrwv_ */ - -#undef nform - -