home *** CD-ROM | disk | FTP | other *** search
/ Turbo Toolbox / Turbo_Toolbox.iso / 1991 / 07_08 / tricks / dettool.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1991-04-19  |  3.5 KB  |  138 lines

  1. (* ------------------------------------------------------ *)
  2. (*                     DETTOOL.PAS                        *)
  3. (*          (c) 1991 Jens Burmeister & TOOLBOX            *)
  4. (* ------------------------------------------------------ *)
  5. UNIT DetTool;
  6.  
  7. INTERFACE
  8.  
  9. TYPE
  10.   RealMatrixPtr = ^RealMatrix;
  11.   RealMatrix    = OBJECT
  12.     Dim,
  13.     RowSize,
  14.     ColSize,
  15.     RealSize : WORD;
  16.     Buf      : Pointer;
  17.  
  18.     CONSTRUCTOR Init(N : WORD );
  19.     DESTRUCTOR  Done; VIRTUAL;
  20.     FUNCTION    GetDimension : WORD;
  21.     FUNCTION    GetValue(i, j : WORD) : REAL;
  22.     PROCEDURE   SetValue(i, j : WORD; Value : REAL);
  23.     FUNCTION    Determinant : REAL;
  24.   END;
  25.  
  26. IMPLEMENTATION
  27.  
  28. CONST
  29.   MaxNrOfPtrs  = 65520 DIV SizeOf(Pointer);
  30.   MaxNrOfReals = 65520 DIV SizeOf(REAL);
  31.  
  32. TYPE
  33.   Pointers = ARRAY [1..MaxNrOfPtrs]  OF Pointer;
  34.   Reals    = ARRAY [1..MaxNrOfReals] OF REAL;
  35.  
  36.   CONSTRUCTOR RealMatrix.Init(N : WORD);
  37.   VAR
  38.     i : WORD;
  39.   BEGIN
  40.     Dim := N;
  41.     RealSize := SizeOf(REAL);
  42.     ColSize  := Dim * SizeOf(Pointer);
  43.     RowSize  := Dim * RealSize;
  44.  
  45.     GetMem(Buf, ColSize);
  46.     FOR i := 1 TO GetDimension DO BEGIN
  47.       GetMem(Pointers(Buf^)[i], RowSize);
  48.       FillChar(Pointers(Buf^)[i]^, RowSize, 0);
  49.     END;
  50.   END;
  51.  
  52.   DESTRUCTOR RealMatrix.Done;
  53.   VAR
  54.     i : WORD;
  55.   BEGIN
  56.     FOR i := GetDimension DOWNTO 1 DO
  57.       FreeMem(Pointers(Buf^)[i], RowSize);
  58.     FreeMem(Buf, ColSize);
  59.   END;
  60.  
  61.   FUNCTION RealMatrix.GetDimension : WORD;
  62.   BEGIN
  63.    GetDimension := Dim;
  64.   END;
  65.  
  66.   FUNCTION RealMatrix.GetValue(i, j : WORD) : REAL;
  67.   VAR
  68.     Value : REAL;
  69.   BEGIN
  70.     Move(Reals(Pointers(Buf^)[i]^)[j], Value, RealSize);
  71.     GetValue := Value;
  72.   END;
  73.  
  74.   PROCEDURE RealMatrix.SetValue(i,j : WORD; Value : REAL);
  75.   BEGIN
  76.     Move(Value, Reals(Pointers(Buf^)[i]^)[j], RealSize);
  77.   END;
  78.  
  79.   FUNCTION RealMatrix.Determinant : REAL;
  80.   CONST
  81.     Tol = 1.0e-16;
  82.   VAR
  83.     i, j, k, N   : WORD;
  84.     Beta, Diag,
  85.     r, s, Sum,
  86.     Sigma, Det   : REAL;
  87.     Singular     : BOOLEAN;
  88.  
  89.   { Die Matrix wird mit HOUSEHOLDER - Transformationen }
  90.   { auf Dreiecksgestalt gebracht. Die Originalelemente }
  91.   { werden überschrieben. Die Determinante ist danach  }
  92.   { das Produkt der Diagonalelemente multipliziert mit }
  93.   { dem Faktor (-1)**(n-1). (n=Dimension der Matrix)   }
  94.  
  95.   BEGIN
  96.     N        := GetDimension;
  97.     Det      := 1.0;
  98.     Singular := FALSE;
  99.     i        := 0;
  100.     REPEAT
  101.       Inc(i);
  102.       Sum := 0.0;
  103.       FOR j := i TO N DO
  104.         Sum := Sum + Sqr(GetValue(i, j));
  105.       Sigma := Sqrt(Sum);
  106.       IF Sigma > Tol THEN BEGIN
  107.         IF GetValue(i, i) < 0 THEN
  108.           Sigma := -Sigma;
  109.         r := Sigma + GetValue(i, i);
  110.         Beta := 1.0 / (Sigma * r);
  111.         Diag := -Sigma;
  112.         SetValue(i, i, r);
  113.         FOR j := i+1 TO N DO BEGIN
  114.           Sum := 0.0;
  115.           FOR k := i TO N DO
  116.             Sum := Sum + GetValue(j, k) * GetValue(i, k);
  117.           Sum := Sum * Beta;
  118.           FOR k := i TO N DO
  119.             SetValue(j, k, GetValue(j, k) -
  120.                      Sum * GetValue(i, k));
  121.         END;
  122.         SetValue(i, i, Diag);
  123.         Det := Det * (-Diag);
  124.       END ELSE
  125.         Singular := TRUE;
  126.     UNTIL (i = N) OR Singular;
  127.     IF Singular THEN
  128.       Determinant := 0.0
  129.     ELSE
  130.       Determinant := Det;
  131.   END;
  132.  
  133. BEGIN
  134.   (* keine Initialisierung *)
  135. END.
  136. (* ------------------------------------------------------ *)
  137. (*               Ende von DETTOOL.PAS                     *)
  138.