// File RTCG8.cs --- sparse matrices
// sestoft@itu.dk * 2002-09-23

// A matrix A is an array of arrays (rows) of doubles, so A[i][j] is
// row i column j.

// We consider multiplications A*B where B is a sparse matrix, and
// for given B generate code that computes A*B when given a matrix A.

using System;
using System.Collections;
using System.Reflection;
using System.Reflection.Emit;

public class RTCG8 {
  public static void Main(string[] args) {
    int dim1 = int.Parse(args[0]),
      dim2 = int.Parse(args[1]),
      dim3 = int.Parse(args[2]),
      nonzero = int.Parse(args[3]),
      count = int.Parse(args[4]);
    
    // Generate some random matrices 
    double[][] A = randomSparse(dim1, dim2, dim1*dim2),
      B = randomSparse(dim2, dim3, nonzero),
      R1 = newMatrix(dim1, dim3),
      R2 = newMatrix(dim1, dim3),
      R3 = newMatrix(dim1, dim3),
      R4 = newMatrix(dim1, dim3);

    // Generate a specialized sparse multiplication method
    AssemblyName assemblyName = new AssemblyName();
    assemblyName.Name = "myassembly";
    // Build: run-only assembly
    AssemblyBuilder assemblyBuilder = 
      AppDomain.CurrentDomain.DefineDynamicAssembly(assemblyName,
                                                    AssemblyBuilderAccess.Run);
    // Build: module mymodule
    ModuleBuilder moduleBuilder = 
      assemblyBuilder.DefineDynamicModule("mymodule");
    // Build: public class MyClass { ... }
    TypeBuilder typeBuilder = 
      moduleBuilder.DefineType("MyClass", 
                               TypeAttributes.Class | TypeAttributes.Public, 
                               typeof(Object));
    Type ty;

    // Build: public static void sparseMultB(double[][] A, double[][] R)
    // and possibly additional copies (for timing)
    {
      Timer t = new Timer();
      SparseMatrix sparseB = new SparseMatrix(B);
      sparseMultGen(typeBuilder, "sparseMultB", dim1, dim2, sparseB);
      // If there are 6 or more arguments, build the method repeatedly
      int methodCount = args.Length >= 6 ? int.Parse(args[5]) : 1;
      for (int i=methodCount-1; i>0; i--)
	sparseMultGen(typeBuilder, "sparseMultB"+i, dim1, dim2, sparseB);
      // Load the class file into the VM
      ty = typeBuilder.CreateType();
      // Probably we need to call the method to convince CLR to JIT it?
      Console.WriteLine("Generating " + methodCount 
			+ " copies of sparseMultB: " + t.Check() + " sec");
    }
    { 
      Timer t = new Timer();
      for (int i=count; i>0; i--)
        matrixMult(A, B, R1);
      Console.WriteLine("matrixMult: " + t.Check() + " sec");
    }
    { 
      Timer t = new Timer();
      for (int i=count; i>0; i--)
        sparseMult(A, B, R2);
      Console.WriteLine("sparseMult: " + t.Check() + " sec");
    }
    { 
      Timer t = new Timer();
      SparseMatrix sparseB = new SparseMatrix(B);
      for (int i=count; i>0; i--)
        sparseMult(A, sparseB, R3);
      Console.WriteLine("sparseMult, two-phase: " + t.Check() + " sec");
    }
    {
      MethodInfo m = ty.GetMethod("sparseMultB");
      Timer t = new Timer();
      for (int i=count; i>0; i--)
        m.Invoke(null, new Object[] { A, R4 });
      Console.WriteLine("Generated sparseMultB: " + t.Check() + " sec");
    }
    Console.WriteLine("R1 == R2 is " + equal(R1, R2));
    Console.WriteLine("R2 == R3 is " + equal(R2, R3));
    Console.WriteLine("R3 == R4 is " + equal(R3, R4));
  }

