MODULE SPMultiSolve EXPORTS Main;

(* Solves the Helmholtz differential equation by multiscale iteration. *)

(*
  This program solves the Helmholz differential equation on the sphere
  
    SLapl(u)(x) + c*u(x) = F(x, u(x))       (1)
    
  where "SLapl" is the spherical Laplacian operator, "c" is a given
  constant, "F" is a given function, and "u" is the function to be
  determined. 
  
  The program uses a previously built hierarchy of meshes (spherical
  triangulations) "tri[minScale..maxScale]", at various scales of resolution
  ("minScale" is the finest, "maxScale" is the coarsest), and associated
  function approximation spaces $V[minScale..maxScale]", with 
  bases "bas[minScale..maxScale]".  
  
  The linear matrices needed to transfer the solutions between meshes
  and to solve the system at each mesh must have been computed and
  factired previously; see "SPGGMaCho" and "SPFGMatrix".  
  
  It is assumed that, at each level of the hierarchy, an approximate
  solution "u" to (1) in the corresponding function space "V" can be
  found by solving the non-linear system
  
    H a = b         (2)
    
  where 
  
    "a[i]" is coefficient "i" of the solution 
           "u" (to be determined) in the basis "bas"; 
  
    "b[i] = <F | bas[i]>, where "F" is the spherical function
            "x -> F(x,u(x))" and "<|>" is the scalar product.
           
    "H[i,j] = <SGrad(bas[i]) | SGrad(bas[j]> + c*<bas[i]|bas[j]>
            where "SGrad" is the spherical gradient operator.
        
  Note that this system is not trivial since the right-hand side "b"
  depends on the function "u", and therefore on the coefficient vector
  "a", usually in a non-linear way.  So (2) must be solved
  iteratively.
  
  The multiscale method accelerates the resolution of (2) by
  performing the iterations simultaneously on all meshes.  After
  performing a few iterations at some scale "r", we either raise the
  current solution "u[r]" from space "V[r]" to the coarser space
  "V[r+1]", or lower it onto the finer space "V[r-1]".
  
  The raising and lowering operations are defined in terms of the
  spaces "W[r]" and "Z[r]", where "W[r]" is the orthogonal projection
  of the coarser space "V[s] = V[r+1]" onto "V[r]", and "Z[r]" is
  the orthogonal complement of "W[r]" in "V[r]".  It is assumed that
  "W[r]" and "V[s]" have the same dimension, and that, moreover,
  
     inf { |v'|/|v| : v in V[s] } >> 0     (3)
     
  where "v'" is the orthogonal projection "v" onto "V[r]".
  
  The raising operation takes a vector "u[r]" in "V[r]" and produces a
  closest approximation "u[s]" in "V[s]", which may be simply the
  orthogonal projection of "u[r]" onto "V[s]".  Note that this
  projection can be viewed as a projection onto "W[r]" followed by a
  dimension-preserving projection from "W[r]" to "V[s]".
  
  The lowering operation takes two vectors "u[s]" in "V[s]" and
  "v[r]" in "V[r]", and produces a new vector "u[r]" in "V[r]" that
  combines all the information that can be taken from "u[s]", whith
  watever information is present in "v[r]" but cannot be represented
  in "V[s]".  That is, we decompose "v[r]" into a "fine-scale"
  component "z[r]" in "Z[r]" and a "coarse-scale" component in "W[r]";
  we discard the latter, and define "u[r]" as "z[r]" plus the best
  possible approximation "w[r]" of "u[s]" in "W[r]".
  
  For the raising and lowering operations we need three matrices:
  "L[r]" is a bijection from "V[s]" to "W[r]", "R[r]" is a
  projection from "V[r]" onto "V[s]", and "T[r]" is the orthogonal
  projection of "V[r]" onto "W[r]".  Then raising a vector "u[r]" from
  "V[r]" to "V[s]" means computing "u[s] = R[r]u[r]", and lowering a
  vector "u[s]" from "V[s]" to "V[r]" means computing "u[r] = (I - T[r])v[r]
  + L[r]u[s]", where "v[r]" is an "a priori" solution in "V[r]".
  
  It turns out that (check!) 
  
    L[r] = G[r]^{-1} F[r]
    R[r] = G[s]^{-1} F[r]^*
    
  where "A^*" is the transpose of "A", and
   
    G[r][i,j] = <bas[r][i] | bas[r][j]>
    F[r][i,j] = <bas[s][i] | bas[r][j]>
    
  Note that "G[r]" is square and non-singular, with "d[r]" rows and
  columns.  Note also that "F[r]" is rectangular, with "d[s]" rows
  and "d[r]" columns; and has full rank because of condition (3) above.
  
*)

IMPORT Math, Wr, Thread, ParseParams, Fmt, Process;
IMPORT SPTriang, SPMatrix, SPNLFunction, SPPWFunction, LR3, LRN, SPRange; 
FROM SPFile IMPORT OpenRead, OpenWrite;
IMPORT SPAnalyticPWFunction;
FROM SPPWFunction IMPORT Id, Point;

FROM Stdio IMPORT stderr;

TYPE 
  Options = RECORD
      triName: TEXT;       (* Template for triangulation file names. *)
      baseName: TEXT;      (* Template for basis file names. *)
      rhsName: TEXT;       (* Name of right-hand-side function. *)
      solName: TEXT;       (* Name of correct solution function. *)
      matName: TEXT;       (* Template for names of Choleski matrix files. *)
      minScale: CARDINAL;  (* Minimum resolution scale *)
      maxScale: CARDINAL;  (* Maximum resolution scale *)
      outName: TEXT;       (* Template for output file names. *)
      maxIter: CARDINAL;   (* Maximum iterations *)
      absTol: REF LONGS;   (* Max error compared to true sol, per scale. *)
      relTol: REF LONGS;   (* Max change in "u" coeffs between iterations, per scale. *)
      order: CARDINAL;     (* Will use "4^order" sub-triangles per triangle. *)
    END;

TYPE 
  LONG = LONGREAL;
  LONGS = ARRAY OF LONG;
  NAT = CARDINAL;
  Basis = SPPWFunction.Basis;
  
  
TYPE
  ScaleData = RECORD
      (* Input data: *)
      tri: REF SPTriang.T;      (* The triangular mesh for some scale. *)
      bas: REF Basis;           (* A basis for the spline space "V" on "tri". *)
      GL: SPMatrix.T;           (* Choleski factor of scalar product matrix for "bas". *)
      HL: SPMatrix.T;           (* Choleski factor of Helmoltz matrix for "bas". *)
      F: SPMatrix.T;            (* Matrix connecting this "bas" with coarser "bas". *)
      ufCorr: SPPWFunction.T;   (* True solution (for comparisons) *)
      (*N: SPMatrix.T;*)        (* Projection matrix of "V" onto "W" subspace. *)
      (* Working vectors: *)  
      x: REF ARRAY OF LONG;     
      y: REF ARRAY OF LONG;     
      z: REF ARRAY OF LONG;     
      (* Current solution: *)   
      uf: SPPWFunction.T;       (* Current solution as a spherical function. *)
      u: REF ARRAY OF LONG;     (* Coordinates of "uf" relative to "bas". *)
    END;
    (*
      The matrices "GL" and "HL" are stored in Choleski factored form.
      *)
  
  MultiScaleData = ARRAY OF ScaleData;
  
<* FATAL Wr.Failure, Thread.Alerted *>

PROCEDURE Main() =
  VAR d: LONG;
      iter: NAT;
      fMax, eMax: LONG;
  BEGIN
    WITH 
      o = GetOptions(),
      m = ReadMultiScaleData(
          o.minScale, o.maxScale, o.triName, o.baseName, o.matName, o.solName
        )^,
      f = SPNLFunction.FromName(o.rhsName) (* Right-hand side *)
      
    DO 
      FOR r := o.maxScale TO o.minScale BY -1 DO
        WITH
          ur = m[r].u^,
          urf = m[r].uf,
          basr = m[r].bas^,
          ufCorr = m[r].ufCorr,
          dim = NUMBER(basr),
          errWr = OpenWrite(ErrorFileName(o.outName, r) & ".erp")
        DO
          WITH solFileName = SolutionFileName(o.outName, r, 0) DO
            IF r < o.maxScale THEN
              FOR i:= 0 TO dim-1 DO ur[i] := 0.0d0 END;
              urf := SPPWFunction.LinComb(ur, basr);
              LowerSol(
                m[r+1].u^, m[r].bas^, m[r].GL, m[r].F, (* m[r].N,*) solFileName, 
                (*WORK*) (*m[r].x^,*) m[r].y^, m[r].z^,
                (*OUT*)  urf, ur
              )
            ELSE
              GuessSol(m[r].bas^, solFileName, (*OUT*) urf, ur)
            END
          END;
          iter := 0;
          WHILE iter < o.maxIter DO
            INC(iter);
            WITH solFileName = SolutionFileName(o.outName, r, iter) DO
              d := Iteration(
                m[r].tri, m[r].bas^, m[r].HL, f, o.order, 
                solFileName,
                (*I/O*) urf, ur, 
                (*WORK*) m[r].x^, m[r].y^, m[r].z^
              )
            END;
            PrintMaxErrorValues(urf, ufCorr, fMax, eMax);
            Wr.PutText(errWr, FI(iter, 6) & " " & FLR(eMax, 12, 8) & "\n");
            Wr.Flush(errWr);
            IF (d <= o.relTol[r - o.minScale]) 
            OR (eMax <= o.absTol[r - o.minScale]) THEN EXIT END;
          END;
          Wr.Close(errWr);
        END;
      END
    END;  
  END Main;
  
PROCEDURE Iteration(
    tri: REF SPTriang.T;     (* Triangulation to use for integration. *)
    READONLY bas: Basis;     (* Basis for the approximation space. *)
    READONLY HL: SPMatrix.T; (* Choleski factor of Helmhotz matrix. *)          
    F: SPNLFunction.T;       (* Right-hand side of differential eqn. *)
    sOrder: NAT;             (* Order of sampling grid in each triangle *)
    solFileName: TEXT;       (* Name template for output files *)
    (* Input/output: *)
    VAR uf: SPPWFunction.T;  (* Current solution. *)
    VAR u: ARRAY OF LONG;    (* Coordinates of "uf" relative to "bas". *)
    (* Work: *)
    VAR b: ARRAY OF LONG;    (* Right-hand-side of (2). *) 
    VAR y: ARRAY OF LONG;    (* Intermediate vector for Choleski solution. *)
    VAR v: ARRAY OF LONG;    (* Coordinates of new solution relative to "bas". *)
  ): LONG =
  (*
    Performs one iteration of the uniscale integration algorithm.
    Upon entry, "uf" should be a guess for the solution, and "u"
    its coordinates relative to "bas".  Upon exit, "uf" and "u"
    will have been updated by solving the linear system (2), 
    where the right-hand side "b' is computed from the given "uf". 
    Returns the distance between the old and new solutions. *)
    VAR sb2: LONG;
  BEGIN
    WITH
      dim = NUMBER(bas)
    DO
      Wr.PutText(stderr, "computing right-hand side \"b\"\n");
      sb2 := 0.0d0;
      FOR i := 0 TO dim-1 DO
        Wr.PutText(stderr, ".");
        WITH 
          bNew = SPPWFunction.DotProduct(uf, bas[i], sOrder, F, Id, tri),
          db = b[i] - bNew
        DO
          sb2 := sb2 + db*db;
          b[i] := bNew;
        END;
        Wr.PutText(stderr, "  b[" & Fmt.Pad(Fmt.Int(i),3) & "] = ");
        Wr.PutText(stderr, 
          Fmt.Pad(Fmt.LongReal(b[i], style := Fmt.Style.Fix, prec := 6), 9)
          & "\n"
        );
      END;
      Wr.PutText(stderr, "\n");
      Wr.PutText(stderr, " change in b =" & 
         Fmt.Pad(Fmt.LongReal(Math.sqrt(sb2), style := Fmt.Style.Fix, prec := 6), 9)
          & "\n"
      );
      Wr.PutText(stderr, "solving system\n");
      SPMatrix.DivCol(HL, b, y);
      SPMatrix.DivRow(y, HL, v);
      WITH d = LRN.Dist(u, v) DO
        Wr.PutText(stderr, "change = " & Fmt.LongReal(d) & "\n");
        u := v;
        uf := SPPWFunction.LinComb(u, bas);
        Wr.PutText(stderr, "writing " & solFileName & "\n");
        WITH wr = OpenWrite(solFileName & "-i.sfn") DO 
          uf.write(wr); Wr.Close(wr)
        END;
        PrintSomeValues(uf);
        RETURN d
      END;
    END;
  END Iteration;
  
PROCEDURE GuessSol(
    READONLY basr: Basis;        (* Basis for space "V[r]". *)
    solFileName: TEXT;           (* Name template for output file. *)
    (* Input/Output: *)
    VAR urf: SPPWFunction.T;     (* Approximation to "usf" in "V[r]" *)
    VAR ur: ARRAY OF LONG;       (* Coordinates of "urf" relative to  "basr" *)
  ) =
  (*
    Stores in "urf" some arbitrary function in "V[r]". *)
  BEGIN
    FOR i:= 0 TO LAST(ur) DO ur[i] := 0.0d0 END;
    Wr.PutText(stderr, "guessing initial solution\n");
    urf := SPPWFunction.LinComb(ur, basr);
    Wr.PutText(stderr, "writing " & solFileName & "\n");
    WITH wr = OpenWrite(solFileName & "-g.sfn") DO
      urf.write(wr); Wr.Close(wr)
    END;
    PrintSomeValues(urf);
  END GuessSol;

PROCEDURE PrintMaxErrorValues(
    u: SPPWFunction.T; 
    ufCorr: SPPWFunction.T;
    VAR fMax: LONG;     (* OUT: max abs function value. *)
    VAR eMax: LONG;     (* OUT: max abs error value. *)
  ) =

  PROCEDURE UValue(READONLY p: Point): LONGREAL =
    BEGIN
      RETURN u.eval(p)
    END UValue;

  PROCEDURE UError(READONLY p: Point): LONGREAL =
    BEGIN
      RETURN u.eval(p) - ufCorr.eval(p)
    END UError;

  BEGIN
    fMax := SPRange.OnSphere(UValue, 40);
    Wr.PutText(stderr, "MAX(ABS(fCorr))   = " & FLR(fMax, 9,6) & "\n");
    eMax := SPRange.OnSphere(UError, 40);
    Wr.PutText(stderr, "MAX(ABS(fAppr-fCorr)) = " & FLR(eMax, 9,6) & "\n");
  END PrintMaxErrorValues;

PROCEDURE MakeCorrectSolution(
    name: TEXT; 
    tri: REF SPTriang.T;
  ): SPPWFunction.T =
  BEGIN
    WITH f = SPAnalyticPWFunction.FromName(name, tri) DO
      IF f = NIL THEN 
        Wr.PutText(stderr, "unrecognized function name\"" & name & "\"\n");
        Process.Exit(1)
      END;
      RETURN f
    END;
  END MakeCorrectSolution;
  
<*UNUSED*>
PROCEDURE RaiseSol(
    READONLY ur: ARRAY OF LONG;  (* Coordinates of a function "urf" in "V[r]" *)
    READONLY bass: Basis;        (* Basis for space "V[s]". *)
    READONLY GLs: SPMatrix.T;    (* Choleski factor of scalar prod mat for "bas[s]". *)
    READONLY Frs: SPMatrix.T;    (* Matrix of products for "bas[r]" and "bas[s]". *)
    solFileName: TEXT;           (* Name template for output file. *)
    (* Work: *)
    VAR v: ARRAY OF LONG;        (* Dimension "d[s]". *)
    VAR y: ARRAY OF LONG;        (* Dimension "d[s]". *)
    VAR w: ARRAY OF LONG;        (* Dimension "d[s]". *)
    (* Output: *)
    VAR usf: SPPWFunction.T;     (* Approximation to "urf" in "V[s]" *)
    VAR us: ARRAY OF LONG;       (* Coordinates of "usf" relative to "bas[s]" *)
  ) =
  BEGIN
    SPMatrix.MulCol(Frs, ur, w);
    SPMatrix.DivCol(GLs, w, y);
    SPMatrix.DivRow( y, GLs, v);
    WITH d = LRN.Dist(us, v) DO
      Wr.PutText(stderr, "change = " & Fmt.LongReal(d) & "\n");
      us := v;
      usf := SPPWFunction.LinComb(us, bass);
      Wr.PutText(stderr, "writing " & solFileName & "\n");
      WITH wr = OpenWrite(solFileName & "-r.sfn") DO
        usf.write(wr); Wr.Close(wr)
      END;
      PrintSomeValues(usf);
    END
  END RaiseSol;

PROCEDURE LowerSol(
    READONLY us: ARRAY OF LONG;  (* Coordinates of a function "usf" in "V[s]" *)
    READONLY basr: Basis;        (* Basis for space "V[r]". *)
    READONLY GLr: SPMatrix.T;    (* Choleski factor of scalar prod mat for "bas[r]". *)
    READONLY Fsr: SPMatrix.T;    (* Matrix of products for "bas[s]" and "bas[r]". *)
    (*READONLY Nr: SPMatrix.T;*) (* Projection of "V[r]" onto "Z[r]". *)
    solFileName: TEXT;           (* Name template for output file. *)
    (* Work: *)
    (*VAR n: ARRAY OF LONG; *)   (* Dimension "d[r]". *)
    VAR y: ARRAY OF LONG;        (* Dimension "d[r]". *)
    VAR v: ARRAY OF LONG;        (* Dimension "d[r]". *)
    (* Input/Output: *)
    VAR urf: SPPWFunction.T;     (* Approximation to "usf" in "V[r]" *)
    VAR ur: ARRAY OF LONG;       (* Coordinates of "urf" relative to  "basr" *)
  ) =
  (*
    Given functions "usf" in "V[s]" and "urf" in "V[r]", replaces the 
    latter by a function "urf" in "V[r]" that combines a good approximation
    of "usf" in "V[r]" with the component of "urf" that is orthogonal
    to "V[s]". *)
  BEGIN
    (*SPMatrix.MulCol(Nr, ur, n); *)
    SPMatrix.MulRow(us, Fsr, v); 
    SPMatrix.DivCol(GLr, v, y); 
    SPMatrix.DivRow(y, GLr, v); 
    (* LRN.Add(n, v);*)
    WITH d = LRN.Dist(ur, v) DO
      Wr.PutText(stderr, "change = " & Fmt.LongReal(d) & "\n");
      ur := v;
      urf := SPPWFunction.LinComb(ur, basr);
      Wr.PutText(stderr, "writing " & solFileName & "\n");
      WITH wr = OpenWrite(solFileName & "-l.sfn") DO
        urf.write(wr); Wr.Close(wr)
      END;
      PrintSomeValues(urf);
    END
  END LowerSol;

PROCEDURE ReadMultiScaleData(
    minScale: NAT;  
    maxScale: NAT;  
    triName: TEXT;
    basName: TEXT; 
    matName: TEXT;
    solName: TEXT;
  ): REF MultiScaleData =
  BEGIN
    WITH
      mRef = NEW(REF MultiScaleData, maxScale+1), m = mRef^
    DO
      FOR r := minScale TO maxScale DO
        m[r] := ReadScaleData(triName, basName, matName, solName, r);
      END;
      FOR s := minScale+1 TO maxScale DO
        WITH
          r = s - 1
        DO
          m[r].F := ReadTransferMatrix(matName, s, r)
        END
      END;
      RETURN mRef
    END
  END ReadMultiScaleData;

PROCEDURE ReadScaleData(
    triName: TEXT;  (* Name prefix for triangulation *)
    baseName: TEXT; (* Name prefix for basis (incl. triang.) *)
    matName: TEXT;  (* Name prefix for basis matrices *)
    solName: TEXT;  (* Name of correct solution function. *) 
    r: NAT;         (* Scale index (0 = finest) *)
  ): ScaleData =
  BEGIN
    WITH
      tri = ReadTriang(TriangFileName(triName, r)),
      bas = ReadBasis(BasisFileName(baseName, r), tri),
      dim = NUMBER(bas^),
      ufCorr = MakeCorrectSolution(solName, tri),
      GL = ReadMatrix(MatrixFileName(matName, r) & "-GL"),
      HL = ReadMatrix(MatrixFileName(matName, r) & "-HL"),
      (* N = ReadMatrix(MatrixFileName(matName, r) & "-N"),*)
      xRef = NEW(REF ARRAY OF LONG, dim),
      yRef = NEW(REF ARRAY OF LONG, dim),
      zRef = NEW(REF ARRAY OF LONG, dim),
      uRef = NEW(REF ARRAY OF LONG, dim)
    DO
      RETURN ScaleData{
        tri := tri, bas := bas, GL := GL, HL := HL, 
        F := SPMatrix.T{rows := 0, cols := 0, r := NIL},
        ufCorr := ufCorr,
        (*N := N,*)
        x := xRef, y := yRef, z := zRef,
        uf := NIL, u := uRef
      }
    END
  END ReadScaleData;
    
PROCEDURE ReadTransferMatrix(
    matName: TEXT;  (* Name prefix for basis matrices *)
    s: NAT;         (* Row scale *)
    r: NAT;         (* Column scale *)
  ): SPMatrix.T =
  (*
    Reads the matrix "Fsr" whose element "[i,j]" is the 
    dot product of basis element "i" for space "Vs" 
    with basis element "j" from space "Vr". *)
  BEGIN
    RETURN ReadMatrix(TransferMatrixName(matName, s, r) & "-F")
  END ReadTransferMatrix;

PROCEDURE ReadTriang(name: TEXT): REF SPTriang.T =
  BEGIN
    WITH
      rd = OpenRead(name & ".tri"),
      tri = MakeTriRef(SPTriang.Read(rd))
    DO
      RETURN tri
    END
  END ReadTriang;

PROCEDURE ReadBasis(name: TEXT; tri: REF SPTriang.T): REF SPPWFunction.Basis =
  BEGIN
    WITH
      rd = OpenRead(name & ".bas"),
      basRef = SPPWFunction.ReadBasis(rd, tri)
    DO
      RETURN basRef
    END 
  END ReadBasis;
      
PROCEDURE ReadMatrix(name: TEXT): SPMatrix.T =
  VAR M: SPMatrix.T;
  BEGIN
    WITH
      rd = OpenRead(name & ".mat")
    DO
      SPMatrix.Read(rd, M);
      RETURN M
    END
  END ReadMatrix;
  
PROCEDURE ErrorFileName(outName: TEXT; scale: NAT): TEXT =
  BEGIN
    WITH
      scaleT = Fmt.Pad(Fmt.Int(scale), 2, '0')
    DO
      RETURN Fmt.F(outName, scaleT)
    END
  END ErrorFileName;
  
PROCEDURE SolutionFileName(outName: TEXT; scale: NAT; iter: NAT): TEXT =
  BEGIN
    WITH
      scaleT = Fmt.Pad(Fmt.Int(scale), 2, '0'),
      iterT = Fmt.Pad(Fmt.Int(iter), 4, '0')
    DO
      RETURN Fmt.F(outName, scaleT) & "-" & iterT
    END
  END SolutionFileName;

PROCEDURE TriangFileName(triName: TEXT; scale: NAT): TEXT =
  BEGIN
    WITH
      scaleT = Fmt.Pad(Fmt.Int(scale), 2, '0')
    DO
      RETURN Fmt.F(triName, scaleT)
    END
  END TriangFileName;

PROCEDURE BasisFileName(basName: TEXT; scale: NAT): TEXT =
  BEGIN
    WITH
      scaleT = Fmt.Pad(Fmt.Int(scale), 2, '0')
    DO
      RETURN Fmt.F(basName, scaleT)
    END
  END BasisFileName;

PROCEDURE MatrixFileName(matName: TEXT; scale: NAT): TEXT =
  BEGIN
    WITH
      scaleT = Fmt.Pad(Fmt.Int(scale), 2, '0')
    DO
      RETURN Fmt.F(matName, scaleT)
    END
  END MatrixFileName;

PROCEDURE TransferMatrixName(matName: TEXT; scale1, scale2: NAT): TEXT =
  BEGIN
    WITH
      scale1T = Fmt.Pad(Fmt.Int(scale1), 2, '0'),
      scale2T = Fmt.Pad(Fmt.Int(scale2), 2, '0')
    DO
      RETURN Fmt.F(matName, scale1T & "-" & scale2T)
    END
  END TransferMatrixName;

PROCEDURE PrintSomeValues(u: SPPWFunction.T) =
  BEGIN
    FOR i := -1 TO +1 DO
      FOR j := -1 TO +1 DO 
        FOR k := -1 TO +1 DO
          IF i # 0 OR j # 0 OR k # 0 THEN
            WITH
              x = FLOAT(i, LONG),
              y = FLOAT(j, LONG),
              z = FLOAT(k, LONG),
              p = LR3.Dir(LR3.T{x, y, z})
            DO
              Wr.PutText(stderr, "u(");
              Wr.PutText(stderr, FI(i,2) & "," & FI(j,2) & "," & FI(k,2));
              Wr.PutText(stderr, ") = " & FLR(u.eval(p), 9,6) & "\n");
            END
          END
        END
      END
    END
  END PrintSomeValues;
  
PROCEDURE FLR(x: LONGREAL; wid, prec: CARDINAL): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.LongReal(x, style:= Fmt.Style.Sci, prec:= prec), wid)
  END FLR;
  
