Preparations

Content

Introduction

For this session we have a new version of our lecture library which is to be found on our hosts at the public directory /home/numerik/pub/hpc/ws18/session15/hpc. You can use this library on the computers at our pool in E.44 and at Theon by adding the compiler option -I/home/numerik/pub/hpc/ws18/session15. Please do not forget to add the option -std=c++17 as support of C++17 is required. As the library is header-only this is all you need.

This version of the lecture library comes already with some extensions. Take for example the iterators for vectors and matrices. They allow to access all elements without worrying whether the matrix is stored row- or column-wise:

#include <printf.hpp>
#include <hpc/matvec/gematrix.hpp>
#include <hpc/matvec/iterators.hpp>
#include <hpc/matvec/print.hpp>

int main() {
   using namespace hpc::matvec;
   DenseVector<double> x(10);
   for (auto [i, xi]: x) {
      xi = i + 1;
   }
   fmt::printf("x:\n"); print(x);

   fmt::printf("A:\n");
   GeMatrix<double> A(7, 8);
   for (auto [i, j, Aij]: A) {
      Aij = j * A.numRows() + i + 1;
   }
   print(A);
}
theon$ g++ -O3 -std=c++17 -o init_matvec -I/home/numerik/pub/hpc/ws18/session15 init_matvec.cpp
theon$ ./init_matvec
x:
   1.00   2.00   3.00   4.00   5.00   6.00   7.00   8.00   9.00  10.00
A:
     1.00     8.00    15.00    22.00    29.00    36.00    43.00    50.00
     2.00     9.00    16.00    23.00    30.00    37.00    44.00    51.00
     3.00    10.00    17.00    24.00    31.00    38.00    45.00    52.00
     4.00    11.00    18.00    25.00    32.00    39.00    46.00    53.00
     5.00    12.00    19.00    26.00    33.00    40.00    47.00    54.00
     6.00    13.00    20.00    27.00    34.00    41.00    48.00    55.00
     7.00    14.00    21.00    28.00    35.00    42.00    49.00    56.00
theon$ 

This interesting variant of a for loop is possible thanks to the range-based for loops which came with C++11 which work on base of iterators and to the introduction of so-called structured bindings in C++17 which provides an elegant approach to let a function return tuples which can be assigned to multiple auto declared variables.

Exercise

In the attached source code, randomInit does not use yet the apply template function. Instead it decides itself how to run through a matrix. Simplify this function by using the iterators from <hpc/matvec/iterators.hpp>. You can compile this program using following command:

theon$ g++ -O3 -std=c++17 -o random_init1 -I/home/numerik/pub/hpc/ws18/session15 random_init1.cpp
theon$ ./random_init1
    18.48   -73.57   -62.91    92.17    58.05    47.30    47.19    62.48
    29.85    56.81   -94.93   -12.35   -77.04   -94.86    -6.88   -33.12
   -15.89   -75.60   -12.59   -32.40    37.45   -97.21   -37.84    84.18
   -52.74   -63.05   -27.98    57.85   -92.37    13.37   -24.76   -62.48
    -5.13    31.35    -8.28    62.95   -43.17    54.05    61.95   -92.05
   -23.81   -66.08   -33.04    93.21   -12.73   -53.19   -10.25   -87.29
   -45.72    67.48    79.02    27.15    91.88   -59.00    78.56    48.90
theon$ 

On our computers in E.44 you will need a newer g++ version named g++-7.2:

heim$ g++-7.2 -O3 -std=c++17 -o random_init1 -I/home/numerik/pub/hpc/ws18/session15 random_init1.cpp
heim$ ./random_init1
   -86.91   -70.92    56.54   -28.64    68.26   -68.98   -47.15    40.06
   -75.42    -3.44   -35.81   -16.99   -29.42    65.44    94.91   -67.03
    90.86     7.27    69.03   -94.92    -2.09   -97.08    33.94   -67.51
    -8.85    69.06    71.07   -19.22   -75.32   -93.85   -97.73    73.56
   -29.47    63.58    12.41   -93.35    44.03   -43.22    11.24   -15.29
   -29.14    20.59    92.87    48.07    89.08    42.70   -96.33   -19.92
   -22.99    21.64    39.79    59.28    55.86   -79.37    98.47   -82.21
heim$ 

Source code

#include <random>
#include <hpc/matvec/gematrix.hpp>
#include <hpc/matvec/print.hpp>
#include <hpc/matvec/traits.hpp>

using namespace hpc;
using namespace hpc::matvec;

template <
   template<typename> class MatrixA, typename T,
   Require< Ge<MatrixA<T>> > = true
>
void randomInit(MatrixA<T>& A) {
   std::random_device random;
   std::mt19937 mt{random()};
   std::uniform_real_distribution<T> uniform(-100,100);
   
   if (A.incRow() < A.incCol()) {
      for (std::size_t j = 0; j < A.numCols(); ++j) {
         for (std::size_t i = 0; i < A.numRows(); ++i) {
	    A(i, j) = uniform(mt);
         }
      }
   } else {
      for (std::size_t i = 0; i < A.numRows(); ++i) {
	 for (std::size_t j = 0; j < A.numCols(); ++j) {
	    A(i, j) = uniform(mt);
         }
      }
   }
}

int main() {
   GeMatrix<double> A(7, 8);
   randomInit(A);
   print(A);
}