home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 1.ddi / CHAP6.ARC / DIRFACT.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1987-12-30  |  11.3 KB  |  362 lines
  1. program Direct_Factorization_Prog;
  2.  
  3. {--------------------------------------------------------------------------}
  4. {-                                                                        -}
  5. {-     Turbo Pascal Numerical Methods Toolbox                             -}
  6. {-     Copyright (c) 1986, 87 by Borland International, Inc.              -}
  7. {-                                                                        -}
  8. {-       Purpose : This program demonstrates how to solve a system of     -}
  9. {-                 linear equations with direct factorization (LU         -}
  10. {-                 decomposition).                                        -}
  11. {-                                                                        -}
  12. {-       Unit    : Matrix    procedure LU_Decompose                       -}
  13. {-                           procedure LU_Solve                           -}
  14. {-                                                                        -}
  15. {--------------------------------------------------------------------------}
  16.  
  17. {$R+}                  { Enable range checking }
  18. {$I-}                  { Disable I/O checking }
  19. {$M 32768, 0, 655360}
  20.  
  21. uses
  22.   Matrix, Dos, Crt, Common;
  23.  
  24. type
  25.   Ptr = ^StackItem;
  26.  
  27.   QueueItem = record
  28.                 Front, Back : Ptr;
  29.               end;
  30.  
  31.   StackItem = record
  32.                 Info : TNvector;
  33.                 Next : Ptr;
  34.               end;
  35.  
  36.  
  37. var
  38.   Dimen : integer;          { Size of the square matrix }
  39.   Coefficients : TNmatrix;  { The matrix }
  40.   Queue : QueueItem;        { Pointer to Queue of constant vectors }
  41.   Decomp : TNmatrix;        { LU decompostion of Coefficients matrix }
  42.   Permute : TNmatrix;       { Permutation matrix from partial pivoting }
  43.   Solution : TNvector;      { Solution to the set of equations }
  44.   Error : byte;             { Flags if something went wrong }
  45.  
  46. procedure InitializeStack(var Queue : QueueItem);
  47.  
  48. begin
  49.   Queue.Front := nil;
  50.   Queue.Back := nil;
  51. end; { procedure InitializeStack }
  52.  
  53. procedure InsertQueue(Data  : TNvector;
  54.                   var Queue : QueueItem);
  55.  
  56. {--------------------------------}
  57. {- Input: Data, Queue           -}
  58. {- Output: Queue                -}
  59. {-                              -}
  60. {- Insert Data onto the Queue   -}
  61. {--------------------------------}
  62.  
  63. var
  64.   NewNode : Ptr;
  65.  
  66. begin
  67.   New(NewNode);
  68.   NewNode^.Info := Data;
  69.   NewNode^.Next := nil;
  70.   if Queue.Back = nil then
  71.     Queue.Front := NewNode
  72.   else
  73.     Queue.Back^.Next := NewNode;
  74.   Queue.Back := NewNode;
  75. end; { procedure InsertQueue }
  76.  
  77. procedure RemoveQueue(var Data  : TNvector;
  78.                       var Queue : QueueItem);
  79.  
  80. {---------------------------------}
  81. {- Input: Queue                  -}
  82. {- Output: Data, Queue           -}
  83. {-                               -}
  84. {- Remove Data from the Queue    -}
  85. {---------------------------------}
  86.  
  87. var
  88.   OldNode : Ptr;
  89.  
  90. begin
  91.   OldNode := Queue.Front;
  92.   Queue.Front := OldNode^.Next;
  93.   Data := OldNode^.Info;
  94.   Dispose(OldNode);
  95. end; { procedure RemoveQueue }
  96.  
  97.  
  98. procedure Initial(var Dimen        : integer;
  99.                   var Coefficients : TNmatrix;
  100.                   var Queue        : QueueItem);
  101.  
  102. {----------------------------------------------------------}
  103. {- Output: Dimen, Coefficients, Queue                     -}
  104. {-                                                        -}
  105. {- This procedure intializes the above variables to zero. -}
  106. {----------------------------------------------------------}
  107.  
  108. begin
  109.   Dimen := 0;
  110.   FillChar(Coefficients, SizeOf(Coefficients), 0);
  111.   InitializeStack(Queue);
  112. end; { procedure Initial }
  113.  
  114. procedure GetData(var Dimen        : integer;
  115.                   var Coefficients : TNmatrix;
  116.                   var Queue        : QueueItem);
  117.  
  118. {-----------------------------------------------------------}
  119. {- Output: Dimen, Coefficients, Queue                      -}
  120. {-                                                         -}
  121. {- This procedure sets the value of Dimen and Coefficients -}
  122. {- keyboard input or file input. This procedure also sets  -}
  123. {- the Queue pointer, Queue, so that it points to a Queue  -}
  124. {- of constant vectors.                                    -}
  125. {-----------------------------------------------------------}
  126.  
  127. var
  128.   Ch : char;
  129.  
  130. procedure GetDataFromKeyboard(var Dimen        : integer;
  131.                               var Coefficients : TNmatrix;
  132.                               var Queue        : QueueItem);
  133.  
  134. {----------------------------------------------------}
  135. {- Output: Dimen, Coefficients, Queue               -}
  136. {-                                                  -}
  137. {- This procedure sets the value of Dimen,          -}
  138. {- Coefficients, and the Queue from keyboard input. -}
  139. {----------------------------------------------------}
  140.  
  141. var
  142.   Row, Column : integer;
  143.   Constants : TNvector;
  144.   Done : boolean;
  145.   Ch : char;
  146.  
  147. begin
  148.   Writeln;
  149.   repeat
  150.     Write('Dimension of the coefficient matrix (1-', TNArraySize,')? ');
  151.     Readln(Dimen);
  152.     IOCheck;
  153.   until not IOerr and (Dimen >= 1) and (Dimen <= TNArraySize);
  154.   Writeln;
  155.   for Row := 1 to Dimen do
  156.     for Column := 1 to Dimen do
  157.       repeat
  158.         Write('Matrix[', Row, ', ', Column, ']: ');
  159.         Readln(Coefficients[Row, Column]);
  160.         IOCheck;
  161.       until not IOerr;
  162.   Done := false;
  163.   while not Done do
  164.   begin
  165.     Writeln;
  166.     Writeln('Elements of the constant vector: ');
  167.     for Row := 1 to Dimen do
  168.       repeat
  169.         Write('Vector[', Row, ']: ');
  170.         Readln(Constants[Row]);
  171.         IOCheck;
  172.       until not IOerr;
  173.     InsertQueue(Constants, Queue);
  174.     Writeln;
  175.     Write('Another constant vector (Y/N)? ');
  176.     repeat
  177.       Ch := ReadKey;
  178.       Ch := UpCase(Ch);
  179.     until Ch in ['Y', 'N'];
  180.     Writeln(Ch);
  181.     if Ch = 'N' then
  182.       Done := true;
  183.   end; { while }
  184. end; { procedure GetDataFromKeyboard }
  185.  
  186. procedure GetDataFromFile(var Dimen        : integer;
  187.                           var Coefficients : TNmatrix;
  188.                           var Queue        : QueueItem);
  189.  
  190. {--------------------------------------------------}
  191. {- Output: Dimen, Coefficients, Queue             -}
  192. {-                                                -}
  193. {- This procedure sets the value of Dimen,        -}
  194. {- Coefficients and the Queue from file input.    -}
  195. {--------------------------------------------------}
  196.  
  197. var
  198.   FileName    : string[255];
  199.   Constants   : TNvector;
  200.   InFile      : text;
  201.   Row, Column : integer;
  202.  
  203. begin
  204.   Writeln;
  205.   repeat
  206.     Writeln;
  207.     repeat
  208.       Write('File name? ');
  209.       Readln(FileName);
  210.       Assign(InFile, FileName);
  211.       Reset(InFile);
  212.       IOCheck;
  213.     until not IOerr;
  214.     Read(InFile, Dimen);
  215.     IOCheck;
  216.     Row := 0;
  217.     while (not IOerr) and (Row < Dimen) do
  218.     begin
  219.       Row := Succ(Row);
  220.       Column := 0;
  221.       while (not IOerr) and (Column < Dimen) do
  222.       begin
  223.         Column := Succ(Column);
  224.         Read(InFile, Coefficients[Row, Column]);
  225.         IOCheck;
  226.       end;
  227.     end;
  228.     while (not EOF(InFile)) and (not IOerr) do
  229.     begin
  230.       Row := 0;
  231.       while (not IOerr) and (Row < Dimen) do
  232.       begin
  233.         Row := Succ(Row);
  234.         Readln(InFile, Constants[Row]);
  235.         IOCheck;
  236.       end;
  237.       if not IOerr then
  238.         InsertQueue(Constants, Queue);
  239.     end;
  240.   until not IOerr;
  241.   Close(InFile);
  242. end; { procedure GetDataFromFile }
  243.  
  244. begin { procedure GetData }
  245.   case InputChannel('Input Data From') of
  246.     'K' : GetDataFromKeyboard(Dimen, Coefficients, Queue);
  247.     'F' : GetDataFromFile(Dimen, Coefficients, Queue);
  248.   end;
  249.   GetOutputFile(OutFile);
  250. end; { procedure GetData }
  251.  
  252. procedure PrintCoefMatrix(Dimen        : integer;
  253.                       var Coefficients : TNmatrix;
  254.                           Error        : byte);
  255.  
  256. {------------------------------------------------------------}
  257. {- Input: Dimen, Coefficients, Error                        -}
  258. {-                                                          -}
  259. {- This procedure prints the coefficients matrix to         -}
  260. {- the device OutFile.                                      -}
  261. {------------------------------------------------------------}
  262.  
  263. var
  264.   Column, Row : integer;
  265.  
  266. begin
  267.   Writeln(OutFile);
  268.   Writeln(OutFile, 'The coefficients: ');
  269.   for Row := 1 to Dimen do
  270.   begin
  271.     for Column := 1 to Dimen do
  272.       Write(OutFile, Coefficients[Row, Column]:13:9);
  273.     Writeln(OutFile);
  274.   end;
  275.   Writeln(OutFile);
  276.   if Error >= 1 then
  277.     DisplayError;
  278.  
  279.   case Error of
  280.     1 : Writeln(OutFile, 'The dimension of the matrix must be greater than 0.');
  281.  
  282.     2 : Writeln(OutFile, 'There is no solution to this set of equations.');
  283.  
  284.   end; { case }
  285. end; { PrintCoefMatrix }
  286.  
  287. procedure Results(Dimen     : integer;
  288.               var Constants : TNvector;
  289.               var Solution  : TNvector;
  290.                   Error     : byte);
  291.  
  292. {------------------------------------------------------------}
  293. {- Input: Dimen, Constants, Solution, Error                 -}
  294. {-                                                          -}
  295. {- This procedure outputs the solution to the system of     -}
  296. {- equations with the constant terms contained in the       -}
  297. {- vector Constants.                                        -}
  298. {------------------------------------------------------------}
  299.  
  300. var
  301.   Row : integer;
  302.  
  303. begin
  304.   Writeln(OutFile);
  305.   Writeln(OutFile);
  306.   Writeln(OutFile, 'The constants:');
  307.   for Row := 1 to Dimen do
  308.     Writeln(OutFile, Constants[Row]);
  309.   Writeln(OutFile);
  310.   if Error >= 1 then
  311.     DisplayError;
  312.  
  313.   case Error of
  314.     0 : begin
  315.           Writeln(OutFile, 'The solution:');
  316.           for Row := 1 to Dimen do
  317.             Writeln(OutFile, Solution[Row]);
  318.         end;
  319.  
  320.     1 : Writeln(OutFile, 'The dimension of the matrix must be greater than 0.');
  321.   end; { case }
  322. end; { procedure Results }
  323.  
  324. procedure FindSolutions(Dimen   : integer;
  325.                     var Decomp  : TNmatrix;
  326.                     var Permute : TNmatrix;
  327.                         Queue   : QueueItem);
  328.  
  329. {--------------------------------------------------------------}
  330. {- Input: Dimen, Decomp, Permute, Queue                       -}
  331. {-                                                            -}
  332. {- This procedure solves each of the system of equations      -}
  333. {- represented by the Constants vectors stored on the Queue.  -}
  334. {- It pulls each Constants vector off the Queue and sends it  -}
  335. {- to the procedure LU_Solve where backwards and forwards     -}
  336. {- substitution is used to solve the matrix equation.  Each   -}
  337. {- of the Solutions returned from this procedure are output   -}
  338. {- in the procedure Results.                                  -}
  339. {--------------------------------------------------------------}
  340.  
  341. var
  342.   Constants : TNvector;
  343.  
  344. begin
  345.   while Queue.Front <> nil do
  346.   begin
  347.     RemoveQueue(Constants, Queue);
  348.     LU_Solve(Dimen, Decomp, Constants, Permute, Solution, Error);
  349.     Results(Dimen, Constants, Solution, Error);
  350.   end;
  351. end; { procedure FindSolutions }
  352.  
  353. begin { program Direct_Factorization }
  354.   ClrScr;
  355.   Initial(Dimen, Coefficients, Queue);
  356.   GetData(Dimen, Coefficients, Queue);
  357.   LU_Decompose(Dimen, Coefficients, Decomp, Permute, Error);
  358.   PrintCoefMatrix(Dimen, Coefficients, Error);
  359.   if Error = 0 then
  360.     FindSolutions(Dimen, Decomp, Permute, Queue);
  361.   Close(OutFile);
  362. end. { program Direct_Factorization }