/* Routines to manipulate PE array matrices
 * represented by the 'dmat' structure
 * HM 11/92
 *
 * last arg of mat_*() routines is output
 */

#include "mat.h"
extern void describe_mat();
extern void exit();

/* create space on PE array for a new matrix of type dmat
 * and return a pointer to the new dmat.
 */
struct dmat * new_mat(m, n)
   int m, n;
{
   extern char * malloc();
   struct dmat *p;
   extern plural char * p_malloc();
   int nxtiles, nytiles, ntiles;

   nytiles = (m-1) / nyproc + 1;
   nxtiles = (n-1) / nxproc + 1;
   ntiles = nxtiles * nytiles;

   p = (struct dmat *) malloc(sizeof(struct dmat));
   p->base = (plural float *) p_malloc(ntiles * sizeof(float));
   if (p == (struct dmat *)0 || p->base == (plural float *)0) {
      fprintf(stderr, "new_mat(%d, %d): out of memory\n", m, n);
      exit(127);
      }
   p->nrows = m;
   p->ncols = n;
   p->nxtiles = nxtiles;
   p->ntiles = ntiles;
   { 
     int i;
     for (i=0; i < p->ntiles; i++) 
	*(p->base + i) = (plural float) 0;
   }
   return p;
   }

/* read and write from ACU to PE matrix element A[i][j]
 * NOTE: can be improved by replacing % with bit mask!
 */
void set_mat_elem(A, i, j, v)
   struct dmat *A; int i, j; float v;
{
   int xpos, ypos, tile;
   int tilerow, tilecol, nxtiles;
   ypos = i % nyproc;	/* row in processor array */
   xpos = j % nxproc;	/* col in processor array */
   tilerow = i/nyproc;	/* tile row */
   tilecol = j/nxproc;	/* tile col */
   nxtiles = (A->ncols - 1) / nxproc + 1;
   tile  = tilerow * nxtiles + tilecol;
   proc[ypos][xpos].*(A->base + tile) = v;
   }

float get_mat_elem(A, i, j)
   struct dmat *A; int i, j;
{
   int xpos, ypos, tile;
   int tilerow, tilecol, nxtiles;
   ypos = i % nyproc;	/* row */
   xpos = j % nxproc;	/* col */
   tilerow = i/nyproc;	/* tile row */
   tilecol = j/nxproc;	/* tile col */
   nxtiles = (A->ncols - 1) / nxproc + 1;
   tile  = tilerow * nxtiles + tilecol;
   return proc[ypos][xpos].*(A->base + tile);
   }

/* copy matrix A onto matrix B
 */
void mat_copy(A, B)
   struct dmat *A, *B;
{
   int i;
   for (i=0; i < A->ntiles; i++)
      *(B->base + i) = *(A->base + i);
   }

/* add two matrices:  mat_add(A, B, C) gives C := A + B
 */
void mat_add(A, B, C) 
   struct dmat *A, *B, *C;
{
   int i;

   if (A->nrows != B->nrows || A->ncols != B->ncols) {
      fprintf(stderr, "mat_add(): matrices do not conform\n");
      exit(127);
      }
   for (i=0; i < A->ntiles; i++)
      *(C->base + i) = *(A->base + i) + *(B->base + i);
   }

/* subtract two matrices:  mat_sub(A, B, C) gives C := A - B
 */
void mat_sub(A, B, C) 
   struct dmat *A, *B, *C;
{
   int i;

   if (A->nrows != B->nrows || A->ncols != B->ncols) {
      fprintf(stderr, "mat_sub(): matrices do not conform\n");
      exit(127);
      }
   for (i=0; i < A->ntiles; i++)
      *(C->base + i) = *(A->base + i) - *(B->base + i);
   }

/* element-wise multiply matrices:  
 * mat_elem_mult(A, B, C) gives C := A .* B
 */
void mat_elem_mult(A, B, C) 
   struct dmat *A, *B, *C;
{
   int i;

   if (A->nrows != B->nrows || A->ncols != B->ncols) {
      fprintf(stderr, "mat_elem_mult(): matrices do not conform\n");
      exit(127);
      }
   for (i=0; i < A->ntiles; i++)
      *(C->base + i) = *(A->base + i) * *(B->base + i);
   }

