Task

Solve   Ax=b   using Gaussian elimination then backwards substitution.

A   being an   n by n   matrix.

Also,   x and b   are   n by 1   vectors.

To improve accuracy, please use partial pivoting and scaling.



See also

   the Wikipedia entry:   Gaussian elimination


GENERIC INTERFACE Matrix(RingElem);
 
(*
"RingElem" must export the following:
- a type T;
- procedures
  + "Nonzero(a: T): BOOLEAN", which indicates whether "a" is nonzero;
  + "Minus(a, b: T): T" and "Times(a, b: T): T",
    which return the results you'd guess from the procedures' names; and
  + "Print(a: T)", which does what the name implies.
*)
 
TYPE
 
  T <: Public;
 
  Public = OBJECT
  METHODS
    init(READONLY data: ARRAY OF ARRAY OF RingElem.T): T;
    (* use this to copy the entries in "data"; returns "self" *)
    initDimensions(m, n: CARDINAL): T;
    (* use this for an mxn matrix of random entries *)
    num_rows(): CARDINAL;
    (* returns the number of rows in "self" *)
    num_cols(): CARDINAL;
    (* returns the number of columns in "self" *)
    entries(): REF ARRAY OF ARRAY OF RingElem.T;
    (* returns the entries in "self" *)
    triangularize();
    (*
      Performs Gaussian elimination in the context of a ring.
      We can add scalar multiples of rows,
      and we can swap rows, but we may lack multiplicative inverses,
      so we cannot necessarily obtain 1 as a row's first entry.
    *)
  END;
 
  PROCEDURE PrintMatrix(m: T);
  (* prints the matrix row-by-row; sorry, no special padding to line up columns *)
 
END Matrix.