tc_test_numerics-A100-tf32.cu

/*
 * Copyright (c) 2020, Massimiliano Fasi and Mantas Mikaitis
 *
 * This program is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License as published by the Free Software
 * Foundation, version 2.
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
 * details.
 *
 *  You should have received a copy of the GNU General Public License along with
 *  this program. If not, see <http://www.gnu.org/licenses/>.
 */

#include <assert.h>
#include <unistd.h>
#include <cstdint>
#include <chrono>
#include <iostream>
#include <mma.h>
#include <iomanip>

using namespace nvcuda;

/*******************
 * Debug functions *
 *******************/
/* Print the elements of the m x n matrix A. The elements are assumed to be
   stored by columns if `bycols` is `true` and by rows if `bycols` is false. */
template <typename floattype>
void print_matrix (float *a,
                   size_t m, size_t n,
                   bool bycols) {
  int i, j;
  if (bycols) {
    for (i=0; i<m; i++) {
      for (j=0; j<n; j++)
        std::cout << a[j*n+i] << " ";
      std::cout << std::endl;
    }
    std::cout << std::endl;
  } else {
    for (i=0; i<m; i++ ) {
      for (j=0; j<n; j++)
        std::cout << a[i*m+j] << " ";
      std::cout  << std::endl;
    }
    std::cout << std::endl;
   }
}


/****************************************************
 * Memory management and wmma::mma_sync() interface *
 ****************************************************/

/* Set the entries of host arrays to zero. */
void host_reset(float *a, float *b, float *c) {
  memset(a, 0, 16*16*sizeof(float));
  memset(b, 0, 16*16*sizeof(float));
  memset(c, 0, 16*16*sizeof(float));
}

/* Compute C += A*B, where A, B, and C are 16x16x16 matrices.
   The matrix C is initialized to 0 when `init` is true. */
__global__ void wmma_ker(float *a, float *b, float *c, bool init) {

  // Declare fragments.
  wmma::fragment<wmma::matrix_a, 16, 16, 8, wmma::precision::tf32,
    wmma::row_major> a_fragment;
  wmma::fragment<wmma::matrix_b, 16, 16, 8, wmma::precision::tf32,
    wmma::col_major> b_fragment;
  wmma::fragment<wmma::accumulator, 16, 16, 8, float> c_fragment;

  // Load input matrices and initialize output (if required).
  wmma::load_matrix_sync(a_fragment, a, 16);
  wmma::load_matrix_sync(b_fragment, b, 16);
  if (init)
    wmma::fill_fragment(c_fragment, 0.0f);
  else
    wmma::load_matrix_sync(c_fragment, c, 16, wmma::mem_col_major);

  // Multiply
  wmma::mma_sync(c_fragment, a_fragment, b_fragment, c_fragment);

  // Store the output
  wmma::store_matrix_sync(c, c_fragment, 16, wmma::mem_col_major);
}

/* Copy data from host to device, perform the operation, and copy result back to
   host. */
void wmma_init_run (float *h_a, float *h_b, float *h_c,
                    float *d_a, float *d_b, float *d_c,
                    bool init) {

  // Copy input from host to device.
  cudaMemcpy(d_a, h_a, 16*16*sizeof(float), cudaMemcpyHostToDevice);
  cudaMemcpy(d_b, h_b, 16*16*sizeof(float), cudaMemcpyHostToDevice);
  cudaMemcpy(d_c, h_c, 16*16*sizeof(float), cudaMemcpyHostToDevice);

  // Perform matrix multiplication.
  wmma_ker<<<1,32>>>(d_a, d_b, d_c, init);

  // Copy result from device to host.
  cudaMemcpy(h_c, d_c, 16*16*sizeof(float), cudaMemcpyDeviceToHost);
}


/**********************
 * Printing functions *
 **********************/
void printheader(FILE *outfile, const char *string) {
  fprintf(outfile,
          "+--------------------------------------------------------------+\n");
  fprintf(outfile, "| %-60s |\n", string);
  fprintf(outfile,
          "+--------------------------------------------------------------+\n");
}
void printitem(FILE *outfile, const char *string) {
  fprintf(outfile, "  | %-49s", string);
}

void printpass(FILE *outfile, bool status) {
  if (status)
    fprintf(outfile, " [PASS] |\n");
  else
    fprintf(outfile, " [FAIL] |\n");
}
void printfooter(FILE *outfile) {
  fprintf(outfile,
          "  +----------------------------------------------------------+\n\n");
}


/***************
 * EXPERIMENTS *
 ***************/