  // Compute A*B in R 
  // Assume A, B and R are rectangular and non-empty

  public static void matrixMult(double[][] A, double[][] B, double[][] R) {
    int aRows = A.Length, aCols = A[0].Length,
      bRows = B.Length, bCols = B[0].Length,
      rRows = R.Length, rCols = R[0].Length;
    if (aCols != bRows || aRows != rRows || bCols != rCols)
      throw new Exception("Matrix dimension mismatch");
    else {
      for (int i=0; i<rRows; i++)
        for (int j=0; j<rCols; j++) {
          double sum = 0.0;
          for (int k=0; k<aCols; k++)
            sum += A[i][k] * B[k][j];
          R[i][j] = sum;
        }
    }
  }

  // Compute A*B in R in two steps, first building a sparse
  // representation of B (an array of lists of nonzero column elements)

  // Assume A, B and R are rectangular and non-empty

  public static void sparseMult(double[][] A, double[][] B, double[][] R) {
    int aRows = A.Length, aCols = A[0].Length,
      bRows = B.Length, bCols = B[0].Length,
      rRows = R.Length, rCols = R[0].Length;
    if (aCols != bRows || aRows != rRows || bCols != rCols)
      throw new Exception("Matrix dimension mismatch");
    else {
      SparseMatrix sparseB = new SparseMatrix(B);
      sparseMult(A, sparseB, R);
    }
  }

  // Compute A*B in R, where B is represented as a SparseMatrix
  // Assume A, B and R conform, and B and R are rectangular and non-empty

  private static void sparseMult(double[][] A, SparseMatrix B, double[][] R) {
    int rRows = R.Length, rCols = R[0].Length;
    for (int i=0; i<rRows; i++) {
      double[] Ai = A[i];
      double[] Ri = R[i];
      for (int j=0; j<rCols; j++) {
        double sum = 0.0;
	foreach (NonZero nz in B.getCol(j))
          sum += Ai[nz.k] * nz.Bkj;
        Ri[j] = sum;
      }
    }
  }

  // Given a SparseMatrix B, generate code to compute A*B in R for any A
  // Build: public static void sparseMultB(double[][] A, double[][] R)
  // in class object co

  private static void sparseMultGen(TypeBuilder typeBuilder, String name,
				    int aRows, int aCols, SparseMatrix B) {
    Timer t = new Timer();
    MethodBuilder methodBuilder =
      typeBuilder.DefineMethod(name,
			       MethodAttributes.Static | MethodAttributes.Public,
			       typeof(void),
			       new Type[] { typeof(double[][]), 
					    typeof(double[][]) });
    ILGenerator ilg = methodBuilder.GetILGenerator();
    // The generated method's parameters: double[][] A, double[][] R
    // arg 0 is A of type double[][]
    // arg 1 is R of type double[][]
    // Variable varR = scope.addVariable(jvmg, double2D_type, "R");
    ilg.DeclareLocal(typeof(double[]));		   // local 0 = Ai
    ilg.DeclareLocal(typeof(double[]));		   // local 1 = Ri
    ilg.DeclareLocal(typeof(int));		   // local 2 = i
    ilg.Emit(OpCodes.Ldc_I4_0);
    ilg.Emit(OpCodes.Stloc_S, 2);		   // i = 0
    
    Label loop = ilg.DefineLabel();
    ilg.MarkLabel(loop);                           // do {
    ilg.Emit(OpCodes.Ldarg_0);
    ilg.Emit(OpCodes.Ldloc_2);
    ilg.Emit(OpCodes.Ldelem_Ref);
    ilg.Emit(OpCodes.Stloc_0);                     // Ai = A[i]
    ilg.Emit(OpCodes.Ldarg_1);
    ilg.Emit(OpCodes.Ldloc_2);
    ilg.Emit(OpCodes.Ldelem_Ref);
    ilg.Emit(OpCodes.Stloc_1);                     // Ri = R[i]
    for (int j=0; j<B.cols; j++) {
      ilg.Emit(OpCodes.Ldloc_1);                   // Load Ri
      ilg.Emit(OpCodes.Ldc_I4, j); 
      ilg.Emit(OpCodes.Ldc_R8, 0.0);               // sum = 0.0
      foreach (NonZero nz in B.getCol(j)) {
	ilg.Emit(OpCodes.Ldc_R8, nz.Bkj);          // load B[k][j]
	ilg.Emit(OpCodes.Ldloc_0);                 // load A[i]
	ilg.Emit(OpCodes.Ldc_I4, nz.k);
	ilg.Emit(OpCodes.Ldelem_R8);               // load A[i][k]
	ilg.Emit(OpCodes.Mul);                     // prod = A[i][k]*B[k][j]
        ilg.Emit(OpCodes.Add);                     // sum += prod
      }
      ilg.Emit(OpCodes.Stelem_R8);                 // R[i][j] = sum
    }
    ilg.Emit(OpCodes.Ldloc_2);
    ilg.Emit(OpCodes.Ldc_I4_1);
    ilg.Emit(OpCodes.Add);
    ilg.Emit(OpCodes.Stloc_2);                     // i++
    ilg.Emit(OpCodes.Ldloc_2);
    ilg.Emit(OpCodes.Ldc_I4, aRows); 
    ilg.Emit(OpCodes.Blt, loop);                   // } while (i<aRows);
    ilg.Emit(OpCodes.Ret);
  }

