(* * * * * * * * * * *
 * Resource Aware ML *
 * * * * * * * * * * *
 *
 * File:
 *   matrix.raml
 *
 * Author:
 *   Jan Hoffmann (2010,2011)
 *
 * Resource Bound:
 *   O(n*m*k*l)
 *
 * Description:
 *   Matrix multiplication. 
 *
 *)

matrixMultList : (L(L(int)),L(L(L(int)))) -> L(L(int))

matrixMultList (acc,mm) = let _=tick(1) in
                       match mm with
                       | nil -> acc
                       | (m::ms) -> matrixMultList( matrixMult (acc,m),ms);



(* Multiply three matrices.  Needs maximal degree 3. *)

matrixMult3 : (L(L(int)),L(L(int)),L(L(int))) -> L(L(int))
matrixMult3 (m1,m2,m3) =  let _=tick(1) in
                           matrixMult (matrixMult (m1,m2), m3);

(*
eq(matrixmult3(A,B,C,D,E,F),1,[matrixmult(A,B,C,D),],[]).
*)


(* Multiply two matrices.  Needs maximal degree 3. *)

matrixMultOld : (L(L(int)),L(L(int))) -> L(L(int))

matrixMultOld (m1,m2) = let _=tick(1) in
                        matrixMult'(m1,transpose m2);


(* A new version that uses an alternate implementation
   of transpose.  Allows for a better bound.*)
matrixMult : (L(L(int)),L(L(int))) -> L(L(int))

matrixMult (m1,m2) =  let _=tick(1) in
                      matrixMult'(m1,transAcc(m2,makeBase m2));
(*

eq(matrixmult(A,B,C,D),1,[transAcc(C,D,C,0),matrixmultp(A,B,C,L),makeBase(C,D)],[L=A+D,A>=0,B>=0,C>=0,D>=0]).

*)

matrixMult' : (L(L(int)),L(L(int))) -> L(L(int))

matrixMult' (m1,m2) =  let _=tick(1) in
                       match m1 with
                        | nil -> nil
                        | (l::ls) -> (lineMult (l,m2))::(matrixMult' (ls,m2));
(*
eq(matrixmultp(A,B,C,D),1,[],[A=0]).
eq(matrixmultp(A,B,C,D),1,[linemult(B,C,D),matrixmultp(Ap,B,C,D)],[A>=1,Ap=A-1]).
*)

lineMult : (L(int),L(L(int))) -> L(int)

lineMult (l,m2) = let _=tick(1) in
                   match m2 with
                    | nil -> nil
                    | (x::xs) -> (mult (l,x))::(lineMult (l,xs));

(*
eq(linemult(A,B,C),1,[],[B=0]).
eq(linemult(A,B,C),1,[mult(A,C),linemult(A,Bp,C)],[B>=1,Bp=B-1]).
*)

mult : (L(int),L(int)) -> int

mult (l1,l2) = let _=tick(1) in
  match l1 with
    | nil -> +0
    | (x::xs) -> match l2 with
                   | nil -> +0
                   | (y::ys) -> x*y + (mult (xs,ys));

(*
eq(mult(A,B),1,[],[A=0]).
eq(mult(A,B),1,[],[A>=1,B=0]).
eq(mult(A,B),1,[mult(Ap,Bp)],[A>=1,B>=1, Ap=A-1,Bp=B-1]).
*) 

(*Transposing matrices*)

transpose : L(L(int)) -> L(L(int)) 
transpose m = let _=tick(1) in
   match m with
     | nil -> nil
     | (xs::xss) -> let (l,m') = split m in
                    match m' with
                      | nil -> nil
                      | (y::ys) -> l::transpose (y::ys);
(*
eq(transpose(R,C),1,[],[R=0]).
eq(transpose(R,C),1,[split(R,C)],[R>=1,C=1]).
eq(transpose(R,C),1,[split(R,C),transpose(R,Cp)],[R>=1,C>=2,Cp=C-1]).
*)
 
split : L(L(int)) -> (L(int),L(L(int)))

split m = let _=tick(1) in
   match m with
     | nil -> (nil,nil)
     | (l::ls) -> match l with
                    | nil -> (nil,nil)   (* This case won't happen! *)
                    | (x::xs) -> let (ys,m') = split ls 
                                 in (x::ys,xs::m');
(*

eq(out1_split(R,C),0,[],[R=0]).
eq(out1_split(R,C),1,[out1_split(Rp,C)],[R>=1,Rp=R-1]).

eq(out2_split(R,C),0,[],[R=0]).
eq(out2_split(R,C),1,[out2_split(Rp,C)],[R>=1,Rp=R-1]).

eq(out3_split(R,C),nat(C-1),[],[R=0]).
eq(out3_split(R,C),0,[out3_split(Rp,C)],[R>=1,Rp=R-1]).

eq(split(R,C),1,[],[R=0]).
eq(split(R,C),1,[split(Rp,C)],[R>=1, Rp=R-1]).
*)

transAcc : (L(L(int)), L(L(int))) -> L(L(int))

transAcc (m,base) =  let _=tick(1) in
                     match m with
                     | nil -> base
                     | (l::m') -> attach(l,transAcc(m',base));

(*

eq(transAcc(A,B,C,D),1,[],[A=0]).
eq(transAcc(A,B,C,D),1,[attach(B,Ap,L),transAcc(Ap,B,C,D)],[A>=1,Ap=A-1, L=D+Ap]).

eq(out1_transAcc(A,B,C,D),nat(C),[],[A=0]).
eq(out1_transAcc(A,B,C,D),0,[out1_transAcc(Ap,B,C,D)],[A>=1,Ap=A-1]).

eq(out2_transAcc(A,B,C,D),nat(D),[],[A=0]).
eq(out2_transAcc(A,B,C,D),1,[out2_transAcc(Ap,B,C,D)],[A>=1,Ap=A-1]).

*)


makeBase : L(L(int)) -> L(L(int))

makeBase m =  let _=tick(1) in
               match m with
               | nil -> nil
               | (l::m') -> mkBase l; 


mkBase : L(int) -> L(L(int))

mkBase m =    let _=tick(1) in
               match m with
               | nil -> nil
               | (l::m') -> nil::mkBase m';

(*

eq(out1_makeBase(A,B),0,[],[A=0]).    %row
eq(out1_makeBase(A,B),0,[out_mkBase(B)],[A=0]).
eq(out2_makeBase(A,B),0,[],[A>=0]).    %col
eq(out_mkBase(A),0,[],[A=0]).
eq(out_mkBase(A),1,[out_mkBase(Ap)],[A>=1,Ap=A-1]).

eq(makeBase(A,B),1,[],[A=0]).
eq(makeBase(A,B),1,[mkBase(B)],[A>=1]).
eq(mkBase(A),1,[],[A=0]).
eq(mkBase(A),1,[mkBase(Ap)],[A>=1,Ap=A-1]).

*)

attach : (L(int),L(L(int))) -> L(L(int))

attach (line,m) = let _=tick(1) in
                  match line with
                 | nil -> nil
                 | (x::xs) -> match m with
                                | nil -> nil
                                | (l::ls) -> (x::l)::attach(xs,ls);


(*

eq(attach(A,B,C),1,[],[A=0,B>=0]).
eq(attach(A,B,C),1,[],[A>=1,B=0]).
eq(attach(A,B,C),1,[attach(Ap,Bp,C)],[A>=1,Ap=A-1,Bp=B-1]).

*)

main = ()