home *** CD-ROM | disk | FTP | other *** search
/ CP/M / CPM_CDROM.iso / simtel / sigm / vols000 / vol019 / linearp.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1984-04-29  |  10.4 KB  |  461 lines
  1. Program LINEAR;
  2. (*  PROGRAM TITLE:    Linear Programming
  3. **
  4. **  WRITTEN BY:        W.M. Yarnall
  5. **            19 Angus Lane
  6. **            Warren, N.J. 07060
  7. **  DATE WRITTEN:    March 1980
  8. **
  9. **  WRITTEN FOR:    S100 MICROSYSTEMS
  10. **            MAR 1980
  11. **
  12. **  SUMMARY:        Minimize a cost function to constraints.
  13. **            Maximize negative of 'profit' function.
  14. **            This program uses the Revised Simplex Algorithm.
  15. **
  16. **  MODIFICATION RECORD:
  17. **    25 MAY 1980    -MODIFIED FOR PASCAL/Z BY RAYMOND E. PENLEY
  18. **
  19. **        ---NOTE---
  20. **
  21. ** The first logical record in Pascal/Z is No.1, NOT record
  22. ** No. 0 as in Pascal/M or UCSD Pascal. This can be rectified
  23. ** very eaisly by adding a "BIAS" to each record number.
  24. **    Pascal/Z : bias = 1    |    Pascal/M : bias = 0
  25. **
  26. *)
  27. LABEL    99;      { File not found exit }
  28.  
  29. CONST
  30.   maxrow = 32;
  31.   maxcol = 64;
  32.   bias   =  1;       (* Bias added to each record *)
  33.   FID_LENGTH = 14; (* MAXIMUM LENGTH ALLOWED FOR A FILE NAME *)
  34.  
  35. TYPE
  36.   FID  = STRING FID_LENGTH;
  37.   ROW = array [1..maxrow] of real;
  38.   COL = array [1..maxcol] of real;
  39.   Frec = record
  40.        CASE TAG : integer of
  41.         0: (name : STRING 6; num1, num2 : integer);
  42.         1: (header : STRING 64);
  43.         2: (Rname : STRING 6; Rindex : integer; RHS : real);
  44.         4: (Cname : STRING 6; Cindex : integer; OBJ : real);
  45.         6: (R, S : integer; T : real);
  46.        99: () {End_Of_File}
  47.      end;
  48.  
  49.   STRING80 = STRING 80;
  50.  
  51. VAR
  52.   ABAR         : array [1..maxrow, 1..maxcol] of real;
  53.   Colname     : array [1..maxcol] of STRING 6;
  54.   fa        : FILE of Frec;    (*---File descriptor <FCB>---*)
  55.   File_ID    : FID;        (*---File Identifier <FID>---*)
  56.   F        : Frec;
  57.   heading    : STRING 64;
  58.   hdrflag    : boolean;
  59.   list        : array [1..maxrow] of integer;
  60.   M, N,
  61.   MP, M1    : integer;
  62.   PNAME        : STRING 6;
  63.   Result    : integer;
  64.   Rowname     : array [1..maxrow] of STRING 6;
  65.   U         : array [1..maxrow, 1..maxrow] of real;
  66.   X,XIK        : ROW;
  67.  
  68. PROCEDURE GETID( MESSAGE : STRING80; VAR ID: FID );
  69. (**
  70.     FID_LENGTH = 14;
  71.     STRING80 = STRING 80;
  72.     FID      = STRING FID_LENGTH;
  73. **)
  74. CONST    SPACE = ' ';
  75. TYPE
  76. (*----Required for PASCAL/Z supplied functions----*)
  77.   STR0 = STRING 0;
  78.   STR255 = STRING 255;
  79.         (*----required by PASCAL/Z----*)
  80.     FUNCTION  LENGTH(X: STR255): INTEGER; EXTERNAL;
  81.     PROCEDURE SETLENGTH(VAR X: STR0; Y: INTEGER); EXTERNAL;
  82.  
  83. begin{GetID}
  84.   SETLENGTH(ID,0);
  85.   writeln;
  86.   write(message);
  87.   READLN(ID);
  88.   While Length(ID)<FID_LENGTH DO APPEND(ID,SPACE)
  89. End{---of GETID---};
  90.  
  91. Procedure PRINTH;
  92. begin
  93.   writeln;
  94.   writeln(' PROG. NAME = ', Pname);
  95.   writeln(' No. Rows   = ', M:6);
  96.   writeln(' No. Cols   = ', N:6);
  97.   writeln
  98. end(*---of PRINTH---*);
  99.  
  100. Procedure PRINTC( B : row ; C : col );
  101. VAR    I : integer;
  102. begin
  103.   writeln;
  104.   writeln('  ---Initial Data---');
  105.   writeln;
  106.   writeln(' Objective Vector');
  107.   writeln;
  108.   For I:=1 to N do
  109.     writeln( Colname[I]:8, C[I]:14:8 );
  110.   writeln;
  111.   writeln(' RHS Vector');
  112.   WRITELN;
  113.   For I:=1 to M do
  114.     writeln( Rowname[I]:8, B[I]:14:8 );
  115.   writeln
  116. end(*---of PRINTC---*);
  117.  
  118. Procedure PRINTD;
  119. VAR    I, J : integer;
  120. begin
  121.   writeln;
  122.   writeln(' ABAR array');
  123.   writeln;
  124.   For J:=1 to N do
  125.     For I:=1 to M do
  126.       Writeln(I:6, J:6,'  ', Rowname[I]:8, Colname[J]:8, ABAR[I,J]:14:8);
  127.   writeln;
  128.   writeln(' ABAR(M+1), ABAR(M+2)');
  129.   writeln;
  130.   For I:=1 to N do
  131.     writeln( Colname[I]:8, ABAR[M1,I]:14:8, ABAR[MP,I]:14:8 );
  132.   writeln
  133. end(*---of PRINTD---*);
  134.  
  135. Procedure PRINTX;
  136. VAR    I : integer;
  137.     S : STRING 6;
  138. begin
  139.   writeln;
  140.   writeln(' List and X Arrays');
  141.   writeln;
  142.   For I:=1 to MP do
  143.     begin
  144.     S := '      ';
  145.     If (list[I]<=N) then S := Colname[ List[I] ];
  146.     If I>M then S := Rowname[I] ;
  147.     writeln( I:8,'  ', S:8, list[I]:7, X[I]:18:8 )
  148.     end(*FOR*);
  149.   writeln
  150. end(*---of PRINX---*);
  151.  
  152. Procedure EXITER(exitcode, X : integer);
  153. begin
  154.   CASE exitcode of
  155.    1:    begin
  156.     Result := 1; (* Normal exit *)
  157.     Writeln(' End of Phase 1 for ', Pname, ' after', X:3,
  158.         ' Iterations');
  159.     PRINTX
  160.     end;
  161.    2:    begin
  162.     Result := 2; (* Error exit *)
  163.     Writeln(' Error in Iteration', X:3);
  164.     PRINTX
  165.     end;
  166.    3:    begin
  167.     Result := 3; (* No feasible solution *)
  168.     Writeln(' No feasible solution after', X:3, ' Iterations');
  169.     PRINTX
  170.     end;
  171.    4:    begin
  172.     Result := 1; (* Normal exit *)
  173.     Writeln(' End of Phase 2 for ', Pname, ' after', X:3,
  174.         ' Iterations');
  175.     PRINTX
  176.     end;
  177.    5:    begin
  178.     Result := 2; (* Unbounded solution *)
  179.     Writeln(' Unbounded solution for ', Pname);
  180.     PRINTX
  181.     end
  182.    end(* CASE exitcode of *)
  183. end(*---of EXITER---*);
  184.  
  185. Procedure INITIAL;
  186. VAR    Rcd,        {Record counter}
  187.     I,J : integer;
  188.     sum : real;
  189.     XEOF,        {End of File flag for a NON text File}
  190.     firstin : boolean;
  191.     B : ROW;
  192.     C : COL;
  193. begin
  194.   For I:=1 to maxrow do
  195.     For J:=1 to maxcol do ABAR[I,J] := 0.0 ;
  196.   firstin := false;
  197.   Rcd := 0;{start at the beginning}
  198.   READ(fa:Rcd+bias, F);
  199.   XEOF := (F.tag=99);
  200.   If F.tag=0 then
  201.     begin
  202.       firstin    := true;
  203.       Pname    := F.name;
  204.       M        := F.num1; {No. Rows}
  205.       N        := F.num2; {No. Columns}
  206.       MP    := M + 2;
  207.       M1    := M + 1;
  208.       PRINTH
  209.     end(* IF *)
  210.   Else
  211.     begin
  212.       writeln;
  213.       writeln(' Bad file format');
  214.       writeln;
  215.       Result := 2
  216.     end(* ELSE *);
  217.   While (firstin) AND (NOT XEOF) do
  218.     begin
  219.     With F do
  220.       CASE TAG of
  221.     1: begin(* heading *)
  222.        heading := header;
  223.        hdrflag := true
  224.        end;
  225.     2: begin(* row_name & RHS *)
  226.        Rowname[Rindex] := Rname;
  227.        B[Rindex] := RHS
  228.        end;
  229.     4: begin(* col_name & OBJ *)
  230.        Colname[Cindex] := Cname;
  231.        C[Cindex] := OBJ
  232.        end;
  233.     6: ABAR[R,S] := T;
  234.        99: (* NULL *)
  235.       End{With/Case};
  236.     Rcd := Rcd + 1;
  237.     READ(fa:Rcd+bias, F);
  238.     XEOF := (F.tag=99);
  239.     end(* While *);
  240.   If firstin then
  241.     begin
  242.     PRINTC(B,C);
  243.     For J:=1 to N do ABAR[M1,J] := C[J];
  244.     For I:=1 to M do
  245.       If B[I]<0.0 then
  246.     begin
  247.     B[I] := -B[I];
  248.     For J:=1 to N do ABAR[I,J] := -ABAR[I,J]
  249.     end;
  250.     For J:=1 to N do
  251.       begin
  252.     SUM := 0.0;
  253.     For I:=1 to M do SUM := SUM - ABAR[I,J];
  254.     ABAR[MP,J] := SUM
  255.       end;
  256.     B[M1] := 0.0;
  257.     SUM := 0.0;
  258.     For I:=1 to M do SUM := SUM - B[I];
  259.     B[MP] := SUM;
  260.     For I:=1 to MP do
  261.       begin
  262.     X[I] := B[I];
  263.     list[I] := N +I;
  264.     For J:=1 to MP do U[I,J] := 0.0
  265.       end;
  266.     For I:=1 to MP do U[I,I] := 1.0;
  267.     PRINTD;
  268.     Rowname[M1] := 'M+1   ';
  269.     Rowname[MP] := 'M+2   ';
  270.     PRINTX
  271.     end(* If firstin *);
  272.   Writeln
  273. end(*---of INITIAL---*);
  274.  
  275. Procedure PHASE1;
  276. LABEL    304; (* Exit point *)
  277. CONST    TOL = 1.0E-5;
  278. VAR    iter, I, J, L, ksave : integer;
  279.     sum, temp, theta, Z  : real;
  280.     XL, XLK             : real;
  281.     DEL, V, W         : ROW;
  282.     test             : boolean;
  283. begin
  284.   writeln(' Start Phase 1');
  285.   writeln;
  286.   iter := 0;
  287.   While true do
  288.     begin
  289.     If ABS(X[MP])<tol then {normal exit}
  290.     begin EXITER(1,iter); goto 304 end;
  291.     If X[MP]>tol then {error exit}
  292.     begin EXITER(2,iter); goto 304 end;
  293.     iter := iter +1;
  294.     For J:=1 to N do
  295.       begin
  296.     SUM := 0.0;
  297.     For I:=1 to MP do
  298.       SUM := SUM + U[MP,I] * ABAR[I,J];
  299.     DEL[J] := SUM
  300.       end;
  301.     test := true;
  302.     For J:=1 to N do
  303.       If DEL[J]<0.0 then test := false;
  304.     If test then {no feasible solution exit}
  305.       begin EXITER(3,iter); goto 304 end;
  306.     temp := 1.0E+36;
  307.     ksave := 0;
  308.     For J:=1 to N do
  309.       If DEL[J]<temp then
  310.     begin temp := DEL[J]; ksave := J end;
  311.     For I:=1 to MP do
  312.       begin
  313.     SUM := 0.0;
  314.     For J:=1 to MP do
  315.       SUM := SUM + U[I,J] * ABAR[J,ksave];
  316.     XIK[I] := SUM
  317.       end;
  318.     theta := 1.0E+36;
  319.     L := 0;
  320.     For I:=1 to M do
  321.       If XIK[I]>0.0 then
  322.     begin
  323.     Z := X[I] / XIK[I];
  324.     If (Z=theta) AND (list[I]>N) then
  325.       L := I
  326.     Else
  327.       If Z<theta then
  328.         begin theta := Z; L := I end
  329.     end;
  330.     If L=0 then
  331.       begin EXITER(2,iter); goto 304 end;
  332.     list[L] := ksave;
  333.     For I:=1 to MP do
  334.       begin
  335.     V[I] := XIK[I] / XIK[L];
  336.     W[I] := U[L,I]
  337.       end;
  338.     XL := X[L];
  339.     XLK := XIK[L];
  340.     For I:=1 to MP do
  341.       begin
  342.     Z := theta;
  343.     If (list[I]<>ksave) then Z := X[I] - XL * V[I];
  344.     X[I] := Z;
  345.     For J:=1 to M do
  346.       begin
  347.         Z := W[J] / XLK;
  348.         If I<>L then Z := U[I,J] - W[J] * V[I];
  349.         U[I,J] := Z
  350.       end
  351.       end;
  352.     writeln(' Iteration', iter:3, ' of ', Pname);
  353.     PRINTX
  354.     end(* While true *);
  355. 304: (* Exit point *)
  356. end(*---of PHASE1---*);
  357.  
  358. Procedure PHASE2;
  359. LABEL    403; (* Exit point *)
  360. CONST    TOL = -1.0E-5;
  361. VAR    I, J, L, iter, ksave : integer;
  362.     SUM, temp, theta, Z  : real;
  363.     XL, XLK             : real;
  364.     DEL, V, W         : ROW;
  365.     test             : boolean;
  366. begin
  367.   iter := 0;
  368.   writeln(' Start Phase 2');
  369.   writeln;
  370.   While true do
  371.     begin
  372.     For J:=1 to N do
  373.       begin
  374.     SUM := 0.0;
  375.     For I:=1 to MP do
  376.       SUM := SUM + U[M1,I] * ABAR[I,J];
  377.     DEL[J] := SUM
  378.       end;
  379.     test := true;
  380.     For J:=1 to N do
  381.       If DEL[J]<tol then test := false;
  382.     If test then
  383.       begin EXITER(4,iter); goto 403 end;
  384.     iter := iter +1;
  385.     temp := 1.0E+36;
  386.     ksave := 0;
  387.     For J:=1 to N do
  388.       If DEL[J]<temp then
  389.     begin temp := DEL[J]; ksave := J end;
  390.     For I:=1 to MP do
  391.       begin
  392.     SUM := 0.0;
  393.     For J:=1 to MP do
  394.       SUM := SUM + U[I,J] * ABAR[J,ksave];
  395.     XIK[I] := SUM
  396.       end;
  397.     test := true;
  398.     For I:=1 to MP do
  399.       If XIK[I]>0.0 then test := false;
  400.     If test then
  401.       begin EXITER(5,iter); goto 403 end;
  402.     theta := 1.0E+36;
  403.     L := 0;
  404.     For I:=1 to M do
  405.       If XIK[I]>0.0 then
  406.     begin
  407.       Z := X[I] / XIK[I];
  408.       If Z<theta then
  409.         begin theta := Z; L := I end
  410.     end;
  411.     List[L] := ksave;
  412.     For I:=1 to MP do
  413.       begin
  414.     V[I] := XIK[I] / XIK[L];
  415.     W[I] := U[L,I];
  416.       end;
  417.     XL := X[L];
  418.     XLK := XIK[L];
  419.     For I:=1 to MP do
  420.       begin
  421.     Z := theta;
  422.     If (list[I]<>ksave) then Z := X[I] - XL * V[I];
  423.     X[I] := Z;
  424.     For J:=1 to M do
  425.       begin
  426.         Z := W[J] / XLK;
  427.         If I<>L then Z := U[I,J] - W[J] * V[I];
  428.         U[I,J] := Z
  429.       end
  430.       end;
  431.     writeln(' Iteration', iter:3, ' of ', Pname);
  432.     PRINTX;
  433.     end(* While true *);
  434. 403: (* Exit point *)
  435. end(*---of PHASE2---*);
  436.  
  437. Procedure CLEAR;
  438. (* simple screen clear routine *)
  439. VAR    ix : 1..25;
  440. begin
  441.   for ix:=1 to 25 do writeln
  442. end;
  443.  
  444. BEGIN (***   MAIN PROGRAM   ***)
  445.   CLEAR;
  446.   GETID(' Enter data File Name ---> ', File_ID);
  447.   RESET(File_ID, fa);    (*---RESET( <FID> , <FCB> )---*)
  448.   If EOF(fa) then
  449.     begin
  450.     Writeln(CHR(7),'File ',File_ID,'not found');
  451.     {exit}goto 99
  452.     end;
  453.   Writeln;
  454.   INITIAL;
  455.   If Result<>2 then PHASE1;
  456.   If Result=1 then PHASE2;
  457.   If hdrflag then Writeln(' ', heading);
  458. 99: {File not found exit};
  459.   Writeln;Writeln;Writeln;Writeln;Writeln
  460. end(*---of Linear---*).
  461.