INTERFACE SPMatrix;

IMPORT Rd, Wr;

(* Matrix*)

TYPE
  LONG = LONGREAL;
  Vector = ARRAY OF LONG;
  
  
  T = RECORD
    rows: CARDINAL;
    cols: CARDINAL;
    r: REF ARRAY OF Element;
  END;
  
  Element = RECORD  
     row: CARDINAL;
     col: CARDINAL;
     va: LONG;
  END;

PROCEDURE Add(READONLY A, B: T): T ;  
  (*
    Returns the matrix  "A+B". *)
  
PROCEDURE MulConst(READONLY A : T; c: LONG): T;
  (*
    Returns the product OF "c" by "A" *)
    
PROCEDURE Mul(READONLY A, B: T): T;
  (*
    Returns the matrix product "A*B". *)
  
PROCEDURE MulRow(
    READONLY x: Vector;
    READONLY A: T;
    VAR y: Vector; 
  );
  (*
    Stores in "y" the product of the row vector "x" by the matrix "A". *)
    
PROCEDURE MulCol(
    READONLY A: T;
    READONLY x: Vector;
    VAR y: Vector;
  );
  (*
    Stores in "y" the product of the matrix "A" by the column vector "x". *)
   
PROCEDURE DivRow(
    READONLY y: Vector;
    READONLY L: T;
    VAR x: Vector;
  ) ;
  (*
    Stores in "x" the product of the vector "y" by the inverse of the
    matrix "L".  The matrix must be lower-triangular, with nonzero
    diagonal elements. *)
    
PROCEDURE DivCol(
    READONLY L: T;
    READONLY y: Vector;
    VAR x: Vector;
  ) ;
  
  (*
    Stores in "x" the product of the inverse of matrix "L" by 
    the column vector "y".  The matrix "L" must be 
    lower-trinagular with nonzero diagonal. *)
    
PROCEDURE Choleski(READONLY A: T; VAR L: T);
 (* 
   Factors the positive definite nxn matrix "A" into 
   "L*L^t" where "L" is lower triangular and "L^t" is
   the transpose of "L". *)

PROCEDURE GaussSeidel(READONLY A: T; READONLY b: Vector; omega: LONG; VAR x: Vector);
  (*
    Performs one iteration of the Gauss-Seidel algorithm for solving
    the linear system "A x = b". That is,
    recomputes "x[i] = x[i] + omega*(b[i] - SUM{a[i,j]*x[j]})/a[i,i]"
    for each "i", increasing. Note that "x[j]" in the summation
    is partly old and partly new.  To understand the role of "omega",
    read the textbooks. The matrix "A" must be square.
  *)

PROCEDURE ConjugateGradient(
   READONLY A: T; 
   READONLY b: Vector; 
   VAR x: Vector; 
  ): LONG;
  (* 
    Solves the linear system "A x = b" using the Conjugate Gradient algorithm 
    [Numerical Recipes section 2.10, subroutine SPARSE].
    Returns the squared norm of the residual "A x - b". A large
    value means that "A" is singular and therefore "x" is only a 
    least squares solution.
  *)  
  
PROCEDURE Find(READONLY r: ARRAY OF Element; i, j: CARDINAL): CARDINAL;
  (*
    Finds an index "p" such that "r[p].row = i" and "r[p].col = j".
    Assumes the elements in "r" are stored in row-major order.
    Returns "LAST(CARDINAL)" if not found. *)

PROCEDURE Write(wr: Wr.T; READONLY A: T);
  (* 
    Write  the  matrix "A" *)
  
PROCEDURE Read(rd: Rd.T; VAR A: T);
  (* 
    Reads the matrix "A" from "rd<
    assuming it was written by "Write". *)
  
END SPMatrix.