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
36 // DBCSR example 3
37 // This example shows how to multiply two DBCSR matrices
main(int argc,char * argv[])38 int main(int argc, char* argv[])
39 {
40 // initialize MPI
41 MPI_Init(&argc, &argv);
42
43 // setup the mpi environment
44 int mpi_size, mpi_rank;
45 MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
46 MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
47
48 // make 2D grid
49 int dims[2] = {0};
50 MPI_Dims_create(mpi_size, 2, dims);
51 int periods[2] = {1};
52 int reorder = 0;
53 MPI_Comm group;
54 MPI_Cart_create(MPI_COMM_WORLD, 2, dims, periods, reorder, &group);
55
56 int coord[2];
57 MPI_Cart_coords(group, mpi_rank, 2, coord);
58
59 std::cout
60 << "I'm processor " << mpi_rank
61 << " over " << mpi_size << " proc"
62 << ", (" << coord[0] << ", " << coord[1] << ") in the 2D grid"
63 << std::endl;
64
65 // initialize the DBCSR library
66 c_dbcsr_init_lib(MPI_COMM_WORLD, nullptr);
67
68 // the matrix will contain nblkrows_total row blocks and nblkcols_total column blocks
69 int nblkrows_total = 4;
70 int nblkcols_total = 4;
71
72 // set the block size for each row and column
73 std::vector<int> row_blk_sizes(nblkrows_total, 2), col_blk_sizes(nblkcols_total, 2);
74
75 // set the row and column distributions (here the distribution is set randomly)
76 auto row_dist = random_dist(nblkrows_total, dims[0]);
77 auto col_dist = random_dist(nblkcols_total, dims[1]);
78
79 // set the DBCSR distribution object
80 void* dist = nullptr;
81 c_dbcsr_distribution_new(&dist, group,
82 row_dist.data(), row_dist.size(),
83 col_dist.data(), col_dist.size());
84
85 // Fill all blocks, i.e. dense matrices
86 auto fill_matrix = [&](void*& matrix)
87 {
88 int max_row_size = *std::max_element(row_blk_sizes.begin(),row_blk_sizes.end());
89 int max_col_size = *std::max_element(col_blk_sizes.begin(),col_blk_sizes.end());
90 int max_nze = max_row_size * max_col_size;
91
92 std::vector<double> block;
93 block.reserve(max_nze);
94
95 for(int i = 0; i < nblkrows_total; i++)
96 {
97 for(int j = 0; j < nblkcols_total; j++)
98 {
99 int blk_proc = -1;
100 c_dbcsr_get_stored_coordinates(matrix, i, j, &blk_proc);
101
102 if(blk_proc == mpi_rank)
103 {
104 block.resize(row_blk_sizes[i] * col_blk_sizes[j]);
105 std::generate(block.begin(), block.end(), [&](){ return static_cast<double>(std::rand())/RAND_MAX; });
106 c_dbcsr_put_block_d(matrix, i, j, block.data(), block.size());
107 }
108 }
109 }
110 };
111
112 // create the DBCSR matrices, i.e. a double precision non symmetric matrix
113 // with nblkrows_total x nblkcols_total blocks and
114 // sizes "sum(row_blk_sizes)" x "sum(col_blk_sizes)", distributed as
115 // specified by the dist object
116
117 // create, fill and finalize matrix a
118 void* matrix_a = nullptr;
119 c_dbcsr_create_new_d(&matrix_a, "this is my matrix a", dist, 'N',
120 row_blk_sizes.data(), row_blk_sizes.size(),
121 col_blk_sizes.data(), col_blk_sizes.size());
122 fill_matrix(matrix_a);
123 c_dbcsr_finalize(matrix_a);
124
125 // create, fill and finalize matrix b
126 void* matrix_b = nullptr;
127 c_dbcsr_create_new_d(&matrix_b, "this is my matrix b", dist, 'N',
128 row_blk_sizes.data(), row_blk_sizes.size(),
129 col_blk_sizes.data(), col_blk_sizes.size());
130 fill_matrix(matrix_b);
131 c_dbcsr_finalize(matrix_b);
132
133 // create and finalize matrix c (empty)
134 void* matrix_c = nullptr;
135 c_dbcsr_create_new_d(&matrix_c, "matrix c", dist, 'N',
136 row_blk_sizes.data(), row_blk_sizes.size(),
137 col_blk_sizes.data(), col_blk_sizes.size());
138 c_dbcsr_finalize(matrix_c);
139
140 // multiply the matrices
141 c_dbcsr_multiply_d('N', 'N', 1.0, &matrix_a, &matrix_b, 0.0, &matrix_c, nullptr);
142
143 // print the matrices
144 c_dbcsr_print(matrix_a);
145 c_dbcsr_print(matrix_b);
146 c_dbcsr_print(matrix_c);
147
148 // release the matrices
149 c_dbcsr_release(&matrix_a);
150 c_dbcsr_release(&matrix_b);
151 c_dbcsr_release(&matrix_c);
152
153 c_dbcsr_distribution_release(&dist);
154
155 // free comm
156 MPI_Comm_free(&group);
157
158 // finalize the DBCSR library
159 c_dbcsr_finalize_lib();
160
161 // finalize MPI
162 MPI_Finalize();
163
164 return 0;
165 }
166