1 /*------------------------------------------------------------------------------------------------*
2  * Copyright (C) by the DBCSR developers group - All rights reserved                              *
3  * This file is part of the DBCSR library.                                                        *
4  *                                                                                                *
5  * For information on the license, see the LICENSE file.                                          *
6  * For further information please visit https://dbcsr.cp2k.org                                    *
7  * SPDX-License-Identifier: GPL-2.0+                                                              *
8  *------------------------------------------------------------------------------------------------*/
9 
10 #include <vector>
11 #include <iostream>
12 #include <algorithm>
13 #include <cstdlib>
14 #include <cstdio>
15 #include <cstdint>
16 #include <random>
17 
18 #include <mpi.h>
19 
20 #include <dbcsr.h>
21 
22 
23 // Random distribution by using round-robin assignment
24 // of blocks to processors
random_dist(int dist_size,int nbins)25 std::vector<int> random_dist(int dist_size, int nbins)
26 {
27     std::vector<int> dist(dist_size);
28 
29     for(int i=0; i < dist_size; i++)
30         dist[i] = (nbins-i+1) % nbins;
31 
32     return dist;
33 }
34 
35 
main(int argc,char * argv[])36 int main(int argc, char* argv[])
37 {
38     MPI_Init(&argc, &argv);
39 
40     int mpi_size, mpi_rank;
41     MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
42     MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
43 
44     // Make 2D grid
45     int dims[2] = {0};
46     MPI_Dims_create(mpi_size, 2, dims);
47     int periods[2] = {1};
48     int reorder = 0;
49     MPI_Comm group;
50     MPI_Cart_create(MPI_COMM_WORLD, 2, dims, periods, reorder, &group);
51 
52     int coord[2];
53     MPI_Cart_coords(group, mpi_rank, 2, coord);
54 
55     std::cout
56         << "I'm processor " << mpi_rank
57         << " over " << mpi_size << " proc"
58         << ", (" << coord[0] << ", " << coord[1] << ") in the 2D grid"
59         << std::endl;
60 
61     c_dbcsr_init_lib(MPI_COMM_WORLD, nullptr);
62 
63     // Total number of blocks
64     int nblkrows_total = 4;
65     int nblkcols_total = 4;
66 
67     // Block sizes
68     std::vector<int> row_blk_sizes(nblkrows_total, 2), col_blk_sizes(nblkcols_total, 2);
69 
70     auto row_dist = random_dist(nblkrows_total, dims[0]);
71     auto col_dist = random_dist(nblkcols_total, dims[1]);
72 
73     void* dist = nullptr;
74 
75     c_dbcsr_distribution_new(&dist, group,
76         row_dist.data(), row_dist.size(),
77         col_dist.data(), col_dist.size());
78 
79     // Fill all blocks, i.e. dense matrices
80     auto fill_matrix = [&](void*& matrix)
81     {
82         int max_row_size = *std::max_element(row_blk_sizes.begin(),row_blk_sizes.end());
83         int max_col_size = *std::max_element(col_blk_sizes.begin(),col_blk_sizes.end());
84         int max_nze = max_row_size * max_col_size;
85 
86         std::vector<double> block;
87         block.reserve(max_nze);
88 
89         for(int i = 0; i < nblkrows_total; i++)
90         {
91             for(int j = 0; j < nblkcols_total; j++)
92             {
93                 int blk_proc = -1;
94                 c_dbcsr_get_stored_coordinates(matrix, i, j, &blk_proc);
95 
96                 if(blk_proc == mpi_rank)
97                 {
98                     block.resize(row_blk_sizes[i] * col_blk_sizes[j]);
99                     std::generate(block.begin(), block.end(), [&](){ return static_cast<double>(std::rand())/RAND_MAX; });
100                     c_dbcsr_put_block_d(matrix, i, j, block.data(), block.size());
101                 }
102             }
103         }
104     };
105 
106     // create and fill matrix a
107     void* matrix_a = nullptr;
108     c_dbcsr_create_new_d(&matrix_a, "this is my matrix a", dist, 'N',
109         row_blk_sizes.data(), row_blk_sizes.size(),
110         col_blk_sizes.data(), col_blk_sizes.size());
111     fill_matrix(matrix_a);
112     c_dbcsr_finalize(matrix_a);
113 
114     // create and fill matrix b
115     void* matrix_b = nullptr;
116     c_dbcsr_create_new_d(&matrix_b, "this is my matrix b", dist, 'N',
117         row_blk_sizes.data(), row_blk_sizes.size(),
118         col_blk_sizes.data(), col_blk_sizes.size());
119     fill_matrix(matrix_b);
120     c_dbcsr_finalize(matrix_b);
121 
122     // create matrix c, empty
123     void* matrix_c = nullptr;
124     c_dbcsr_create_new_d(&matrix_c, "matrix c", dist, 'N',
125         row_blk_sizes.data(), row_blk_sizes.size(),
126         col_blk_sizes.data(), col_blk_sizes.size());
127     c_dbcsr_finalize(matrix_c);
128 
129     // multiply the matrices
130     c_dbcsr_multiply_d('N', 'N', 1.0, &matrix_a, &matrix_b, 0.0, &matrix_c, nullptr);
131 
132     c_dbcsr_print(matrix_a);
133     c_dbcsr_print(matrix_b);
134     c_dbcsr_print(matrix_c);
135 
136     c_dbcsr_release(&matrix_a);
137     c_dbcsr_release(&matrix_b);
138     c_dbcsr_release(&matrix_c);
139 
140     c_dbcsr_distribution_release(&dist);
141 
142 
143     MPI_Comm_free(&group);
144 
145     c_dbcsr_finalize_lib();
146 
147     MPI_Finalize();
148 
149     return 0;
150 }
151