Source Code 1992 March/Source_Code_CD-ROM_Walnut_Creek_March_1992.iso/usenet/compsrcs/misc/volume07/minit

Internet Message Format | 1991-08-27 | 24.1 KB

From decwrl!ucbvax!tut.cis.ohio-state.edu!brutus.cs.uiuc.edu!wuarchive!wugate!uunet!allbery Thu Aug  3 08:51:40 PDT 1989
Article 999 of comp.sources.misc:
Path: decwrl!ucbvax!tut.cis.ohio-state.edu!brutus.cs.uiuc.edu!wuarchive!wugate!uunet!allbery
From: [e-mail] (Brandon S. Allbery - comp.sources.misc)
Newsgroups: comp.sources.misc
Subject: v07i107: A Linear Programming Package for comp.sources.misc
Message-ID: <[e-mail]>
Date: 28 Jul 89 01:23:23 GMT
Sender: [e-mail]
Reply-To: [e-mail] (Badri Lokanathan )
Lines: 1042
Approved: [e-mail] (Brandon S. Allbery - comp.sources.misc)

Posting-number: Volume 7, Issue 107
Submitted-by: [e-mail] (Badri Lokanathan )
Archive-name: minit

I am enclosing the following Linear Programming package for submission
to comp.sources.misc. I keep seeing postings to comp.sources.wanted
for such a package, so I figure some folks out there will have use for it!

