Programmer 7500/MAX_PROGRAMMERS.iso/PASCAL/ECO30603.ZIP/ECO30603.LZH/ECOLIBII/ECO_MTRX.PAS

Pascal/Delphi Source File | 1992-08-28 | 10.4 KB | 375 lines

(*
    ▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓·──                                              ──·▓▓▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓│                                                  │░░▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓   ECO_MTRX was NOT Conceived, Designed and         ░░▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓   Written by Floor A.C. Naaijkens for              ░░▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓   UltiHouse Software / The ECO Group.              ░░▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓                                                    ░░▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓   (C) MCMXCII by EUROCON PANATIONAL CORPORATION.   ░░▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓   All Rights Reserved for The ECO Group.           ░░▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓                                                    ░░▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓   Credit where credit is due.                      ░░▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓│                                                  │░░▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓·──                                              ──·░░▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓▓▓░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓
    ▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓
*)
{

Msg # 639
Date: 19 Aug 92 15:41:00
From: Greg Smith
  To: Kevin Lamonte
Subj: [1/3] Matrix Problem
____________________________________________________________________________

AREA:PASCAL
EID:ac30 19137d20
In a message to All <18 Aug 92  14:07> Kevin Lamonte wrote:

 KL> Help!  I'm trying to write a general matrices-math unit.  One of
the
 KL> functions is to calculate a determinant on an NxN matrice (limit is
set

Here's a copy of a Matrix unit that I wrote earlier this year.
Unfortunately, as you'll find out, it'd dead SLOW!  (Part due to the OOP
design, and part due to me not wanting to get into piles of math to
figure out the best algorithms to use).  It should work for matrices up
to 80x80, however I've never tested it on anything greater than 9x9 due
to speed problems.
}



(*----------------------------------------------------------*)
(* unit for doing everything you'll ever want to do with    *)
(* matrices... and more.                                    *)
(* requirements: turbo pascal 5.5 or better                 *)
(* no rights reserved, no copyright, released into the pd   *)
(* by the author: gregory p. smith of boulder, co. usa      *)
(* credit not required, but is greatly appreciated          *)
(*----------------------------------------------------------*)
(* last update: april  4th, 1992   -- gregory p. smith      *)
(*----------------------------------------------------------*)
unit eco_mtrx;
{$N+,E+,G+  -- support high precision numbers }


interface


type
{$IFDEF HIPREC}  Real = Extended; {$ENDIF}
  MXArray = Array[0..6399] of Real;              { allows 80x80 matrix }
  mxarrptr = ^mxarray;
  mxptr = ^matrix;
  matrix = object
    {------------------- variables -----------------------}
    m     : mxarrptr;                       { matrix data }
    cols,                             { columns in matrix }
    rows  : byte;                        { rows in matrix }
    {------------------- procedures ----------------------}
    constructor init(c, r: byte);     { create the matrix }
    procedure  setelement(c, r: byte; e: real);
    function   getelement(c, r: byte): real;
    procedure  multbyconst(con: real); { mult. by a const. }
    function   determinant: real;
    function   minor(c, r: byte): real;     { matrix minor }
    function   cofactor(c, r: byte): real;  { signed minor }
    function   adjoint: mxptr; { returns ptr to adj matrix }
    function   inverse: mxptr; { returns ptr to inv matrix }
    destructor done;                     { free up memory }
  end; { matrix }


const
  needmemory = -1; { need memory for temporary matrix. }
  notsquare  = -2; { square matrix required for this operation. }
  noinverse  = -3; { no inverse could be computed. }
  badsizes   = -4; { matrices were not comptiable for multiplication. }
  matherror  : shortint = 0; { holds one of the above error codes. }



  procedure writemx(var out: text; mx: mxptr);
  procedure showmatherror(err: integer);
  function  multmx(a, b: mxptr): mxptr;





implementation






(*----------------- matrix implementation --------------------------*)

  constructor matrix.init(c, r: byte);
  begin
    if c*r*sizeof(real) > sizeof(mxarray) then fail;
    getmem(m, c*r*sizeof(real));
  {$IFDEF HEAPMGR}
    if m = nil then fail;
  {$ENDIF}
    cols := c;
    rows := r;
  end;



  procedure matrix.setelement(c, r: byte; e: real);
  begin  { warning: no bounds checking to see if you're sane... }
    m^[(pred(r)*cols)+pred(c)] := e;
  end;



  function matrix.getelement(c, r: byte): real;
  begin  { warning: no bounds checking to see if you're sane... }
    getelement := m^[(pred(r)*cols)+pred(c)];
  end;



  procedure matrix.multbyconst(con: real);
  var
    t : byte;
  begin
    t := pred(rows * cols);
    for t := pred(rows * cols) downto 0 do
     m^[t] := m^[t] * con;
  end;



  function matrix.determinant: real;
  var
    s   : byte;
    tmp : real;

  begin
    if rows <> cols then begin
      matherror := notsquare;
      determinant := 0;
      exit;
    end;
    case rows of
      3..255 : begin { 3x3 and up }
        tmp := 0;
        for s := 1 to cols do tmp := tmp + (m^[pred(s)] * cofactor(s, 1));
        determinant := tmp;
      end; { 3x3 and up }
      2 : begin { 2x2 }  { ((1,1) * (2,2)) - ((1,2) * (2,1)) }
        determinant := (m^[0] * m^[3]) - (m^[2] * m^[1]);
      end; { 2x2 }
      1 : begin
        determinant := m^[0]; { (1,1) }
      end;
    end; { case rows ... }
    matherror := 0;
  end;





  function matrix.minor(c, r: byte): real;
  var
    tm : matrix;
    x, y,               { indexes to this matrix }
    tx, ty : byte;      { indexes to temp matrix }
  begin
    if (cols <> rows) or (cols < 2) { or (rows < 2) } then begin
      matherror := notsquare;
      minor := 0;
      exit;
    end;
    tm.init(pred(cols), pred(rows));
    { the minor is the determinant of the lower order matrix formed by }
    { excluding the row and column of the element we're computing for. }
    ty := 1;
    for y := 1 to rows do begin        { create lower order matrix }
      tx := 1;
      if y <> r then begin
        for x := 1 to cols do if x <> c then begin        { getelement(x, y) }
          tm.m^[(pred(ty)*tm.cols)+pred(tx)] := m^[(pred(y)*cols)+pred(x)];
          (* tm.setelement(tx, ty, m^[(pred(y)*cols)+pred(x)]); *)
          inc(tx);
        end;
        inc(ty);
      end;
    end;
    minor := tm.determinant;
    tm.done;
    matherror := 0;
  end;




  function matrix.cofactor(c, r: byte): real;
  begin
    { cofactor is the signed minor, depending upon whether the sum of }
    { the row and column is even or odd. }
    if odd(c+r) then cofactor := minor(c, r) * -1
     else cofactor := minor(c, r);
  end;



  function matrix.adjoint: mxptr;
  var
    tm   : mxptr;
    x, y : byte;    { index into both matrices }

  begin
    if (cols < 2) or (rows < 2) then begin
      matherror := notsquare;
      adjoint := nil;
      exit;
    end;
    new(tm, init(rows, cols));
    {$IFDEF HEAPMGR}
    if tm = nil then begin
      matherror := needmemory;
      adjoint := nil;
      exit;
    end;
    {$ENDIF}

    { adjoint matrix'es row elements correspond to the cofactors of the }
    { original matrices column elements.  not fun by hand... ;-}
    for y := 1 to rows do for x := 1 to cols do
      tm^.setelement(x, y, cofactor(y, x));
    adjoint := tm;
    matherror := 0;
  end;




  function matrix.inverse: mxptr;
  var    { warning: returns nil for undefined inverses }
    tm  : mxptr;
    det :  real;

  begin
    det := determinant;
    if det = 0 then begin
      inverse := nil;
      matherror := noinverse;
      exit;
    end;
    { inverse a is 1/deta * adja }
    tm := adjoint;
    if tm = nil then begin
      inverse := nil;
      exit;
    end;
    tm^.multbyconst(1 / det);
    inverse := tm;
  end;



  destructor matrix.done;
  begin
    freemem(m, cols*rows*sizeof(real));
  end;




(*------------- end matrix implementation --------------------------*)




  procedure showmatherror(err: integer);
  begin
    write('ERROR: ');
    case err of
    {$IFDEF HEAPMGR}
      needmemory : writeln('Out of memory for temporary matrix!');
    {$ENDIF}
      notsquare  : writeln(
        'Attempted to perform square operation on a rectangular matrix!'
      );
      noinverse  : writeln('Inverse undefined.');
      badsizes   : writeln(
        'Matrices not of compatible sizes to be multiplied!'
      );
    end;
    halt(abs(err));
  end;




  procedure writemx(var out: text; mx: mxptr);
  var x, y : byte;
  begin
    with mx^ do for x := 1 to rows do begin
      for y := 1 to cols do write(out,getelement(y, x):10:3);
      writeln(out);
    end;
  end;





  function multmx(a, b: mxptr): mxptr;
  var
    tm       : mxptr;
    x, y, tx :  byte;
    tmp      :  real;

  begin
    if (a^.cols <> b^.rows) then begin
      matherror := badsizes;
      multmx := nil;
      exit;
    end;
    new(tm, init(a^.rows, b^.cols));
    if tm = nil then begin
      matherror := needmemory;
      multmx := nil;
      exit;
    end;
    for x := 1 to b^.cols do for y := 1 to a^.rows do begin
      tmp := 0;
      for tx := 1 to a^.cols { b^.rows } do
      tmp := tmp + a^.getelement(y, tx) * b^.getelement(tx, x);
      tm^.setelement(y, x, tmp);
    end;
    multmx := tm;
  end;



end. { unit }


===  hope that works for you.
   << greg smith >>


---
 * origin: interconnect - littleton, co. [phone] (1:104/60.0)
seen-by: 2/777 13/13 28/777 29/777 105/27 42 138/112 234/50 270/18
seen-by: 280/0 283/1 500/0 1 2 3 4 5 7 8 9 13 222 512/0 1 4 8 23 32
seen-by: 512/37 44 45 47 112 133 146 148 149 167 168 172 1007
path: 518 512 1 13/13 138/112 105/27 42 500/1 9 512/0 167 172 172.1