Turbo Toolbox/Turbo_Toolbox.iso/sonderh1/matrix.inc

Text File | 1979-12-31 | 13.9 KB | 470 lines

(****************************************************************************)
(*                                                                          *)
(*      MATRIX.INC  -  Basisroutinen zur Behandlung reeller Matrizen        *)
(*                                                                          *)
(*                         ---  Vers 4.0  ---                               *)
(*                                                                          *)
(****************************************************************************)

(* folgende Konstanten-, Typen- und Variablen-Vereinbarungen muessen je
   nach verwendetem Compiler ev. in das Hauptprogramm uebernommen werden !  *)

CONST Zero       = 1.0e-11;                         (* Maschinengenauigkeit *)

TYPE  Range      = 1..MaxIndex;                            (* Indexbereich  *)
      Vector     = ARRAY [Range] OF REAL;
      Matrix     = RECORD
                      Rows,                                (* Zeilenzahl    *)
                      Columns: Range;                      (* Spaltenzahl   *)
                      Coeff  : ARRAY [Range] OF Vector;    (* Koeffizienten *)
                   END;
      MedString  = STRING [14];
      RealArray  = ARRAY [0..MaxIndex] OF REAL;

VAR   Sign       : INTEGER;                  (* Vorzeichen bei Determinante *)


(* ------------------------------------------------------------------------ *)
(*                     Fehlermeldungen ausgeben:                            *)

PROCEDURE Error (n: INTEGER);

BEGIN
  Write('Fehler bei ');
  CASE n OF
    1: WriteLn('Systemloesung:  mehr Gleichungen als Variablen!');
    2: WriteLn('Systemloesung:  System ist nicht loesbar!');
    3: WriteLn('Systemloesung:  erweiterte Matrix ausserhalb Indexbereich!');
    4: WriteLn('Addition/Subtraktion:  verschiedene Zeilengroessen !');
    5: WriteLn('Addition/Subtraktion:  verschiedene Spaltengroessen !');
    6: WriteLn('Multiplikation:  Spaltenzahl A  ungleich  Zeilenzahl B!');
    7: WriteLn('Determinante:  Matrix ist nicht quadratisch!');
    8: WriteLn('Eingabe:  zulaessiger Indexbereich ueberschritten!');
    9: WriteLn('Invertierung:  Matrix ist nicht quadratisch!');
   10: WriteLn('Invertierung:  Matrix ist singulaer!');
  END;
  WriteLn;
  Halt;                                             (* Programm abbrechen ! *)
END;

(* ------------------------------------------------------------------------ *)
(* Die Anweisungsfolge 'Assign (Medium, MedStr); ReWrite (Medium);' bzw.
   'ReSet (Medium)' ist bei anderen Compilern ev. zu 'ReWrite (Medium,
   MedStr);' bzw. 'ReSet (Medium, MedStr);' in den folgenden zwei Pro-
   zeduren zusammen zu fassen!                                              *)
(* ------------------------------------------------------------------------ *)
(*                          Matrizen-Ausgabe:                               *)

PROCEDURE WriteMat (MedStr: MedString;       (* Ausgabemedium (z.B. 'LST:') *)
                    A     : Matrix;          (* auszugebende Matrix         *)
                    m,n   : Byte);           (* Formatparameter bei Ausgabe *)

VAR i,j   : Range;
    k     : INTEGER;
    Medium: TEXT;

BEGIN
  FOR k := 1 TO Length(MedStr) DO   (* 'MedStr' in Grossbuchstaben wandeln. *)
    MedStr[k] := UpCase(MedStr[k]);
  Assign(Medium, MedStr);          (* Zuordnung des Dateinamens zu 'Medium' *)
  With A DO
  BEGIN
    IF NOT ((MedStr = 'LST:') OR (MedStr = 'CON:')) THEN
      BEGIN
        ReWrite(Medium);
        WriteLn(Medium, Rows);
        WriteLn(Medium, Columns);
      END
    ELSE
      ReSet(Medium);
    FOR i := 1 TO Rows DO
    BEGIN
      FOR j := 1 TO Columns DO
        Write(Medium,Coeff[i,j]:m:n);
      WriteLn(Medium);
    END;
  END;
  Close(Medium);
