home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 1.ddi / CHAP7.ARC / EIGENVAL.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1987-12-30  |  18.2 KB  |  286 lines
  1. unit EigenVal;
  2.  
  3. {----------------------------------------------------------------------------}
  4. {-                                                                          -}
  5. {-     Turbo Pascal Numerical Methods Toolbox                               -}
  6. {-     Copyright (c) 1986, 87 by Borland International, Inc.                -}
  7. {-                                                                          -}
  8. {-  This unit provides procedures for finding eigenvalues and eigenvectors. -}
  9. {-                                                                          -}
  10. {----------------------------------------------------------------------------}
  11.  
  12. {$I Float.inc} { Determines the setting of the $N compiler directive }
  13.  
  14. interface
  15.  
  16. {$IFOPT N+}
  17. type
  18.   Float = Double; { 8 byte real, requires 8087 math chip }
  19.  
  20. const
  21.   TNNearlyZero = 1E-015;
  22. {$ELSE}
  23. type
  24.   Float = real;   { 6 byte real, no math chip required }
  25.  
  26. const
  27.   TNNearlyZero = 1E-07;
  28. {$ENDIF}
  29.  
  30.   TNArraySize = 30;               { Maximum size of matrix }
  31.  
  32. type
  33.   TNvector    = array[1..TNArraySize] of Float;
  34.   TNmatrix    = array[1..TNArraySize] of TNvector;
  35.   TNIntVector = array[1..TNArraySize] of integer;
  36.  
  37. procedure Power(Dimen       : integer;
  38.             var Mat         : TNmatrix;
  39.             var GuessVector : TNvector;
  40.                 MaxIter     : integer;
  41.                 Tolerance   : Float;
  42.             var Eigenvalue  : Float;
  43.             var Eigenvector : TNvector;
  44.             var Iter        : integer;
  45.             var Error       : byte);
  46.  
  47. {----------------------------------------------------------------------------}
  48. {-                                                                          -}
  49. {-              Input:  Dimen, Mat, GuessVector, MaxIter, Tolerance         -}
  50. {-             Output:  Eigenvalue, Eigenvector, Iter, Error                -}
  51. {-                                                                          -}
  52. {-            Purpose:  The power method approximates the dominant          -}
  53. {-                      eigenvalue of a matrix.  The dominant               -}
  54. {-                      eigenvalue is the eigenvalue of largest             -}
  55. {-                      absolute magnitude. Given a square matrix Mat       -}
  56. {-                      and an arbitrary vector OldApprox, the vector       -}
  57. {-                      NewApprox is constructed by the matrix              -}
  58. {-                      operation NewApprox = Mat - OldApprox .             -}
  59. {-                      NewApprox is divided by its largest element         -}
  60. {-                      ApproxEigenval, thereby normalizing                 -}
  61. {-                      NewApprox. If NewApprox is the same as              -}
  62. {-                      OldApprox then ApproxEigenval is the dominant       -}
  63. {-                      eigenvalue and NewApprox is the associated          -}
  64. {-                      eigenvector of the matrix Mat. If NewApprox         -}
  65. {-                      is not the same as OldApprox then OldApprox         -}
  66. {-                      is set equal to NewApprox and the operation         -}
  67. {-                      repeats until a solution is reached. Aitken's       -}
  68. {-                      delta-squared acceleration is used to speed         -}
  69. {-                      the convergence from linear to quadratic.           -}
  70. {-                                                                          -}
  71. {- User-defined Types:  TNvector = array[1..TNArraySize] of real;           -}
  72. {-                      TNmatrix = array[1..TNArraySize] of TNvector;       -}
  73. {-                                                                          -}
  74. {-   Global Variables:  Dimen       : integer;  Dimension of the matrix     -}
  75. {-                      Mat         : TNmatrix; The matrix                  -}
  76. {-                      GuessVector : TNvector; An initial guess of an      -}
  77. {-                                              eigenvector                 -}
  78. {-                      MaxIter     : integer;  Max. number of iterations   -}
  79. {-                      Tolerance   : real;     Tolerance in answer         -}
  80. {-                      Eigenvalue  : real;     Eigenvalue of the matrix    -}
  81. {-                      Eigenvector : TNvector  Eigenvector of the matrix   -}
  82. {-                      Iter        : integer;  Number of iterations        -}
  83. {-                      Error       : byte;     Flags if something goes     -}
  84. {-                                              wrong                       -}
  85. {-             Errors:  0: No Errors                                        -}
  86. {-                      1: Dimen < 2                                        -}
  87. {-                      2: Tolerance <= 0                                   -}
  88. {-                      3: MaxIter < 1                                      -}
  89. {-                      4: Iter >= MaxIter                                  -}
  90. {-                                                                          -}
  91. {----------------------------------------------------------------------------}
  92.  
  93. procedure InversePower(Dimen       : integer;
  94.                        Mat         : TNmatrix;
  95.                    var GuessVector : TNvector;
  96.                        ClosestVal  : Float;
  97.                        MaxIter     : integer;
  98.                        Tolerance   : Float;
  99.                    var Eigenvalue  : Float;
  100.                    var Eigenvector : TNvector;
  101.                    var Iter        : integer;
  102.                    var Error       : byte);
  103.  
  104. {----------------------------------------------------------------------------}
  105. {-                                                                          -}
  106. {-          Input:  Dimen, Mat, GuessVector, ClosestVal, MaxIter, Tolerance -}
  107. {-         Output:  Eigenvalue, Eigenvector, Iter, Error                    -}
  108. {-                                                                          -}
  109. {-            Purpose:  Whereas the power method converges onto the         -}
  110. {-                      dominant eigenvalue of a matrix (see                -}
  111. {-                      POWER.INC), the inverse power method                -}
  112. {-                      converges onto the eigenvalue closest to a          -}
  113. {-                      user-supplied value.  The user supplies a           -}
  114. {-                      square matrix Mat, an initial approximation         -}
  115. {-                      ClosestVal to the eigenvalue and an initial         -}
  116. {-                      vector OldApprox.  The linear system                -}
  117. {-                      (Mat - ClosestVal - I)OldApprox = NewApprox is      -}
  118. {-                      solved via LU decomposition (see                    -}
  119. {-                      DECMP-LU.INC\SOLVE-LU.INC).  The vector             -}
  120. {-                      NewApprox is divided by its largest element         -}
  121. {-                      ApproxEigenval (thereby normalizing the             -}
  122. {-                      vector).  If NewApprox is identical to              -}
  123. {-                      OldApprox then (1/ApproxEigenval + ClosestVal)      -}
  124. {-                      is the eigenvalue of A closest to ClosestVal        -}
  125. {-                      and OldApprox is the associated eigenvector.        -}
  126. {-                      If NewApprox is not identical to OldApprox          -}
  127. {-                      then OldApprox is set equal to NewApprox and        -}
  128. {-                      the process repeats until a solution is             -}
  129. {-                      reached.                                            -}
  130. {-                                                                          -}
  131. {- User-defined Types:  TNvector = array[1..TNArraySize] of real;           -}
  132. {-                      TNmatrix = array[1..TNArraySize] of TNvector;       -}
  133. {-                                                                          -}
  134. {-   Global Variables:  Dimen       : integer;  Dimension of the matrix     -}
  135. {-                      Mat         : TNmatrix; The matrix                  -}
  136. {-                      GuessVector : TNvector; Initial guess of an         -}
  137. {-                                              Eigenvector                 -}
  138. {-                      ClosestVal  : real;     Converge to eigenvalue      -}
  139. {-                                              Closest to this             -}
  140. {-                      MaxIter     : integer;  Max. number of iterations   -}
  141. {-                      Tolerance   : real;     Tolerance in answer         -}
  142. {-                      Eigenvalue  : real;     Eigenvalue of the matrix    -}
  143. {-                      Eigenvector : TNvector  Eigenvector of the matrix   -}
  144. {-                      Iter        : integer;  Number of iterations        -}
  145. {-                      Error       : byte;     Flags if something goes     -}
  146. {-                                              wrong                       -}
  147. {-             Errors:  0: No Errors                                        -}
  148. {-                      1: Dimen < 2                                        -}
  149. {-                      2: Tolerance <= 0                                   -}
  150. {-                      3: MaxIter < 1                                      -}
  151. {-                      4: Iter >= MaxIter                                  -}
  152. {-                      5: eigenvalue/eigenvector not calculated            -}
  153. {-                         See note below.                                  -}
  154. {-                                                                          -}
  155. {-               Note:  If the matrix Mat - EigenValue-I, where I           -}
  156. {-                      is the identity matrix, is singular, then           -}
  157. {-                      the inverse power method may not be used            -}
  158. {-                      to approximate an eigenvalue and eigenvector        -}
  159. {-                      of Mat and Error = 5 will be returned.              -}
  160. {-                                                                          -}
  161. {----------------------------------------------------------------------------}
  162.  
  163. procedure Wielandt(Dimen        : integer;
  164.                    Mat          : TNmatrix;
  165.                var GuessVector  : TNvector;
  166.                    MaxEigens    : integer;
  167.                    MaxIter      : integer;
  168.                    Tolerance    : Float;
  169.                var NumEigens    : integer;
  170.                var Eigenvalues  : TNvector;
  171.                var Eigenvectors : TNmatrix;
  172.                var Iter         : TNIntVector;
  173.                var Error        : byte);
  174.  
  175. {----------------------------------------------------------------------------}
  176. {-                                                                          -}
  177. {-        Input: Dimen, Mat, GuessVector, MaxEigens, MaxIter, Tolerance     -}
  178. {-       Output: NumEigens, Eigenvalues, Eigenvectors, Iter, Error          -}
  179. {-                                                                          -}
  180. {-            Purpose:  This procedure attempts to approximate some (or     -}
  181. {-                      all, depending on MaxEigens)  of the                -}
  182. {-                      eigenvectors and eigenvalues of a matrix.  The      -}
  183. {-                      power method is used in conjunction with Wielandt's -}
  184. {-                      deflation.                                          -}
  185. {-                                                                          -}
  186. {- User-defined Types:  TNvector = array[1..RowSize] of real;               -}
  187. {-                      TNmatrix = array[1..ColumnSize] of TNvector;        -}
  188. {-                      TNIntVector = array[1..RowSize] of integer;         -}
  189. {-                                                                          -}
  190. {-   Global Variables:  Dimen        : integer;  Dimension of the matrix    -}
  191. {-                      Mat          : TNmatrix; The matrix                 -}
  192. {-                      GuessVector  : TNvector; An initial guess of an     -}
  193. {-                                               eigenvector                -}
  194. {-                      MaxEigens    : integer;  Maximum number of eigens   -}
  195. {-                                               to find                    -}
  196. {-                      MaxIter      : integer;  Max. number of iterations  -}
  197. {-                      Tolerance    : real;     Tolerance in answer        -}
  198. {-                      NumEigens    : integer;  Number of eigenvalues      -}
  199. {-                                               calculated                 -}
  200. {-                      Eigenvalues  : TNvector; Eigenvalues of the matrix  -}
  201. {-                      Eigenvectors : TNmatrix  Eigenvectors of the matrix -}
  202. {-                      Iter         : TNInTVector; Number of iterations    -}
  203. {-                      Error        : byte;     Flags if something goes    -}
  204. {-                                               wrong                      -}
  205. {-             Errors:  0: No Errors                                        -}
  206. {-                      1: Dimen < 2                                        -}
  207. {-                      2: Tolerance <= 0                                   -}
  208. {-                      3: MaxIter < 1                                      -}
  209. {-                      4: MaxEigens < 1                                    -}
  210. {-                      5: Iter >= MaxIter                                  -}
  211. {-                      6: Last two roots aren't real                       -}
  212. {-                                                                          -}
  213. {----------------------------------------------------------------------------}
  214.  
  215. procedure Jacobi(Dimen        : integer;
  216.                  Mat          : TNmatrix;
  217.                  MaxIter      : integer;
  218.                  Tolerance    : Float;
  219.              var Eigenvalues  : TNvector;
  220.              var Eigenvectors : TNmatrix;
  221.              var Iter         : integer;
  222.              var Error        : byte);
  223.  
  224. {----------------------------------------------------------------------------}
  225. {-                                                                          -}
  226. {-     Input: Dimen, Mat, MaxIter, Tolerance                                -}
  227. {-    Output: Eigenvalues, Eigenvector, Iter, Error                         -}
  228. {-                                                                          -}
  229. {-   Purpose: The eigensystem of a symmetric matrix can be                  -}
  230. {-            computed much more simply than the                            -}
  231. {-            eigensystem of an arbitrary matrix.  The                      -}
  232. {-            cyclic Jacobi method is an iterative                          -}
  233. {-            technique for approximating the complete                      -}
  234. {-            eigensystem of a symmetric matrix to a given                  -}
  235. {-            tolerance. The method consists of multiplying                 -}
  236. {-            the matrix, Mat, by a series of rotation                      -}
  237. {-            matrices, r@-[i].  The rotation matrices are                  -}
  238. {-            chosen so that the elements of the upper                      -}
  239. {-            triangular part of Mat are systematically                     -}
  240. {-            annihilated.  That is, r@-[1] is chosen so                    -}
  241. {-            that Mat[1, 1] is identically zero; r@-[2] is                 -}
  242. {-            chosen so that Mat[1, 2] is identically zero;                 -}
  243. {-            etc.  Since each operation will probably                      -}
  244. {-            change the value of elements annihilated in                   -}
  245. {-            previous operations, the method is iterative.                 -}
  246. {-            Eventually, the matrix will be diagonal. The                  -}
  247. {-            eigenvalues will be the elements along the                    -}
  248. {-            diagonal of the matrix and the eigenvectors                   -}
  249. {-            will be the rows of the matrix created by the                 -}
  250. {-            product of all the rotation matrices r@-[i].                  -}
  251. {-                                                                          -}
  252. {-            The user inputs the matrix, tolerance and maximum             -}
  253. {-            number of iterations. The procedure returns the               -}
  254. {-            eigenvalues and eigenvectors (or error code) of the           -}
  255. {-            matrix.                                                       -}
  256. {-                                                                          -}
  257. {-   User-Defined Types: TNvector = array[1..TNArraySize] of real;          -}
  258. {-                       TNmatrix = array[1..TNArraySize] of TNvector;      -}
  259. {-                                                                          -}
  260. {-   Global Variables:  Dimen        : integer   Dimension of square matrix -}
  261. {-                      Mat          : TNmatrix  Square matrix              -}
  262. {-                      MaxIter      : integer   Maximum number of          -}
  263. {-                                               Iterations                 -}
  264. {-                      Tolerance    : real      Tolerance in answer        -}
  265. {-                      Eigenvalues  : TNvector  Eigenvalues of Mat         -}
  266. {-                      Eigenvectors : TNmatrix  Eigenvectors of Mat        -}
  267. {-                      Iter         : integer   Number of iterations       -}
  268. {-                      Error        : byte      Error code                 -}
  269. {-                                                                          -}
  270. {-             Errors:  0: No error                                         -}
  271. {-                      1: Dimen < 1                                        -}
  272. {-                      2: Tolerance < TNNearlyZero                         -}
  273. {-                      3: MaxIter < 1                                      -}
  274. {-                      4: Mat not symmetric                                -}
  275. {-                      5: Iter > MaxIter                                   -}
  276. {-                                                                          -}
  277. {----------------------------------------------------------------------------}
  278.  
  279. implementation
  280.  
  281. {$I Eigen1.inc}   { Include procedure code }
  282.  
  283. {$I Eigen2.inc}
  284.  
  285. end. { EigenVal }
  286.