home *** CD-ROM | disk | FTP | other *** search
/ Amiga MA Magazine 1998 #6 / amigamamagazinepolishissue1998.iso / coders / matrix / matrix.mod < prev    next >
Text File  |  1995-04-21  |  9KB  |  437 lines
  1. (*--------------------------------------------------------------------
  2.  
  3. :Program.       MATRIX.mod
  4. :Author.        Michael Meyer
  5. :Address.       Würmseestrasse 28, 81476 München
  6. :EMail.         ZNet: M.MEYER@AMC.ZER.SUB.ORG
  7. :EMail.         Fido: Michael Meyer 2:2480/21.100
  8. :Version.       1.0
  9. :Date.          6 Feb 1994
  10. :Copyright.     1994 Michael Meyer
  11. :Contents.      Dieses Modul unterstüzt das Rechnen mit Matrizen
  12. :Language.      Oberon-2
  13. :Translator.    Amiga Oberon 3.00d
  14.  
  15. --------------------------------------------------------------------*)
  16.  
  17. (* Matrix = elements( [0,  0] [0,  1] [0,  2] ... [0,  n-1]
  18.  *                    [1,  0] [1,  1] [1,  2] ... [1,  n-1]
  19.  *                    [2,  0] [2,  1] [2,  2] ... [2,  n-1]
  20.  *                       .       .       .           .
  21.  *                       .       .       .           .
  22.  *                       .       .       .           .
  23.  *                    [m-1,0] [m-1,1] [m-1,2] ... [m-1,n-1] )
  24.  *)
  25.  
  26. MODULE MATRIX;
  27.  
  28. IMPORT BT * := BasicTypes;
  29.  
  30. TYPE
  31.   MATRIX * = POINTER TO MATRIXDesc;
  32.   MATRIXDesc * = RECORD (BT.GROUPDesc)
  33.     m, n: LONGINT;
  34.     elements: POINTER TO ARRAY OF ARRAY OF LONGREAL;
  35.   END;
  36.  
  37. PROCEDURE Create * (m, n: LONGINT): MATRIX;
  38. (*
  39.  * require
  40.  *   m>0 & n>0
  41.  *)
  42. VAR
  43.   new: MATRIX;
  44.   i, j: LONGINT;
  45. BEGIN
  46.   NEW(new);
  47.   NEW(new.elements, m, n);
  48.   new.m := m;
  49.   new.n := n;
  50.   FOR i := 0 TO m-1 DO
  51.     FOR j := 0 TO n-1 DO
  52.       new.elements[j, i] := 0.0;
  53.     END;
  54.   END;
  55.   RETURN new;
  56. END Create;
  57.  
  58. PROCEDURE NeutralMatrix * (m: LONGINT): MATRIX;
  59. (*
  60.  * require
  61.  *   m>0 & n>0
  62.  *   a.DimM() = a.DimN()
  63.  *)
  64. VAR
  65.   new: MATRIX;
  66.   i, j: LONGINT;
  67. BEGIN
  68.   NEW(new);
  69.   NEW(new.elements, m, m);
  70.   new.m := m;
  71.   new.n := m;
  72.   FOR i := 0 TO m-1 DO
  73.     FOR j := 0 TO m-1 DO
  74.       new.elements[j, i] := 0.0;
  75.     END;
  76.   END;
  77.   FOR i := 0 TO m-1 DO
  78.     new.elements[i, i] := 1.0;
  79.   END;
  80.   RETURN new;
  81. END NeutralMatrix;
  82.  
  83. PROCEDURE (a: MATRIX) Get * (m, n: LONGINT): LONGREAL;
  84. (*
  85.  * require
  86.  *   a.DimM()>m &
  87.  *   a.DimN()>n;
  88.  *   m>0 &
  89.  *   n>0
  90.  *)
  91. BEGIN
  92.   RETURN a.elements[m, n];
  93. END Get;
  94.  
  95. PROCEDURE (a: MATRIX) Put * (m, n: LONGINT; x: LONGREAL);
  96. (*
  97.  * require
  98.  *   a.DimM()>m &
  99.  *   a.DimN()>n;
  100.  *   m>0 &
  101.  *   n>0
  102.  *)
  103. BEGIN
  104.   a.elements[m, n] := x;
  105. END Put;
  106.  
  107. PROCEDURE (a: MATRIX) DimM * (): LONGINT;
  108. BEGIN
  109.   RETURN a.m;
  110. END DimM;
  111.  
  112. PROCEDURE (a: MATRIX) DimN * (): LONGINT;
  113. BEGIN
  114.   RETURN a.n;
  115. END DimN;
  116.  
  117. PROCEDURE (a: MATRIX) Add * (b: BT.GROUP): BT.GROUP;
  118. (*
  119.  * require
  120.  *   a.DimM() = b(MATRIX).DimM()
  121.  *   a.DimN() = b(MATRIX).DimN()
  122.  *)
  123. VAR
  124.   new: MATRIX;
  125.   i, j: LONGINT;
  126. BEGIN
  127.   WITH b: MATRIX DO
  128.     new := Create(a.m, a.n);
  129.     FOR i := 0 TO a.m-1 DO
  130.       FOR j := 0 TO a.n-1 DO
  131.         new.elements[i, j] := a.elements[i, j] + b.elements[i, j]
  132.       END;
  133.     END;
  134.   END;
  135.   RETURN new;
  136. END Add;
  137.  
  138. PROCEDURE (a: MATRIX) Sub * (b: BT.GROUP): BT.GROUP;
  139. (*
  140.  * require
  141.  *   a.DimM() = b(MATRIX).DimM()
  142.  *   a.DimN() = b(MATRIX).DimN()
  143.  *)
  144. VAR
  145.   new: MATRIX;
  146.   i, j: LONGINT;
  147. BEGIN
  148.   WITH b: MATRIX DO
  149.     new := Create(a.m, a.n);
  150.     FOR i := 0 TO a.m-1 DO
  151.       FOR j := 0 TO a.n-1 DO
  152.         new.elements[i, j] := a.elements[i, j] - b.elements[i, j]
  153.       END;
  154.     END;
  155.   END;
  156.   RETURN new;
  157. END Sub;
  158.  
  159. PROCEDURE (a: MATRIX) Mul * (b: BT.GROUP): BT.GROUP;
  160. (*
  161.  * c = a * b
  162.  *
  163.  * require
  164.  * a.DimN() = b(MATRIX).DimM()
  165.  *
  166.  * => c.DimM() = a.DimM()
  167.  *    c.DimN() = b(MATRIX).DimN()
  168.  *)
  169. VAR
  170.   new: MATRIX;
  171.   i, j, o: LONGINT;
  172.   c: LONGREAL;
  173. BEGIN
  174.   WITH b: MATRIX DO
  175.     new := Create(a.m, b.n);
  176.     FOR i := 0 TO a.m-1 DO
  177.       FOR j := 0 TO b.n-1 DO
  178.         c := 0.0;
  179.         FOR o := 0 TO a.n-1 DO
  180.           c := c + a.elements[i, o] * b.elements[o, j];
  181.         END;
  182.         new.elements[i, j] := c;
  183.       END;
  184.     END;
  185.   END;
  186.   RETURN new;
  187. END Mul;
  188.  
  189. PROCEDURE (a: MATRIX) SkalarMul * (r: LONGREAL): BT.GROUP;
  190. (*
  191.  * multiply each element of a by r.
  192.  *)
  193. VAR
  194.   new: MATRIX;
  195.   i, j: LONGINT;
  196. BEGIN
  197.   WITH a: MATRIX DO
  198.     new := Create(a.m, a.n);
  199.     FOR i := 0 TO a.m-1 DO
  200.       FOR j := 0 TO a.n-1 DO
  201.         new.elements[i, j] := r * a.elements[i, j];
  202.       END;
  203.     END;
  204.   END;
  205.   RETURN new;
  206. END SkalarMul;
  207.  
  208. PROCEDURE (a: MATRIX) isEqual * (b: BT.COMPAREABLE): BOOLEAN;
  209. VAR
  210.   i, j: LONGINT;
  211. BEGIN
  212.   WITH b: MATRIX DO
  213.     IF ~((a.m = b.m) & (a.n = b.n)) THEN RETURN FALSE END;
  214.     FOR i := 0 TO a.m-1 DO
  215.       FOR j := 0 TO a.n-1 DO
  216.         IF a.elements[i, j] # b.elements[i, j] THEN RETURN FALSE END;
  217.       END;
  218.     END;
  219.   END;
  220.   RETURN TRUE;
  221. END isEqual;
  222.  
  223. PROCEDURE (a: MATRIX) isDimEqual * (b: BT.COMPAREABLE): BOOLEAN;
  224. BEGIN
  225.   WITH b: MATRIX DO
  226.     IF ~((a.m = b.m) & (a.n = b.n)) THEN
  227.       RETURN FALSE;
  228.     ELSE
  229.       RETURN TRUE;
  230.     END;
  231.   END;
  232. END isDimEqual;
  233.  
  234. PROCEDURE (a: MATRIX) TransposeMatrix * (): BT.GROUP;
  235. VAR
  236.   i, j: LONGINT;
  237.   new: MATRIX;
  238. BEGIN
  239.   new := Create(a.n, a.m);
  240.   FOR i := 0 TO a.m-1 DO
  241.     FOR j := 0 TO a.n-1 DO
  242.       new.elements[j, i] := a.elements[i, j];
  243.     END;
  244.   END;
  245.   RETURN new;
  246. END TransposeMatrix;
  247.  
  248. PROCEDURE FindPivot (a: MATRIX; c, l : LONGINT): LONGINT;
  249. (*
  250.  * find Pivot-element
  251.  * in column c
  252.  * from line l on
  253.  *
  254.  * returns line of Pivot-element
  255.  *)
  256. VAR
  257.   i : LONGINT;
  258.   pivot, next : LONGREAL;
  259.   line : LONGINT;
  260. BEGIN
  261.   pivot := ABS(a.elements[l, c]);
  262.   line := l;
  263.   FOR i := l+1 TO a.m-1 DO
  264.     next := ABS(a.elements[i, c]);
  265.     IF next > pivot THEN
  266.       pivot := next;
  267.       line := i;
  268.     END;
  269.   END;
  270. RETURN line;
  271. END FindPivot;
  272.  
  273. PROCEDURE ChangeTwoLines (a: MATRIX; l1, l2, c: LONGINT);
  274. (*
  275.  * change two lines
  276.  * from column c to column a.DimN()-1
  277.  *)
  278. VAR
  279.   i : LONGINT;
  280.   x : LONGREAL;
  281. BEGIN
  282.   FOR i := c TO a.n-1 DO
  283.     x := a.elements[l1, i];
  284.     a.elements[l1, i] := a.elements[l2, i];
  285.     a.elements[l2, i] := x;
  286.   END;
  287. END ChangeTwoLines;
  288.  
  289. PROCEDURE MulLineAddLine(a: MATRIX; l1, l2, c: LONGINT; s: LONGREAL);
  290. (*
  291.  * multiplicate each element of line l1 from column c with
  292.  * skalar s and add the result to line l2
  293.  *)
  294. VAR
  295.   i : LONGINT;
  296.   x : LONGREAL;
  297. BEGIN
  298.   FOR i := c TO a.n-1 DO
  299.     a.elements[l2, i] := a.elements[l2, i] + a.elements[l1, i] * s;
  300.   END;
  301. END MulLineAddLine;
  302.  
  303.  
  304. PROCEDURE (a: MATRIX) Neg * (): BT.GROUP;
  305. (*
  306.  * require
  307.  *   a.DimM() = a.DimN()
  308.  *   a.Det() # 0
  309.  *)
  310. VAR
  311.   new, newN, return: MATRIX;
  312.   i, j, k: LONGINT;
  313.   ln: LONGINT;
  314.   t1, t2, skalar, pivot, det: LONGREAL;
  315.   sign: BOOLEAN;
  316. BEGIN
  317.   new := a;
  318.   newN := NeutralMatrix(a.m);
  319.   return := Create(a.m, a.n);
  320. (*------------------ Gauss ------------------*)
  321.   sign := TRUE;
  322.   FOR i := 0 TO new.m-2 DO
  323.     ln := FindPivot(new, i, i);
  324.     pivot := -1.0 * new.elements[ln, i];
  325.     IF ~(ln = i) THEN
  326.       sign := ~sign;
  327.       ChangeTwoLines(new, i, ln, i);
  328.       ChangeTwoLines(newN, i, ln, 0);
  329.     END;
  330.     FOR j := i+1 TO new.m-1 DO
  331.       skalar := new.elements[j, i]/pivot;
  332.       MulLineAddLine(new, i, j, i, skalar);
  333.       MulLineAddLine(newN, i, j, 0, skalar);
  334.     END;
  335.   END;
  336.   IF ~sign THEN
  337.     FOR i := 0 TO new.n-1 DO
  338.       new.elements[new.m-1, i] := -1.0 * new.elements[new.m-1, i];
  339.       newN.elements[new.m-1, i] := -1.0 * newN.elements[new.m-1, i];
  340.     END;
  341.  
  342.   END;
  343. (*------------------ Jordan ------------------*)
  344.   FOR i := new.m-1 TO 1 BY -1 DO
  345.     t1 := -1.0 * new.elements[i, i];
  346.     FOR j := i-1 TO 0 BY -1 DO
  347.       t2 := new.elements[j, i] / t1;
  348.       MulLineAddLine(new, i, j, j, t2);
  349.       MulLineAddLine(newN, i, j, 0, t2);
  350.     END;
  351.   END;
  352.   FOR i := 0 TO new.m-1 DO
  353.     FOR j := 0 TO new.n-1 DO
  354.       return.elements[i, j] := newN.elements[i, j]/new.elements[i, i];
  355.     END;
  356.   END;
  357.   RETURN return;
  358. END Neg;
  359.  
  360. PROCEDURE (a: MATRIX) Gauss * (): BT.GROUP;
  361. (*
  362.  * require
  363.  *   a.DimM() = a.DimN()
  364.  *)
  365. VAR
  366.   new: MATRIX;
  367.   sign: BOOLEAN;  (* sign of the matrix: TRUE=plus and FALSE=minus *)
  368.   i, j: LONGINT;
  369.   ln: LONGINT;  (* line-number in which Pivot-element is *)
  370.   skalar, pivot: LONGREAL;
  371. BEGIN
  372.   new := Create(a.m, a.n);
  373.   new := a;
  374.   sign := TRUE;
  375.   FOR i := 0 TO new.m-2 DO
  376.     ln := FindPivot(new, i, i);
  377.     pivot := -1.0 * new.elements[ln, i];
  378.     IF ~(ln = i) THEN
  379.       sign := ~sign;
  380.       ChangeTwoLines(new, i, ln, i);
  381.     END;
  382.     FOR j := i+1 TO new.m-1 DO
  383.       skalar := new.elements[j, i]/pivot;
  384.       MulLineAddLine(new, i, j, i, skalar);
  385.     END;
  386.   END;
  387.   IF ~sign THEN
  388.     FOR i := 0 TO new.n-1 DO
  389.       new.elements[new.m-1, i] := -1.0 * new.elements[new.m-1, i];
  390.     END;
  391.   END;
  392.   RETURN new;
  393. END Gauss;
  394.  
  395. PROCEDURE (a: MATRIX) Det * (): LONGREAL;
  396. (*
  397.  * require
  398.  *   a.DimM() = a.DimN()
  399.  *)
  400. VAR
  401.   new: MATRIX;
  402.   i: LONGINT;
  403.   result: LONGREAL;
  404. BEGIN
  405.   new := a.Gauss()(MATRIX);
  406.   result := new.elements[0, 0];
  407.   FOR i := 1 TO new.m-1 DO
  408.     result := result * new.elements[i, i];
  409.   END;
  410.   RETURN result;
  411. END Det;
  412.  
  413. PROCEDURE (a: MATRIX) Norm * (): LONGREAL;
  414. BEGIN
  415.   RETURN a.Det();
  416. END Norm;
  417.  
  418. PROCEDURE (a: MATRIX) Compare * (b: BT.COMPAREABLE): LONGINT;
  419. (*
  420.  * require
  421.  *   a.DimM() = a.DimN()
  422.  *   b(MATRIX).DimM() = b(MATRIX).DimN()
  423.  *)
  424. VAR
  425.   i, j: LONGINT;
  426.   an, bn: LONGREAL;
  427. BEGIN
  428.   WITH b: MATRIX DO
  429.     an := a.Norm(); bn := b.Norm();
  430.     IF    an = bn THEN RETURN  0
  431.     ELSIF an > bn THEN RETURN  1
  432.     ELSE               RETURN -1 END;
  433.   END;
  434. END Compare;
  435.  
  436. END MATRIX.
  437.