home *** CD-ROM | disk | FTP | other *** search
/ Turbo Toolbox / Turbo_Toolbox.iso / sonderh1 / matrix.inc < prev    next >
Text File  |  1979-12-31  |  14KB  |  470 lines
  1. (****************************************************************************)
  2. (*                                                                          *)
  3. (*      MATRIX.INC  -  Basisroutinen zur Behandlung reeller Matrizen        *)
  4. (*                                                                          *)
  5. (*                         ---  Vers 4.0  ---                               *)
  6. (*                                                                          *)
  7. (****************************************************************************)
  8.  
  9. (* folgende Konstanten-, Typen- und Variablen-Vereinbarungen muessen je
  10.    nach verwendetem Compiler ev. in das Hauptprogramm uebernommen werden !  *)
  11.  
  12. CONST Zero       = 1.0e-11;                         (* Maschinengenauigkeit *)
  13.  
  14. TYPE  Range      = 1..MaxIndex;                            (* Indexbereich  *)
  15.       Vector     = ARRAY [Range] OF REAL;
  16.       Matrix     = RECORD
  17.                       Rows,                                (* Zeilenzahl    *)
  18.                       Columns: Range;                      (* Spaltenzahl   *)
  19.                       Coeff  : ARRAY [Range] OF Vector;    (* Koeffizienten *)
  20.                    END;
  21.       MedString  = STRING [14];
  22.       RealArray  = ARRAY [0..MaxIndex] OF REAL;
  23.  
  24. VAR   Sign       : INTEGER;                  (* Vorzeichen bei Determinante *)
  25.  
  26.  
  27. (* ------------------------------------------------------------------------ *)
  28. (*                     Fehlermeldungen ausgeben:                            *)
  29.  
  30. PROCEDURE Error (n: INTEGER);
  31.  
  32. BEGIN
  33.   Write('Fehler bei ');
  34.   CASE n OF
  35.     1: WriteLn('Systemloesung:  mehr Gleichungen als Variablen!');
  36.     2: WriteLn('Systemloesung:  System ist nicht loesbar!');
  37.     3: WriteLn('Systemloesung:  erweiterte Matrix ausserhalb Indexbereich!');
  38.     4: WriteLn('Addition/Subtraktion:  verschiedene Zeilengroessen !');
  39.     5: WriteLn('Addition/Subtraktion:  verschiedene Spaltengroessen !');
  40.     6: WriteLn('Multiplikation:  Spaltenzahl A  ungleich  Zeilenzahl B!');
  41.     7: WriteLn('Determinante:  Matrix ist nicht quadratisch!');
  42.     8: WriteLn('Eingabe:  zulaessiger Indexbereich ueberschritten!');
  43.     9: WriteLn('Invertierung:  Matrix ist nicht quadratisch!');
  44.    10: WriteLn('Invertierung:  Matrix ist singulaer!');
  45.   END;
  46.   WriteLn;
  47.   Halt;                                             (* Programm abbrechen ! *)
  48. END;
  49.  
  50. (* ------------------------------------------------------------------------ *)
  51. (* Die Anweisungsfolge 'Assign (Medium, MedStr); ReWrite (Medium);' bzw.
  52.    'ReSet (Medium)' ist bei anderen Compilern ev. zu 'ReWrite (Medium,
  53.    MedStr);' bzw. 'ReSet (Medium, MedStr);' in den folgenden zwei Pro-
  54.    zeduren zusammen zu fassen!                                              *)
  55. (* ------------------------------------------------------------------------ *)
  56. (*                          Matrizen-Ausgabe:                               *)
  57.  
  58. PROCEDURE WriteMat (MedStr: MedString;       (* Ausgabemedium (z.B. 'LST:') *)
  59.                     A     : Matrix;          (* auszugebende Matrix         *)
  60.                     m,n   : Byte);           (* Formatparameter bei Ausgabe *)
  61.  
  62. VAR i,j   : Range;
  63.     k     : INTEGER;
  64.     Medium: TEXT;
  65.  
  66. BEGIN
  67.   FOR k := 1 TO Length(MedStr) DO   (* 'MedStr' in Grossbuchstaben wandeln. *)
  68.     MedStr[k] := UpCase(MedStr[k]);
  69.   Assign(Medium, MedStr);          (* Zuordnung des Dateinamens zu 'Medium' *)
  70.   With A DO
  71.   BEGIN
  72.     IF NOT ((MedStr = 'LST:') OR (MedStr = 'CON:')) THEN
  73.       BEGIN
  74.         ReWrite(Medium);
  75.         WriteLn(Medium, Rows);
  76.         WriteLn(Medium, Columns);
  77.       END
  78.     ELSE
  79.       ReSet(Medium);
  80.     FOR i := 1 TO Rows DO
  81.     BEGIN
  82.       FOR j := 1 TO Columns DO
  83.         Write(Medium,Coeff[i,j]:m:n);
  84.       WriteLn(Medium);
  85.     END;
  86.   END;
  87.   Close(Medium);
  88. END;
  89.  
  90. (* ------------------------------------------------------------------------ *)
  91. (*                            Matrizen-Eingabe:                             *)
  92.  
  93. PROCEDURE ReadMat (    MedStr: MedString;    (* Eingabemedium (z.B. 'CON:') *)
  94.                    VAR A     : Matrix);      (* einzulesende Matrix         *)
  95.  
  96. VAR i,j   : Range;
  97.     k     : INTEGER;
  98.     Medium: TEXT;
  99.  
  100. BEGIN
  101.   FOR k := 1 TO Length(MedStr) DO
  102.     MedStr[k] := UpCase(MedStr[k]);
  103.   Assign(Medium, MedStr);
  104.   Reset(Medium);
  105.   WriteLn;
  106.   WITH A DO
  107.   BEGIN
  108.     IF MedStr = 'CON:' THEN
  109.       Write('Zeilenzahl : ');
  110.     ReadLn(Medium, Rows);
  111.     IF MedStr = 'CON:' THEN
  112.       Write('Spaltenzahl: ');
  113.     ReadLn(Medium, Columns);
  114.     IF MedStr = 'CON:' THEN
  115.       WriteLn('-----------------');
  116.     IF NOT ([Rows, Columns] <= [1..MaxIndex]) THEN
  117.       Error(8);
  118.     FOR i := 1 TO Rows DO
  119.     BEGIN
  120.       IF MedStr = 'CON:' THEN
  121.         Write(i, '. Zeile:  ');
  122.       FOR j := 1 TO Columns DO
  123.         Read(Medium, Coeff[i,j]);
  124.     END;
  125.   END;
  126.   IF MedStr = 'CON:' THEN
  127.     WriteLn;
  128.   Close(Medium);
  129. END;
  130.  
  131. (****************************************************************************)
  132. (*                             Utility-Prozeduren:                          *)
  133. (* ------------------------------------------------------------------------ *)
  134. (*                      Zeilen 'm' und 'n' vertauschen:                     *)
  135.  
  136. PROCEDURE SwapRow (VAR A: Matrix; m, n: Range);
  137.  
  138. VAR i: Range;
  139.     h: REAL;
  140.  
  141. BEGIN
  142.   WITH A DO
  143.     FOR i := 1 TO Columns DO
  144.     BEGIN
  145.       h := Coeff[n,i];
  146.       Coeff[n,i] := Coeff[m,i];
  147.       Coeff[m,i] := h
  148.     END;
  149. END;
  150.  
  151. (* ------------------------------------------------------------------------ *)
  152. (*                    Spalten 'm' und 'n' vertauschen:                      *)
  153.  
  154. PROCEDURE SwapCol (VAR A: Matrix; m, n: Range);
  155.  
  156. VAR i: Range;
  157.     h: REAL;
  158.  
  159. BEGIN
  160.   WITH A DO
  161.     FOR i := 1 TO Rows DO
  162.     BEGIN
  163.       h := Coeff[i,n];
  164.       Coeff[i,n] := Coeff[i,m];
  165.       Coeff[i,m] := h;
  166.     END;
  167. END;
  168.  
  169. (* ------------------------------------------------------------------------ *)
  170. (*                           Pivot-Zeile suchen:                            *)
  171.  
  172. FUNCTION PivotRow (A: Matrix; n: Range): Range;
  173.  
  174. VAR i, r: Range;
  175.     p   : REAL;
  176.  
  177. BEGIN
  178.   WITH A DO
  179.   BEGIN
  180.     p := Coeff[n,n];
  181.     r := n;
  182.     FOR i := Succ(n) TO Rows DO
  183.       IF Abs(Coeff[i,n]) > Abs(p) THEN
  184.       BEGIN
  185.         p := Coeff[i,n];
  186.         r := i;
  187.       END;
  188.   END;
  189.   PivotRow := r;
  190. END;
  191.  
  192. (* ------------------------------------------------------------------------ *)
  193. (*               Matrix auf obere Dreiecksform bringen:                     *)
  194.  
  195. PROCEDURE MakeLU (VAR A: Matrix);
  196.  
  197. VAR i, j, n, p: Range;
  198.  
  199.   (* ---------------------------------------------------------------------- *)
  200.  
  201.   PROCEDURE SubtractRow (m, n: Range);
  202.  
  203.   VAR i      : Range;
  204.       Lambda,
  205.       Pivot  : REAL;
  206.  
  207.   BEGIN
  208.     WITH A DO
  209.     BEGIN
  210.       Pivot  := Coeff[m,m];
  211.       Lambda := Coeff[n,m];
  212.       FOR i := 1 TO Columns DO
  213.         IF Abs(Pivot) >= Zero THEN
  214.           Coeff[n,i] := Coeff[n,i] - Lambda/Pivot * Coeff[m,i]
  215.         ELSE
  216.           Coeff[n,i] := 0;
  217.     END;
  218.   END;
  219.  
  220.   (* ---------------------------------------------------------------------- *)
  221.  
  222. BEGIN (* MakeLU *)
  223.   Sign := 1;
  224.   n := A.Rows;
  225.   FOR i := 1 TO n DO
  226.   BEGIN
  227.     p := PivotRow(A,i);
  228.     IF p <> i THEN
  229.     BEGIN
  230.       SwapRow(A,p,i);
  231.       Sign := -Sign;
  232.     END;
  233.     FOR j := Succ(i) TO n DO
  234.       SubtractRow (i,j);
  235.   END;
  236. END;
  237.  
  238. (****************************************************************************)
  239. (*                   Initialisierung als Einheitsmatrix:                    *)
  240.  
  241. PROCEDURE SetUpId (VAR A: Matrix; n: Range);
  242.  
  243. VAR i, j: Range;
  244.  
  245. BEGIN
  246.   WITH A DO
  247.   BEGIN
  248.     Rows := n;
  249.     Columns := n;
  250.     FOR i := 1 TO n DO
  251.       FOR j := 1 TO n DO
  252.         Coeff[i,j] := Ord(i=j);
  253.   END;
  254. END;
  255.  
  256. (* ------------------------------------------------------------------------ *)
  257. (*                           Matrizen-Addition:                             *)
  258.  
  259. PROCEDURE AddMat (A, B: Matrix; VAR C: Matrix);              (*  A + B = C  *)
  260.  
  261. VAR i, j: Range;
  262.  
  263. BEGIN
  264.   IF A.Rows <> B.Rows THEN
  265.     Error(4);
  266.   IF A.Columns <> B.Columns THEN
  267.     Error(5);
  268.   C.Rows := A.Rows;
  269.   C.Columns := A.Columns;
  270.   FOR i := 1 TO A.Rows DO
  271.     FOR j := 1 TO A.Columns DO
  272.       C.Coeff[i,j] := A.Coeff[i,j] + B.Coeff[i,j];
  273. END;
  274.  
  275. (* ------------------------------------------------------------------------ *)
  276. (*                           Matrizen-Subtraktion:                          *)
  277.  
  278. PROCEDURE SubMat (A, B: Matrix; VAR C: Matrix);              (*  A - B = C  *)
  279.  
  280. VAR i, j: Range;
  281.  
  282. BEGIN
  283.   WITH B DO
  284.     FOR i := 1 TO Rows DO
  285.       FOR j := 1 TO Columns DO
  286.         Coeff[i,j] := -Coeff[i,j];
  287.   AddMat(A,B,C);
  288. END;
  289.  
  290. (* ------------------------------------------------------------------------ *)
  291. (*                         Matrizen-Multiplikation:                         *)
  292.  
  293. PROCEDURE MultMat (A, B: Matrix; VAR C: Matrix);             (*  A * B = C  *)
  294.  
  295. VAR i, j, k: Range;
  296.       Sigma: REAL;
  297.  
  298. BEGIN
  299.   IF A.Columns <> B.Rows THEN
  300.     Error(6);
  301.   C.Rows := A.Rows;
  302.   C.Columns := B.Columns;
  303.   FOR i := 1 TO A.Rows DO
  304.     FOR j := 1 TO B.Columns DO
  305.     BEGIN
  306.       Sigma := 0;
  307.       FOR k := 1 TO A.Columns DO
  308.         Sigma := Sigma + A.Coeff[i,k] * B.Coeff[k,j];
  309.       C.Coeff[i,j] := Sigma;
  310.     END;
  311. END;
  312.  
  313. (* ------------------------------------------------------------------------ *)
  314. (*                   Determinante einer Matrix berechnen:                   *)
  315.  
  316. FUNCTION Det (A: Matrix): REAL;
  317.  
  318. VAR d: REAL;
  319.     i: Range;
  320.  
  321. BEGIN
  322.   IF A.Rows <> A.Columns THEN
  323.     Error(7);
  324.   MakeLU(A);
  325.   d := Sign;
  326.   WITH A DO
  327.     FOR i := 1 TO Rows DO
  328.       d := d * Coeff[i,i];
  329.   Det := d;
  330. END;
  331.  
  332. (* ------------------------------------------------------------------------ *)
  333. (*                       Inverse Matrix berechnen:                          *)
  334.  
  335. PROCEDURE Inverse (VAR A: Matrix);
  336.  
  337. VAR i, j, n, p,
  338.          NoExch: Range;                     (* Anzahl Spaltenvertauschungen *)
  339.     Exch: ARRAY [Range,1..2] OF INTEGER;    (* Protokoll         -          *)
  340.  
  341.   (* ---------------------------------------------------------------------- *)
  342.  
  343.   PROCEDURE Transform (m: Range);           (* Transformationsvorschrift    *)
  344.  
  345.   VAR B    : Matrix;
  346.       i, j : Range;
  347.       Pivot: REAL;
  348.  
  349.   BEGIN
  350.     B := A;
  351.     WITH B DO
  352.     BEGIN
  353.       Pivot := Coeff[m,m];
  354.       IF Abs(Pivot) <= Zero THEN
  355.         Error(10);
  356.       FOR i := 1 TO n DO
  357.         FOR j := 1 TO n DO
  358.           IF i <> m THEN
  359.             IF j <> m THEN
  360.               Coeff[i,j] := B.Coeff[i,j] - A.Coeff[i,m] * A.Coeff[m,j] / Pivot
  361.             ELSE
  362.               Coeff[i,j] := A.Coeff[i,j] / Pivot
  363.           ELSE
  364.             If j <> m THEN
  365.               Coeff[i,j] := -A.Coeff[i,j] / Pivot
  366.             ELSE
  367.               Coeff[i,j] := 1 / Pivot;
  368.     END;
  369.     A := B;
  370.   END;
  371.  
  372.   (* ---------------------------------------------------------------------- *)
  373.  
  374. BEGIN (* Inverse *)
  375.   NoExch := 0;
  376.   IF A.Rows <> A.Columns THEN
  377.     Error(9);
  378.   n := A.Rows;
  379.  
  380.   (* Transformationsvorschrift anwenden *)
  381.  
  382.   FOR i := 1 TO n DO
  383.   BEGIN
  384.     p := PivotRow(A,i);
  385.     IF p <> i THEN
  386.     BEGIN
  387.       NoExch := Succ(NoExch);
  388.       Exch[NoExch,1] := i;
  389.       Exch[NoExch,2] := p;
  390.       SwapRow(A,i,p);
  391.     END;
  392.     Transform(i);
  393.   END;
  394.  
  395.   (* vorgenommene Spaltenvertauschungen rueckgaengig machen *)
  396.  
  397.   FOR i := NoExch DOWNTO 1 DO
  398.     SwapCol(A,Exch[i,2],Exch[i,1]);
  399. END;
  400.  
  401. (* ------------------------------------------------------------------------ *)
  402. (*                   Lineares Gleichungssystem loesen:                      *)
  403.  
  404. PROCEDURE SolveSystem (    A: Matrix;         (* Koeffizientenmatrix        *)
  405.                            b: Vector;         (* rechte Seite               *)
  406.                        VAR x: Vector;         (* Loesung inhomogenes System *)
  407.                        VAR U: Matrix);        (* Loesung homogenes System   *)
  408.  
  409. VAR i, j, k, n, m: Range;
  410.  
  411. BEGIN
  412.   IF A.Rows > A.Columns THEN
  413.     Error(1);                            (* mehr Gleichungen als Variablen? *)
  414.   If A.Columns = MaxIndex THEN
  415.     Error(3);                            (* Matrix ausserhalb Indexbereich? *)
  416.   WITH A DO                                 (* erweiterte Matrix generieren *)
  417.   BEGIN
  418.     Columns := Succ(Columns);
  419.     FOR i := 1 TO Rows DO
  420.       Coeff[i,Columns] := b[i];
  421.     FOR i := 1 TO Columns-Rows DO
  422.     BEGIN
  423.       Rows := Succ(Rows);
  424.       FOR j := 1 TO Columns DO
  425.         Coeff[Rows,j] := 0;
  426.     END;
  427.   END;
  428.   MakeLU(A);                       (* Matrix auf obere Dreiecksform bringen *)
  429.   m := Pred(A.Columns);                       (* m = Anzahl der Variablen   *)
  430.   n := m;
  431.   WITH A DO
  432.     WHILE Abs(Coeff[n,n]) < Zero DO
  433.     BEGIN
  434.       IF Abs(Coeff[n,Columns]) > Zero THEN
  435.         Error(2);
  436.       n := Pred(n);                           (* n = Anzahl der Gleichungen *)
  437.     END;
  438.   FOR i := 1 TO n DO                          (* Vektor x initialisieren    *)
  439.     x[i] := A.Coeff[i,Succ(m)];
  440.   FOR i := Succ(n) TO m DO
  441.     x[i] := 0;
  442.   WITH U DO                                   (* Matrix U initialisieren    *)
  443.   BEGIN
  444.     Rows := m;
  445.     Columns := m-n;
  446.     FOR i := 1 TO Rows DO
  447.       FOR j := 1 TO Columns DO
  448.         Coeff[i,j] := 0;
  449.     FOR i := 1 TO Columns DO                  (* freie Parameter eintragen  *)
  450.       Coeff[Succ(Rows-i),i] := 1;
  451.     FOR i := n DOWNTO 1 DO                              (* Auswertung von x *)
  452.     BEGIN
  453.       FOR j := n DOWNTO Succ(i) DO
  454.         x[i] := x[i] - A.Coeff[i,j] * x[j];
  455.       x[i] := x[i] / A.Coeff[i,i];
  456.     END;
  457.     FOR i := n DOWNTO 1 DO                              (* Auswertung von U *)
  458.       FOR j := 1 TO m-n DO
  459.       BEGIN
  460.         Coeff[i,j] := -A.Coeff[i,Succ(m-j)];
  461.         FOR k := n DOWNTO Succ(i) DO
  462.           Coeff[i,j] := Coeff[i,j] - A.Coeff[i,k] * Coeff[k,j];
  463.         Coeff[i,j] := Coeff[i,j] / A.Coeff[i,i];
  464.       END;
  465.   END;
  466. END;
  467.  
  468. (* ------------------------------------------------------------------------ *)
  469. (*                           Ende von MATRIX.INC                            *)
  470.