home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / misc / volume7 / minit < prev    next >
Text File  |  1989-07-27  |  24KB  |  1,049 lines
  1. Newsgroups: comp.sources.misc
  2. subject: v07i107: A Linear Programming Package for comp.sources.misc
  3. From: allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc)
  4. Reply-To: badri@ee.rochester.edu (Badri Lokanathan )
  5.  
  6. Posting-number: Volume 7, Issue 107
  7. Submitted-by: badri@ee.rochester.edu (Badri Lokanathan )
  8. Archive-name: minit
  9.  
  10. I am enclosing the following Linear Programming package for submission
  11. to comp.sources.misc. I keep seeing postings to comp.sources.wanted
  12. for such a package, so I figure some folks out there will have use for it!
  13.  
  14. Badri Lokanathan
  15. # Shell archive; cut here and unpack with /bin/sh
  16. echo x - README
  17. sed 's/^X//' >README <<'*-*-END-of-README-*-*'
  18. XThis bundle contains a linear programming package based on the dual
  19. Xsimplex method. The original code was given in Collected Algorithms
  20. Xfrom CACM (1968) in algol 60 by Rudolfo Salazar and Subrata Sen. The
  21. Xevaluators of this algorithm suggested some fixes that have also been
  22. Xincorporated.
  23. XI translated it into C and have been using it for sometime now without
  24. Xany problems. The code consists of two parts:
  25. X
  26. X(1) A function called minit_space_handler to initialize arrays by
  27. X    mallocing the appropriate amount of space, and minit, the actual
  28. X    LP solver. These routines can be invoked directly by another
  29. X    package. Both are in file minit.c
  30. X
  31. X(2) A front end to use the minit routines stand alone. Usage is
  32. X    documented in the man page.
  33. X    
  34. XIf a laser printer is available, the man page may be printed as
  35. X    eqn minit.1 | troff -t -man | lpr -t -Pwhatever
  36. XOtherwise
  37. X    nroff -man minit.1
  38. X
  39. XThe implementation has been tested on a VAX 750 under BSD 4.2/4.3
  40. Xand Sun-3's with OS 3.5 and 4.0.3.
  41. X
  42. XBadri Lokanathan
  43. XDept. of Electrical Engineering
  44. XUniversity of Rochester
  45. *-*-END-of-README-*-*
  46. echo x - Makefile
  47. sed 's/^X//' >Makefile <<'*-*-END-of-Makefile-*-*'
  48. XCFLAGS    = -O -c
  49. X#DEFS    = -DDEBUG
  50. XDEFS    =
  51. X
  52. X.c.o:
  53. X    cc ${CFLAGS} ${DEFS} $*.c
  54. X
  55. Xminit: minit.o minmain.o
  56. X    cc -O -o minit minit.o minmain.o
  57. *-*-END-of-Makefile-*-*
  58. echo x - minit.1
  59. sed 's/^X//' >minit.1 <<'*-*-END-of-minit.1-*-*'
  60. X.TH MINIT 1 "14 Dec 1988" "UR VLSI"
  61. X.SH NAME
  62. Xminit \- an efficient Linear Programming package
  63. X.SH SYNOPSIS
  64. X\fBminit\fP [\fB-p\fPprec] [\fB-i\fP[filename]]
  65. X.SH DESCRIPTION
  66. X.I Minit
  67. Xsolves a Linear Programming problem of n variables and m constraints,
  68. Xthe last p of which are equality constraints by the dual simplex method.
  69. XThe problem statement is
  70. X.sp 1
  71. X.if t \{
  72. X.EQ I
  73. XMaximize ~~z ~=~ c bar sup T ~ x bar ~~~~subject ~~"to"
  74. X.EN
  75. X.EQ I
  76. Xbold A x bar ~<=~ b, bar
  77. X~~~~x bar ~>=~ 0 bar
  78. X.EN
  79. X\}
  80. X.if n \{
  81. X.RS
  82. X.nf
  83. XMaximize z = cX
  84. Xsubject to
  85. XAX <= b
  86. XX >= 0,
  87. Xwhere
  88. X\(bu c is a (1*n) row vector
  89. X\(bu X is a (n*1) column vector
  90. X\(bu A is a (m*n) matrix
  91. X\(bu b is a (m*1) column vector.
  92. X.RE
  93. X.fi
  94. X\}
  95. X.SH OPTIONS
  96. X.IP "-p"
  97. XThe following integer specifies the precision (number of decimal places)
  98. Xof the output.
  99. XThe default is free-format.
  100. X.IP "-i"
  101. XThe following argument is the name of the input file.
  102. XIf no filename is specified,
  103. Xminit goes into interactive mode and prompts for input.
  104. X.LP
  105. XIf the -i option is omitted,
  106. Xminit takes input (in the appropriate format) from standard input. 
  107. X.SH EXAMPLE
  108. XAn example with 4 variables, 2 inequality constraints and 1
  109. Xequality constraint follows. The input corresponds to the problem
  110. X.sp 1
  111. X.if t \{
  112. X.EQ I
  113. XMaximize ~~6 x sub 1 ~+~ x sub 2 ~-~ x sub 3 ~-~ x sub 4
  114. X~~~~subj. ~~"to"
  115. X.EN
  116. X.br
  117. X.EQ I
  118. Xx sub 1 ~+~ 2 x sub 2 ~+~ x sub 3 ~+~ x sub 4 mark ~<=~ 5
  119. X.EN
  120. X.br
  121. X.EQ I
  122. X3 x sub 1 ~+~ x sub 2 ~-~ x sub 3 lineup ~<=~ 8
  123. X.EN
  124. X.br
  125. X.EQ I
  126. Xx sub 2 ~+~ x sub 3 ~+~ x sub 4  lineup ~=~ 1
  127. X.EN
  128. X.br
  129. X.EQ I
  130. Xx sub i ~>=~ 0, ~~i ~=~ 1 ,..., 4
  131. X.EN
  132. X\}
  133. X.if n \{
  134. X.RS
  135. X.nf
  136. XMaximize 6x1 + x2 - x3 - x4 subj. to
  137. X
  138. Xx1 + 2x2 + x3 + x4 <= 5
  139. X3x1 + x2 - x3      <= 8
  140. X      x2 + x3 + x4  = 1
  141. X
  142. Xx1, x2, x3, x4 >= 0
  143. X.RE
  144. X.fi
  145. X\}
  146. X.SH DATA
  147. X.nf
  148. X4 2 1
  149. X
  150. X1 2  1 1
  151. X3 1 -1 0
  152. X0 1  1 1
  153. X
  154. X5 8 1
  155. X6 1 -1 -1
  156. X.fi
  157. X.SH CAVEATS
  158. XInput can have any amount of white space and integer or floating
  159. Xpoint format. No missing values, however. Also, all equality
  160. Xconstraints are to be specified after the inequality constraints.
  161. X.SH REFERENCE
  162. XRudolfo C. Salazar and Subrata K. Sen,
  163. X``MINIT (MINimum ITerations) algorithm for Linear Programming,''
  164. X\fICollected Algorithms from CACM,\fP
  165. Xalgorithm #333, 1968
  166. X.SH AUTHOR
  167. XTranslated into C from Algol 60 by Badri Lokanathan,
  168. XDept. of EE, University of Rochester
  169. *-*-END-of-minit.1-*-*
  170. echo x - minit.c
  171. sed 's/^X//' >minit.c <<'*-*-END-of-minit.c-*-*'
  172. X/****************************** MINIT ******************************************
  173. X * This file contains procedures to solve a Linear Programming
  174. X * problem of n variables and m constraints, the last p of which
  175. X * are equality constraints by the dual simplex method.
  176. X *
  177. X * The code was originally in Algol 60 and was published in Collected
  178. X * Algorithms from CACM (alg #333) by Rudolfo C. Salazar and Subrata
  179. X * K. Sen in 1968 under the title MINIT (MINimum ITerations) algorithm
  180. X * for Linear Programming.
  181. X * It was directly translated into C by Badri Lokanathan, Dept. of EE,
  182. X * University of Rochester, with no modification to code structure. 
  183. X *
  184. X * The problem statement is
  185. X * Maximise z = cX
  186. X *
  187. X * subject to
  188. X * AX <= b
  189. X * X >=0
  190. X *
  191. X * c is a (1*n) row vector
  192. X * A is a (m*n) matrix
  193. X * b is a (m*1) column vector.
  194. X * e is a matrix with with (m+1) rows and lcol columns (lcol = m+n-p+1)
  195. X *   and forms the initial tableau of the algorithm.
  196. X * td is read into the procedure and should be a very small number,
  197. X *   say 10**-8
  198. X * 
  199. X * The condition of optimality is the non-negativity of
  200. X * e[1,j] for j = 1 ... lcol-1 and of e[1,lcol] = 2 ... m+1.
  201. X * If the e[i,j] values are greater than or equal to -td they are
  202. X * considered to be non-negative. The value of td should reflect the 
  203. X * magnitude of the co-efficient matrix.
  204. X *
  205. X ******************************************************************************/
  206. X
  207. X/***************************** INCLUDES ***************************************/
  208. X#include <stdio.h>
  209. X/************************ CONSTANT DEFINITIONS ********************************/
  210. X#ifndef LARGENUMBER
  211. X#define LARGENUMBER 1e6
  212. X#endif  LARGENUMBER
  213. X#ifndef VLARGENUMBER
  214. X#define VLARGENUMBER 1e8
  215. X#endif  VLARGENUMBER
  216. X#ifndef VSMALLNUMBER
  217. X#define VSMALLNUMBER 1e-8
  218. X#endif  VSMALLNUMBER
  219. X/************************** MACRO DEFINITIONS  ********************************/
  220. X#define ABS(A) (A > 0 ? A: -A)
  221. X#define MAX(A,B) (A > B ? A : B)
  222. X#define MIN(A,B) (A < B ? A : B)
  223. X/************************ FILE-LIMITED VARIABLES ******************************/
  224. Xstatic int k, m, n, p, lcol, L, im, jmin, jm, imax, debug;
  225. Xstatic float td=VSMALLNUMBER;    /* This is application specific. */
  226. Xstatic float gmin, phimax;
  227. X/************************** FILE-LIMITED ARRAYS *******************************/
  228. Xstatic int *imin, *jmax, *ind, *ind1, *chk;
  229. Xstatic float **e, *thmin, *delmax;
  230. X
  231. X/********************************* ZX3LP ***************************************
  232. X * This is a user-specific front end to minit.
  233. X * In its present form, it has been written for compatibility with ifplan. 
  234. X * Note that zx3lp was originally a fortran subroutine. The weird return
  235. X * values and flags (such as IER) are inherited from there.
  236. X * Variable IA and arrays RW, IW are not necessary.
  237. X ******************************************************************************/
  238. Xzx3lp(A,IA,B,C,N,M1,M2,S,PSOL,DSOL,RW,IW,IER)
  239. Xfloat **A, *B, *C, *PSOL, *DSOL, *RW, *S;
  240. Xint *IA, *N, *M1, *M2, *IW, *IER;
  241. X{
  242. X    int err;
  243. X
  244. X    /* Allocate space for arrays and initialize globals. */
  245. X    if(!minit_space_handler(*N,*M1+*M2,*M2))
  246. X    {
  247. X        (void) fprintf(stderr,"Space problems in LP handler.\n");
  248. X        return;
  249. X    }
  250. X
  251. X    /* Invoke minit here. */
  252. X    err = minit(A,B,C,PSOL,DSOL,S);
  253. X    switch(err)
  254. X    {
  255. X        case 1:
  256. X        case 2:
  257. X        *IER = 133;
  258. X        break;
  259. X
  260. X        case 3:
  261. X        *IER = 131;
  262. X        break;
  263. X
  264. X        default:
  265. X        *IER = 0;
  266. X    }
  267. X}
  268. X
  269. X/************************ minit_space_handler **********************************
  270. X * This routine performs space maintenance and initializes global arrays.     
  271. X * It looks at n, m and p, which are static in this file, to look for
  272. X * previously allocated space.
  273. X ******************************************************************************/
  274. Xint
  275. Xminit_space_handler(N,M,P)
  276. Xint N, M, P;
  277. X{
  278. X    register int i, j;
  279. X    float **tmp;
  280. X    static int MS, LS;        /* To keep track of array sizes. */
  281. X    char *calloc(), *malloc(), *realloc();
  282. X
  283. X    /* Initialize globals. */
  284. X    m = M;
  285. X    n = N;
  286. X    p = P;
  287. X    lcol = m + n - p + 1;
  288. X
  289. X    if(LS == 0)
  290. X    {
  291. X        /* First pass. Allocate space for every array. */
  292. X
  293. X        MS = m + 1;
  294. X        LS = m + n - p + 1;
  295. X        if(!(jmax = (int *) calloc((unsigned) MS,sizeof(int))))
  296. X            return(0);
  297. X
  298. X        if(!(ind1 = (int *) calloc((unsigned) MS,sizeof(int))))
  299. X            return(0);
  300. X
  301. X        if(!(chk = (int *) calloc((unsigned) MS,sizeof(int))))
  302. X            return(0);
  303. X
  304. X        if(!(delmax = (float *) calloc((unsigned) MS,sizeof(float))))
  305. X            return(0);
  306. X
  307. X        if(!(imin = (int *) calloc((unsigned) LS,sizeof(int))))
  308. X            return(0);
  309. X
  310. X        if(!(ind = (int *) calloc((unsigned) LS,sizeof(int))))
  311. X            return(0);
  312. X
  313. X        if(!(thmin = (float *) calloc((unsigned) LS,sizeof(float))))
  314. X            return(0);
  315. X
  316. X        if(!(tmp=e=(float **) calloc((unsigned) MS,sizeof(float *))))
  317. X            return(0);
  318. X
  319. X        for(i = 0; i < MS; i++)
  320. X            if(!(*(tmp++)=(float *) calloc((unsigned) LS,
  321. X                    sizeof(float))))
  322. X                return(0);
  323. X
  324. X        return(1);
  325. X    }
  326. X
  327. X    if(M + 1 > MS)
  328. X    {
  329. X        /* Need to reallocate space. */
  330. X        MS = M + 1;
  331. X        if(!(jmax = (int *) realloc((char *)jmax,
  332. X                (unsigned)(MS*sizeof(int)))))
  333. X            return(0);
  334. X
  335. X        if(!(ind1 = (int *) realloc((char *) ind1,
  336. X                (unsigned)(MS*sizeof(int)))))
  337. X            return(0);
  338. X
  339. X        if(!(chk = (int *) realloc((char *) chk,
  340. X                (unsigned)(MS*sizeof(int)))))
  341. X            return(0);
  342. X
  343. X        if(!(delmax = (float *) realloc((char *) delmax,
  344. X                (unsigned)(MS*sizeof(float)))))
  345. X            return(0);
  346. X
  347. X        if(!(e=(float **) realloc((char *) e,
  348. X                (unsigned)(MS*sizeof(float *)))))
  349. X            return(0);
  350. X
  351. X        for(tmp = e + m, i = 0; i < MS; i++)
  352. X            if(!(*(tmp++)=(float *)
  353. X                    malloc((unsigned)(LS*sizeof(float)))))
  354. X                return(0);
  355. X    }
  356. X
  357. X    if(lcol > LS)
  358. X    {
  359. X        LS = lcol;
  360. X        if(!(ind = (int *) realloc((char *) ind,
  361. X                (unsigned)(LS*sizeof(int)))))
  362. X            return(0);
  363. X
  364. X        if(!(imin = (int *) realloc((char *) imin,
  365. X                (unsigned)(LS*sizeof(int)))))
  366. X            return(0);
  367. X
  368. X        if(!(thmin = (float *) realloc((char *) thmin,
  369. X                (unsigned)(LS*sizeof(float)))))
  370. X            return(0);
  371. X
  372. X        for(tmp = e, i = 0; i < MS; tmp++, i++)
  373. X            if(!(*tmp=(float *) realloc((char *) tmp,
  374. X                    (unsigned)(LS*sizeof(float)))))
  375. X                return(0);
  376. X    }
  377. X
  378. X    /* Finally, initialize all arrays to 0. */
  379. X    for(i = 0; i <= m; i++)
  380. X    {
  381. X        jmax[i] = ind1[i] = chk[i] = 0;
  382. X        delmax[i] = 0.0;
  383. X        for(j = 0; j < lcol; j++)
  384. X            e[i][j] = 0.0;
  385. X    }
  386. X
  387. X    for(j = 0; j < lcol; j++)
  388. X    {
  389. X        imin[j] = ind[j] = 0;
  390. X        thmin[j] = 0.0;
  391. X    }
  392. X
  393. X    return(1);
  394. X}
  395. X
  396. X/********************************* MINIT ***************************************
  397. X * This is the main procedure. It is invoked with the various matrices,
  398. X * after space for other arrays has been generated elsewhere (in zx3lp.)
  399. X * It returns
  400. X * 0 if everything went OK and x, w, z as the solutions.
  401. X * 1 if no solution existed.
  402. X * 2 if primal had no feasible solutions.
  403. X * 3 if primal had unbounded solutions.
  404. X ******************************************************************************/
  405. Xint
  406. Xminit(A,b,c,x,w,z)
  407. Xfloat **A, *b, *c, *x, *w, *z;
  408. X{
  409. X    register int i, j;
  410. X    void tab();
  411. X
  412. X    /* Generate initial tableau. */
  413. X    for(j = 0; j < n; j++)
  414. X        e[0][j] = -c[j];
  415. X
  416. X    for(i = 0; i < m; i++)
  417. X    {
  418. X        for(j = 0; j < n; j++)
  419. X            e[i+1][j] = A[i][j];
  420. X    }
  421. X
  422. X    for(j = 0; j < m; j++)
  423. X        e[j+1][lcol-1] = b[j];
  424. X
  425. X    for(i = 0; i < m - p; i++)
  426. X        e[1+i][n+i] = 1.0;
  427. X
  428. X    /* Now begins the actual algorithm. */
  429. X    for(i = 1; i < m+1; i++)
  430. X        chk[i] = -1;    /* Indexing problem; Algol
  431. X                 * version treates 0 as out
  432. X                 * of bounds, in C we prefer -1.
  433. X                 */
  434. X
  435. X    if(!p)    /* There are no equality constraints */
  436. X        goto RCS;
  437. X    else
  438. X    {
  439. X    if(phase1())
  440. X        return(1);
  441. X    }
  442. X        
  443. X    RCS: L = 1; k = 1;
  444. X
  445. X    if(debug)
  446. X        tab();    /* Print the tableau. */
  447. X
  448. X    for(j = 0; j < lcol - 1; j++)
  449. X    {
  450. X        if(e[0][j] < -td)
  451. X            ind[(L++)-1] = j; /* ind[L] keeps track of the
  452. X                       * columns in which e[0][j]
  453. X                       * is negative.
  454. X                       */
  455. X    }
  456. X
  457. X    for(i = 1; i < m + 1; i++)
  458. X    {
  459. X        if(e[i][lcol-1] < -td)
  460. X            ind1[(k++)-1] = i; /* ind1[k] keeps track of the
  461. X                        * rows in which e[i][lcol]
  462. X                        * is negative.
  463. X                        */
  464. X    }
  465. X
  466. X    if(L == 1)
  467. X    {
  468. X        if(k == 1) /* Results */
  469. X            goto RESULTS;
  470. X        else
  471. X        {
  472. X            if(k == 2)
  473. X            {
  474. X                for(j = 0; j < lcol - 1; j++)
  475. X                {
  476. X                    if(e[ind1[0]][j] < 0)
  477. X                        goto R;
  478. X                }
  479. X                /* Primal problem has no feasible solutions. */
  480. X                return(2);
  481. X            }
  482. X            else
  483. X                goto R;
  484. X        }
  485. X    }
  486. X    else
  487. X    {
  488. X        if(L == 2)
  489. X        {
  490. X            if(k == 1)
  491. X            {
  492. X                for(i = 1; i < m + 1; i++)
  493. X                {
  494. X                    if(e[i][ind[0]] > 0)
  495. X                        goto C;
  496. X                }
  497. X
  498. X                /* Primal problem is unbounded. */
  499. X                return(3);
  500. X            }
  501. X            else
  502. X                goto S;
  503. X        }
  504. X
  505. X        if(k == 1)
  506. X            goto C;
  507. X        else
  508. X            goto S;
  509. X    }
  510. X
  511. X    R:
  512. X    prophi();
  513. X    if(rowtrans(imax,jm))
  514. X        return(1);
  515. X    goto RCS;
  516. X
  517. X    C:
  518. X    progamma();
  519. X    if(rowtrans(im,jmin))
  520. X        return(1);
  521. X    goto RCS;
  522. X
  523. X    S:
  524. X    progamma();
  525. X    prophi();
  526. X    if(gmin == LARGENUMBER)
  527. X    {
  528. X        if(rowtrans(imax,jm))
  529. X            return(1);
  530. X        goto RCS;
  531. X    }
  532. X
  533. X    if(phimax == - LARGENUMBER)
  534. X    {
  535. X        if(rowtrans(im,jmin))
  536. X            return(1);
  537. X        goto RCS;
  538. X    }
  539. X
  540. X    if(ABS(phimax) > ABS(gmin))
  541. X    {
  542. X        if(rowtrans(imax,jm))
  543. X            return(1);
  544. X    }
  545. X    else
  546. X    {
  547. X        if(rowtrans(im,jmin))
  548. X            return(1);
  549. X    }
  550. X    goto RCS;
  551. X
  552. X    RESULTS:
  553. X    /* Output results here. */
  554. X    *z = e[0][lcol-1];
  555. X    for(i = 0; i < n; i++)
  556. X        x[i] = 0.0;
  557. X
  558. X    for(j = 0; j < m; j++)
  559. X        x[j] = 0.0;
  560. X
  561. X    for(i = 1; i < m + 1; i++)
  562. X    {
  563. X        if(chk[i] >= n)
  564. X            chk[i] = -1;
  565. X
  566. X        if(chk[i] > -1)
  567. X            x[chk[i]] = e[i][lcol-1];
  568. X    }
  569. X
  570. X    for(j = n; j < lcol - 1; j++)
  571. X        w[j-n] = e[0][j];
  572. X    
  573. X    return(0);
  574. X}
  575. X
  576. X/****************************** rowtrans ***************************************
  577. X * Performs the usual tableau transformations in a linear programming
  578. X * problem, (p,q) being the pivotal element.
  579. X * Returns the following error codes:
  580. X * 0: Everything was OK.
  581. X * 1: No solution.
  582. X ******************************************************************************/
  583. Xint
  584. Xrowtrans(p,q)
  585. Xint p, q;
  586. X{
  587. X    register int i, j;
  588. X    float dummy;
  589. X
  590. X    if(p == -1 || q == -1) /* No solution. */
  591. X        return(1);
  592. X
  593. X    dummy = e[p][q];
  594. X
  595. X    if(debug)
  596. X        (void) printf("--\nPivot element is e[%d][%d] = %f\n",p,q,dummy);
  597. X
  598. X    for(j = 0; j < lcol; j++)
  599. X        e[p][j] /= dummy; 
  600. X
  601. X    for(i = 0; i < m + 1; i++)
  602. X    {
  603. X        if(i != p)
  604. X        {
  605. X            if(e[i][q] != 0.0)
  606. X            {
  607. X                dummy = e[i][q];
  608. X                for(j = 0; j < lcol; j++)
  609. X                    e[i][j] -= e[p][j] * dummy;
  610. X            }
  611. X        }
  612. X    }
  613. X
  614. X    chk[p] = q;
  615. X    return(0);
  616. X} /* rowtrans */
  617. X
  618. X/****************************** progamma ***************************************
  619. X * Performs calculations over columns to determine the pivot element.
  620. X ******************************************************************************/
  621. Xprogamma()
  622. X{
  623. X    float theta, gamma;
  624. X    register int i, L1;
  625. X
  626. X    gmin = LARGENUMBER;    /* gmin is set to a large no. for
  627. X                 * initialization purposes.
  628. X                 */
  629. X    jmin = -1;        /* Array indices in C start from 0 */
  630. X
  631. X    for(L1 = 0; L1 < L - 1; L1++)
  632. X    {
  633. X        imin[ind[L1]] = 0;
  634. X        thmin[ind[L1]] = LARGENUMBER;
  635. X        for(i = 1; i < m + 1; i++)
  636. X        {
  637. X            if(e[i][ind[L1]] > td && e[i][lcol-1] >= -td)
  638. X            {
  639. X                theta = e[i][lcol-1] / e[i][ind[L1]];
  640. X                if(theta < thmin[ind[L1]])
  641. X                {
  642. X                    thmin[ind[L1]] = theta;
  643. X                    imin[ind[L1]] = i;
  644. X                }
  645. X            }
  646. X        }
  647. X
  648. X        if(thmin[ind[L1]] == LARGENUMBER)
  649. X            gamma = VLARGENUMBER;
  650. X        else
  651. X            gamma = thmin[ind[L1]] * e[0][ind[L1]];
  652. X
  653. X        if(gamma < gmin)
  654. X        {
  655. X            gmin = gamma;
  656. X            jmin = ind[L1];
  657. X        }
  658. X    }
  659. X    if(jmin > -1)
  660. X        im = imin[jmin];
  661. X} /* progamma */
  662. X
  663. X/****************************** prophi *****************************************
  664. X * Performs calculations over rows to determine the pivot element.
  665. X ******************************************************************************/
  666. Xprophi()
  667. X{
  668. X    float delta, phi;
  669. X    register int j, k1;
  670. X
  671. X    phimax = - LARGENUMBER;    /* phimax is set to a small no. for
  672. X                 * initialization purposes.
  673. X                 */
  674. X    imax = -1;        /* Array indices in C start from 0 */
  675. X    
  676. X    for(k1 = 0; k1 < k - 1; k1++)
  677. X    {
  678. X        jmax[ind1[k1]] = 0;
  679. X        delmax[ind1[k1]] = - LARGENUMBER;
  680. X        for(j = 0; j < lcol - 1; j++)
  681. X        {
  682. X            if(e[ind1[k1]][j] < -td && e[0][j] >= -td)
  683. X            {
  684. X                delta = e[0][j] / e[ind1[k1]][j];
  685. X                if(delta > delmax[ind1[k1]])
  686. X                {
  687. X                    delmax[ind1[k1]] = delta;
  688. X                    jmax[ind1[k1]] = j;
  689. X                }
  690. X            }
  691. X        }
  692. X
  693. X        if(delmax[ind1[k1]] == - LARGENUMBER)
  694. X            phi = - VLARGENUMBER;
  695. X        else
  696. X            phi = delmax[ind1[k1]] * e[ind1[k1]][lcol-1];
  697. X        
  698. X        if(phi > phimax)
  699. X        {
  700. X            phimax = phi;
  701. X            imax = ind1[k1];
  702. X        }
  703. X    }
  704. X    if(imax > -1)
  705. X        jm = jmax[imax];
  706. X} /* prophi */
  707. X
  708. X/****************************** phase1 *****************************************
  709. X * Applied only to equality constraints if any.
  710. X ******************************************************************************/
  711. Xphase1()
  712. X{
  713. X    float theta, gamma;
  714. X    register int i, j, r;
  715. X    /* Fix suggested by Holmgren, Obradovic, Kolm. */
  716. X    register int im1, jmin1, first;
  717. X    void tab();
  718. X
  719. X    im1 = jmin1 = -1;
  720. X    /* Fix suggested by Messham to allow negative coeffs. in
  721. X     * equality constraints.
  722. X     */
  723. X    for(i = m - p + 1; i < m + 1; i++)
  724. X    {
  725. X        if(e[i][lcol - 1] < 0)
  726. X        {
  727. X            for(j = 0; j < lcol; j++)
  728. X                e[i][j] = -e[i][j];
  729. X        }
  730. X    }
  731. X
  732. X    for(r = 0; r < p; r++)
  733. X    {
  734. X        gmin = LARGENUMBER;    /* gmin is set to a large no. for
  735. X                     * initialization purposes.
  736. X                     */
  737. X        L = 1;
  738. X        jmin = -1;
  739. X        first = 1;
  740. X        for(j = 0; j < n; j++)
  741. X        {
  742. X            thmin[j] = LARGENUMBER;
  743. X            /* Fix suggested by Kolm and Dahlstrand */
  744. X            /* if(e[0,j] < 0) */
  745. X            if(e[0][j] < -td)
  746. X                ind[(L++)-1] = j;
  747. X        }
  748. X
  749. X    L1:    if(L == 1)
  750. X        {
  751. X            for(j = 0; j < n; j++)
  752. X                ind[j] = j;
  753. X
  754. X            L = n + 1;
  755. X        }
  756. X
  757. X        for(k = 0; k < L - 1; k++)
  758. X        {
  759. X            for(i = m - p + 1; i < m + 1; i++)
  760. X                if(chk[i] == -1)
  761. X                {
  762. X                    /* Fix suggested by Kolm
  763. X                     * and Dahlstrand
  764. X                     */
  765. X                    /* if(e[i][ind[k]] > 0.0) */
  766. X                    if(e[i][ind[k]] > td)
  767. X                    {
  768. X                        if((theta = e[i][lcol-1] /
  769. X                          e[i][ind[k]]) < thmin[ind[k]])
  770. X                        {
  771. X                            thmin[ind[k]] = theta;
  772. X                            imin[ind[k]] = i;
  773. X                        }
  774. X                    }
  775. X                }
  776. X
  777. X            /* Fix suggested by Obradovic overrides
  778. X             * fixes suggested by Kolm and Dahstrand
  779. X             * as well as Messham.
  780. X             */
  781. X            if(thmin[ind[k]] < LARGENUMBER)
  782. X            {
  783. X                gamma = thmin[ind[k]] * e[0][ind[k]];
  784. X                if(gamma < gmin)
  785. X                {
  786. X                    gmin = gamma;
  787. X                    jmin = ind[k];
  788. X                }
  789. X            }
  790. X        }
  791. X        if(jmin == -1)
  792. X        {
  793. X            if(first)
  794. X            {
  795. X                first = 0;
  796. X                L = 1;
  797. X                goto L1;
  798. X            }
  799. X            else
  800. X                im = -1;
  801. X        }
  802. X        else
  803. X            im = imin[jmin];
  804. X
  805. X        if(im == im1 && jmin == jmin1)
  806. X        {
  807. X            L = 1;
  808. X            goto L1;
  809. X        }
  810. X
  811. X        if(debug)
  812. X            tab();    /* Print the tableau. */
  813. X
  814. X        if(rowtrans(im,jmin))
  815. X            return(1);
  816. X        
  817. X        im1 = im;
  818. X        jmin1 = jmin;
  819. X    }
  820. X    return(0);
  821. X} /* phase1 */
  822. X
  823. X/****************************** tab *****************************************
  824. X * The following procedure is for debugging. It simply prints the
  825. X * current tableau.
  826. X ******************************************************************************/
  827. Xstatic void
  828. Xtab()
  829. X{
  830. X    register int i, j;
  831. X
  832. X    (void) printf("\n");
  833. X
  834. X    for(i = 0; i < 35; i++)
  835. X        (void) printf("-");
  836. X
  837. X    (void) printf(" TABLEAU ");
  838. X
  839. X    for(i = 0; i < 35; i++)
  840. X        (void) printf("-");
  841. X
  842. X    (void) printf("\n");
  843. X
  844. X    for(i = 0; i < m+1; i++)
  845. X    {
  846. X        for(j = 0; j < lcol; j++)
  847. X            (void) printf("%6.3f ",e[i][j]);
  848. X
  849. X        (void) printf("\n");
  850. X    }
  851. X}
  852. *-*-END-of-minit.c-*-*
  853. echo x - minmain.c
  854. sed 's/^X//' >minmain.c <<'*-*-END-of-minmain.c-*-*'
  855. X/*
  856. X * This is a front end to test the linear programming package
  857. X * minit. This algorithm was written by Das and Salazar in algol.
  858. X *
  859. X * Solve the LP      T
  860. X *        max C X  subj. to
  861. X *              AX <= B
  862. X *               X >= 0
  863. X */
  864. X
  865. X#include <stdio.h>
  866. X#include <string.h>
  867. Xchar *calloc();
  868. X
  869. Xmain(argc,argv)
  870. Xint argc;
  871. Xchar **argv;
  872. X{
  873. X    int  ier,m,n,m1,m2;    /* Dimensions of various matrices */
  874. X    float **A,**A1,*b,*c,*x,*w,z,tmp;
  875. X    char format[64], *usage, *sprintf();
  876. X    int suppress;
  877. X    int int_part, frac_part;
  878. X    float precision;
  879. X    register int i,j,arg;
  880. X    FILE *fopen(), *in;
  881. X
  882. X    n = m1 = m2 = 0;
  883. X    in = stdin;
  884. X    suppress = 1;    /* Suppress printing interactive messages. */
  885. X    precision = 0.0;
  886. X
  887. X    usage = "\
  888. XUsage:\tminit [-pprec] [-i[filename]]\n\
  889. Xor:\tminit [-pprec] < filename\n\
  890. XFor detailed help: minit -help\n";
  891. X
  892. X    for (arg = 1; arg < argc; arg++)
  893. X    {
  894. X        if (argv[arg][0] != '-' ||
  895. X            (argv[arg][1] != 'p' && argv[arg][1] != 'i'))
  896. X        {
  897. X            if (!strcmp(argv[arg],"-help") || !strcmp(argv[arg],"help"))
  898. X            (void) system("man minit");
  899. X            else
  900. X            (void) fputs(usage, stderr);
  901. X            exit(1);
  902. X        }
  903. X
  904. X        if (argv[arg][1] == 'p')
  905. X            (void) sscanf(&(argv[arg][2]),"%f",&precision);
  906. X        else    /* if argv[arg][1] == 'i') */
  907. X        {
  908. X            if (in != stdin)    /* in already assigned. */
  909. X            {
  910. X                (void) fputs("Only one input file permitted.\n",
  911. X                    stderr);
  912. X                exit(1);
  913. X            }
  914. X            if((in = fopen(&(argv[arg][2]),"r")) == NULL)
  915. X            {
  916. X                (void) fprintf(stderr,
  917. X                    "No input file %s; interactive mode\n",
  918. X                    &(argv[arg][2]));
  919. X
  920. X                in = stdin;
  921. X                suppress = 0;
  922. X            }
  923. X        }
  924. X    }
  925. X
  926. X    if (!suppress)
  927. X        (void) fputs("Enter value of n \(No. of unknowns\): ", stderr);
  928. X
  929. X    (void) fscanf(in,"%d",&n);
  930. X
  931. X    if (!suppress)
  932. X        (void) fputs("Enter no. of inequality constraints: ", stderr);
  933. X
  934. X    (void) fscanf(in,"%d",&m1);
  935. X
  936. X    if (!suppress)
  937. X        (void) fputs("Enter no. of equality constraints: ", stderr);
  938. X
  939. X    (void) fscanf(in,"%d",&m2);
  940. X
  941. X    m = m1+m2;
  942. X
  943. X    /* Allocate space for the arrays */
  944. X        A1 = A     = (float **) calloc((unsigned)m,sizeof(float *));
  945. X    for(i=0;i<m;i++)
  946. X        *(A1++) = (float *) calloc((unsigned)n,sizeof(float));
  947. X
  948. X    b     = (float *) calloc((unsigned)m,sizeof(float));
  949. X    c     = (float *) calloc((unsigned)n,sizeof(float));
  950. X    x     = (float *) calloc((unsigned)n,sizeof(float));
  951. X    w     = (float *) calloc((unsigned)m,sizeof(float));
  952. X
  953. X    if (!suppress)
  954. X        (void) fprintf(stderr,"Enter \(%d by %d\) elements of A:\n",
  955. X            m,n);
  956. X    for(i=0;i<m;i++)
  957. X        for(j=0;j<n;j++)
  958. X        {
  959. X            (void) fscanf(in,"%f",&tmp);
  960. X            A[i][j] = tmp;;
  961. X        }
  962. X
  963. X    if (!suppress)
  964. X        (void) fprintf(stderr,
  965. X            "Enter %d elements of b, entered as a linear array:\n ",
  966. X                m);
  967. X    for(i=0;i<m;i++)
  968. X    {
  969. X        (void) fscanf(in,"%f",&tmp);
  970. X        b[i] = tmp;
  971. X    }
  972. X
  973. X    if (!suppress)
  974. X        (void) fprintf(stderr,
  975. X            "Enter %d elements of c, entered as a linear array:\n",
  976. X                n);
  977. X    for(j=0;j<n;j++)
  978. X    {
  979. X        (void) fscanf(in,"%f",&tmp);
  980. X        c[j] = tmp;
  981. X    }
  982. X
  983. X    if (suppress)
  984. X    {
  985. X        (void) fputs("INPUT DATA:\n",stderr);
  986. X        (void) fprintf(stderr,"No. of variables \(n\) = %d\n", n);
  987. X        (void) fprintf(stderr,
  988. X            "No. of inequality constraints \(m1\) = %d\n", m1);
  989. X        (void) fprintf(stderr,
  990. X            "No. of equality constraints \(m2\) = %d\n", m2);
  991. X    }
  992. X
  993. X    (void) fputs("------------------------------",stderr);
  994. X    (void) fputs(" EXECUTING MINIT ",stderr);
  995. X    (void) fputs("------------------------------",stderr);
  996. X    (void) fputc('\n',stderr);
  997. X
  998. X    zx3lp(A,(int *)NULL,b,c,&n,&m1,&m2,&z,x,w,(float *)NULL,
  999. X        (int *)NULL,&ier);
  1000. X    switch(ier)
  1001. X    {
  1002. X        case 133:
  1003. X        break;
  1004. X
  1005. X        case 1:
  1006. X        (void) fprintf(stderr,"No solution.\n");
  1007. X        exit(1);
  1008. X
  1009. X        case 131:
  1010. X        (void) fprintf(stderr,"Primal has unbounded solution.\n");
  1011. X        exit(3);
  1012. X    }
  1013. X
  1014. X    /* Make sure precision is acceptible. */
  1015. X    int_part = precision;
  1016. X    frac_part = precision - int_part;
  1017. X    if (int_part <= frac_part || int_part > 25)
  1018. X        precision = 10.3;    /* The default precision. */
  1019. X
  1020. X    (void) sprintf(format,"\nValue of z: %%%gg\n",precision);
  1021. X    (void) printf(format,z);
  1022. X    (void) puts("Primal solution:");
  1023. X
  1024. X    (void) sprintf(format,"%%%gg ",precision);
  1025. X    for(j=0;j<n;j++)
  1026. X        (void) printf(format,x[j]);
  1027. X
  1028. X    (void) puts("\nDual solution:");
  1029. X    for(i=0;i<m;i++)
  1030. X        (void) printf(format,w[i]);
  1031. X
  1032. X    (void) putchar('\n');
  1033. X    exit(0);
  1034. X}
  1035. *-*-END-of-minmain.c-*-*
  1036. echo x - example
  1037. sed 's/^X//' >example <<'*-*-END-of-example-*-*'
  1038. X4 2 1
  1039. X
  1040. X1 2  1 1
  1041. X3 1 -1 0
  1042. X0 1  1 1
  1043. X
  1044. X5 8 1
  1045. X6 1 -1 -1
  1046. *-*-END-of-example-*-*
  1047. exit
  1048.  
  1049.