research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
cantera/Cantera/src/equil/vcs_dbolsm.c
Harry Moffat c96226700e got the Linear programming routines to work within the Cantera 
directory. I had forgot to include several files. You can now turn
these on via a Makefile flag
2008-01-03 21:38:25 +00:00

1574 lines
48 KiB
C
Raw Blame History

/* 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 integer c_n4 = -4;
static integer c__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_(integer *);
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 = bl[j];
rdum = 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_(&c__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_(&c__4));
tolsze = sqrt(d1mach_(&c__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_(&c__4)) {
nerr = 25;
nlevel = 0;
nchar = 88;
rdum2 = tolind;
rdum = d1mach_(&c__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 = 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 = 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 = 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_ */
Reference in a new issue View git blame Copy permalink