INTERFACE SPLinearSystem;

(* Procedures for solving linear systems *)

IMPORT  SPMatrix, SPPWFunction, SPTriang, HLR3; 
FROM SPPWFunction IMPORT Basis;


TYPE
  LONG = LONGREAL;
  NAT = CARDINAL;

PROCEDURE CholeskiSolve(
    READONLY HL: SPMatrix.T;   (* Choleski factor of Helmholtz matrix. *)          
    READONLY b: ARRAY OF LONG; (* Right-hand-side of (2). *) 
    (* Work: *)
    VAR y: ARRAY OF LONG;      (* Intermediate vector for Choleski solution. *)
    (* Output: *)
    VAR v: ARRAY OF LONG;      (* Coordinates of new solution relative to "bas". *)
  );
  (*
    Solves the system "H u = b" by two-step substitution, given a
    Choleski decomposition of the matrix "H".  Upon entry, "HL" is the
    lower trinagular Choleski factor of a square matrix H.  Upon exit,
    "v" is the solution of the linear system "H v = b". *)
  
PROCEDURE GaussSeidelIteration(
    READONLY H: SPMatrix.T;    (* Helmholtz matrix. *)          
    READONLY b: ARRAY OF LONG; (* Right-hand-side of (2). *) 
    omega: LONG;               (* Overrelaxation factor. *)
    (* Input: *)
    READONLY u: ARRAY OF LONG; (* Coordinates of current sol relative to "bas". *)
    (* Out: *) 
    VAR v: ARRAY OF LONG;      (* Coordinates of new solution relative to "bas". *)
  );
  (*
    Performs one iteration of the Gauss-Seidel algorithm for solving
    the system "A u = b".  Upon entry, "u" should be the coordinates
    of a guess for the solution, relative to some basis.  Upon exit,
    "v" will have been updated by performing one iteration of the
    Gauss-Seidel algorithm. *)
  
PROCEDURE ConjugateGradientSolve(
    READONLY H: SPMatrix.T;    (* Helmholtz matrix. *)          
    READONLY b: ARRAY OF LONG; (* Right-hand-side of (2). *) 
    (* Input: *)
    READONLY u: ARRAY OF LONG; (* Coordinates of current sol relative to "bas". *)
    (* Out: *) 
    VAR v: ARRAY OF LONG;      (* Coordinates of new solution relative to "bas". *)
  );
  (*
    Attempts to solve the linear system "H u = b" by conjugate
    gradient minimization of the squared norm of the residual vector.
    Upon entry, "u" should be the coordinates of a guess for the
    solution, relative to some basis.  Upon exit, "v" will be a
    tentative solution for the system, computed by the Conjugate
    Gradient algorithm. *)
  
PROCEDURE GuessSol(VAR u: ARRAY OF LONG);
  (*
    Fills "u" with "random" values.
  *)

PROCEDURE BuildFunction(
    READONLY bas: Basis;         (* Basis for space "V[r]". *)
    READONLY u: ARRAY OF LONG;   (* Coordinates of "g" relative to  "bas" *)
    solName: TEXT;               (* Name prefix for output file. *)
    f: SPPWFunction.T;           (* True solution (for comparisons) *)
    (* Output: *)
    VAR g: SPPWFunction.T;       (* The function "g" *)
  );
  (*
    Stores in "g" the function with coefficients "u" relative 
    to the basis "bas". Also writes "g" to disk and compares it
    to the given function "f". *)

PROCEDURE UpdateSolution(
    READONLY v: ARRAY OF LONG; (* New solution vector. *)
    (* Input/output: *)
    VAR u: ARRAY OF LONG;      (* Current solution vector. *)
  ): LONG;
  (*
    This procedure is handy when using iterative methods like
    GaussSeidel for solving functional approximation or integration problems.
    It replaces the previous solution vector "u"
    by the newly computed solution vector "v", comparing the two.
    Returns the L2 norm of the difference. *)

PROCEDURE PlotFunctionAndError(
    outName: TEXT;
    g: SPPWFunction.T;
    f: SPPWFunction.T;
    tri: REF SPTriang.T;
    fMax: LONG;
    eMax: LONG;
    obs: HLR3.Point;
    gridSize: NAT;
  );
  (*
    Plots the function "g" comparing it to the function "f".
    Writes the plots to "NAME-val-{f,v}.eps" and "NAME-err-{f,v}.eps"
    where "f" means "front" and "v" means "back". *)

PROCEDURE PrintSomeValues(u: SPPWFunction.T; f: SPPWFunction.T);
  (*
    Prints values of "f(p)", "g(p)", and their difference 
    at a few points (mostly the octahedron vertices, edge
    centers, and face centers). *)

PROCEDURE PrintMaxErrorValues(
    g: SPPWFunction.T; 
    f: SPPWFunction.T;
    VAR fMax: LONG;     (* OUT: max abs function value. *)
    VAR eMax: LONG;     (* OUT: max abs error VALUE. *)
  );
  (*
    Computes and prints the maximum absolute value "fMax" of the
    function "g", and the maximum absolute error "eMax" when
    compared to "f". *)

END SPLinearSystem.