UNSAFE MODULE TestVecLoops 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 = 10000 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],
      (*
      x2 = x[2],
      x3 = x[3],
      x4 = x[4],
      x5 = x[5],
      *)
      
      s = x0*x0 + x1*x1 (* + x2*x2 + x3*x3 + x4*x4 + x5*x5 *)
    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 (* , s, s , 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;
      (*
      r[2] := s;
      r[3] := s;
      r[4] := s;
      r[5] := 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[2]*m[2,0] + x[3]*m[3,0] + x[4]*m[4,0] + x[5]*m[5,0] *),
      x[0]*m[0,1] + x[1]*m[1,1] (* + x[2]*m[2,1] + x[3]*m[3,1] + x[4]*m[4,1] + x[5]*m[5,1] *)
      (*,
      x[0]*m[0,2] + x[1]*m[1,2] (* + x[2]*m[2,2] + x[3]*m[3,2] + x[4]*m[4,2] + x[5]*m[5,2] *),
      x[0]*m[0,3] + x[1]*m[1,3] (* + x[2]*m[2,3] + x[3]*m[3,3] + x[4]*m[4,3] + x[5]*m[5,3] *),
      x[0]*m[0,4] + x[1]*m[1,4] (* + x[2]*m[2,4] + x[3]*m[3,4] + x[4]*m[4,4] + x[5]*m[5,4] *),
      x[0]*m[0,5] + x[1]*m[1,5] (* + x[2]*m[2,5] + x[3]*m[3,5] + x[4]*m[4,5] + x[5]*m[5,5] *),
      *)
    }
  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] (*, x2 = x[2], x3 = x[3], x4 = x[4], x5 = x[5] *),
      
      m00 = m[0,0], m01 = m[0,1] (*, m02 = m[0,2], m03 = m[0,3], m04 = m[0,4], m05 = m[0,5] *),
      m10 = m[1,0], m11 = m[1,1] (*, m12 = m[1,2], m13 = m[1,3], m14 = m[1,4], m15 = m[1,5] *)
      (*,                           
      m20 = m[2,0], m21 = m[2,1] (*, m22 = m[2,2], m23 = m[2,3], m24 = m[2,4], m25 = m[2,5] *),
      m30 = m[3,0], m31 = m[3,1] (*, m32 = m[3,2], m33 = m[3,3], m34 = m[3,4], m35 = m[3,5] *),
      m40 = m[4,0], m41 = m[4,1] (*, m42 = m[4,2], m43 = m[4,3], m44 = m[4,4], m45 = m[4,5] *),
      m50 = m[5,0], m51 = m[5,1] (*, m52 = m[5,2], m53 = m[5,3], m54 = m[5,4], m55 = m[5,5] *),
      *)
    DO
      RETURN Vec{
        x0*m00 + x1*m10 (* + x2*m20 + x3*m30 + x4*m40 + x5*m50 *),
        x0*m01 + x1*m11 (* + x2*m21 + x3*m31 + x4*m41 + x5*m51 *)
        (*,
        x0*m02 + x1*m12 (* + x2*m22 + x3*m32 + x4*m42 + x5*m52 *),
        x0*m03 + x1*m13 (* + x2*m23 + x3*m33 + x4*m43 + x5*m53 *),
        x0*m04 + x1*m14 (* + x2*m24 + x3*m34 + x4*m44 + x5*m54 *),
        x0*m05 + x1*m15 (* + x2*m25 + x3*m35 + x4*m45 + x5*m55 *),
        *)
      }
    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 TestVecLoops.