END;

(* ------------------------------------------------------------------------ *)
(*                            Matrizen-Eingabe:                             *)

PROCEDURE ReadMat (    MedStr: MedString;    (* Eingabemedium (z.B. 'CON:') *)
                   VAR A     : Matrix);      (* einzulesende Matrix         *)

VAR i,j   : Range;
    k     : INTEGER;
    Medium: TEXT;

BEGIN
  FOR k := 1 TO Length(MedStr) DO
    MedStr[k] := UpCase(MedStr[k]);
  Assign(Medium, MedStr);
  Reset(Medium);
  WriteLn;
  WITH A DO
  BEGIN
    IF MedStr = 'CON:' THEN
      Write('Zeilenzahl : ');
    ReadLn(Medium, Rows);
    IF MedStr = 'CON:' THEN
      Write('Spaltenzahl: ');
    ReadLn(Medium, Columns);
    IF MedStr = 'CON:' THEN
      WriteLn('-----------------');
    IF NOT ([Rows, Columns] <= [1..MaxIndex]) THEN
      Error(8);
    FOR i := 1 TO Rows DO
    BEGIN
      IF MedStr = 'CON:' THEN
        Write(i, '. Zeile:  ');
      FOR j := 1 TO Columns DO
        Read(Medium, Coeff[i,j]);
    END;
  END;
  IF MedStr = 'CON:' THEN
    WriteLn;
  Close(Medium);
END;

(****************************************************************************)
(*                             Utility-Prozeduren:                          *)
(* ------------------------------------------------------------------------ *)
(*                      Zeilen 'm' und 'n' vertauschen:                     *)

PROCEDURE SwapRow (VAR A: Matrix; m, n: Range);

VAR i: Range;
    h: REAL;

BEGIN
  WITH A DO
    FOR i := 1 TO Columns DO
    BEGIN
      h := Coeff[n,i];
      Coeff[n,i] := Coeff[m,i];
      Coeff[m,i] := h
    END;
END;

(* ------------------------------------------------------------------------ *)
(*                    Spalten 'm' und 'n' vertauschen:                      *)

PROCEDURE SwapCol (VAR A: Matrix; m, n: Range);

VAR i: Range;
    h: REAL;

BEGIN
  WITH A DO
    FOR i := 1 TO Rows DO
    BEGIN
      h := Coeff[i,n];
      Coeff[i,n] := Coeff[i,m];
      Coeff[i,m] := h;
    END;
END;

(* ------------------------------------------------------------------------ *)
(*                           Pivot-Zeile suchen:                            *)

FUNCTION PivotRow (A: Matrix; n: Range): Range;

VAR i, r: Range;
    p   : REAL;

BEGIN
  WITH A DO
  BEGIN
    p := Coeff[n,n];
    r := n;
    FOR i := Succ(n) TO Rows DO
      IF Abs(Coeff[i,n]) > Abs(p) THEN
      BEGIN
        p := Coeff[i,n];
        r := i;
      END;
  END;
  PivotRow := r;
END;

(* ------------------------------------------------------------------------ *)
(*               Matrix auf obere Dreiecksform bringen:                     *)

PROCEDURE MakeLU (VAR A: Matrix);

VAR i, j, n, p: Range;

  (* ---------------------------------------------------------------------- *)

  PROCEDURE SubtractRow (m, n: Range);

  VAR i      : Range;
      Lambda,
      Pivot  : REAL;

  BEGIN
    WITH A DO
    BEGIN
      Pivot  := Coeff[m,m];
      Lambda := Coeff[n,m];
      FOR i := 1 TO Columns DO
        IF Abs(Pivot) >= Zero THEN
          Coeff[n,i] := Coeff[n,i] - Lambda/Pivot * Coeff[m,i]
        ELSE
          Coeff[n,i] := 0;
    END;
  END;

  (* ---------------------------------------------------------------------- *)

