// Level 3 cblas - Workshop 4
 // w4_cblas.cpp

 #include <iostream>
 #include <iomanip>
 #include <cstdlib>
 #include <chrono>
 // add cblas header file
 using namespace std::chrono;

 // indexing function (column major order)
 //
 inline int idx(int r, int c, int n)
 {
     // ... add indexing formula
 }

 // display matrix M, which is stored in column-major order
 //
 void display(const char* str, const float* M, int nr, int nc)
 {
     std::cout << str << std::endl;
     std::cout << std::fixed << std::setprecision(4);
     for (int i = 0; i < nr; i++) {
         for (int j = 0; j < nc; j++)
             std::cout << std::setw(10)
              << // ... access in column-major order;
         std::cout << std::endl;
     }
     std::cout << std::endl;
 }

 // report system time
 //
 void reportTime(const char* msg, steady_clock::duration span) {
     auto ms = duration_cast<milliseconds>(span);
     std::cout << msg << " - took - " <<
      ms.count() << " millisecs" << std::endl;
 }


 // matrix multiply
 //
 void sgemm(const float* A, const float* B, float* C, int n) {
     steady_clock::time_point ts, te;

     // level 3 calculation: C = alpha * A * B + beta * C

     // add any preliminaries

     ts = steady_clock::now();
     // ... add call to cblas sgemm
     te = steady_clock::now();
     reportTime("matrix-matrix multiplication", te - ts);
     
 }

 int main(int argc, char* argv[]) {
     if (argc != 2) {
         std::cerr << argv[0] << ": invalid number of arguments\n"; 
         std::cerr << "Usage: " << argv[0] << "  size_of_matrices\n"; 
         return 1;
     }
     int n = std::atoi(argv[1]); // no of rows/columns in A, B, C 

     // allocate host memory
     float* h_A = new float[n * n];
     float* h_B = new float[n * n];
     float* h_C = new float[n * n];

     // populate host matrices a and b
     for (int i = 0, kk = 0; i < n; i++)
         for (int j = 0; j < n; j++, kk++)
             h_A[kk] = h_B[kk] = (float)kk;

     // C = A * B
     sgemm(h_A, h_B, h_C, n);

     // display results
     if (n <= 5) {
         display("A :", h_A, n, n);
         display("B :", h_B, n, n);
         display("C = A B :", h_C, n, n);
     }

     // deallocate host memory
     delete [] h_A;
     delete [] h_B;
     delete [] h_C;
 }

# Makefile for w4
#
GCC_VERSION = 7.2.0
PREFIX = /usr/local/gcc/${GCC_VERSION}/bin/
CC = ${PREFIX}gcc
CPP = ${PREFIX}g++

w4_cblas: w4_cblas.o
        $(CPP) -ow4_cblas w4_cblas.o -lgslcblas 

w4_cblas.o: w4_cblas.cpp
        $(CPP) -c -O2 -std=c++17 w4_cblas.cpp

clean:
        rm *.o
 

 > make

 > w4_cblas 4

 matrix-matrix multiplication - took 0 secs 
 A :
     0.0000    4.0000    8.0000   12.0000
     1.0000    5.0000    9.0000   13.0000
     2.0000    6.0000   10.0000   14.0000
     3.0000    7.0000   11.0000   15.0000

 B :
     0.0000    4.0000    8.0000   12.0000
     1.0000    5.0000    9.0000   13.0000
     2.0000    6.0000   10.0000   14.0000
     3.0000    7.0000   11.0000   15.0000

 C = A B :
    56.0000  152.0000  248.0000  344.0000
    62.0000  174.0000  286.0000  398.0000
    68.0000  196.0000  324.0000  452.0000
    74.0000  218.0000  362.0000  506.0000

 // Level 3 cuBLAS - Workshop 4
 // w4_cublas.cu

 // ...

 // matrix multiply
 //
 void sgemm(const float* h_A, const float* h_B, float* h_C, int n) {
     steady_clock::time_point ts, te;

     // level 3 calculation: C = alpha * A * B + beta * C
     
     ts = steady_clock::now();
     // ... allocate memory on the device
     te = steady_clock::now();
     reportTime("allocation of device memory for matrices d_A, d_B and d_C",
      te - ts);

     // ... create cuBLAS context

     ts = steady_clock::now();
     // ... copy host matrices to the device
     te = steady_clock::now();
     reportTime("copying of matrices h_A and h_B to device memory", te - ts);

     ts = steady_clock::now();
     // ... calculate matrix-matrix product
     te = steady_clock::now();
     reportTime("matrix-matrix multiplication", te - ts);

     // ... copy result matrix from the device to the host
     te = steady_clock::now();
     reportTime("copying of matrix d_C from device", te - ts);

     // ... destroy cuBLAS context

     ts = steady_clock::now();
     // ... deallocate device memory
     te = steady_clock::now();
     reportTime("deallocation of device memory for matrices A, B and C",
      te - ts);
     
 }