/* element-wise divide matrices:  
 * mat_elem_div(A, B, C) gives C := A ./ B
 */
void mat_elem_div(A, B, C) 
   struct dmat *A, *B, *C;
{
   int i;

   if (A->nrows != B->nrows || A->ncols != B->ncols) {
      fprintf(stderr, "mat_elem_div(): matrices do not conform\n");
      exit(127);
      }
   for (i=0; i < A->ntiles; i++)
      *(C->base + i) = *(A->base + i) / *(B->base + i);
   }

/* matrix multiply:  mat_mult(A, B, C) gives C := A * B
 */
void mat_mult(A, B, C) 
   struct dmat *A, *B, *C;
{
   void fp_matmul();
   if (A->ncols != B->nrows || A->nrows != C->nrows || B->ncols != C->ncols) {
      fprintf(stderr, "mat_mult(): matrices do not conform\n");
      describe_mat(A); describe_mat(B); describe_mat(C);
      exit(127);
      }
   fp_matmul(C->base, A->base, B->base, A->nrows, A->ncols, B->ncols);
   }

/* multiply each element of a matrix by a scalar:  B := s .* A
 */
void mat_scalar_mult(s, A, B)
   float s;
   struct dmat *A, *B;
{
   int i;
   for (i=0; i < A->ntiles; i++)
      *(B->base + i) = *(A->base + i) * (plural float) s;
   }

/* apply a unary dp function elementwise: B := f(A)
 */
void mat_apply_unary(f, A, B)
   plural float (*f)();
   struct dmat *A, *B;
{
   int i;
   for (i=0; i < A->ntiles; i++)
      *(B->base + i) = (*f)( *(A->base + i) );
   }

/* apply a binary dp function elementwise: C := f(A, B)
 */
void mat_apply_binary(f, A, B, C)
   plural float (*f)();
   struct dmat *A, *B, *C;
{
   int i;
   for (i=0; i < A->ntiles; i++)
      *(C->base + i) = (*f)( *(A->base + i), *(B->base + i) );
   }

/* transpose a matrix: B is the transpose of A
 */
void mat_transp(A, B)
   struct dmat *A, *B;
{
   if (A->nrows != B->ncols || B->nrows != A->ncols) {
      fprintf(stderr, "mat_transp(): matrices do not conform\n");
      exit(127);
      }
   fp_transpose(A->nrows, A->ncols, A->base, B->base);
   }

/* Make A a diagonal matrix
 * note: this routine uses an ACU loop & so is slower than
 * other mat routines, and should not be used inside a loop.
 */
void diag_mat(A)
   struct dmat *A;
{
   int i, j, k;
   k = A->nrows < A->ncols ? A->nrows : A->ncols;
   for (i=0; i<k; i++)
      for (j=0; j<k; j++) 
	 if (i==j)
	    set_mat_elem(A,i,j, (float) 1);
         else
	    set_mat_elem(A,i,j, (float) 0);
   }

/* set all elements of A to 1
 */
void ones_mat(A)
   struct dmat * A;
{
   int i;
   for (i=0; i < A->ntiles; i++)
      *(A->base + i) = (plural float) 1;
   }


/* return a plural random number in the interval [0,1]
 */
#define PF (plural float)
#define P1 8191
#define P2 190843
plural int seed = 42;  
plural float myrand()
{
   seed = (seed * (plural) P1) % (plural) P2;
   return PF (seed & (plural) 1023) / PF 1022;
   }

/* add s * r to each element of A, where median value of
 * r is 1 and max value of r is 4.   (e^4-1 = 53.5982)
 */
void mat_rand_add(s, A, B)
   float s;
   struct dmat *A, *B;
{
   int i;
   plural float r;
   for (i=0; i < A->ntiles; i++) {
      r = PF 4.0 - fp_log(myrand() * PF 53.5982 + PF 1.0);
      if (seed & 8) 
	 *(B->base + i) = *(A->base + i) + r * s;
      else
	 *(B->base + i) = *(A->base + i) - r * s;
      }
   }

void describe_mat(A) 
   struct dmat * A; 
{
   printf("base=%x, rows=%d, cols=%d, tiles=%d\n", 
	  A->base, A->nrows, A->ncols, A->ntiles);
   }

