diff --git a/Cantera/src/equil/Makefile.in b/Cantera/src/equil/Makefile.in index 1723a95c0..afdb336d3 100644 --- a/Cantera/src/equil/Makefile.in +++ b/Cantera/src/equil/Makefile.in @@ -9,7 +9,7 @@ ############################################################### .SUFFIXES : -.SUFFIXES : .cpp .d .o .h +.SUFFIXES : .cpp .c .d .o .h INCDIR = ../../../build/include/cantera/kernel INSTALL_TSC = ../../../bin/install_tsc @@ -17,6 +17,10 @@ do_ranlib = @DO_RANLIB@ do_VCSnonideal = @COMPILE_VCSNONIDEAL@ +# +# Decide whether ot use the linear programmaing module lookalike +# (from Goodwin's code). Both should produce the same results. +do_AltLinProg = 1 debug_mode = @CANTERA_DEBUG_MODE@ ifeq ($(debug_mode), 1) @@ -35,7 +39,13 @@ endif # #LOCAL_DEFS=-DDEBUG_BASISOPTIMIZE # -LOCAL_DEFS=-DDEBUG_BASISOPTIMIZE -DDEBUG_CHEMEQUIL -DCANTERA_SRC_TREE -DALTLINPROG -DDEBUG + +ifeq ($(do_AltLinProg), 1) +DALT_STR=-DALTLINPROG +endif +LOCAL_DEFS=-DDEBUG_BASISOPTIMIZE -DDEBUG_CHEMEQUIL \ + -DDEBUG $(DALT_STR) +# # PIC_FLAG=@PIC@ @@ -46,6 +56,15 @@ EQUIL_OBJ = BasisOptimize.o ChemEquil.o MultiPhase.o MultiPhaseEquil.o equilibra EQUIL_H = ChemEquil.h MultiPhase.h MultiPhaseEquil.h equil.h PropertyCalculator.h +# depending on the option, different object files are used +ifeq ($(do_AltLinProg), 1) +DALT_OBJ=vcs_setMolesLinProg.o +else +DALT_OBJ=vcs_linmaxc.o vcs_dbocls.o vcs_dbols.o vcs_dbolsm.o vcs_dmout.o \ + vcs_dvout.o vcs_ivout.o vcs_xerror.o +endif + + ifeq ($(do_VCSnonideal), 1) VCSNONIDEAL_OBJ = vcs_solve_TP.o vcs_VolPhase.o vcs_solve.o vcs_prob.o \ vcs_TP.o vcs_TV.o vcs_report.o vcs_util.o \ @@ -54,8 +73,10 @@ VCSNONIDEAL_OBJ = vcs_solve_TP.o vcs_VolPhase.o vcs_solve.o vcs_prob.o \ vcs_nasa_poly.o vcs_nondim.o vcs_Exception.o \ vcs_funcVtot.o vcs_inest.o vcs_rearrange.o \ vcs_root1d.o vcs_rxnadj.o vcs_timer_generic.o \ - vcs_SpeciesProperties.o vcs_setMolesLinProg.o \ - vcs_prep.o vcs_species_thermo.o vcs_Gibbs.o + vcs_SpeciesProperties.o \ + vcs_prep.o vcs_species_thermo.o vcs_Gibbs.o \ + $(DALT_OBJ) + VCSNONIDEAL_H = vcs_internal.h vcs_VolPhase.h vcs_solve.h vcs_prob.h \ vcs_IntStarStar.h vcs_DoubleStarStar.h vcs_defs.h \ vcs_MultiPhaseEquil.h vcs_nasa_poly.h vcs_Exception.h \ @@ -77,12 +98,18 @@ all: $(LIB) $(VLIB) .depends $(INSTALL_TSC) "$${lh}" $(INCDIR) ; \ done) -%.d: Makefile %.o - @CXX_DEPENDS@ $(CXX_FLAGS) $(CXX_INCLUDES) $*.cpp > $*.d +.cpp.d: Makefile %.o + g++ -MM $(CXX_FLAGS) $(CXX_INCLUDES) $*.cpp > $*.d + +.c.d: Makefile %.o + g++ -MM $(CXX_FLAGS) $(CXX_INCLUDES) -I../../../ext/f2c_libs $*.c > $*.d .cpp.o: @CXX@ -c $< $(CXX_FLAGS) $(CXX_INCLUDES) +.c.o: + @CC@ -c $< $(CXX_FLAGS) $(CXX_INCLUDES) -I../../../ext/f2c_libs + $(LIB): $(EQUIL_OBJ) $(EQUIL_H) $(VCSNONIDEAL_OBJ) $(VCSNONIDEAL_H) @ARCHIVE@ $(LIB) $(EQUIL_OBJ) $(VCSNONIDEAL_OBJ) > /dev/null ifeq ($(do_ranlib),1) @@ -96,9 +123,8 @@ ifeq ($(do_ranlib),1) endif - clean: - @(for lh in dummy.h $(EQUIL_H) ; do \ + @(for lh in dummy.h $(EQUIL_H) $(VCSNONIDEAL_H); do \ th=$(INCDIR)/"$${lh}" ; \ if test -f "$${th}" ; then \ $(RM) "$${th}" ; \ diff --git a/Cantera/src/equil/vcs_MultiPhaseEquil.h b/Cantera/src/equil/vcs_MultiPhaseEquil.h index 481717de6..3063e4373 100644 --- a/Cantera/src/equil/vcs_MultiPhaseEquil.h +++ b/Cantera/src/equil/vcs_MultiPhaseEquil.h @@ -13,8 +13,7 @@ #ifndef VCS_MULTIPHASEEQUIL_H #define VCS_MULTIPHASEEQUIL_H - -#ifdef CANTERA_NOTIN_SRC_TREE +#ifdef CANTERA_APP #include "cantera/kernel/ct_defs.h" #include "cantera/kernel/MultiPhase.h" #else diff --git a/Cantera/src/equil/vcs_dbocls.c b/Cantera/src/equil/vcs_dbocls.c index 2ab550e45..65765ebc5 100644 --- a/Cantera/src/equil/vcs_dbocls.c +++ b/Cantera/src/equil/vcs_dbocls.c @@ -719,7 +719,8 @@ static integer c__4 = 4; drelpr = d1mach_(&c__4); checkl = FALSE_; filter = TRUE_; - lenx = (*ncols + *mcon << 1) + 2; + /*lenx = (*ncols + *mcon << 1) + 2;*/ + lenx = ((*ncols + *mcon) << 1) + 2; iscale = 1; igo = 1; accum = FALSE_; @@ -1302,11 +1303,13 @@ L230: liw = 1; /* SET THE NEW TRIANGULARIZATION FACTOR. */ - x[(*ncols + *mcon << 1) + 1] = zero; + /*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; */ + 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[ diff --git a/Cantera/src/equil/vcs_dbolsm.c b/Cantera/src/equil/vcs_dbolsm.c new file mode 100644 index 000000000..f1898d298 --- /dev/null +++ b/Cantera/src/equil/vcs_dbolsm.c @@ -0,0 +1,1574 @@ +/* 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_ */ + diff --git a/Cantera/src/equil/vcs_dmout.c b/Cantera/src/equil/vcs_dmout.c new file mode 100644 index 000000000..574f8ea61 --- /dev/null +++ b/Cantera/src/equil/vcs_dmout.c @@ -0,0 +1,476 @@ +/* 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/Cantera/src/equil/vcs_inest.cpp b/Cantera/src/equil/vcs_inest.cpp index ec43793e7..31d2c2d7d 100644 --- a/Cantera/src/equil/vcs_inest.cpp +++ b/Cantera/src/equil/vcs_inest.cpp @@ -24,11 +24,6 @@ namespace VCSnonideal { static char pprefix[20] = " --- vcs_inest: "; - -#ifdef ALTLINPROG -#else - extern int linprogmax(double *, double *, double *, double *, int, int, int); -#endif /*****************************************************************************/ /*****************************************************************************/ /*****************************************************************************/ diff --git a/Cantera/src/equil/vcs_ivout.c b/Cantera/src/equil/vcs_ivout.c new file mode 100644 index 000000000..ef4cbbcd1 --- /dev/null +++ b/Cantera/src/equil/vcs_ivout.c @@ -0,0 +1,298 @@ +/* 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/Cantera/src/equil/vcs_linmaxc.cpp b/Cantera/src/equil/vcs_linmaxc.cpp index 1512b317d..a209eb56e 100644 --- a/Cantera/src/equil/vcs_linmaxc.cpp +++ b/Cantera/src/equil/vcs_linmaxc.cpp @@ -21,6 +21,7 @@ #endif #include "vcs_internal.h" +#include "vcs_solve.h" extern "C" void dbocls_(double *W, int *MDW, int *MCON, int *MROWS, int *NCOLS, @@ -34,6 +35,8 @@ extern "C" void dbocls_(double *W, int *MDW, int *MCON, int *MROWS, ******************************************************************************/ namespace VCSnonideal { +#ifdef ALTLINPROG +#else int linprogmax(double *XMOLES, double *CC, double *AX, double *BB, int NE, int M, int NE0) @@ -146,5 +149,5 @@ int linprogmax(double *XMOLES, double *CC, double *AX, double *BB, return 0; } - +#endif } diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 04c1ebc06..c6fd4d963 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -802,6 +802,11 @@ public: }; +#ifdef ALTLINPROG +#else +int linprogmax(double *, double *, double *, double *, int, int, int); +#endif + } #endif diff --git a/Cantera/src/equil/vcs_xerror.c b/Cantera/src/equil/vcs_xerror.c new file mode 100644 index 000000000..a7d98d18a --- /dev/null +++ b/Cantera/src/equil/vcs_xerror.c @@ -0,0 +1,261 @@ +/* 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; + +/* ======================================================================= */ +/* -------------------------------------------------- */ +/* | CVS Head Information | */ +/* -------------------------------------------------- */ +/* $Author$ */ +/* $Date$ */ +/* $Revision$ */ +/* ======================================================================= */ + +/* 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 + +