int main(int argc, char** argv){

  // Declare pointers and allocate memory.
  float *h_a, *h_b, *h_c, *d_a, *d_b, *d_c,
    minsubnormaltf32 = ldexp(1., -136), // smallest subnormal tf32
    minsubnormal32 = ldexp(1., -149), // smallest subnormal binary32
    belowone = nextafterf(1., 0.) ,   // largest float smaller than 1.0
    gapbelowone = 1. - belowone,
    aboveone = nextafterf(1., 2.),    // smallest float larger than 1.0
    belowtwo = 2. - ldexp(1., -23);   // largest float smaller than 2.0

  assert(belowone == 1. - ldexp(1., -24));
  assert(aboveone == 1. + ldexp(1., -23));

  h_a = new float[16*16];
  h_b = new float[16*16];
  h_c = new float[16*16];
 
  cudaMalloc(&d_a, 16*16*sizeof(float));
  cudaMalloc(&d_b, 16*16*sizeof(float));
  cudaMalloc(&d_c, 16*16*sizeof(float));

  FILE *outfile = stdout;
  bool pass;

  printheader(outfile, "A. Support for subnormal numbers");// ;

   printitem(outfile, "*) TF32 subnormals in input (binary32 mode)");
  host_reset(h_a, h_b, h_c);
  h_a[0] = minsubnormaltf32;
  h_b[0] = ldexp(1, 10);
  wmma_init_run(h_a, h_b, h_c, d_a, d_b, d_c, false);
  printpass(outfile, h_c[0] == ldexp(1., -126));

  printitem(outfile, "*) Binary32 subnormals in input");
  host_reset(h_a, h_b, h_c);
  h_c[0] = minsubnormal32;
  wmma_init_run (h_a, h_b, h_c, d_a, d_b, d_c, false);
  printpass(outfile, h_c[0] == minsubnormal32);

   printitem(outfile, "*) TF32/Binary32 subnormals in output (binary32 mode)");
  host_reset(h_a, h_b, h_c);
  h_a[0] = ldexp(1., -126);
  h_b[0] = ldexp(1., -1);
  wmma_init_run (h_a, h_b, h_c, d_a, d_b, d_c, false);
  pass = h_c[0] == ldexp(1, -127);
  h_a[0] = ldexp(1., -126);
  h_b[0] = 1.0;
  h_c[0] = ldexp(-1., -127);
  wmma_init_run (h_a, h_b, h_c, d_a, d_b, d_c, false);
  pass = pass && (h_c[0] == ldexp(1, -127));
  printpass(outfile, pass);

  printfooter(outfile);

  printheader(outfile, "B. Accuracy of the dot products ");// ;

  printitem(outfile, "*) Products are computed exactly ");
  host_reset(h_a, h_b, h_c);
  h_a[0] = belowone;
  h_b[0] = belowone;
  wmma_init_run (h_a, h_b, h_c, d_a, d_b, d_c, false);
  pass = (h_c[0] == 1-ldexp(1, -10) + ldexp(1, -22));
  size_t i,j;
  for (i=0; i<4; i++) {
    h_a[i] = belowone;
    h_b[i] = belowone;
  }
  h_c[0] = 0.0;
  wmma_init_run (h_a, h_b, h_c, d_a, d_b, d_c, false);
  pass = (h_c[0] == (4 * (1-ldexp(1, -10) + ldexp(1, -22))));
  printpass(outfile, pass);

  printitem(outfile, "*) Products are accumulated in binary32 ");
  host_reset(h_a, h_b, h_c);
  pass = true;
  for (i=0; i<4; i++) {
    h_a[i] = 0.5;
    h_b[i] = ldexp(1, -24);
  }
  h_c[0] = 1.;
  wmma_init_run (h_a, h_b, h_c, d_a, d_b, d_c, false);
  pass = pass && h_c[0] == 1;
  printpass(outfile, pass);

  printitem(outfile, "*) Sum starts from largest element");
  host_reset(h_a, h_b, h_c);
  pass = true;
  for (i=0; i<4; i++) {
    h_a[i] = 0.5;
    h_b[i] = ldexp(1, -24);
  }
  for (j=0; j<4; j++) {
    h_c[0] = ldexp(1, -24);
    if (j>0)
      h_a[j-1] = 0.5;
    h_b[j-1] = ldexp(1, -24);
    h_a[j] = 1.0;
    h_b[j] = 1.0;
    wmma_init_run (h_a, h_b, h_c, d_a, d_b, d_c, false);
    pass = pass && h_c[0] == 1;
  }
  printpass(outfile, pass);

  printfooter(outfile);

  printheader(outfile, "C. Rounding modes in tensor core computations ");

  printitem(outfile, "*) Round-down for positive values ");
  host_reset(h_a, h_b, h_c);
  for (i=0; i<4; i++) {
    h_a[i] = 1.0;
  }
  h_b[0] = 2.;
  h_b[1] = ldexp(1., -23) + ldexp(1., -24);
  wmma_init_run (h_a, h_b, h_c, d_a, d_b, d_c, false);
  printpass(outfile, h_c[0] == 2.);

  printitem(outfile, "*) Round-up for negative values ");
  host_reset(h_a, h_b, h_c);
  for (i=0; i<4; i++) {
    h_a[i] = 1.0;
  }
  h_b[0] = -2.;
  h_b[1] = -ldexp(1., -23) - ldexp(1., -24);
  wmma_init_run (h_a, h_b, h_c, d_a, d_b, d_c, false);
  printpass(outfile, h_c[0] == -2.);

  printfooter(outfile);

  printheader(outfile, "D. Features of the accumulator");

  printitem(outfile, "1) Extra bits in the significand alignment");
  host_reset(h_a, h_b, h_c);
  h_a[0] = 1.0;
  h_b[0] = 1.0;
  h_c[0] = -belowone;
  wmma_init_run (h_a, h_b, h_c, d_a, d_b, d_c, false);
  assert(1 - belowone == ldexp(1., -24));
  assert(gapbelowone == ldexp(1., -24));
  printpass(outfile, h_c[0] == ldexp(1., -24));

  printitem(outfile, "2) Normalization in addition");
  host_reset(h_a, h_b, h_c);
  for (i=0; i<4; i++) {
    h_a[i] = 1.0;
    h_b[i] = ldexp(1, -24);
  }
  h_c[0] = 1. - ldexp(1., -24);
  wmma_init_run (h_a, h_b, h_c, d_a, d_b, d_c, false);
  pass = h_c[0] == 1. + ldexp(1., -23);
  printpass(outfile, pass);

  printitem(outfile, "3) Normalization in subtraction");
  host_reset(h_a, h_b, h_c);
  h_a[0] = 1.0;
  h_a[1] = 1.0;
  h_b[0] = 1.0;
  h_b[1] = -ldexp(1.,-24);
  h_c[0] = -1. + ldexp(1.,-24);
  wmma_init_run (h_a, h_b, h_c, d_a, d_b, d_c, false);
  pass = pass && h_c[0] == 0.0;
  printpass(outfile, pass);

  printitem(outfile, "4) Extra bits for carry out");
  host_reset(h_a, h_b, h_c);
  for (i=0; i<4; i++) {
    h_a[i] = 1.0;
    h_b[i] = 1.0;
  }
  pass = true;
  for (i=0; i<4; i++) {
    if (i>0)
      h_b[i-1] = 1.0;
    h_b[i] = ldexp(1.,-23);
    h_c[0] = 1. + ldexp(1.,-22) + ldexp(1.,-23);
    wmma_init_run (h_a, h_b, h_c, d_a, d_b, d_c, false);
    pass = pass && h_c[0] == 4. + ldexp(1., -21);
  }

  // Test for the third bit
  host_reset(h_a, h_b, h_c);
  for (i=0; i<4; i++) {
    h_a[i] = 1.0;
    h_b[i] = 1.0;
  }
  pass = true;
  h_b[0] = 1.0;
  h_b[1] = 1.5;
  h_b[2] = 1.75;
  h_b[3] = 1.875;
  h_c[0] = 1.875;
  wmma_init_run (h_a, h_b, h_c, d_a, d_b, d_c, false);
  pass = pass && h_c[0] == 8.;

  // Round-down in normalization of positive values.
  host_reset(h_a, h_b, h_c);
  for (i=0; i<4; i++) {
    h_a[i] = 1.0;
    h_b[i] = 1.0;
  }
  h_c[0] = 1. + ldexp(1., -22) + ldexp(1., -23);
  wmma_init_run (h_a, h_b, h_c, d_a, d_b, d_c, false);
  pass = pass && (h_c[0] == 5.);

  // Round-up in normalization of negative values.
  host_reset(h_a, h_b, h_c);
  for (i=0; i<4; i++) {
    h_a[i] = 1.0;
    h_b[i] = -1;
  }
  h_c[0] = -1. - ldexp(1., -22) - ldexp(1., -23);
  wmma_init_run (h_a, h_b, h_c, d_a, d_b, d_c, false);
  pass = pass && (h_c[0] == -5.);

  printpass(outfile, pass);

  printitem(outfile, "5) Monotonicity of dot product");
  host_reset(h_a, h_b, h_c);
  for (i=0; i<3; i++) {
    h_a[i] = 0.5;
    h_b[i] = ldexp(1., -24);
  }
  for (i=3; i<4; i++) {
    h_a[i] = 0.5;
    h_b[i] = ldexp(1., -23) + ldexp(1., -24);
  }
  h_c[0] = 1. - ldexp(1., -24);
  wmma_init_run (h_a, h_b, h_c, d_a, d_b, d_c, false);
  float partial = h_c[0];
  h_c[0] = 1.0;
  wmma_init_run (h_a, h_b, h_c, d_a, d_b, d_c, false);
  printpass(outfile, h_c[0] < partial);

  printfooter(outfile);

  // Free dynamically allocated memory.
  //  free(h_a);
  //  free(h_b);
  free(h_c);
  cudaFree(d_a);
  cudaFree(d_b);
  cudaFree(d_c);
}