# Parallel algorithms#

This program will perform a matrix multiplication in parallel. The output will look something like this:

```Matrix A is :
4 9 6
1 9 8

Matrix B is :
4 9
6 1
9 8

Resultant Matrix is :
124 93
111 127
```

## Setup#

The source code for this example can be found here: `matrix_multiplication.cpp`.

To compile this program, go to your HPX build directory (see Building HPX for information on configuring and building HPX) and enter:

```\$ make examples.quickstart.matrix_multiplication
```

To run the program type:

```\$ ./bin/matrix_multiplication
```

or:

```\$ ./bin/matrix_multiplication --n 2 --m 3 --k 2 --s 100 --l 0 --u 10
```

where the first matrix is n x m and the second m x k, s is the seed for creating the random values of the matrices and the range of these values is [l,u]

This should print:

```Matrix A is :
4 9 6
1 9 8

Matrix B is :
4 9
6 1
9 8

Resultant Matrix is :
124 93
111 127
```

Notice that the numbers may be different because of the random initialization of the matrices.

## Walkthrough#

Now that you have compiled and run the code, let’s look at how the code works.

First, `main()` is used to initialize the runtime system and pass the command line arguments to the program. `hpx::init` calls `hpx_main()` after setting up HPX, which is where our program is implemented.

```int main(int argc, char* argv[])
{
using namespace hpx::program_options;
options_description cmdline("usage: " HPX_APPLICATION_STRING " [options]");
// clang-format off
("n",
hpx::program_options::value<std::size_t>()->default_value(2),
"Number of rows of first matrix")
("m",
hpx::program_options::value<std::size_t>()->default_value(3),
"Number of columns of first matrix (equal to the number of rows of "
"second matrix)")
("k",
hpx::program_options::value<std::size_t>()->default_value(2),
"Number of columns of second matrix")
("seed,s",
hpx::program_options::value<unsigned int>(),
"The random number generator seed to use for this run")
("l",
hpx::program_options::value<int>()->default_value(0),
"Lower limit of range of values")
("u",
hpx::program_options::value<int>()->default_value(10),
"Upper limit of range of values");
// clang-format on
hpx::local::init_params init_args;
init_args.desc_cmdline = cmdline;

return hpx::local::init(hpx_main, argc, argv, init_args);
}
```

Proceeding to the `hpx_main()` function, we can see that matrix multiplication can be done very easily.

```int hpx_main(hpx::program_options::variables_map& vm)
{
using element_type = int;

// Define matrix sizes
std::size_t const rowsA = vm["n"].as<std::size_t>();
std::size_t const colsA = vm["m"].as<std::size_t>();
std::size_t const rowsB = colsA;
std::size_t const colsB = vm["k"].as<std::size_t>();
std::size_t const rowsR = rowsA;
std::size_t const colsR = colsB;

// Initialize matrices A and B
std::vector<int> A(rowsA * colsA);
std::vector<int> B(rowsB * colsB);
std::vector<int> R(rowsR * colsR);

// Define seed
unsigned int seed = std::random_device{}();
if (vm.count("seed"))
seed = vm["seed"].as<unsigned int>();

gen.seed(seed);
std::cout << "using seed: " << seed << std::endl;

// Define range of values
int const lower = vm["l"].as<int>();
int const upper = vm["u"].as<int>();

// Matrices have random values in the range [lower, upper]
std::uniform_int_distribution<element_type> dis(lower, upper);
auto generator = std::bind(dis, gen);
hpx::ranges::generate(A, generator);
hpx::ranges::generate(B, generator);

// Perform matrix multiplication
hpx::experimental::for_loop(hpx::execution::par, 0, rowsA, [&](auto i) {
hpx::experimental::for_loop(0, colsB, [&](auto j) {
R[i * colsR + j] = 0;
hpx::experimental::for_loop(0, rowsB, [&](auto k) {
R[i * colsR + j] += A[i * colsA + k] * B[k * colsB + j];
});
});
});

// Print all 3 matrices
print_matrix(A, rowsA, colsA, "A");
print_matrix(B, rowsB, colsB, "B");
print_matrix(R, rowsR, colsR, "R");

return hpx::local::finalize();
}
```

First, the dimensions of the matrices are defined. If they were not given as command-line arguments, their default values are 2 x 3 for the first matrix and 3 x 2 for the second. We use standard vectors to define the matrices to be multiplied as well as the resultant matrix.

To give some random initial values to our matrices, we use std::uniform_int_distribution. Then, `std::bind()` is used along with `hpx::ranges::generate()` to yield two matrices A and B, which contain values in the range of [0, 10] or in the range defined by the user at the command-line arguments. The seed to generate the values can also be defined by the user.

The next step is to perform the matrix multiplication in parallel. This can be done by just using an `hpx::experimental::for_loop` combined with a parallel execution policy `hpx::execution::par` as the outer loop of the multiplication. Note that the execution of `hpx::experimental::for_loop` without specifying an execution policy is equivalent to specifying `hpx::execution::seq` as the execution policy.

Finally, the matrices A, B that are multiplied as well as the resultant matrix R are printed using the following function.

```void print_matrix(std::vector<int> const& M, std::size_t rows, std::size_t cols,
char const* message)
{
std::cout << "\nMatrix " << message << " is:" << std::endl;
for (std::size_t i = 0; i < rows; i++)
{
for (std::size_t j = 0; j < cols; j++)
std::cout << M[i * cols + j] << " ";
std::cout << "\n";
}
}
```