Functions as parameters

Content

So far, the initialization of a matrix was lacking in flexibility:

template<typename Matrix>
void init_matrix(Matrix& A) {
   for (std::size_t i = 0; i < A.m; ++i) {
      for (std::size_t j = 0; j < A.n; ++j) {
	 A(i, j) = j * A.m + i + 1;
      }
   }
}

In dependence of the test or a benchmark it could be preferable to fill the matrix with pseudo-random numbers, to select other parameters for the distribution, to chose another distribution, or to seed the pseudo-random number sequence with a fixed value to assure reproducibility.

The computation of an individual value for an element could be delegated to a function that is passed as parameter to init_matrix. This function would get the indices i and j along with the dimensions m and n, and would return \(A_{i,j}\) as value. The following template function would support \(A_{i,j} \leftarrow j*m+i+1\):

template<typename T>
T init_value(std::size_t i, std::size_t j, std::size_t m, std::size_t n) {
   return j*m + i + 1;
}

Function types can be specified by the return type, followed by parentheses that enclose the parameter types. Hence, the general type of an initialization function taking dimensions and indices would be T(std::size_t, std::size_t, std::size_t, std::size_t). Functions of this type would expect four parameters of type std::size_t and return a value of type T. Within a declarator, a name is to be added, e.g. T init_value(std::size_t, std::size_t, std::size_t, std::size_t). Following variant of init_matrix expects such a function parameter:

template<typename Matrix>
void init_matrix(Matrix& A,
      typename Matrix::Element
	 init_value(std::size_t, std::size_t, std::size_t, std::size_t)) {
   for (std::size_t i = 0; i < A.m; ++i) {
      for (std::size_t j = 0; j < A.n; ++j) {
	 A(i, j) = init_value(i, j, A.m, A.n);
      }
   }
}

Once we have provided a more general solution, it is worthwile to make it more cache-friendly:

template<typename Matrix>
void init_matrix(Matrix& A,
      typename Matrix::Element
	 init_value(std::size_t, std::size_t, std::size_t, std::size_t)) {
   if (A.incRow > A.incCol) {
      for (std::size_t i = 0; i < A.m; ++i) {
	 for (std::size_t j = 0; j < A.n; ++j) {
	    A(i, j) = init_value(i, j, A.m, A.n);
	 }
      }
   } else {
      for (std::size_t j = 0; j < A.n; ++j) {
	 for (std::size_t i = 0; i < A.m; ++i) {
	    A(i, j) = init_value(i, j, A.m, A.n);
	 }
      }
   }
}

Note that init_value is not the function itself but implicitly a pointer to a function. This is the reason why the parameter could have also declared as T (*init_value)(std::size_t, std::size_t, std::size_t, std::size_t) but the derefering operator * for functions can be omitted since the event of ANSI C.

The initialization function could then be used as follows:

GeMatrix<double> A(3, 7, StorageOrder::ColMajor);
init_matrix(A, init_value<double>);

Exercise

Sources

The matrix class has been moved into the header file gematrix.hpp and is now named GeMatrix for general matrix. The test program provides init_matrix and also two variants of init_value, i.e. unique_init_value (as above) and random_init_value which is based upon a pseudo-random-generator.

#ifndef HPC_GEMATRIX_HPP
#define HPC_GEMATRIX_HPP

#include <cassert> /* needed for assert */
#include <cstddef> /* needed for std::size_t and std::ptrdiff_t */

enum class StorageOrder {ColMajor, RowMajor};

template<typename T>
struct GeMatrix {
   const std::size_t m; /* number of rows */
   const std::size_t n; /* number of columns */
   const std::size_t incRow;
   const std::size_t incCol;
   T* const data;
   using Element = T;

   GeMatrix(std::size_t m, std::size_t n, StorageOrder order) :
	 m(m), n(n),
	 incRow(order == StorageOrder::ColMajor? 1: n),
	 incCol(order == StorageOrder::RowMajor? 1: m),
	 data(new T[m*n]) {
   }

   ~GeMatrix() {
      delete[] data;
   }

   const T& operator()(std::size_t i, std::size_t j) const {
      assert(i < m && j < n);
      return data[i*incRow + j*incCol];
   }
   
   T& operator()(std::size_t i, std::size_t j) {
      assert(i < m && j < n);
      return data[i*incRow + j*incCol];
   }
};

template<typename T>
struct GeMatrixView {
   const std::size_t m; /* number of rows */
   const std::size_t n; /* number of columns */
   const std::size_t incRow;
   const std::size_t incCol;
   T* const data;
   using Element = T;

   GeMatrixView(std::size_t m, std::size_t n,
	    T* data,
	    std::size_t incRow, std::size_t incCol) :
	 m(m), n(n), 
	 incRow(incRow), incCol(incCol),
	 data(data) {
   }

   const T& operator()(std::size_t i, std::size_t j) const {
      assert(i < m && j < n);
      return data[i*incRow + j*incCol];
   }
   
   T& operator()(std::size_t i, std::size_t j) {
      assert(i < m && j < n);
      return data[i*incRow + j*incCol];
   }
};

#endif
#include <cstddef>
#include <random>
#include <printf.hpp>
#include "gematrix.hpp"

template<typename Matrix>
void init_matrix(Matrix& A,
      typename Matrix::Element
	 init_value(std::size_t, std::size_t, std::size_t, std::size_t)) {
   if (A.incRow > A.incCol) {
      for (std::size_t i = 0; i < A.m; ++i) {
	 for (std::size_t j = 0; j < A.n; ++j) {
	    A(i, j) = init_value(i, j, A.m, A.n);
	 }
      }
   } else {
      for (std::size_t j = 0; j < A.n; ++j) {
	 for (std::size_t i = 0; i < A.m; ++i) {
	    A(i, j) = init_value(i, j, A.m, A.n);
	 }
      }
   }
}

template<typename T>
T init_value(std::size_t i, std::size_t j, std::size_t m, std::size_t n) {
   return j*m + i + 1;
}

template<typename Matrix>
void print_matrix(const Matrix& A) {
   for (std::size_t i = 0; i < A.m; ++i) {
      fmt::printf("  ");
      for (std::size_t j = 0; j < A.n; ++j) {
	 fmt::printf(" %7.2lf", A(i, j));
      }
      fmt::printf("\n");
   }
}

int main() {
   GeMatrix<double> A(3, 7, StorageOrder::ColMajor);
   init_matrix(A, init_value<double>);
   fmt::printf("A:\n"); print_matrix(A);
}
theon$ g++ -Wall -std=gnu++11 -o test_initmatrix test_initmatrix.cpp
theon$ ./test_initmatrix
A:
      1.00    4.00    7.00   10.00   13.00   16.00   19.00
      2.00    5.00    8.00   11.00   14.00   17.00   20.00
      3.00    6.00    9.00   12.00   15.00   18.00   21.00
theon$