Badri Lokanathan
# Shell archive; cut here and unpack with /bin/sh
echo x - README
sed 's/^X//' >README <<'*-*-END-of-README-*-*'
XThis bundle contains a linear programming package based on the dual
Xsimplex method. The original code was given in Collected Algorithms
Xfrom CACM (1968) in algol 60 by Rudolfo Salazar and Subrata Sen. The
Xevaluators of this algorithm suggested some fixes that have also been
Xincorporated.
XI translated it into C and have been using it for sometime now without
Xany problems. The code consists of two parts:
X
X(1) A function called minit_space_handler to initialize arrays by
X    mallocing the appropriate amount of space, and minit, the actual
X    LP solver. These routines can be invoked directly by another
X    package. Both are in file minit.c
X
X(2) A front end to use the minit routines stand alone. Usage is
X    documented in the man page.
X
XIf a laser printer is available, the man page may be printed as
X    eqn minit.1 | troff -t -man | lpr -t -Pwhatever
XOtherwise
X    nroff -man minit.1
X
XThe implementation has been tested on a VAX 750 under BSD 4.2/4.3
Xand Sun-3's with OS 3.5 and 4.0.3.
X
XBadri Lokanathan
XDept. of Electrical Engineering
XUniversity of Rochester
*-*-END-of-README-*-*
echo x - Makefile
sed 's/^X//' >Makefile <<'*-*-END-of-Makefile-*-*'
XCFLAGS    = -O -c
X#DEFS    = -DDEBUG
XDEFS    =
X
X.c.o:
X    cc ${CFLAGS} ${DEFS} $*.c
X
Xminit: minit.o minmain.o
X    cc -O -o minit minit.o minmain.o
*-*-END-of-Makefile-*-*
echo x - minit.1
sed 's/^X//' >minit.1 <<'*-*-END-of-minit.1-*-*'
X.TH MINIT 1 "14 Dec 1988" "UR VLSI"
X.SH NAME
Xminit \- an efficient Linear Programming package
X.SH SYNOPSIS
X\fBminit\fP [\fB-p\fPprec] [\fB-i\fP[filename]]
X.SH DESCRIPTION
X.I Minit
Xsolves a Linear Programming problem of n variables and m constraints,
Xthe last p of which are equality constraints by the dual simplex method.
XThe problem statement is
X.sp 1
X.if t \{
X.EQ I
XMaximize ~~z ~=~ c bar sup T ~ x bar ~~~~subject ~~"to"
X.EN
X.EQ I
Xbold A x bar ~<=~ b, bar
X~~~~x bar ~>=~ 0 bar
X.EN
X\}
X.if n \{
X.RS
X.nf
XMaximize z = cX
Xsubject to
XAX <= b
XX >= 0,
Xwhere
X\(bu c is a (1*n) row vector
X\(bu X is a (n*1) column vector
X\(bu A is a (m*n) matrix
X\(bu b is a (m*1) column vector.
X.RE
X.fi
X\}
X.SH OPTIONS
X.IP "-p"
XThe following integer specifies the precision (number of decimal places)
Xof the output.
XThe default is free-format.
X.IP "-i"
XThe following argument is the name of the input file.
XIf no filename is specified,
Xminit goes into interactive mode and prompts for input.
X.LP
XIf the -i option is omitted,
Xminit takes input (in the appropriate format) from standard input.
X.SH EXAMPLE
XAn example with 4 variables, 2 inequality constraints and 1
Xequality constraint follows. The input corresponds to the problem
X.sp 1
X.if t \{
X.EQ I
XMaximize ~~6 x sub 1 ~+~ x sub 2 ~-~ x sub 3 ~-~ x sub 4
X~~~~subj. ~~"to"
X.EN
X.br
X.EQ I
Xx sub 1 ~+~ 2 x sub 2 ~+~ x sub 3 ~+~ x sub 4 mark ~<=~ 5
X.EN
X.br
X.EQ I
X3 x sub 1 ~+~ x sub 2 ~-~ x sub 3 lineup ~<=~ 8
X.EN
X.br
X.EQ I
Xx sub 2 ~+~ x sub 3 ~+~ x sub 4  lineup ~=~ 1
X.EN
X.br
X.EQ I
Xx sub i ~>=~ 0, ~~i ~=~ 1 ,..., 4
X.EN
X\}
X.if n \{
X.RS
X.nf
XMaximize 6x1 + x2 - x3 - x4 subj. to
X
Xx1 + 2x2 + x3 + x4 <= 5
X3x1 + x2 - x3      <= 8
X      x2 + x3 + x4  = 1
X
Xx1, x2, x3, x4 >= 0
X.RE
X.fi
X\}
X.SH DATA
X.nf
X4 2 1
X
X1 2  1 1
X3 1 -1 0
X0 1  1 1
X
X5 8 1
X6 1 -1 -1
X.fi
X.SH CAVEATS
XInput can have any amount of white space and integer or floating
Xpoint format. No missing values, however. Also, all equality
Xconstraints are to be specified after the inequality constraints.
X.SH REFERENCE
XRudolfo C. Salazar and Subrata K. Sen,
X``MINIT (MINimum ITerations) algorithm for Linear Programming,''
X\fICollected Algorithms from CACM,\fP
Xalgorithm #333, 1968
X.SH AUTHOR
XTranslated into C from Algol 60 by Badri Lokanathan,
XDept. of EE, University of Rochester
*-*-END-of-minit.1-*-*
echo x - minit.c
sed 's/^X//' >minit.c <<'*-*-END-of-minit.c-*-*'
X/****************************** MINIT ******************************************
X * This file contains procedures to solve a Linear Programming
X * problem of n variables and m constraints, the last p of which
X * are equality constraints by the dual simplex method.
X *
X * The code was originally in Algol 60 and was published in Collected
X * Algorithms from CACM (alg #333) by Rudolfo C. Salazar and Subrata
X * K. Sen in 1968 under the title MINIT (MINimum ITerations) algorithm
X * for Linear Programming.
X * It was directly translated into C by Badri Lokanathan, Dept. of EE,
X * University of Rochester, with no modification to code structure.
X *
X * The problem statement is
X * Maximise z = cX
X *
X * subject to
X * AX <= b
X * X >=0
X *
X * c is a (1*n) row vector
X * A is a (m*n) matrix
X * b is a (m*1) column vector.
X * e is a matrix with with (m+1) rows and lcol columns (lcol = m+n-p+1)
X *   and forms the initial tableau of the algorithm.
X * td is read into the procedure and should be a very small number,
X *   say 10**-8
X *
X * The condition of optimality is the non-negativity of
X * e[1,j] for j = 1 ... lcol-1 and of e[1,lcol] = 2 ... m+1.
X * If the e[i,j] values are greater than or equal to -td they are
X * considered to be non-negative. The value of td should reflect the
X * magnitude of the co-efficient matrix.
X *
X ******************************************************************************/
X
X/***************************** INCLUDES ***************************************/
X#include <stdio.h>
X/************************ CONSTANT DEFINITIONS ********************************/
X#ifndef LARGENUMBER
X#define LARGENUMBER 1e6
X#endif  LARGENUMBER
X#ifndef VLARGENUMBER
X#define VLARGENUMBER 1e8
X#endif  VLARGENUMBER
X#ifndef VSMALLNUMBER
X#define VSMALLNUMBER 1e-8
X#endif  VSMALLNUMBER
X/************************** MACRO DEFINITIONS  ********************************/
X#define ABS(A) (A > 0 ? A: -A)
X#define MAX(A,B) (A > B ? A : B)
X#define MIN(A,B) (A < B ? A : B)
X/************************ FILE-LIMITED VARIABLES ******************************/
Xstatic int k, m, n, p, lcol, L, im, jmin, jm, imax, debug;
Xstatic float td=VSMALLNUMBER;    /* This is application specific. */
Xstatic float gmin, phimax;
X/************************** FILE-LIMITED ARRAYS *******************************/
Xstatic int *imin, *jmax, *ind, *ind1, *chk;
Xstatic float **e, *thmin, *delmax;
X
X/********************************* ZX3LP ***************************************
X * This is a user-specific front end to minit.
X * In its present form, it has been written for compatibility with ifplan.
X * Note that zx3lp was originally a fortran subroutine. The weird return
X * values and flags (such as IER) are inherited from there.
X * Variable IA and arrays RW, IW are not necessary.
X ******************************************************************************/
Xzx3lp(A,IA,B,C,N,M1,M2,S,PSOL,DSOL,RW,IW,IER)
Xfloat **A, *B, *C, *PSOL, *DSOL, *RW, *S;
Xint *IA, *N, *M1, *M2, *IW, *IER;
X{
X    int err;
X
X    /* Allocate space for arrays and initialize globals. */
X    if(!minit_space_handler(*N,*M1+*M2,*M2))
X    {
X        (void) fprintf(stderr,"Space problems in LP handler.\n");
X        return;
X    }
X
X    /* Invoke minit here. */
X    err = minit(A,B,C,PSOL,DSOL,S);
X    switch(err)
X    {
X        case 1:
X        case 2:
X        *IER = 133;
X        break;
X
X        case 3:
X        *IER = 131;
X        break;
X
X        default:
X        *IER = 0;
X    }
X}
X
X/************************ minit_space_handler **********************************
X * This routine performs space maintenance and initializes global arrays.
X * It looks at n, m and p, which are static in this file, to look for
X * previously allocated space.
X ******************************************************************************/
Xint
Xminit_space_handler(N,M,P)
Xint N, M, P;
X{
X    register int i, j;
X    float **tmp;
X    static int MS, LS;        /* To keep track of array sizes. */
X    char *calloc(), *malloc(), *realloc();
X
X    /* Initialize globals. */
X    m = M;
X    n = N;
X    p = P;
X    lcol = m + n - p + 1;
X
X    if(LS == 0)
X    {
X        /* First pass. Allocate space for every array. */
X
X        MS = m + 1;
X        LS = m + n - p + 1;
X        if(!(jmax = (int *) calloc((unsigned) MS,sizeof(int))))
X            return(0);
X
X        if(!(ind1 = (int *) calloc((unsigned) MS,sizeof(int))))
X            return(0);
X
X        if(!(chk = (int *) calloc((unsigned) MS,sizeof(int))))
X            return(0);
X
X        if(!(delmax = (float *) calloc((unsigned) MS,sizeof(float))))
X            return(0);
X
X        if(!(imin = (int *) calloc((unsigned) LS,sizeof(int))))
X            return(0);
X
X        if(!(ind = (int *) calloc((unsigned) LS,sizeof(int))))
X            return(0);
X
X        if(!(thmin = (float *) calloc((unsigned) LS,sizeof(float))))
X            return(0);
X
X        if(!(tmp=e=(float **) calloc((unsigned) MS,sizeof(float *))))
X            return(0);
X
X        for(i = 0; i < MS; i++)
X            if(!(*(tmp++)=(float *) calloc((unsigned) LS,
X                    sizeof(float))))
X                return(0);
X
X        return(1);
X    }
X
X    if(M + 1 > MS)
X    {
X        /* Need to reallocate space. */
X        MS = M + 1;
X        if(!(jmax = (int *) realloc((char *)jmax,
X                (unsigned)(MS*sizeof(int)))))
X            return(0);
X
X        if(!(ind1 = (int *) realloc((char *) ind1,
X                (unsigned)(MS*sizeof(int)))))
X            return(0);
X
X        if(!(chk = (int *) realloc((char *) chk,
X                (unsigned)(MS*sizeof(int)))))
X            return(0);
X
X        if(!(delmax = (float *) realloc((char *) delmax,
X                (unsigned)(MS*sizeof(float)))))
X            return(0);
X
X        if(!(e=(float **) realloc((char *) e,
X                (unsigned)(MS*sizeof(float *)))))
X            return(0);
X
X        for(tmp = e + m, i = 0; i < MS; i++)
X            if(!(*(tmp++)=(float *)
X                    malloc((unsigned)(LS*sizeof(float)))))
X                return(0);
X    }
X
X    if(lcol > LS)
X    {
X        LS = lcol;
X        if(!(ind = (int *) realloc((char *) ind,
X                (unsigned)(LS*sizeof(int)))))
X            return(0);
X
X        if(!(imin = (int *) realloc((char *) imin,
X                (unsigned)(LS*sizeof(int)))))
X            return(0);
X
X        if(!(thmin = (float *) realloc((char *) thmin,
X                (unsigned)(LS*sizeof(float)))))
X            return(0);
X
X        for(tmp = e, i = 0; i < MS; tmp++, i++)
X            if(!(*tmp=(float *) realloc((char *) tmp,
X                    (unsigned)(LS*sizeof(float)))))
X                return(0);
X    }
X
X    /* Finally, initialize all arrays to 0. */
X    for(i = 0; i <= m; i++)
X    {
X        jmax[i] = ind1[i] = chk[i] = 0;
X        delmax[i] = 0.0;
X        for(j = 0; j < lcol; j++)
X            e[i][j] = 0.0;
X    }
X
X    for(j = 0; j < lcol; j++)
X    {
X        imin[j] = ind[j] = 0;
X        thmin[j] = 0.0;
X    }
X
X    return(1);
X}
X
X/********************************* MINIT ***************************************
X * This is the main procedure. It is invoked with the various matrices,
X * after space for other arrays has been generated elsewhere (in zx3lp.)
X * It returns
X * 0 if everything went OK and x, w, z as the solutions.
X * 1 if no solution existed.
X * 2 if primal had no feasible solutions.
X * 3 if primal had unbounded solutions.
X ******************************************************************************/
Xint
Xminit(A,b,c,x,w,z)
Xfloat **A, *b, *c, *x, *w, *z;
X{
X    register int i, j;
X    void tab();
X
X    /* Generate initial tableau. */
X    for(j = 0; j < n; j++)
X        e[0][j] = -c[j];
X
X    for(i = 0; i < m; i++)
X    {
X        for(j = 0; j < n; j++)
X            e[i+1][j] = A[i][j];
X    }
X
X    for(j = 0; j < m; j++)
X        e[j+1][lcol-1] = b[j];
X
X    for(i = 0; i < m - p; i++)
X        e[1+i][n+i] = 1.0;
X
X    /* Now begins the actual algorithm. */
X    for(i = 1; i < m+1; i++)
X        chk[i] = -1;    /* Indexing problem; Algol
X                 * version treates 0 as out
X                 * of bounds, in C we prefer -1.
X                 */
X
X    if(!p)    /* There are no equality constraints */
X        goto RCS;
X    else
X    {
X    if(phase1())
X        return(1);
X    }
X
X    RCS: L = 1; k = 1;
X
X    if(debug)
X        tab();    /* Print the tableau. */
X
X    for(j = 0; j < lcol - 1; j++)
X    {
X        if(e[0][j] < -td)
X            ind[(L++)-1] = j; /* ind[L] keeps track of the
X                       * columns in which e[0][j]
X                       * is negative.
X                       */
X    }
X
X    for(i = 1; i < m + 1; i++)
X    {
X        if(e[i][lcol-1] < -td)
X            ind1[(k++)-1] = i; /* ind1[k] keeps track of the
X                        * rows in which e[i][lcol]
X                        * is negative.
X                        */
X    }
X
X    if(L == 1)
X    {
X        if(k == 1) /* Results */
X            goto RESULTS;
X        else
X        {
X            if(k == 2)
X            {
X                for(j = 0; j < lcol - 1; j++)
X                {
X                    if(e[ind1[0]][j] < 0)
X                        goto R;
X                }
X                /* Primal problem has no feasible solutions. */
X                return(2);
X            }
X            else
X                goto R;
X        }
X    }
X    else
X    {
X        if(L == 2)
X        {
X            if(k == 1)
X            {
X                for(i = 1; i < m + 1; i++)
X                {
X                    if(e[i][ind[0]] > 0)
X                        goto C;
X                }
X
X                /* Primal problem is unbounded. */
X                return(3);
X            }
X            else
X                goto S;
X        }
X
X        if(k == 1)
X            goto C;
X        else
X            goto S;
X    }
X
X    R:
X    prophi();
X    if(rowtrans(imax,jm))
X        return(1);
X    goto RCS;
X
X    C:
X    progamma();
X    if(rowtrans(im,jmin))
X        return(1);
X    goto RCS;
X
X    S:
X    progamma();
X    prophi();
X    if(gmin == LARGENUMBER)
X    {
X        if(rowtrans(imax,jm))
X            return(1);
X        goto RCS;
X    }
X
X    if(phimax == - LARGENUMBER)
X    {
X        if(rowtrans(im,jmin))
X            return(1);
X        goto RCS;
X    }
X
X    if(ABS(phimax) > ABS(gmin))
X    {
X        if(rowtrans(imax,jm))
X            return(1);
X    }
X    else
X    {
X        if(rowtrans(im,jmin))
X            return(1);
X    }
X    goto RCS;
X
X    RESULTS:
X    /* Output results here. */
X    *z = e[0][lcol-1];
X    for(i = 0; i < n; i++)
X        x[i] = 0.0;
X
X    for(j = 0; j < m; j++)
X        x[j] = 0.0;
X
X    for(i = 1; i < m + 1; i++)
X    {
X        if(chk[i] >= n)
X            chk[i] = -1;
X
X        if(chk[i] > -1)
X            x[chk[i]] = e[i][lcol-1];
X    }
X
X    for(j = n; j < lcol - 1; j++)
X        w[j-n] = e[0][j];
X
X    return(0);
X}
X
X/****************************** rowtrans ***************************************
X * Performs the usual tableau transformations in a linear programming
X * problem, (p,q) being the pivotal element.
X * Returns the following error codes:
X * 0: Everything was OK.
X * 1: No solution.
X ******************************************************************************/
Xint
Xrowtrans(p,q)
Xint p, q;
X{
X    register int i, j;
X    float dummy;
X
X    if(p == -1 || q == -1) /* No solution. */
X        return(1);
X
X    dummy = e[p][q];
X
X    if(debug)
X        (void) printf("--\nPivot element is e[%d][%d] = %f\n",p,q,dummy);
X
X    for(j = 0; j < lcol; j++)
X        e[p][j] /= dummy;
X
X    for(i = 0; i < m + 1; i++)
X    {
X        if(i != p)
X        {
X            if(e[i][q] != 0.0)
X            {
X                dummy = e[i][q];
X                for(j = 0; j < lcol; j++)
X                    e[i][j] -= e[p][j] * dummy;
X            }
X        }
X    }
X
X    chk[p] = q;
X    return(0);
X} /* rowtrans */
X
X/****************************** progamma ***************************************
X * Performs calculations over columns to determine the pivot element.
X ******************************************************************************/
Xprogamma()
X{
X    float theta, gamma;
X    register int i, L1;
X
X    gmin = LARGENUMBER;    /* gmin is set to a large no. for
X                 * initialization purposes.
X                 */
X    jmin = -1;        /* Array indices in C start from 0 */
X
X    for(L1 = 0; L1 < L - 1; L1++)
X    {
X        imin[ind[L1]] = 0;
X        thmin[ind[L1]] = LARGENUMBER;
X        for(i = 1; i < m + 1; i++)
X        {
X            if(e[i][ind[L1]] > td && e[i][lcol-1] >= -td)
X            {
X                theta = e[i][lcol-1] / e[i][ind[L1]];
X                if(theta < thmin[ind[L1]])
X                {
X                    thmin[ind[L1]] = theta;
X                    imin[ind[L1]] = i;
X                }
X            }
X        }
X
X        if(thmin[ind[L1]] == LARGENUMBER)
X            gamma = VLARGENUMBER;
X        else
X            gamma = thmin[ind[L1]] * e[0][ind[L1]];
X
X        if(gamma < gmin)
X        {
X            gmin = gamma;
X            jmin = ind[L1];
X        }
X    }
X    if(jmin > -1)
X        im = imin[jmin];
X} /* progamma */
X
X/****************************** prophi *****************************************
X * Performs calculations over rows to determine the pivot element.
X ******************************************************************************/
Xprophi()
X{
X    float delta, phi;
X    register int j, k1;
X
X    phimax = - LARGENUMBER;    /* phimax is set to a small no. for
X                 * initialization purposes.
X                 */
X    imax = -1;        /* Array indices in C start from 0 */
X
X    for(k1 = 0; k1 < k - 1; k1++)
X    {
X        jmax[ind1[k1]] = 0;
X        delmax[ind1[k1]] = - LARGENUMBER;
X        for(j = 0; j < lcol - 1; j++)
X        {
X            if(e[ind1[k1]][j] < -td && e[0][j] >= -td)
X            {
X                delta = e[0][j] / e[ind1[k1]][j];
X                if(delta > delmax[ind1[k1]])
X                {
X                    delmax[ind1[k1]] = delta;
X                    jmax[ind1[k1]] = j;
X                }
X            }
X        }
X
X        if(delmax[ind1[k1]] == - LARGENUMBER)
X            phi = - VLARGENUMBER;
X        else
X            phi = delmax[ind1[k1]] * e[ind1[k1]][lcol-1];
X
X        if(phi > phimax)
X        {
X            phimax = phi;
X            imax = ind1[k1];
X        }
X    }
X    if(imax > -1)
X        jm = jmax[imax];
X} /* prophi */
X
X/****************************** phase1 *****************************************
X * Applied only to equality constraints if any.
X ******************************************************************************/
Xphase1()
X{
X    float theta, gamma;
X    register int i, j, r;
X    /* Fix suggested by Holmgren, Obradovic, Kolm. */
X    register int im1, jmin1, first;
X    void tab();
X
X    im1 = jmin1 = -1;
X    /* Fix suggested by Messham to allow negative coeffs. in
X     * equality constraints.
X     */
X    for(i = m - p + 1; i < m + 1; i++)
X    {
X        if(e[i][lcol - 1] < 0)
X        {
X            for(j = 0; j < lcol; j++)
X                e[i][j] = -e[i][j];
X        }
X    }
X
X    for(r = 0; r < p; r++)
X    {
X        gmin = LARGENUMBER;    /* gmin is set to a large no. for
X                     * initialization purposes.
X                     */
X        L = 1;
X        jmin = -1;
X        first = 1;
X        for(j = 0; j < n; j++)
X        {
X            thmin[j] = LARGENUMBER;
X            /* Fix suggested by Kolm and Dahlstrand */
X            /* if(e[0,j] < 0) */
X            if(e[0][j] < -td)
X                ind[(L++)-1] = j;
X        }
X
X    L1:    if(L == 1)
X        {
X            for(j = 0; j < n; j++)
X                ind[j] = j;
X
X            L = n + 1;
X        }
X
X        for(k = 0; k < L - 1; k++)
X        {
X            for(i = m - p + 1; i < m + 1; i++)
X                if(chk[i] == -1)
X                {
X                    /* Fix suggested by Kolm
X                     * and Dahlstrand
X                     */
X                    /* if(e[i][ind[k]] > 0.0) */
X                    if(e[i][ind[k]] > td)
X                    {
X                        if((theta = e[i][lcol-1] /
X                          e[i][ind[k]]) < thmin[ind[k]])
X                        {
X                            thmin[ind[k]] = theta;
X                            imin[ind[k]] = i;
X                        }
X                    }
X                }
X
X            /* Fix suggested by Obradovic overrides
X             * fixes suggested by Kolm and Dahstrand
X             * as well as Messham.
X             */
X            if(thmin[ind[k]] < LARGENUMBER)
X            {
X                gamma = thmin[ind[k]] * e[0][ind[k]];
X                if(gamma < gmin)
X                {
X                    gmin = gamma;
X                    jmin = ind[k];
X                }
X            }
X        }
X        if(jmin == -1)
X        {
X            if(first)
X            {
X                first = 0;
X                L = 1;
X                goto L1;
X            }
X            else
X                im = -1;
X        }
X        else
X            im = imin[jmin];
X
X        if(im == im1 && jmin == jmin1)
X        {
X            L = 1;
X            goto L1;
X        }
X
X        if(debug)
X            tab();    /* Print the tableau. */
X
X        if(rowtrans(im,jmin))
X            return(1);
X
X        im1 = im;
X        jmin1 = jmin;
X    }
X    return(0);
X} /* phase1 */
X
X/****************************** tab *****************************************
X * The following procedure is for debugging. It simply prints the
X * current tableau.
X ******************************************************************************/
Xstatic void
Xtab()
X{
X    register int i, j;
X
X    (void) printf("\n");
X
X    for(i = 0; i < 35; i++)
X        (void) printf("-");
X
X    (void) printf(" TABLEAU ");
X
X    for(i = 0; i < 35; i++)
X        (void) printf("-");
X
X    (void) printf("\n");
X
X    for(i = 0; i < m+1; i++)
X    {
X        for(j = 0; j < lcol; j++)
X            (void) printf("%6.3f ",e[i][j]);
X
X        (void) printf("\n");
X    }
X}
*-*-END-of-minit.c-*-*
echo x - minmain.c
sed 's/^X//' >minmain.c <<'*-*-END-of-minmain.c-*-*'
X/*
X * This is a front end to test the linear programming package
X * minit. This algorithm was written by Das and Salazar in algol.
X *
X * Solve the LP      T
X *        max C X  subj. to
X *              AX <= B
X *               X >= 0
X */
X
X#include <stdio.h>
X#include <string.h>
Xchar *calloc();
X
Xmain(argc,argv)
Xint argc;
Xchar **argv;
X{
X    int  ier,m,n,m1,m2;    /* Dimensions of various matrices */
X    float **A,**A1,*b,*c,*x,*w,z,tmp;
X    char format[64], *usage, *sprintf();
X    int suppress;
X    int int_part, frac_part;
X    float precision;
X    register int i,j,arg;
X    FILE *fopen(), *in;
X
X    n = m1 = m2 = 0;
X    in = stdin;
X    suppress = 1;    /* Suppress printing interactive messages. */
X    precision = 0.0;
X
X    usage = "\
XUsage:\tminit [-pprec] [-i[filename]]\n\
Xor:\tminit [-pprec] < filename\n\
XFor detailed help: minit -help\n";
X
X    for (arg = 1; arg < argc; arg++)
X    {
X        if (argv[arg][0] != '-' ||
X            (argv[arg][1] != 'p' && argv[arg][1] != 'i'))
X        {
X            if (!strcmp(argv[arg],"-help") || !strcmp(argv[arg],"help"))
X            (void) system("man minit");
X            else
X            (void) fputs(usage, stderr);
X            exit(1);
X        }
X
X        if (argv[arg][1] == 'p')
X            (void) sscanf(&(argv[arg][2]),"%f",&precision);
X        else    /* if argv[arg][1] == 'i') */
X        {
X            if (in != stdin)    /* in already assigned. */
X            {
X                (void) fputs("Only one input file permitted.\n",
X                    stderr);
X                exit(1);
X            }
X            if((in = fopen(&(argv[arg][2]),"r")) == NULL)
X            {
X                (void) fprintf(stderr,
X                    "No input file %s; interactive mode\n",
X                    &(argv[arg][2]));
X
X                in = stdin;
X                suppress = 0;
X            }
X        }
X    }
X
X    if (!suppress)
X        (void) fputs("Enter value of n \(No. of unknowns\): ", stderr);
X
X    (void) fscanf(in,"%d",&n);
X
X    if (!suppress)
X        (void) fputs("Enter no. of inequality constraints: ", stderr);
X
X    (void) fscanf(in,"%d",&m1);
X
X    if (!suppress)
X        (void) fputs("Enter no. of equality constraints: ", stderr);
X
X    (void) fscanf(in,"%d",&m2);
X
X    m = m1+m2;
X
X    /* Allocate space for the arrays */
X        A1 = A     = (float **) calloc((unsigned)m,sizeof(float *));
X    for(i=0;i<m;i++)
X        *(A1++) = (float *) calloc((unsigned)n,sizeof(float));
X
X    b     = (float *) calloc((unsigned)m,sizeof(float));
X    c     = (float *) calloc((unsigned)n,sizeof(float));
X    x     = (float *) calloc((unsigned)n,sizeof(float));
X    w     = (float *) calloc((unsigned)m,sizeof(float));
X
X    if (!suppress)
X        (void) fprintf(stderr,"Enter \(%d by %d\) elements of A:\n",
X            m,n);
X    for(i=0;i<m;i++)
X        for(j=0;j<n;j++)
X        {
X            (void) fscanf(in,"%f",&tmp);
X            A[i][j] = tmp;;
X        }
X
X    if (!suppress)
X        (void) fprintf(stderr,
X            "Enter %d elements of b, entered as a linear array:\n ",
X                m);
X    for(i=0;i<m;i++)
X    {
X        (void) fscanf(in,"%f",&tmp);
X        b[i] = tmp;
X    }
X
X    if (!suppress)
X        (void) fprintf(stderr,
X            "Enter %d elements of c, entered as a linear array:\n",
X                n);
X    for(j=0;j<n;j++)
X    {
X        (void) fscanf(in,"%f",&tmp);
X        c[j] = tmp;
X    }
X
X    if (suppress)
X    {
X        (void) fputs("INPUT DATA:\n",stderr);
X        (void) fprintf(stderr,"No. of variables \(n\) = %d\n", n);
X        (void) fprintf(stderr,
X            "No. of inequality constraints \(m1\) = %d\n", m1);
X        (void) fprintf(stderr,
X            "No. of equality constraints \(m2\) = %d\n", m2);
X    }
X
X    (void) fputs("------------------------------",stderr);
X    (void) fputs(" EXECUTING MINIT ",stderr);
X    (void) fputs("------------------------------",stderr);
X    (void) fputc('\n',stderr);
X
X    zx3lp(A,(int *)NULL,b,c,&n,&m1,&m2,&z,x,w,(float *)NULL,
X        (int *)NULL,&ier);
X    switch(ier)
X    {
X        case 133:
X        break;
X
X        case 1:
X        (void) fprintf(stderr,"No solution.\n");
X        exit(1);
X
X        case 131:
X        (void) fprintf(stderr,"Primal has unbounded solution.\n");
X        exit(3);
X    }
X
X    /* Make sure precision is acceptible. */
X    int_part = precision;
X    frac_part = precision - int_part;
X    if (int_part <= frac_part || int_part > 25)
X        precision = 10.3;    /* The default precision. */
X
X    (void) sprintf(format,"\nValue of z: %%%gg\n",precision);
X    (void) printf(format,z);
X    (void) puts("Primal solution:");
X
X    (void) sprintf(format,"%%%gg ",precision);
X    for(j=0;j<n;j++)
X        (void) printf(format,x[j]);
X
X    (void) puts("\nDual solution:");
X    for(i=0;i<m;i++)
X        (void) printf(format,w[i]);
X
X    (void) putchar('\n');
X    exit(0);
X}
*-*-END-of-minmain.c-*-*
echo x - example
sed 's/^X//' >example <<'*-*-END-of-example-*-*'
X4 2 1
X
X1 2  1 1
X3 1 -1 0
X0 1  1 1
X
X5 8 1
X6 1 -1 -1
*-*-END-of-example-*-*
exit