BEGIN (* MakeLU *)
  Sign := 1;
  n := A.Rows;
  FOR i := 1 TO n DO
  BEGIN
    p := PivotRow(A,i);
    IF p <> i THEN
    BEGIN
      SwapRow(A,p,i);
      Sign := -Sign;
    END;
    FOR j := Succ(i) TO n DO
      SubtractRow (i,j);
  END;
END;

(****************************************************************************)
(*                   Initialisierung als Einheitsmatrix:                    *)

PROCEDURE SetUpId (VAR A: Matrix; n: Range);

VAR i, j: Range;

BEGIN
  WITH A DO
  BEGIN
    Rows := n;
    Columns := n;
    FOR i := 1 TO n DO
      FOR j := 1 TO n DO
        Coeff[i,j] := Ord(i=j);
  END;
END;

(* ------------------------------------------------------------------------ *)
(*                           Matrizen-Addition:                             *)

PROCEDURE AddMat (A, B: Matrix; VAR C: Matrix);              (*  A + B = C  *)

VAR i, j: Range;

BEGIN
  IF A.Rows <> B.Rows THEN
    Error(4);
  IF A.Columns <> B.Columns THEN
    Error(5);
  C.Rows := A.Rows;
  C.Columns := A.Columns;
  FOR i := 1 TO A.Rows DO
    FOR j := 1 TO A.Columns DO
      C.Coeff[i,j] := A.Coeff[i,j] + B.Coeff[i,j];
END;

(* ------------------------------------------------------------------------ *)
(*                           Matrizen-Subtraktion:                          *)

PROCEDURE SubMat (A, B: Matrix; VAR C: Matrix);              (*  A - B = C  *)

VAR i, j: Range;

BEGIN
  WITH B DO
    FOR i := 1 TO Rows DO
      FOR j := 1 TO Columns DO
        Coeff[i,j] := -Coeff[i,j];
  AddMat(A,B,C);
END;

(* ------------------------------------------------------------------------ *)
(*                         Matrizen-Multiplikation:                         *)

PROCEDURE MultMat (A, B: Matrix; VAR C: Matrix);             (*  A * B = C  *)

VAR i, j, k: Range;
      Sigma: REAL;

BEGIN
  IF A.Columns <> B.Rows THEN
    Error(6);
  C.Rows := A.Rows;
  C.Columns := B.Columns;
  FOR i := 1 TO A.Rows DO
    FOR j := 1 TO B.Columns DO
    BEGIN
      Sigma := 0;
      FOR k := 1 TO A.Columns DO
        Sigma := Sigma + A.Coeff[i,k] * B.Coeff[k,j];
      C.Coeff[i,j] := Sigma;
    END;
END;

(* ------------------------------------------------------------------------ *)
(*                   Determinante einer Matrix berechnen:                   *)

FUNCTION Det (A: Matrix): REAL;

VAR d: REAL;
    i: Range;

BEGIN
  IF A.Rows <> A.Columns THEN
    Error(7);
  MakeLU(A);
  d := Sign;
  WITH A DO
    FOR i := 1 TO Rows DO
      d := d * Coeff[i,i];
  Det := d;
END;

(* ------------------------------------------------------------------------ *)
(*                       Inverse Matrix berechnen:                          *)

PROCEDURE Inverse (VAR A: Matrix);

VAR i, j, n, p,
         NoExch: Range;                     (* Anzahl Spaltenvertauschungen *)
    Exch: ARRAY [Range,1..2] OF INTEGER;    (* Protokoll         -          *)

  (* ---------------------------------------------------------------------- *)

  PROCEDURE Transform (m: Range);           (* Transformationsvorschrift    *)

  VAR B    : Matrix;
      i, j : Range;
      Pivot: REAL;

  BEGIN
    B := A;
    WITH B DO
    BEGIN
      Pivot := Coeff[m,m];
      IF Abs(Pivot) <= Zero THEN
        Error(10);
      FOR i := 1 TO n DO
        FOR j := 1 TO n DO
          IF i <> m THEN
            IF j <> m THEN
              Coeff[i,j] := B.Coeff[i,j] - A.Coeff[i,m] * A.Coeff[m,j] / Pivot
            ELSE
              Coeff[i,j] := A.Coeff[i,j] / Pivot
          ELSE
            If j <> m THEN
              Coeff[i,j] := -A.Coeff[i,j] / Pivot
            ELSE
              Coeff[i,j] := 1 / Pivot;
    END;
    A := B;
  END;

  (* ---------------------------------------------------------------------- *)