  // A SparseMatrix has a dimension, and an array of the NonZeros in
  // each of B's columns

  private class SparseMatrix {
    public readonly int rows, cols;
    readonly ArrayList[] /* of NonZero */ nonzeros;

    public SparseMatrix(double[][] B) {
      rows = B.Length; cols = B[0].Length;
      nonzeros = new ArrayList[cols];
      for (int j=0; j<cols; j++) {
        nonzeros[j] = new ArrayList();
        for (int k=0; k<rows; k++) 
          if (B[k][j] != 0.0)
            nonzeros[j].Add(new NonZero(k, B[k][j]));
      }
    }
    
    public ArrayList /* of NonZero */ getCol(int j) {
      return nonzeros[j];
    }
  }

  // A pair of a row number k and a non-zero element B[k][-]

  private class NonZero {
    public readonly int k;
    public readonly double Bkj;

    public NonZero(int k, double Bkj) {
      this.k = k; this.Bkj = Bkj;
    }
  }
    
  // Generate a rectangular matrix with at most n non-zero elements,

  public static double[][] randomSparse(int rows, int cols, int n) {
    Random rnd = new Random();
    double[][] R = newMatrix(rows, cols);      // All zeros initially
    for (int k=0; k<n; k++) {
      int i = (int)(rows * rnd.NextDouble()), 
	j = (int)(cols * rnd.NextDouble());
      R[i][j] = rnd.NextDouble();
    }
    return R;
  }

  // Allocate matrices

  public static double[][] newMatrix(int rows, int cols) {
    double[][] res = new double[rows][];
    for (int i=0; i<rows; i++)
      res[i] = new double[cols];
    return res;
  }

  // Compare matrices

  public static bool equal(double[][] A, double[][] B) {
    int aRows = A.Length, aCols = A[0].Length,
      bRows = B.Length, bCols = B[0].Length;
    if (aCols != bCols || aRows != bRows)
      return false;
    else {
      for (int i=0; i<aRows; i++)
        for (int j=0; j<aCols; j++) 
          if (A[i][j] != B[i][j])
            return false;
      return true;
    }
  }
}

// Crude timing utility ----------------------------------------
   
public class Timer {
  private DateTime start;

  public Timer() {
    start = DateTime.Now;
  }

  public double Check() {
    TimeSpan dur = DateTime.Now - start;
    return dur.TotalSeconds;
  }
}