(* A program to compare the accuracy of the five pivoting strategies: *)
(* direct diagonal, partial, partial scaled, total, and total scaled  *)
(* methods to solve a system AX = B up to 10^p cases.  *)

PROGRAM accuracy_comparison(OUTPUT);

CONST n = 5;                                      (* size of matrix A *)
      m = 1;                               (* column size of matrix B *)

TYPE indextypen = ARRAY[1..n] OF INTEGER;
     indextypem = ARRAY[1..m] OF INTEGER;
     matrixn = ARRAY[1..n,1..n] OF DOUBLE;
     matrixm = ARRAY[1..m,1..m] OF DOUBLE;
     matrixns = ARRAY[1..n,1..(n + m)] OF DOUBLE;
     matrixms = ARRAY[1..m,1..(m + n)] OF DOUBLE;
     matrixnm = ARRAY[1..n,1..m] OF DOUBLE;
     matrixmn = ARRAY[1..m,1..n] OF DOUBLE;
     row_max_typen = ARRAY[1..n] OF DOUBLE;
     row_max_typem = ARRAY[1..m] OF DOUBLE;

VAR index, rindex : indextypen;
    indext, rindext : indextypem;
    anbm,                          (* the solution of the system AX = B *)
    xprime, anbmprime, yprime, b, p, pp : matrixnm;
    xprimet, anbmt, yprimet, bt, at : matrixmn;
    aprime : matrixns;
    bprime : matrixms;
    aa : matrixn;
    bb : matrixm;
    ok : BOOLEAN;
    c, ans : CHAR;
    row_max : row_max_typen;
    row_maxt : row_max_typem;
    po,         (* the exponent of the total number of comparison cases *)
    seed, i, j, k, l, ll, en, em, casenum1, casenum2 : INTEGER;
    count : ARRAY['0'..'4'] OF INTEGER;
    h1, hmin, tol : DOUBLE;
    h : ARRAY['0'..'4'] OF DOUBLE;

[EXTERNAL,ASYNCHRONOUS] FUNCTION MTH$RANDOM(VAR seed:INTEGER):REAL;EXTERN;

[EXTERNAL] PROCEDURE product
             (VAR xy:ARRAY[xyr1..xyrh:INTEGER;xycl..xych:INTEGER] OF DOUBLE;
              x:ARRAY[xr1..xrh:INTEGER;xcl..xch:INTEGER] OF DOUBLE;
              y:ARRAY[yr1..yrh:INTEGER;ycl..ych:INTEGER] OF DOUBLE);EXTERN;

[EXTERNAL] PROCEDURE A_n_B_m
             (n,m : INTEGER;
              pivoting_strategy:CHAR;
              tol : DOUBLE;
              a : ARRAY[arl..arh:INTEGER;acl..ach:INTEGER] OF DOUBLE;
              b : ARRAY[brl..brh:INTEGER;bcl..bch:INTEGER] OF DOUBLE;
              VAR at: ARRAY[atrl..atrh:INTEGER;atcl..atch:INTEGER] OF DOUBLE;
              VAR bt: ARRAY[btrl..btrh:INTEGER;btcl..btch:INTEGER] OF DOUBLE;
              VAR aprime:ARRAY[aprl..aprh:INTEGER;apcl..apch:INTEGER] OF DOUBLE;
              VAR bprime:ARRAY[bprl..bprh:INTEGER;bpcl..bpch:INTEGER] OF DOUBLE;
              VAR ok :BOOLEAN;
              VAR anbm, anbmprime, xprime, yprime:
                    ARRAY[xrl..xrh:INTEGER;xcl..xch:INTEGER] OF DOUBLE;
              VAR anbmt, xprimet, yprimet:
                    ARRAY[trl..trh:INTEGER;tcl..tch:INTEGER] OF DOUBLE;
              VAR row_max:ARRAY[rl..rh:INTEGER] OF DOUBLE;
              VAR row_maxt:ARRAY[rtl..rth:INTEGER] OF DOUBLE;
              VAR index:ARRAY[il..ih:INTEGER] OF INTEGER;
              VAR indext:ARRAY[itl..ith:INTEGER] OF INTEGER;
              VAR rindex:ARRAY[ril..rih:INTEGER] OF INTEGER;
              VAR rindext:ARRAY[ritl..rith:INTEGER] OF INTEGER);EXTERN;

BEGIN
po := 5;                                  (* May be read from the keybord *)
en := -1;                                 (* if input is declared.        *)
em := 1;
tol := 1.0E-20;
casenum1 := 1;               (* To be used to control the total number of *)
casenum2 := 1;               (* random cases. See the next next comment.  *)
FOR ll := 1 TO po DO
  BEGIN 
  casenum2 := casenum2 * 10;
  FOR c := '0' TO '4' DO count[c] := 0;
  FOR k := casenum1 TO casenum2 DO
    BEGIN
    seed := k;               (* May be read from keybord at the beginning *)
    FOR i := 1 TO n DO       (* or at present, or any appropriate time.   *)
      FOR j := 1 TO n DO aa[i,j] := MTH$RANDOM(seed);
    FOR l := 1 to m DO
      FOR i := 1 TO n DO
        BEGIN
        b[i,l] := 0;
        FOR j := 1 TO n DO b[i,l] := aa[i,j] + b[i,l]
        END;
    FOR c := '0' TO '4' DO
      BEGIN
      ok := true;
      a_n_b_m(en,em,c,tol,aa,b,at,bt,aprime,bprime,ok,anbm,anbmprime,
              xprime,yprime,anbmt,xprimet,yprimet,
              row_max,row_maxt,index,indext,rindex,rindext);
      IF ok THEN
        BEGIN
        h[c] := 0;
        FOR l := 1 TO m DO
          FOR i := 1 TO n DO
            BEGIN 
            h1 := anbm[i,l] - 1; 
            h[c] := h[c] + h1 * h1 
            END
        END
      ELSE h[c] := 10
      END;
    hmin := 1;
    FOR c := '0' TO '4' DO IF h[c] < hmin THEN hmin := h[c];
    FOR c := '0' TO '4' DO IF h[c] = hmin THEN count[c] := count[c] + 1
    END;
  writeln('In',casenum2 - casenum1 + 1,           (* The total number of *)
          ' random cases there are');             (* random cases.       *)
  FOR c := '0' TO '4' DO BEGIN
                         write(count[c],' cases in which ');
                         CASE c OF
                           '0': write('Diagonal');
                           '1': write('Partial');
                           '2': write('Partial Scaled');
                           '3': write('Total');
                           '4': write('Total Scaled');
                         END;
                         writeln(' strategy is the most accurate.')
                         END;
  writeln; writeln; writeln
  END
END.