BEGIN (* Inverse *)
  NoExch := 0;
  IF A.Rows <> A.Columns THEN
    Error(9);
  n := A.Rows;

  (* Transformationsvorschrift anwenden *)

  FOR i := 1 TO n DO
  BEGIN
    p := PivotRow(A,i);
    IF p <> i THEN
    BEGIN
      NoExch := Succ(NoExch);
      Exch[NoExch,1] := i;
      Exch[NoExch,2] := p;
      SwapRow(A,i,p);
    END;
    Transform(i);
  END;

  (* vorgenommene Spaltenvertauschungen rueckgaengig machen *)

  FOR i := NoExch DOWNTO 1 DO
    SwapCol(A,Exch[i,2],Exch[i,1]);
END;

(* ------------------------------------------------------------------------ *)
(*                   Lineares Gleichungssystem loesen:                      *)

PROCEDURE SolveSystem (    A: Matrix;         (* Koeffizientenmatrix        *)
                           b: Vector;         (* rechte Seite               *)
                       VAR x: Vector;         (* Loesung inhomogenes System *)
                       VAR U: Matrix);        (* Loesung homogenes System   *)

VAR i, j, k, n, m: Range;

BEGIN
  IF A.Rows > A.Columns THEN
    Error(1);                            (* mehr Gleichungen als Variablen? *)
  If A.Columns = MaxIndex THEN
    Error(3);                            (* Matrix ausserhalb Indexbereich? *)
  WITH A DO                                 (* erweiterte Matrix generieren *)
  BEGIN
    Columns := Succ(Columns);
    FOR i := 1 TO Rows DO
      Coeff[i,Columns] := b[i];
    FOR i := 1 TO Columns-Rows DO
    BEGIN
      Rows := Succ(Rows);
      FOR j := 1 TO Columns DO
        Coeff[Rows,j] := 0;
    END;
  END;
  MakeLU(A);                       (* Matrix auf obere Dreiecksform bringen *)
  m := Pred(A.Columns);                       (* m = Anzahl der Variablen   *)
  n := m;
  WITH A DO
    WHILE Abs(Coeff[n,n]) < Zero DO
    BEGIN
      IF Abs(Coeff[n,Columns]) > Zero THEN
        Error(2);
      n := Pred(n);                           (* n = Anzahl der Gleichungen *)
    END;
  FOR i := 1 TO n DO                          (* Vektor x initialisieren    *)
    x[i] := A.Coeff[i,Succ(m)];
  FOR i := Succ(n) TO m DO
    x[i] := 0;
  WITH U DO                                   (* Matrix U initialisieren    *)
  BEGIN
    Rows := m;
    Columns := m-n;
    FOR i := 1 TO Rows DO
      FOR j := 1 TO Columns DO
        Coeff[i,j] := 0;
    FOR i := 1 TO Columns DO                  (* freie Parameter eintragen  *)
      Coeff[Succ(Rows-i),i] := 1;
    FOR i := n DOWNTO 1 DO                              (* Auswertung von x *)
    BEGIN
      FOR j := n DOWNTO Succ(i) DO
        x[i] := x[i] - A.Coeff[i,j] * x[j];
      x[i] := x[i] / A.Coeff[i,i];
    END;
    FOR i := n DOWNTO 1 DO                              (* Auswertung von U *)
      FOR j := 1 TO m-n DO
      BEGIN
        Coeff[i,j] := -A.Coeff[i,Succ(m-j)];
        FOR k := n DOWNTO Succ(i) DO
          Coeff[i,j] := Coeff[i,j] - A.Coeff[i,k] * Coeff[k,j];
        Coeff[i,j] := Coeff[i,j] / A.Coeff[i,i];
      END;
  END;
END;

(* ------------------------------------------------------------------------ *)
(*                           Ende von MATRIX.INC                            *)