PROCEDURE FI(x: INTEGER; wid: CARDINAL): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.Int(x), wid)
  END FI;

PROCEDURE MakeTriRef(READONLY tri: SPTriang.T): REF SPTriang.T = 
  (*
    Droga de Modula-3... *)
  BEGIN
    WITH t = NEW(REF SPTriang.T) DO t^ := tri; RETURN t END
  END MakeTriRef;

PROCEDURE GetOptions (): Options =
  VAR o: Options;
  BEGIN
    TRY
      WITH  pp = NEW(ParseParams.T).init(stderr) DO 
          
        pp.getKeyword("-triName");                               
        o.triName := pp.getNext();  
       
        pp.getKeyword("-baseName");                               
        o.baseName := pp.getNext();  
       
        pp.getKeyword("-rhsName");                               
        o.rhsName := pp.getNext();  
        
        pp.getKeyword("-solName");                               
        o.solName := pp.getNext();
        
        pp.getKeyword("-matName");                               
        o.matName := pp.getNext();
        
        pp.getKeyword("-minScale");
        o.minScale := pp.getNextInt(0, 99);
       
        pp.getKeyword("-maxScale");
        o.maxScale := pp.getNextInt(o.minScale+1, 99);
       
        pp.getKeyword("-outName");                               
        o.outName := pp.getNext();  

        o.relTol := NEW(REF LONGS, o.maxScale - o.minScale + 1);
        IF pp.keywordPresent("-relTol") THEN
          FOR scale := o.minScale TO o.maxScale DO
            o.relTol[scale - o.minScale] := pp.getNextLongReal(0.0d0,LAST(LONGREAL))
          END
        ELSE
          FOR scale := o.minScale TO o.maxScale DO
            o.relTol[scale - o.minScale] := 0.0d0
          END
        END;
      
        o.absTol := NEW(REF LONGS, o.maxScale - o.minScale + 1);
        IF pp.keywordPresent("-absTol") THEN
          FOR scale := o.minScale TO o.maxScale DO
            o.absTol[scale - o.minScale] := pp.getNextLongReal(0.0d0,LAST(LONGREAL))
          END
        ELSE
          FOR scale := o.minScale TO o.maxScale DO
            o.absTol[scale - o.minScale] := 0.0d0
          END
        END;
                  
        IF pp.keywordPresent("-maxIter") THEN
          o.maxIter := pp.getNextInt(0, LAST(CARDINAL))
        ELSE
          o.maxIter := 10
        END;
      
        pp.getKeyword("-order");
        o.order := pp.getNextInt(0,LAST(CARDINAL));
        
        pp.finish();
      END
    EXCEPT                                                            
    | ParseParams.Error =>                                              
        Wr.PutText(stderr, "Usage: SPMultiSolve \\\n");
        Wr.PutText(stderr, "  -triName FILEPREFIX \\\n");
        Wr.PutText(stderr, "  -baseName FILEPREFIX \\\n");
        Wr.PutText(stderr, "  -rhsName <functionName> \\\n");
        Wr.PutText(stderr, "  -solName <functionName> \\\n");
        Wr.PutText(stderr, "  -matName FILEPREFIX \\\n");
        Wr.PutText(stderr, "  -outName FILEPREFIX\\\n");
        Wr.PutText(stderr, "  -minScale NUMBER -maxScale NUMBER \\\n");
        Wr.PutText(stderr, "  [ -relTol NUMBER NUMBER ... ] \\\n");
        Wr.PutText(stderr, "  [ -absTol NUMBER NUMBER ... ] \\\n");
        Wr.PutText(stderr, "  [ -maxIter NUMBER ] \\\n");
        Wr.PutText(stderr, "  -order NUMBER \n");
        Process.Exit (1);
    END;
    RETURN o
  END GetOptions;

BEGIN
  Main()
END SPMultiSolve.