gen_matrix.c

C++编写的高性能矩阵乘法的Stranssen算法

/*
 * ---------------------------------------------------------------- *
 * Functions to generate a(i,j), b(i,j) and c(i,j). These functions *
 * use global variables l_lambda and l_gammma as parameters.        *
 * ---------------------------------------------------------------- *
 */
#include "matrix.h"
#include "matrix_test.h"

double l_lambda = .5, l_gamma = -.5;

double aij(int i, int j)
{
  double temp;

#if STRAS_RANDOM
  srand48(rand());
  temp = drand48();
#else
  temp = (1.0*(i+1) + l_lambda*(j+1));
#endif

  return temp;
}

double bij(int i, int j)
{
  double temp;

#if STRAS_RANDOM
  srand48(rand());
  temp = drand48();
#else
  temp = (1.0*(i+1) + l_gamma*(j+1));
#endif

  return temp;
}

double cij(int i, int j)
{
  double temp;

#if STRAS_RANDOM
  temp = drand48();
#else
  temp = 1.0*i + l_lambda*j;
#endif

  return temp;
}

double cij_result(char c_transa, char c_transb, int k, int i, int j, double alpha, double beta)
{

  double temp;
  if ((c_transa == 'T') && (c_transb == 'T'))
    temp = (3.0*l_lambda*l_gamma*k*(k+1)*(i+1)
            + 3.0*k*(k+1)*(j+1)
            + 6.0*l_lambda*k*(i+1)*(j+1)
            + l_gamma*k*(2.0*k*k + 3.0*k+1.0));
  else if (c_transa == 'T')
    temp = (1.0*k*(k+1)*(2.0*k+1.0)
            + 3.0*k*(k+1)*(l_lambda*(i+1)+l_gamma*(j+1))
            + 6.0*l_lambda*l_gamma*k*(i+1)*(j+1));
  else if (c_transb == 'T')
    temp = (3.0*l_gamma*k*(k+1)*(i+1)
            + 3.0*l_lambda*k*(k+1)*(j+1)
            + 6.0*k*(i+1)*(j+1)
            + l_lambda*l_gamma*k*(2.0*k*k + 3.0*k+1.0));
  else 
    temp = (3.0*k*(k+1)*(i+1)
            + 3.0*l_lambda*l_gamma*k*(k+1.0)*(j+1)
            + 6.0*l_gamma*k*(i+1)*(j+1)
            + l_lambda*k*(2.0*k*k + 3.0*k+1.0));

  /* To reduce # of divides, above equations were multiplies by GCD = 6, so
   * now need to divide by 6.0.  First part of equation is the original value
   * for C, from cij() above.
   */
  temp = (beta*(1.0*i + l_lambda*j) + alpha*temp/6.0);

  return temp;
}