Line data Source code
1 : /*----------------------------------------------------------------------------*/
2 : /* CP2K: A general program to perform molecular dynamics simulations */
3 : /* Copyright 2000-2026 CP2K developers group <https://cp2k.org> */
4 : /* */
5 : /* SPDX-License-Identifier: BSD-3-Clause */
6 : /*----------------------------------------------------------------------------*/
7 : #include "dbm_matrix.h"
8 : #include "dbm_hyperparams.h"
9 :
10 : #include <assert.h>
11 : #include <math.h>
12 : #include <omp.h>
13 : #include <stdbool.h>
14 : #include <stddef.h>
15 : #include <stdlib.h>
16 : #include <string.h>
17 :
18 : /*******************************************************************************
19 : * \brief Creates a new matrix.
20 : * \author Ole Schuett
21 : ******************************************************************************/
22 2692336 : void dbm_create(dbm_matrix_t **matrix_out, dbm_distribution_t *dist,
23 : const char name[], const int nrows, const int ncols,
24 : const int row_sizes[nrows], const int col_sizes[ncols]) {
25 2692336 : assert(omp_get_num_threads() == 1);
26 :
27 2692336 : dbm_matrix_t *const matrix = calloc(1, sizeof(dbm_matrix_t));
28 :
29 2692336 : assert(dist->rows.length == nrows);
30 2692336 : assert(dist->cols.length == ncols);
31 2692336 : dbm_distribution_hold(dist);
32 2692336 : matrix->dist = dist;
33 :
34 2692336 : size_t size = (strlen(name) + 1) * sizeof(char);
35 2692336 : matrix->name = malloc(size);
36 2692336 : assert(matrix->name != NULL && name != NULL);
37 2692336 : memcpy(matrix->name, name, size);
38 :
39 2692336 : matrix->nrows = nrows;
40 2692336 : matrix->ncols = ncols;
41 :
42 2692336 : size = nrows * sizeof(int);
43 2692336 : matrix->row_sizes = malloc(size);
44 2692336 : assert(matrix->row_sizes != NULL || size == 0);
45 2692336 : if (size != 0) {
46 2691394 : memcpy(matrix->row_sizes, row_sizes, size);
47 : }
48 :
49 2692336 : size = ncols * sizeof(int);
50 2692336 : matrix->col_sizes = malloc(size);
51 2692336 : assert(matrix->col_sizes != NULL || size == 0);
52 2692336 : if (size != 0) {
53 2686091 : memcpy(matrix->col_sizes, col_sizes, size);
54 : }
55 :
56 2692336 : const int num_shards = dbm_get_num_shards(matrix);
57 2692336 : matrix->shards = malloc(num_shards * sizeof(dbm_shard_t));
58 2692336 : assert(matrix->shards != NULL || num_shards == 0);
59 2692336 : #pragma omp parallel for
60 : for (int ishard = 0; ishard < num_shards; ishard++) {
61 : dbm_shard_init(&matrix->shards[ishard]);
62 : }
63 :
64 2692336 : assert(*matrix_out == NULL);
65 2692336 : *matrix_out = matrix;
66 2692336 : }
67 :
68 : /*******************************************************************************
69 : * \brief Releases a matrix and all its ressources.
70 : * \author Ole Schuett
71 : ******************************************************************************/
72 2692336 : void dbm_release(dbm_matrix_t *matrix) {
73 2692336 : assert(omp_get_num_threads() == 1);
74 5384672 : for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
75 2692336 : dbm_shard_release(&matrix->shards[ishard]);
76 : }
77 2692336 : dbm_distribution_release(matrix->dist);
78 2692336 : free(matrix->row_sizes);
79 2692336 : free(matrix->col_sizes);
80 2692336 : free(matrix->shards);
81 2692336 : free(matrix->name);
82 2692336 : free(matrix);
83 2692336 : }
84 :
85 : /*******************************************************************************
86 : * \brief Copies content of matrix_b into matrix_a.
87 : * Matrices must have the same row/col block sizes and distribution.
88 : * \author Ole Schuett
89 : ******************************************************************************/
90 704694 : void dbm_copy(dbm_matrix_t *matrix_a, const dbm_matrix_t *matrix_b) {
91 704694 : assert(omp_get_num_threads() == 1);
92 :
93 704694 : assert(matrix_b->nrows == matrix_a->nrows);
94 6779872 : for (int i = 0; i < matrix_b->nrows; i++) {
95 6075178 : assert(matrix_b->row_sizes[i] == matrix_a->row_sizes[i]);
96 : }
97 704694 : assert(matrix_b->ncols == matrix_a->ncols);
98 32343650 : for (int i = 0; i < matrix_b->ncols; i++) {
99 31638956 : assert(matrix_b->col_sizes[i] == matrix_a->col_sizes[i]);
100 : }
101 :
102 704694 : assert(matrix_a->dist == matrix_b->dist);
103 :
104 704694 : #pragma omp parallel for DBM_OMP_SCHEDULE
105 : for (int ishard = 0; ishard < dbm_get_num_shards(matrix_a); ishard++) {
106 : dbm_shard_copy(&matrix_a->shards[ishard], &matrix_b->shards[ishard]);
107 : }
108 704694 : }
109 :
110 : /*******************************************************************************
111 : * \brief Copies content of matrix_b into matrix_a.
112 : * Matrices may have different distributions.
113 : * \author Ole Schuett
114 : ******************************************************************************/
115 144 : void dbm_redistribute(const dbm_matrix_t *matrix, dbm_matrix_t *redist) {
116 144 : assert(omp_get_num_threads() == 1);
117 144 : assert(matrix->nrows == redist->nrows);
118 6384 : for (int i = 0; i < matrix->nrows; i++) {
119 6240 : assert(matrix->row_sizes[i] == redist->row_sizes[i]);
120 : }
121 144 : assert(matrix->ncols == redist->ncols);
122 6384 : for (int i = 0; i < matrix->ncols; i++) {
123 6240 : assert(matrix->col_sizes[i] == redist->col_sizes[i]);
124 : }
125 :
126 144 : assert(cp_mpi_comms_are_similar(matrix->dist->comm, redist->dist->comm));
127 144 : const cp_mpi_comm_t comm = redist->dist->comm;
128 144 : const int nranks = cp_mpi_comm_size(comm);
129 :
130 : // 1st pass: Compute send_count.
131 144 : int send_count[nranks];
132 144 : memset(send_count, 0, nranks * sizeof(int));
133 288 : for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
134 144 : const dbm_shard_t *const shard = &matrix->shards[ishard];
135 8936 : for (int iblock = 0; iblock < shard->nblocks; iblock++) {
136 8792 : const dbm_block_t *const blk = &shard->blocks[iblock];
137 8792 : const int row_size = matrix->row_sizes[blk->row];
138 8792 : const int col_size = matrix->col_sizes[blk->col];
139 8792 : const int block_size = row_size * col_size;
140 8792 : const int rank = dbm_get_stored_coordinates(redist, blk->row, blk->col);
141 8792 : assert(0 <= rank && rank < nranks);
142 8792 : send_count[rank] += 2 + block_size;
143 : }
144 : }
145 :
146 : // 1st communication: Exchange counts.
147 144 : int recv_count[nranks];
148 144 : cp_mpi_alltoall_int(send_count, 1, recv_count, 1, comm);
149 :
150 : // Compute displacements and allocate data buffers.
151 144 : int send_displ[nranks + 1], recv_displ[nranks + 1];
152 144 : send_displ[0] = recv_displ[0] = 0;
153 432 : for (int irank = 1; irank <= nranks; irank++) {
154 288 : send_displ[irank] = send_displ[irank - 1] + send_count[irank - 1];
155 288 : recv_displ[irank] = recv_displ[irank - 1] + recv_count[irank - 1];
156 : }
157 144 : const int total_send_count = send_displ[nranks];
158 144 : const int total_recv_count = recv_displ[nranks];
159 144 : double *const data_send = cp_mpi_alloc_mem(total_send_count * sizeof(double));
160 144 : double *const data_recv = cp_mpi_alloc_mem(total_recv_count * sizeof(double));
161 :
162 : // 2nd pass: Fill send_data.
163 144 : int send_data_positions[nranks];
164 144 : memcpy(send_data_positions, send_displ, nranks * sizeof(int));
165 288 : for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
166 144 : const dbm_shard_t *const shard = &matrix->shards[ishard];
167 8936 : for (int iblock = 0; iblock < shard->nblocks; iblock++) {
168 8792 : const dbm_block_t *const blk = &shard->blocks[iblock];
169 8792 : const double *const blk_data = &shard->data[blk->offset];
170 8792 : const int row_size = matrix->row_sizes[blk->row];
171 8792 : const int col_size = matrix->col_sizes[blk->col];
172 8792 : const int block_size = row_size * col_size;
173 8792 : const int rank = dbm_get_stored_coordinates(redist, blk->row, blk->col);
174 8792 : const int pos = send_data_positions[rank];
175 8792 : data_send[pos + 0] = blk->row; // send integers as doubles
176 8792 : data_send[pos + 1] = blk->col;
177 8792 : memcpy(&data_send[pos + 2], blk_data, block_size * sizeof(double));
178 8792 : send_data_positions[rank] += 2 + block_size;
179 : }
180 : }
181 432 : for (int irank = 0; irank < nranks; irank++) {
182 288 : assert(send_data_positions[irank] == send_displ[irank + 1]);
183 : }
184 :
185 : // 2nd communication: Exchange data.
186 144 : cp_mpi_alltoallv_double(data_send, send_count, send_displ, data_recv,
187 : recv_count, recv_displ, comm);
188 144 : cp_mpi_free_mem(data_send);
189 :
190 : // 3rd pass: Unpack data.
191 144 : dbm_clear(redist);
192 144 : int recv_data_pos = 0;
193 8936 : while (recv_data_pos < total_recv_count) {
194 8792 : const int row = (int)data_recv[recv_data_pos + 0];
195 8792 : const int col = (int)data_recv[recv_data_pos + 1];
196 8792 : assert(data_recv[recv_data_pos + 0] == (double)row);
197 8792 : assert(data_recv[recv_data_pos + 1] == (double)col);
198 8792 : dbm_put_block(redist, row, col, false, &data_recv[recv_data_pos + 2]);
199 8792 : const int row_size = matrix->row_sizes[row];
200 8792 : const int col_size = matrix->col_sizes[col];
201 8792 : const int block_size = row_size * col_size;
202 8792 : recv_data_pos += 2 + block_size;
203 : }
204 144 : assert(recv_data_pos == total_recv_count);
205 144 : cp_mpi_free_mem(data_recv);
206 144 : }
207 :
208 : /*******************************************************************************
209 : * \brief Looks up a block from given matrics. This routine is thread-safe.
210 : * If the block is not found then a null pointer is returned.
211 : * \author Ole Schuett
212 : ******************************************************************************/
213 19352608 : void dbm_get_block_p(dbm_matrix_t *matrix, const int row, const int col,
214 : double **block, int *row_size, int *col_size) {
215 19352608 : assert(0 <= row && row < matrix->nrows);
216 19352608 : assert(0 <= col && col < matrix->ncols);
217 19352608 : assert(dbm_get_stored_coordinates(matrix, row, col) == matrix->dist->my_rank);
218 19352608 : *row_size = matrix->row_sizes[row];
219 19352608 : *col_size = matrix->col_sizes[col];
220 19352608 : *block = NULL;
221 :
222 19352608 : const int ishard = dbm_get_shard_index(matrix, row, col);
223 19352608 : const dbm_shard_t *const shard = &matrix->shards[ishard];
224 19352608 : const dbm_block_t *const blk = dbm_shard_lookup(shard, row, col);
225 19352608 : if (blk != NULL) {
226 17727682 : *block = &shard->data[blk->offset];
227 : }
228 19352608 : }
229 :
230 : /*******************************************************************************
231 : * \brief Adds a block to given matrix. This routine is thread-safe.
232 : * If block already exist then it gets overwritten (or summed).
233 : * \author Ole Schuett
234 : ******************************************************************************/
235 36328463 : void dbm_put_block(dbm_matrix_t *matrix, const int row, const int col,
236 : const bool summation, const double *block) {
237 36328463 : assert(0 <= row && row < matrix->nrows);
238 36328463 : assert(0 <= col && col < matrix->ncols);
239 36328463 : assert(dbm_get_stored_coordinates(matrix, row, col) == matrix->dist->my_rank);
240 36328463 : const int row_size = matrix->row_sizes[row];
241 36328463 : const int col_size = matrix->col_sizes[col];
242 36328463 : const int block_size = row_size * col_size;
243 :
244 36328463 : const int ishard = dbm_get_shard_index(matrix, row, col);
245 36328463 : dbm_shard_t *const shard = &matrix->shards[ishard];
246 36328463 : omp_set_lock(&shard->lock);
247 36328463 : const dbm_block_t *const blk =
248 36328463 : dbm_shard_get_or_allocate_block(shard, row, col, block_size);
249 36328463 : double *const blk_data = &shard->data[blk->offset];
250 36328463 : if (summation) {
251 4531254 : assert(blk_data != NULL || 0 == block_size);
252 1615089585 : for (int i = 0; i < block_size; i++) {
253 1610558331 : blk_data[i] += block[i];
254 : }
255 31797209 : } else if (block_size != 0) {
256 31792110 : memcpy(blk_data, block, block_size * sizeof(double));
257 : }
258 36328463 : omp_unset_lock(&shard->lock);
259 36328463 : }
260 :
261 : /*******************************************************************************
262 : * \brief Remove all blocks from matrix, but does not release underlying memory.
263 : * \author Ole Schuett
264 : ******************************************************************************/
265 3063996 : void dbm_clear(dbm_matrix_t *matrix) {
266 3063996 : assert(omp_get_num_threads() == 1);
267 :
268 3063996 : #pragma omp parallel for DBM_OMP_SCHEDULE
269 : for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
270 : dbm_shard_t *const shard = &matrix->shards[ishard];
271 : shard->data_promised = shard->data_size = shard->nblocks = 0;
272 : // Does not deallocate memory, hence data_allocated remains unchanged.
273 : memset(shard->hashtable, 0, shard->hashtable_size * sizeof(int));
274 : }
275 3063996 : }
276 :
277 : /*******************************************************************************
278 : * \brief Removes all blocks from the matrix whose norm is below the threshold.
279 : * Blocks of size zero are always kept.
280 : * \author Ole Schuett
281 : ******************************************************************************/
282 939434 : void dbm_filter(dbm_matrix_t *matrix, const double eps) {
283 939434 : assert(omp_get_num_threads() == 1);
284 :
285 939434 : if (eps == 0.0) {
286 : return;
287 : }
288 862416 : const double eps2 = eps * eps;
289 :
290 862416 : #pragma omp parallel for DBM_OMP_SCHEDULE
291 : for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
292 : dbm_shard_t *const shard = &matrix->shards[ishard];
293 : const int old_nblocks = shard->nblocks;
294 : shard->nblocks = shard->data_promised = 0;
295 : memset(shard->hashtable, 0, shard->hashtable_size * sizeof(int));
296 :
297 : for (int iblock = 0; iblock < old_nblocks; iblock++) {
298 : const dbm_block_t old_blk = shard->blocks[iblock];
299 : const double *old_blk_data = &shard->data[old_blk.offset];
300 : const int row_size = matrix->row_sizes[old_blk.row];
301 : const int col_size = matrix->col_sizes[old_blk.col];
302 : const int block_size = row_size * col_size;
303 : double norm = 0.0;
304 : for (int i = 0; i < block_size; i++) {
305 : const double value = old_blk_data[i];
306 : norm += value * value;
307 : if (eps2 <= norm) {
308 : break;
309 : }
310 : }
311 : // For historic reasons zero-sized blocks are never filtered.
312 : if (block_size > 0 && norm < eps2) {
313 : continue; // filter the block
314 : }
315 : // Re-create block.
316 : dbm_block_t *const new_blk = dbm_shard_promise_new_block(
317 : shard, old_blk.row, old_blk.col, block_size);
318 : assert(new_blk->offset <= old_blk.offset);
319 : if (new_blk->offset != old_blk.offset) {
320 : // Using memmove because it can handle overlapping buffers.
321 : double *new_blk_data = &shard->data[new_blk->offset];
322 : memmove(new_blk_data, old_blk_data, block_size * sizeof(double));
323 : }
324 : }
325 : shard->data_size = shard->data_promised;
326 : // TODO: Could call realloc to release excess memory.
327 : }
328 : }
329 :
330 : /*******************************************************************************
331 : * \brief Adds list of blocks efficiently. The blocks will be filled with zeros.
332 : * This routine must always be called within an OpenMP parallel region.
333 : * \author Ole Schuett and Hans Pabst
334 : ******************************************************************************/
335 2048114 : void dbm_reserve_blocks(dbm_matrix_t *matrix, const int nblocks,
336 : const int rows[], const int cols[]) {
337 2048114 : assert(omp_get_num_threads() == omp_get_max_threads() &&
338 : "Please call dbm_reserve_blocks within an OpenMP parallel region.");
339 2048114 : const int my_rank = matrix->dist->my_rank;
340 2048114 : const int nshards = dbm_get_num_shards(matrix);
341 2048114 : const int nthreads = omp_get_num_threads();
342 2048114 : const int ithread = omp_get_thread_num();
343 : // Shared by the threads: blocks are binned by shard such that a single
344 : // thread owns a shard (no locks).
345 2048114 : static int *reserve_offset = NULL; // per thread and shard
346 2048114 : static int *reserve_start = NULL; // per shard (and end)
347 2048114 : static int *reserve_rowcol = NULL; // row/col pairs grouped by shard
348 :
349 : #pragma omp single
350 : {
351 2048114 : reserve_offset = calloc((size_t)nthreads * nshards, sizeof(int));
352 2048114 : reserve_start = malloc((nshards + 1) * sizeof(int));
353 2048114 : assert(reserve_offset != NULL && reserve_start != NULL);
354 : }
355 2048114 : int *const offset = reserve_offset + (size_t)ithread * nshards;
356 34245834 : for (int i = 0; i < nblocks; i++) {
357 32197720 : const int row = rows[i], col = cols[i];
358 32197720 : assert(0 <= row && row < matrix->nrows);
359 32197720 : assert(0 <= col && col < matrix->ncols);
360 32197720 : assert(dbm_get_stored_coordinates(matrix, row, col) == my_rank);
361 32197720 : ++offset[dbm_get_shard_index(matrix, row, col)];
362 : }
363 2048114 : #pragma omp barrier
364 :
365 : #pragma omp single
366 : { // exclusive prefix sum over shards (major) and threads (minor)
367 : int total = 0;
368 4096228 : for (int ishard = 0; ishard < nshards; ishard++) {
369 2048114 : reserve_start[ishard] = total;
370 4096228 : for (int i = 0; i < nthreads; i++) {
371 2048114 : int *const count = reserve_offset + (size_t)i * nshards + ishard;
372 2048114 : const int n = *count;
373 2048114 : *count = total;
374 2048114 : total += n;
375 : }
376 : }
377 2048114 : reserve_start[nshards] = total;
378 4096228 : reserve_rowcol =
379 2048114 : (0 < total ? malloc(2 * (size_t)total * sizeof(int)) : NULL);
380 2048114 : assert(reserve_rowcol != NULL || 0 == total);
381 : }
382 34245834 : for (int i = 0; i < nblocks; i++) {
383 32197720 : const int j = offset[dbm_get_shard_index(matrix, rows[i], cols[i])]++;
384 32197720 : reserve_rowcol[2 * j + 0] = rows[i];
385 32197720 : reserve_rowcol[2 * j + 1] = cols[i];
386 : }
387 2048114 : #pragma omp barrier
388 :
389 : #pragma omp for DBM_OMP_SCHEDULE
390 2048114 : for (int ishard = 0; ishard < nshards; ishard++) {
391 2048114 : dbm_shard_t *const shard = &matrix->shards[ishard];
392 : // one rebuild of the hashtable rather than one per growth
393 2048114 : dbm_shard_reserve(shard, shard->nblocks + reserve_start[ishard + 1] -
394 2048114 : reserve_start[ishard]);
395 34245834 : for (int j = reserve_start[ishard]; j < reserve_start[ishard + 1]; j++) {
396 32197720 : const int row = reserve_rowcol[2 * j + 0];
397 32197720 : const int col = reserve_rowcol[2 * j + 1];
398 32197720 : const int block_size = matrix->row_sizes[row] * matrix->col_sizes[col];
399 32197720 : dbm_shard_get_or_promise_block(shard, row, col, block_size);
400 : }
401 2048114 : dbm_shard_allocate_promised_blocks(shard);
402 : }
403 :
404 : #pragma omp single
405 : {
406 2048114 : free(reserve_rowcol);
407 2048114 : free(reserve_start);
408 2048114 : free(reserve_offset);
409 2048114 : reserve_rowcol = reserve_start = reserve_offset = NULL;
410 : }
411 2048114 : }
412 :
413 : /*******************************************************************************
414 : * \brief Multiplies all entries in the given matrix by the given factor alpha.
415 : * \author Ole Schuett
416 : ******************************************************************************/
417 768009 : void dbm_scale(dbm_matrix_t *matrix, const double alpha) {
418 768009 : assert(omp_get_num_threads() == 1);
419 768009 : if (alpha == 1.0) {
420 : return;
421 : }
422 581123 : if (alpha == 0.0) {
423 567003 : dbm_zero(matrix);
424 567003 : return;
425 : }
426 :
427 14120 : #pragma omp parallel for DBM_OMP_SCHEDULE
428 : for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
429 : dbm_shard_t *const shard = &matrix->shards[ishard];
430 : for (int i = 0; i < shard->data_size; i++) {
431 : shard->data[i] *= alpha;
432 : }
433 : }
434 : }
435 :
436 : /*******************************************************************************
437 : * \brief Sets all blocks in the given matrix to zero.
438 : * \author Ole Schuett
439 : ******************************************************************************/
440 567003 : void dbm_zero(dbm_matrix_t *matrix) {
441 567003 : assert(omp_get_num_threads() == 1);
442 :
443 567003 : #pragma omp parallel for DBM_OMP_SCHEDULE
444 : for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
445 : dbm_shard_t *const shard = &matrix->shards[ishard];
446 : if (shard->data != NULL) {
447 : memset(shard->data, 0, shard->data_size * sizeof(double));
448 : }
449 : }
450 567003 : }
451 :
452 : /*******************************************************************************
453 : * \brief Adds matrix_b to matrix_a.
454 : * \author Ole Schuett
455 : ******************************************************************************/
456 349908 : void dbm_add(dbm_matrix_t *matrix_a, const dbm_matrix_t *matrix_b) {
457 349908 : assert(omp_get_num_threads() == 1);
458 349908 : assert(matrix_a->dist == matrix_b->dist);
459 :
460 349908 : #pragma omp parallel for DBM_OMP_SCHEDULE
461 : for (int ishard = 0; ishard < dbm_get_num_shards(matrix_b); ishard++) {
462 : dbm_shard_t *const shard_a = &matrix_a->shards[ishard];
463 : const dbm_shard_t *shard_b = &matrix_b->shards[ishard];
464 : for (int iblock = 0; iblock < shard_b->nblocks; iblock++) {
465 : const dbm_block_t blk_b = shard_b->blocks[iblock];
466 : const int row_size = matrix_b->row_sizes[blk_b.row];
467 : const int col_size = matrix_b->col_sizes[blk_b.col];
468 : assert(row_size == matrix_a->row_sizes[blk_b.row]);
469 : assert(col_size == matrix_a->col_sizes[blk_b.col]);
470 : const int block_size = row_size * col_size;
471 : const dbm_block_t *const blk_a = dbm_shard_get_or_allocate_block(
472 : shard_a, blk_b.row, blk_b.col, block_size);
473 : double *const data_a = &shard_a->data[blk_a->offset];
474 : const double *const data_b = &shard_b->data[blk_b.offset];
475 : for (int i = 0; i < block_size; i++) {
476 : data_a[i] += data_b[i];
477 : }
478 : }
479 : }
480 349908 : }
481 :
482 : /*******************************************************************************
483 : * \brief Creates an iterator for the blocks of the given matrix.
484 : * The iteration order is not stable.
485 : * This routine must always be called within an OpenMP parallel region.
486 : * \author Ole Schuett
487 : ******************************************************************************/
488 5044712 : void dbm_iterator_start(dbm_iterator_t **iter_out, const dbm_matrix_t *matrix) {
489 5044712 : assert(omp_get_num_threads() == omp_get_max_threads() &&
490 : "Please call dbm_iterator_start within an OpenMP parallel region.");
491 5044712 : dbm_iterator_t *iter = malloc(sizeof(dbm_iterator_t));
492 5044712 : assert(iter != NULL);
493 5044712 : iter->matrix = matrix;
494 5044712 : iter->next_block = 0;
495 5044712 : iter->next_shard = omp_get_thread_num();
496 5796353 : while (iter->next_shard < dbm_get_num_shards(matrix) &&
497 5044712 : matrix->shards[iter->next_shard].nblocks == 0) {
498 751641 : iter->next_shard += omp_get_num_threads();
499 : }
500 5044712 : assert(*iter_out == NULL);
501 5044712 : *iter_out = iter;
502 5044712 : }
503 :
504 : /*******************************************************************************
505 : * \brief Returns number of blocks the iterator will provide to calling thread.
506 : * \author Ole Schuett
507 : ******************************************************************************/
508 787744 : int dbm_iterator_num_blocks(const dbm_iterator_t *iter) {
509 787744 : int num_blocks = 0;
510 787744 : for (int ishard = omp_get_thread_num();
511 1575488 : ishard < dbm_get_num_shards(iter->matrix);
512 787744 : ishard += omp_get_num_threads()) {
513 787744 : num_blocks += iter->matrix->shards[ishard].nblocks;
514 : }
515 787744 : return num_blocks;
516 : }
517 :
518 : /*******************************************************************************
519 : * \brief Tests whether the given iterator has any block left.
520 : * \author Ole Schuett
521 : ******************************************************************************/
522 85443316 : bool dbm_iterator_blocks_left(const dbm_iterator_t *iter) {
523 85443316 : return iter->next_shard < dbm_get_num_shards(iter->matrix);
524 : }
525 :
526 : /*******************************************************************************
527 : * \brief Returns the next block from the given iterator.
528 : * \author Ole Schuett
529 : ******************************************************************************/
530 91330041 : void dbm_iterator_next_block(dbm_iterator_t *iter, int *row, int *col,
531 : double **block, int *row_size, int *col_size) {
532 91330041 : const dbm_matrix_t *matrix = iter->matrix;
533 91330041 : assert(iter->next_shard < dbm_get_num_shards(matrix));
534 91330041 : const dbm_shard_t *shard = &matrix->shards[iter->next_shard];
535 91330041 : assert(iter->next_block < shard->nblocks);
536 91330041 : dbm_block_t *blk = &shard->blocks[iter->next_block];
537 :
538 91330041 : *row = blk->row;
539 91330041 : *col = blk->col;
540 91330041 : *row_size = matrix->row_sizes[blk->row];
541 91330041 : *col_size = matrix->col_sizes[blk->col];
542 91330041 : *block = &shard->data[blk->offset];
543 :
544 91330041 : iter->next_block++;
545 91330041 : if (iter->next_block >= shard->nblocks) {
546 : // Advance to the next non-empty shard...
547 4293071 : iter->next_shard += omp_get_num_threads();
548 4293071 : while (iter->next_shard < dbm_get_num_shards(matrix) &&
549 0 : matrix->shards[iter->next_shard].nblocks == 0) {
550 0 : iter->next_shard += omp_get_num_threads();
551 : }
552 4293071 : iter->next_block = 0; // ...and continue with its first block.
553 : }
554 91330041 : }
555 :
556 : /*******************************************************************************
557 : * \brief Releases the given iterator.
558 : * \author Ole Schuett
559 : ******************************************************************************/
560 5044712 : void dbm_iterator_stop(dbm_iterator_t *iter) { free(iter); }
561 :
562 : /*******************************************************************************
563 : * \brief Private routine to accumulate using Kahan's summation.
564 : * \author Hans Pabst
565 : ******************************************************************************/
566 14383468 : static double kahan_sum(double value, double *accumulator,
567 : double *compensation) {
568 14383468 : double r, c;
569 14383468 : assert(NULL != accumulator && NULL != compensation);
570 14383468 : c = value - *compensation;
571 14383468 : r = *accumulator + c;
572 14383468 : *compensation = (r - *accumulator) - c;
573 14383468 : *accumulator = r;
574 14383468 : return r;
575 : }
576 :
577 : /*******************************************************************************
578 : * \brief Computes a checksum of the given matrix.
579 : * \author Ole Schuett
580 : ******************************************************************************/
581 190 : double dbm_checksum(const dbm_matrix_t *matrix) {
582 190 : double checksum = 0.0, compensation = 0.0;
583 380 : for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
584 190 : const dbm_shard_t *shard = &matrix->shards[ishard];
585 14383658 : for (int i = 0; i < shard->data_size; i++) {
586 14383468 : const double value = shard->data[i];
587 14383468 : kahan_sum(value * value, &checksum, &compensation);
588 : }
589 : }
590 190 : cp_mpi_sum_double(&checksum, 1, matrix->dist->comm);
591 190 : return checksum;
592 : }
593 :
594 : /*******************************************************************************
595 : * \brief Returns the largest absolute value of all matrix elements.
596 : * \author Ole Schuett
597 : ******************************************************************************/
598 48 : double dbm_maxabs(const dbm_matrix_t *matrix) {
599 48 : double maxabs = 0.0;
600 96 : for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
601 48 : const dbm_shard_t *shard = &matrix->shards[ishard];
602 1599050 : for (int i = 0; i < shard->data_size; i++) {
603 1599002 : maxabs = fmax(maxabs, fabs(shard->data[i]));
604 : }
605 : }
606 48 : cp_mpi_max_double(&maxabs, 1, matrix->dist->comm);
607 48 : return maxabs;
608 : }
609 :
610 : /*******************************************************************************
611 : * \brief Calculates maximum relative difference between matrix_a and matrix_b.
612 : * \author Hans Pabst
613 : ******************************************************************************/
614 0 : double dbm_maxeps(const dbm_matrix_t *matrix_a, const dbm_matrix_t *matrix_b) {
615 0 : const int num_shards = dbm_get_num_shards(matrix_a);
616 0 : double epsilon = 0;
617 :
618 0 : assert(omp_get_num_threads() == 1);
619 0 : assert(matrix_a->dist == matrix_b->dist);
620 0 : assert(num_shards == dbm_get_num_shards(matrix_b));
621 :
622 : // Handles asymmetric sparsity: blocks in A but not B are compared against
623 : // zero, and blocks in B but not A are checked in a second pass. The
624 : // reduction(max) avoids the previous data race on the shared epsilon.
625 0 : #pragma omp parallel for reduction(max : epsilon) DBM_OMP_SCHEDULE
626 : for (int ishard = 0; ishard < num_shards; ++ishard) {
627 : const dbm_shard_t *const shard_a = &matrix_a->shards[ishard];
628 : const dbm_shard_t *const shard_b = &matrix_b->shards[ishard];
629 : for (int iblock = 0; iblock < shard_a->nblocks; ++iblock) {
630 : const dbm_block_t *const blk_a = &shard_a->blocks[iblock];
631 : const dbm_block_t *const blk_b =
632 : dbm_shard_lookup(shard_b, blk_a->row, blk_a->col);
633 : const int row_size = matrix_a->row_sizes[blk_a->row];
634 : const int col_size = matrix_a->col_sizes[blk_a->col];
635 : assert(NULL == blk_b || row_size == matrix_b->row_sizes[blk_b->row]);
636 : assert(NULL == blk_b || col_size == matrix_b->col_sizes[blk_b->col]);
637 : const double *const data_a = &shard_a->data[blk_a->offset];
638 : const double *const data_b =
639 : (NULL == blk_b) ? NULL : &shard_b->data[blk_b->offset];
640 : const int block_size = row_size * col_size;
641 : for (int i = 0; i < block_size; ++i) {
642 : const double d = data_a[i] - (NULL == data_b ? 0.0 : data_b[i]);
643 : const double e = fabs(0 != data_a[i] ? (d / data_a[i]) : d);
644 : if (epsilon < e) {
645 : epsilon = e;
646 : }
647 : }
648 : }
649 : // Second pass: blocks present in B but absent from A (treat A as zero).
650 : for (int iblock = 0; iblock < shard_b->nblocks; ++iblock) {
651 : const dbm_block_t *const blk_b = &shard_b->blocks[iblock];
652 : if (NULL != dbm_shard_lookup(shard_a, blk_b->row, blk_b->col)) {
653 : continue;
654 : }
655 : const int row_size = matrix_b->row_sizes[blk_b->row];
656 : const int col_size = matrix_b->col_sizes[blk_b->col];
657 : const double *const data_b = &shard_b->data[blk_b->offset];
658 : const int block_size = row_size * col_size;
659 : for (int i = 0; i < block_size; ++i) {
660 : const double e = fabs(data_b[i]);
661 : if (epsilon < e) {
662 : epsilon = e;
663 : }
664 : }
665 : }
666 : }
667 :
668 0 : cp_mpi_max_double(&epsilon, 1, matrix_a->dist->comm);
669 0 : return epsilon;
670 : }
671 :
672 : /*******************************************************************************
673 : * \brief Returns the name of the matrix of the given matrix.
674 : * \author Ole Schuett
675 : ******************************************************************************/
676 2899377 : const char *dbm_get_name(const dbm_matrix_t *matrix) { return matrix->name; }
677 :
678 : /*******************************************************************************
679 : * \brief Returns the number of local Non-Zero Elements of the given matrix.
680 : * \author Ole Schuett
681 : ******************************************************************************/
682 2794300 : int dbm_get_nze(const dbm_matrix_t *matrix) {
683 2794300 : int nze = 0;
684 5588600 : for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
685 2794300 : nze += matrix->shards[ishard].data_size;
686 : }
687 2794300 : return nze;
688 : }
689 :
690 : /*******************************************************************************
691 : * \brief Returns the number of local blocks of the given matrix.
692 : * \author Ole Schuett
693 : ******************************************************************************/
694 1083840 : int dbm_get_num_blocks(const dbm_matrix_t *matrix) {
695 1083840 : int nblocks = 0;
696 2167680 : for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
697 1083840 : nblocks += matrix->shards[ishard].nblocks;
698 : }
699 1083840 : return nblocks;
700 : }
701 :
702 : /*******************************************************************************
703 : * \brief Returns the row block sizes of the given matrix.
704 : * \author Ole Schuett
705 : ******************************************************************************/
706 5830177 : void dbm_get_row_sizes(const dbm_matrix_t *matrix, int *nrows,
707 : const int **row_sizes) {
708 5830177 : *nrows = matrix->nrows;
709 5830177 : *row_sizes = matrix->row_sizes;
710 5830177 : }
711 :
712 : /*******************************************************************************
713 : * \brief Returns the column block sizes of the given matrix.
714 : * \author Ole Schuett
715 : ******************************************************************************/
716 3130537 : void dbm_get_col_sizes(const dbm_matrix_t *matrix, int *ncols,
717 : const int **col_sizes) {
718 3130537 : *ncols = matrix->ncols;
719 3130537 : *col_sizes = matrix->col_sizes;
720 3130537 : }
721 :
722 : /*******************************************************************************
723 : * \brief Returns the local row block sizes of the given matrix.
724 : * \author Ole Schuett
725 : ******************************************************************************/
726 468392 : void dbm_get_local_rows(const dbm_matrix_t *matrix, int *nlocal_rows,
727 : const int **local_rows) {
728 468392 : *nlocal_rows = matrix->dist->rows.nlocals;
729 468392 : *local_rows = matrix->dist->rows.local_indicies;
730 468392 : }
731 :
732 : /*******************************************************************************
733 : * \brief Returns the local column block sizes of the given matrix.
734 : * \author Ole Schuett
735 : ******************************************************************************/
736 169160 : void dbm_get_local_cols(const dbm_matrix_t *matrix, int *nlocal_cols,
737 : const int **local_cols) {
738 169160 : *nlocal_cols = matrix->dist->cols.nlocals;
739 169160 : *local_cols = matrix->dist->cols.local_indicies;
740 169160 : }
741 :
742 : /*******************************************************************************
743 : * \brief Returns the MPI rank on which the given block should be stored.
744 : * \author Ole Schuett
745 : ******************************************************************************/
746 96061266 : int dbm_get_stored_coordinates(const dbm_matrix_t *matrix, const int row,
747 : const int col) {
748 96061266 : return dbm_distribution_stored_coords(matrix->dist, row, col);
749 : }
750 :
751 : /*******************************************************************************
752 : * \brief Returns the distribution of the given matrix.
753 : * \author Ole Schuett
754 : ******************************************************************************/
755 1839883 : const dbm_distribution_t *dbm_get_distribution(const dbm_matrix_t *matrix) {
756 1839883 : return matrix->dist;
757 : }
758 :
759 : // EOF
|
