UNSAFE MODULE TestVecLoops2 EXPORTS Main;
(*
  Computes the costs of "Math.VecLoops(x,y)" and "Math.sqrt(x*x + y*y)". *)

IMPORT Wr, Fmt, Uresource, Thread;

FROM Stdio IMPORT stderr;
FROM Fmt IMPORT Style;

CONST
  Microsecond = 1.0d0;
  Second = 1000000.0d0 * Microsecond;
  
CONST N = 2;

TYPE 
  Vec = ARRAY [0..N-1] OF LONGREAL;
  Mat = ARRAY [0..N-1] OF ARRAY [0..N-1] OF LONGREAL;
  VecFn = PROCEDURE (READONLY x: Vec; READONLY m: Mat): Vec;

PROCEDURE DoIt() =
  VAR F: VecFn;
      T: TEXT;
  VAR x, z: Vec;
      m: Mat;
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    FOR k := 0 TO 11 DO
      CASE k OF 
      |  0 => F := VecEmpty; T := "empty";
      |  1 => F := VecLoop;  T := "loop ";
      |  2 => F := VecWith;  T := "with ";
      |  3 => F := VecUnrll; T := "unrll";
      |  4 => F := VecDotsT; T := "dots ";
      |  5 => F := VecFullT; T := "full ";
      |  6 => F := VecAssgn; T := "assgn";
      |  7 => F := VecMLoop; T := "mloop";
      |  8 => F := VecMSwap; T := "mswap";
      |  9 => F := VecMWith; T := "mwith";
      | 10 => F := VecMUnrl; T := "munrl";
      | 11 => F := VecMUnrW; T := "munrw";
      END;
      x := Vec{0.1d0, ..};
      m := Mat{Vec{12.7d0, ..}, ..};
      WITH start = Now(), N = 100000 DO
        FOR i := 0 TO N-1 DO 
          z := F(x, m);
          x[0] := x[0] + 1024.0d0
        END;
        WITH stop = Now(), t = (stop-start)/FLOAT(N, LONGREAL)/Microsecond DO
          Wr.PutText(stderr, T);
          Wr.PutText(stderr, " = ");
          Wr.PutText(stderr, Fmt.LongReal(t, prec := 2, style := Style.Fix));
          Wr.PutText(stderr, "\n");
        END
      END
    END
  END DoIt;
  
PROCEDURE VecEmpty(<*UNUSED*> READONLY x: Vec; <*UNUSED*> READONLY m: Mat): Vec =
  BEGIN
    RETURN Vec{0.0d0, ..}
  END VecEmpty;

PROCEDURE VecLoop(READONLY x: Vec; <*UNUSED*> READONLY m: Mat): Vec =
  VAR s: LONGREAL := 0.0d0;
  BEGIN
    FOR i := 0 TO N-1 DO 
      s := s + x[i] * x[i]
    END;
    RETURN Vec{s, ..}
  END VecLoop;

PROCEDURE VecWith(READONLY x: Vec; <*UNUSED*> READONLY m: Mat): Vec =
  VAR s: LONGREAL := 0.0d0;
  BEGIN
    FOR i := 0 TO N-1 DO 
      WITH xi = x[i] DO
        s := s + xi * xi
      END
    END;
    RETURN Vec{s, ..}
  END VecWith;

PROCEDURE VecUnrll(READONLY x: Vec; <*UNUSED*> READONLY m: Mat): Vec =
  BEGIN
    WITH
      x0 = x[0],
      x1 = x[1],
      
      s = x0*x0 + x1*x1
    DO
      RETURN Vec{s, ..}
    END
  END VecUnrll;

PROCEDURE VecDotsT(READONLY x: Vec; <*UNUSED*> READONLY m: Mat): Vec =
  BEGIN
    WITH s = x[0]*x[0] DO
      RETURN Vec{s, ..}
    END
  END VecDotsT;

PROCEDURE VecFullT(READONLY x: Vec; <*UNUSED*> READONLY m: Mat): Vec =
  BEGIN
    WITH s = x[0]*x[0] DO
      RETURN Vec{s, s}
    END
  END VecFullT;

PROCEDURE VecAssgn(READONLY x: Vec; <*UNUSED*> READONLY m: Mat): Vec =
  VAR r: Vec;
  BEGIN
    WITH
      s = x[0]*x[0]
    DO
      r[0] := s;
      r[1] := s;
    END;
    RETURN r
  END VecAssgn;

PROCEDURE VecMLoop(READONLY x: Vec; READONLY m: Mat): Vec =
  (* Vec-matrix product with vanilla loops: *)
  VAR r: Vec; s: LONGREAL;
  BEGIN
    FOR i := 0 TO N-1 DO
      s := 0.0d0;
      FOR k := 0 TO N-1 DO
        s := s + x[k]*m[k,i]
      END;
      r[i] := s
    END;
    RETURN r
  END VecMLoop;

PROCEDURE VecMSwap(READONLY x: Vec; READONLY m: Mat): Vec =
  (* Vec-matrix product with reversed loops: *)
  VAR r: Vec;
  BEGIN
    r := Vec{0.0d0, ..};
    FOR k := 0 TO N-1 DO
      FOR i := 0 TO N-1 DO
        r[i] := r[i] + x[k]*m[k,i]
      END
    END;
    RETURN r
  END VecMSwap;

PROCEDURE VecMWith(READONLY x: Vec; READONLY m: Mat): Vec =
  (* Vec-matrix product with reversed loops and WITH on each row *)
  VAR r: Vec;
  BEGIN
    r := Vec{0.0d0, ..};
    FOR k := 0 TO N-1 DO
      WITH mk = m[k], xk = x[k] DO
        FOR i := 0 TO N-1 DO
          WITH ri = r[i] DO
            ri := ri + xk*mk[i]
          END
        END;
      END
    END;
    RETURN r
  END VecMWith;

PROCEDURE VecMUnrl(READONLY x: Vec; READONLY m: Mat): Vec =
  (* Unrolled vec-matrix product *)
  BEGIN
    RETURN Vec{
      x[0]*m[0,0] + x[1]*m[1,0],
      x[0]*m[0,1] + x[1]*m[1,1]
    }
  END VecMUnrl;

PROCEDURE VecMUnrW(READONLY x: Vec; READONLY m: Mat): Vec =
  (* Unrolled vec-matrix product with WITH on elements: *)
  BEGIN
    WITH
      x0 = x[0], x1 = x[1],
      
      m00 = m[0,0], m01 = m[0,1],
      m10 = m[1,0], m11 = m[1,1]
    DO
      RETURN Vec{
        x0*m00 + x1*m10,
        x0*m01 + x1*m11
      }
    END
  END VecMUnrW;

PROCEDURE Now(): LONGREAL =
  VAR ru: Uresource.struct_rusage;
  BEGIN
    WITH
      rCode = Uresource.getrusage(Uresource.RUSAGE_SELF, ADR(ru))
    DO
      <* ASSERT rCode = 0 *>
      RETURN 
        FLOAT(ru.ru_utime.tv_sec, LONGREAL) * Second +
        FLOAT(ru.ru_utime.tv_usec, LONGREAL) * Microsecond
    END
  END Now;

BEGIN
  DoIt()
END TestVecLoops2.