home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 1.ddi / CHAP7.ARC / INVPOWER.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1987-12-30  |  12.5 KB  |  379 lines
  1. program InversePower_Prog;
  2.  
  3. {----------------------------------------------------------------------------}
  4. {-                                                                          -}
  5. {-     Turbo Pascal Numerical Methods Toolbox                               -}
  6. {-     Copyright (c) 1986, 87 by Borland International, Inc.                -}
  7. {-                                                                          -}
  8. {-       Purpose : to demonstrate the procedure InversePower for            -}
  9. {-                 approximating an eigenvalue and eigenvector of a         -}
  10. {-                 matrix.                                                  -}
  11. {-                                                                          -}
  12. {-       Unit    : EigenVal    procedure InversePower                       -}
  13. {-                                                                          -}
  14. {----------------------------------------------------------------------------}
  15.  
  16. {$R+}                  { Enable range checking }
  17. {$I-}                  { Disable I/O checking }
  18. {$M 45056, 0, 655360}  { Set MinStack:MinHeap:MaxHeap }
  19.  
  20. uses
  21.   EigenVal, Dos, Crt, Common;
  22.  
  23. var
  24.   Dimen : integer;          { Size of the square matrix }
  25.   Mat : TNmatrix;           { The matrix }
  26.   MaxIter : integer;        { Maximum # iterations allowed }
  27.   Tolerance : Float;        { Tolerance }
  28.   GuessVector : TNvector;   { Initial guess of an eigenvector }
  29.   ClosestVal : Float;       { Initial approximation of the eigenvalues. }
  30.                             { The algorithm will converge to the        }
  31.                             { eigenvalue closest to ClosestVal          }
  32.   Eigenvector : TNvector;   { Eigenvector of the matrix }
  33.   Eigenvalue : Float;       { Associated eigenvector }
  34.   Iter : integer;           { Number of iterations }
  35.   Error : byte;             { Flags if something went wrong }
  36.  
  37. procedure GetData(var Dimen       : integer;
  38.                   var Mat         : TNmatrix;
  39.                   var ClosestVal  : Float;
  40.                   var Tolerance   : Float;
  41.                   var MaxIter     : integer;
  42.                   var GuessVector : TNvector);
  43.  
  44. {-------------------------------------------------}
  45. {- Output: Dimen, Mat, ClosestVal, Tolerance,    -}
  46. {-         MaxIter, GuessVector                  -}
  47. {-                                               -}
  48. {- This procedure assigns values to the above    -}
  49. {- variables from either keyboad or file input.  -}
  50. {-------------------------------------------------}
  51.  
  52. var
  53.   Ch : char;
  54.  
  55. procedure GetDataFromKeyboard(var Dimen       : integer;
  56.                               var Mat         : TNmatrix;
  57.                               var GuessVector : TNvector);
  58.  
  59. {----------------------------------------------}
  60. {- Output: Dimen, Mat, GuessVector            -}
  61. {-                                            -}
  62. {- This procedure assigns values to the above -}
  63. {- variables from keyboard input.             -}
  64. {----------------------------------------------}
  65.  
  66. var
  67.   Row, Column : integer;
  68.  
  69. begin
  70.   Writeln;
  71.   repeat
  72.     Write('Dimension of the matrix (1-', TNArraySize, ')? ');
  73.     Readln(Dimen);
  74.     IOCheck;
  75.   until (not IOerr) and (Dimen >= 1) and (Dimen <= TNArraySize);
  76.   Writeln;
  77.   for Row := 1 to Dimen do
  78.     for Column := 1 to Dimen do
  79.     repeat
  80.       Write('Matrix[', Row, ', ', Column, ']: ');
  81.       Readln(Mat[Row, Column]);
  82.       IOCheck;
  83.     until not IOerr;
  84.   Writeln;
  85.   Writeln('Now input an initial guess for the eigenvector:');
  86.   for Row := 1 to Dimen do
  87.   repeat
  88.     Write('Vector[', Row, ']: ');
  89.     Readln(GuessVector[Row]);
  90.     IOCheck;
  91.   until not IOerr;
  92. end; { procedure GetDataFromKeyboard }
  93.  
  94. procedure GetDataFromFile(var Dimen       : integer;
  95.                           var Mat         : TNmatrix;
  96.                           var GuessVector : TNvector);
  97.  
  98. {----------------------------------------------}
  99. {- Output: Dimen, Mat, GuessVector            -}
  100. {-                                            -}
  101. {- This procedure assigns values to the above -}
  102. {- variables from file input.                 -}
  103. {----------------------------------------------}
  104.  
  105. var
  106.   FileName : string[255];
  107.   InFile : text;
  108.   Row, Column : integer;
  109.  
  110. begin
  111.   Writeln;
  112.   repeat
  113.     Writeln;
  114.     repeat
  115.       Write('File name? ');
  116.       Readln(FileName);
  117.       Assign(InFile,FileName);
  118.       Reset(InFile);
  119.       IOCheck;
  120.     until not IOerr;
  121.     Read(InFile, Dimen);
  122.     IOCheck;
  123.     Row := 0;
  124.     while (not IOerr) and (Row < Dimen) do
  125.     begin
  126.       Row := Succ(Row);
  127.       Column := 0;
  128.       while (not IOerr) and (Column < Dimen) do
  129.       begin
  130.         Column := Succ(Column);
  131.         Read(InFile, Mat[Row, Column]);
  132.         IOCheck;
  133.       end;
  134.     end;
  135.     Row := 0;
  136.     while (not IOerr) and (Row < Dimen) do
  137.     begin
  138.       Row := Succ(Row);
  139.       Read(InFile, GuessVector[Row]);
  140.       IOCheck;
  141.     end;
  142.   until not IOerr;
  143.   Close(InFile);
  144. end; { procedure GetDataFromFile }
  145.  
  146.  
  147. procedure GetClosestVal(Dimen       : integer;
  148.                     var Mat         : TNmatrix;
  149.                     var GuessVector : TNvector;
  150.                     var ClosestVal  : Float);
  151.  
  152. {-------------------------------------------------------------}
  153. {- Input: Dimen, Mat, GuessVector                            -}
  154. {- Output: ClosestVal                                        -}
  155. {-                                                           -}
  156. {- ClosestVal is the user's approximation of the eigenvalue. -}
  157. {- This procedure makes an estimate of the eigenvalue based  -}
  158. {- based on the initial approximation to the eigenvalue      -}
  159. {- (GuessVector) and the matrix (Mat).  The formula used in  -}
  160. {- making this estimate is:                                  -}
  161. {-                                                           -}
  162. {-                (GuessVector - (Mat # GuessVector))        -}
  163. {-   ClosestVal = ------------------------------------       -}
  164. {-                    (GuessVector - GuessVector)            -}
  165. {-                                                           -}
  166. {- where - is a dot product and # is matrix multiplication.  -}
  167. {- If the user wants to input a different value for          -}
  168. {- ClosestVal, she may do so.                                -}
  169. {-------------------------------------------------------------}
  170.  
  171. var
  172.   Numerator, Denominator : Float;
  173.   TempVector : TNvector;
  174.  
  175. procedure Mult_Mat_Vec(Dimen  : integer;
  176.                    var Mat    : TNmatrix;
  177.                    var Vec    : TNvector;
  178.                    var Result : TNvector);
  179.  
  180. {-----------------------------------------}
  181. {- Input: Dimen, Mat, Vec                -}
  182. {- Output: Result                        -}
  183. {-                                       -}
  184. {- Multiply a square matrix and a vector -}
  185. {-----------------------------------------}
  186.  
  187. var
  188.   Row, Column : integer;
  189.   Entry : Float;
  190.  
  191. begin
  192.   for Row := 1 to Dimen do
  193.   begin
  194.     Entry := 0;
  195.     for Column := 1 to Dimen do
  196.       Entry := Entry + Mat[Row, Column] * Vec[Column];
  197.     Result[Row] := Entry;
  198.   end;
  199. end; { procedure Mult_Mat_Vec }
  200.  
  201. procedure Dot_Product(Dimen      : integer;
  202.                   var Vec1, Vec2 : TNvector;
  203.                   var Result     : Float);
  204.  
  205. {--------------------------------------------}
  206. {- Input: Dimen, Vec1, Vec2                 -}
  207. {- Output: Result                           -}
  208. {-                                          -}
  209. {- Calculate the dot product of two vectors -}
  210. {--------------------------------------------}
  211.  
  212. var
  213.   Row : integer;
  214.  
  215. begin
  216.   Result := 0;
  217.   for Row := 1 to Dimen do
  218.     Result := Result + Vec1[Row] * Vec2[Row];
  219. end; { procedure Dot_Product }
  220.  
  221. begin { procedure GetClosestVal }
  222.   Dot_Product(Dimen, GuessVector, GuessVector, Denominator);
  223.   Mult_Mat_Vec(Dimen, Mat, GuessVector, TempVector);
  224.   Dot_Product(Dimen, GuessVector, TempVector, Numerator);
  225.   if Denominator <> 0 then
  226.     ClosestVal := Numerator/Denominator
  227.   else
  228.     ClosestVal := 0;
  229.   Writeln;
  230.   repeat
  231.     Write('Approximate eigenvalue: ');
  232.     ReadFloat(ClosestVal);
  233.     IOCheck;
  234.   until not IOerr;
  235. end; { procedure GetClosestVal }
  236.  
  237. procedure GetTolerance(var Tolerance : Float);
  238.  
  239. {-----------------------------------------------------------}
  240. {- Output: Tolerance                                       -}
  241. {-                                                         -}
  242. {- This procedure reads in the Tolerance from the keyboard -}
  243. {-----------------------------------------------------------}
  244.  
  245. begin
  246.   Writeln;
  247.   Tolerance := 1E-6;
  248.   repeat
  249.     Write('Tolerance (> 0): ');
  250.     ReadFloat(Tolerance);
  251.     IOCheck;
  252.     if Tolerance <= 0 then
  253.     begin
  254.       IOerr := true;
  255.       Tolerance := 1E-6;
  256.     end;
  257.   until not IOerr;
  258. end; { procedure GetTolerance }
  259.  
  260. procedure GetMaxIter(var MaxIter : integer);
  261.  
  262. {-------------------------------------------------}
  263. {- Output: MaxIter                               -}
  264. {-                                               -}
  265. {- This procedure reads in the maximum number of -}
  266. {- iterations from the keyboard                  -}
  267. {-------------------------------------------------}
  268.  
  269. begin
  270.   Writeln;
  271.   MaxIter := 200;
  272.   repeat
  273.     Write('Maximum number of iterations (> 0): ');
  274.     ReadInt(MaxIter);
  275.     IOCheck;
  276.     if MaxIter <= 0 then
  277.     begin
  278.       IOerr := true;
  279.       MaxIter := 200;
  280.     end;
  281.   until not IOerr;
  282. end; { procedure GetMaxIter }
  283.  
  284. begin { procedure GetData }
  285.   Dimen := 0;
  286.   FillChar(Mat, SizeOf(Mat), 0);
  287.   FillChar(GuessVector, SizeOf(GuessVector), 0);
  288.   case InputChannel('Input Data From') of
  289.     'K' : GetDataFromKeyboard(Dimen, Mat, GuessVector);
  290.     'F' : GetDataFromFile(Dimen, Mat, GuessVector);
  291.   end;
  292.   GetClosestVal(Dimen, Mat, GuessVector, ClosestVal);
  293.   GetTolerance(Tolerance);
  294.   GetMaxIter(MaxIter);
  295.   GetOutputFile(OutFile);
  296. end; { procedure GetData }
  297.  
  298. procedure Results(Dimen       : integer;
  299.               var Mat         : TNmatrix;
  300.                   ClosestVal  : Float;
  301.                   Tolerance   : Float;
  302.                   MaxIter     : integer;
  303.                   Eigenvector : TNvector;
  304.                   Eigenvalue  : Float;
  305.                   Iter        : integer;
  306.                   Error       : byte);
  307.  
  308. {---------------------------------------------}
  309. {- Output the results to the device OutFile. -}
  310. {---------------------------------------------}
  311.  
  312. var
  313.   Column, Row : integer;
  314.  
  315. begin
  316.   Writeln(OutFile);
  317.   Writeln(OutFile);
  318.   Writeln(OutFile, 'The matrix: ');
  319.   for Row := 1 to Dimen do
  320.   begin
  321.     for Column := 1 to Dimen do
  322.       Write(OutFile, Mat[Row, Column]);
  323.     Writeln(OutFile);
  324.   end;
  325.   Writeln(OutFile);
  326.   Writeln(OutFile, 'Approximate eigenvalue: ' : 30, ClosestVal);
  327.   Writeln(OutFile, 'Tolerance: ' : 30, Tolerance);
  328.   Writeln(OutFile, 'Maximum number of iterations: ' : 30, MaxIter);
  329.   Writeln(OutFile);
  330.   if Error = 4 then
  331.     DisplayWarning;
  332.   if Error in [1, 2, 3, 5] then
  333.     DisplayError;
  334.   case Error of
  335.     1 : Writeln(OutFile,
  336.                 'The dimension of the matrix must be greater than zero.');
  337.  
  338.     2 : Writeln(OutFile, 'The tolerance must be greater than zero.');
  339.  
  340.     3 : Writeln(OutFile,
  341.                 'Maximum number of iterations must be greater than zero.');
  342.  
  343.     4 : begin
  344.           Writeln(OutFile, 'Convergence did not occur after ',
  345.                             Iter, ' iterations.');
  346.           Writeln(OutFile);
  347.           Writeln(OutFile, 'The results of the last iteration:');
  348.         end;
  349.  
  350.     5 : begin
  351.           Writeln(OutFile, 'The inverse power method may not be able');
  352.           Writeln(OutFile,
  353.                   'to compute an eigenvalue/vector of this matrix.');
  354.           Writeln(OutFile,
  355.                   'Try inputing a different approximate eigenvalue.');
  356.         end;
  357.   end; { case }
  358.  
  359.   if Error in [0, 4] then
  360.   begin
  361.     Writeln(OutFile);
  362.     Writeln(OutFile, 'Number of iterations: ' : 30, Iter : 3);
  363.     Writeln(OutFile, ' The approximate eigenvector:');
  364.     for Row := 1 to Dimen do
  365.       Writeln(OutFile, Eigenvector[Row]);
  366.     Writeln(OutFile);
  367.     Writeln(OutFile, 'The associated eigenvalue: ': 30, Eigenvalue);
  368.   end;
  369. end;  {procedure Results}
  370.  
  371. begin { program InversePower }
  372.   ClrScr;
  373.   GetData(Dimen, Mat, ClosestVal, Tolerance, MaxIter, GuessVector);
  374.   InversePower(Dimen, Mat, GuessVector, ClosestVal, MaxIter, Tolerance,
  375.                Eigenvalue, Eigenvector, Iter, Error);
  376.   Results(Dimen, Mat, ClosestVal, Tolerance, MaxIter,
  377.           Eigenvector, Eigenvalue, Iter, Error);
  378.   Close(OutFile);
  379. end. { program InversePower }