/*****************************************************************************
 *
 * file:   matmul.sac
 *
 * description:
 *
 *   SAC demo program.
 *
 *   This SAC demo program implements matrix multiplication for two matrices
 *   A and B of double precision floating point numbers.
 *
 *   The sizes of A and B are preset as follows:
 *     A is of size SIZE1 x SIZE2, and
 *     B is of size SIZE2 x SIZE3.
 *   These parameters may be set at compile time.
 *
 *   The pragma given in the matmul function forces the compiler to introduce
 *   a blocking of squares containing 32 by 32 elements. Though these numbers
 *   seem to work well on many architectures, it might be worthwhile to vary
 *   them for getting better execution times.
 *
 *   On SUN-SOLARIS systems cc should be used instead of gcc for achieving
 *   much better runtime performance. To do so, specify "-target suncc" as
 *   compiler option!
 *
 *****************************************************************************/

#ifndef SIZE1
#define SIZE1 1000
#endif

#ifndef SIZE2
#define SIZE2 1000
#endif

#ifndef SIZE3
#define SIZE3 1000
#endif


use Array:all;
use SimplePrint:all;

inline double[.,.] transpose( double[.,.] A)
{
  res = with( . <= [i,j] <= .) {
       } genarray( [shape(A)[[1]], shape(A)[[0]]], A[[j,i]], 0d);
  return( res);
}

double[.,.] matmul( double[.,.] A, double[.,.] B)
{

  BT = transpose( B);

  /*  C = #pragma wlcomp BvL0( [32,32], Default) */
  C = with( . <= [i,j] <= .)
      genarray( [ shape(A)[[0]], shape(B)[[1]] ], sum( A[[i]] * BT[[j]]));
  return( C);
}


int main()
{

  A = genarray( [SIZE1,SIZE2], 1.0);
  B = genarray( [SIZE2,SIZE3], 1.0);

  C = matmul( A, B);

#if 0
  printf( "C[0,0] = %f\n", C[0,0]);

  return(0);
#else
  a = print( C[0,0]);
  return( a);
#endif
}