Line data    Source code

       1             : /*----------------------------------------------------------------------------*/
       2             : /*  CP2K: A general program to perform molecular dynamics simulations         */
       3             : /*  Copyright 2000-2024 CP2K developers group <https://cp2k.org>              */
       4             : /*                                                                            */
       5             : /*  SPDX-License-Identifier: BSD-3-Clause                                     */
       6             : /*----------------------------------------------------------------------------*/
       7             : 
       8             : #include <assert.h>
       9             : #include <math.h>
      10             : #include <omp.h>
      11             : #include <stdbool.h>
      12             : #include <stddef.h>
      13             : #include <stdlib.h>
      14             : #include <string.h>
      15             : 
      16             : #include "dbm_hyperparams.h"
      17             : #include "dbm_matrix.h"
      18             : #include "dbm_mpi.h"
      19             : 
      20             : /*******************************************************************************
      21             :  * \brief Creates a new matrix.
      22             :  * \author Ole Schuett
      23             :  ******************************************************************************/
      24     1097007 : void dbm_create(dbm_matrix_t **matrix_out, dbm_distribution_t *dist,
      25             :                 const char name[], const int nrows, const int ncols,
      26             :                 const int row_sizes[nrows], const int col_sizes[ncols]) {
      27     1097007 :   assert(omp_get_num_threads() == 1);
      28             : 
      29     1097007 :   dbm_matrix_t *matrix = calloc(1, sizeof(dbm_matrix_t));
      30             : 
      31     1097007 :   assert(dist->rows.length == nrows);
      32     1097007 :   assert(dist->cols.length == ncols);
      33     1097007 :   dbm_distribution_hold(dist);
      34     1097007 :   matrix->dist = dist;
      35             : 
      36     1097007 :   size_t size = (strlen(name) + 1) * sizeof(char);
      37     1097007 :   matrix->name = malloc(size);
      38     1097007 :   memcpy(matrix->name, name, size);
      39             : 
      40     1097007 :   matrix->nrows = nrows;
      41     1097007 :   matrix->ncols = ncols;
      42             : 
      43     1097007 :   size = nrows * sizeof(int);
      44     1097007 :   matrix->row_sizes = malloc(size);
      45     1097007 :   memcpy(matrix->row_sizes, row_sizes, size);
      46             : 
      47     1097007 :   size = ncols * sizeof(int);
      48     1097007 :   matrix->col_sizes = malloc(size);
      49     1097007 :   memcpy(matrix->col_sizes, col_sizes, size);
      50             : 
      51     1097007 :   matrix->shards = malloc(dbm_get_num_shards(matrix) * sizeof(dbm_shard_t));
      52     1097007 : #pragma omp parallel for
      53             :   for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
      54             :     dbm_shard_init(&matrix->shards[ishard]);
      55             :   }
      56             : 
      57     1097007 :   assert(*matrix_out == NULL);
      58     1097007 :   *matrix_out = matrix;
      59     1097007 : }
      60             : 
      61             : /*******************************************************************************
      62             :  * \brief Releases a matrix and all its ressources.
      63             :  * \author Ole Schuett
      64             :  ******************************************************************************/
      65     1097007 : void dbm_release(dbm_matrix_t *matrix) {
      66     1097007 :   assert(omp_get_num_threads() == 1);
      67     2194014 :   for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
      68     1097007 :     dbm_shard_release(&matrix->shards[ishard]);
      69             :   }
      70     1097007 :   dbm_distribution_release(matrix->dist);
      71     1097007 :   free(matrix->row_sizes);
      72     1097007 :   free(matrix->col_sizes);
      73     1097007 :   free(matrix->shards);
      74     1097007 :   free(matrix->name);
      75     1097007 :   free(matrix);
      76     1097007 : }
      77             : 
      78             : /*******************************************************************************
      79             :  * \brief Copies content of matrix_b into matrix_a.
      80             :  *        Matrices must have the same row/col block sizes and distribution.
      81             :  * \author Ole Schuett
      82             :  ******************************************************************************/
      83      256044 : void dbm_copy(dbm_matrix_t *matrix_a, const dbm_matrix_t *matrix_b) {
      84      256044 :   assert(omp_get_num_threads() == 1);
      85             : 
      86      256044 :   assert(matrix_b->nrows == matrix_a->nrows);
      87     4085745 :   for (int i = 0; i < matrix_b->nrows; i++) {
      88     3829701 :     assert(matrix_b->row_sizes[i] == matrix_a->row_sizes[i]);
      89             :   }
      90      256044 :   assert(matrix_b->ncols == matrix_a->ncols);
      91    34353552 :   for (int i = 0; i < matrix_b->ncols; i++) {
      92    34097508 :     assert(matrix_b->col_sizes[i] == matrix_a->col_sizes[i]);
      93             :   }
      94             : 
      95      256044 :   assert(matrix_a->dist == matrix_b->dist);
      96             : 
      97      256044 : #pragma omp parallel for schedule(dynamic)
      98             :   for (int ishard = 0; ishard < dbm_get_num_shards(matrix_a); ishard++) {
      99             :     dbm_shard_copy(&matrix_a->shards[ishard], &matrix_b->shards[ishard]);
     100             :   }
     101      256044 : }
     102             : 
     103             : /*******************************************************************************
     104             :  * \brief Copies content of matrix_b into matrix_a.
     105             :  *        Matrices may have different distributions.
     106             :  * \author Ole Schuett
     107             :  ******************************************************************************/
     108         144 : void dbm_redistribute(const dbm_matrix_t *matrix, dbm_matrix_t *redist) {
     109         144 :   assert(omp_get_num_threads() == 1);
     110         144 :   assert(matrix->nrows == redist->nrows);
     111        6384 :   for (int i = 0; i < matrix->nrows; i++) {
     112        6240 :     assert(matrix->row_sizes[i] == redist->row_sizes[i]);
     113             :   }
     114         144 :   assert(matrix->ncols == redist->ncols);
     115        6384 :   for (int i = 0; i < matrix->ncols; i++) {
     116        6240 :     assert(matrix->col_sizes[i] == redist->col_sizes[i]);
     117             :   }
     118             : 
     119         144 :   assert(dbm_mpi_comms_are_similar(matrix->dist->comm, redist->dist->comm));
     120         144 :   const dbm_mpi_comm_t comm = redist->dist->comm;
     121         144 :   const int nranks = dbm_mpi_comm_size(comm);
     122             : 
     123             :   // 1st pass: Compute send_count.
     124         144 :   int send_count[nranks];
     125         144 :   memset(send_count, 0, nranks * sizeof(int));
     126         288 :   for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
     127         144 :     dbm_shard_t *shard = &matrix->shards[ishard];
     128        8936 :     for (int iblock = 0; iblock < shard->nblocks; iblock++) {
     129        8792 :       const dbm_block_t *blk = &shard->blocks[iblock];
     130        8792 :       const int row_size = matrix->row_sizes[blk->row];
     131        8792 :       const int col_size = matrix->col_sizes[blk->col];
     132        8792 :       const int block_size = row_size * col_size;
     133        8792 :       const int rank = dbm_get_stored_coordinates(redist, blk->row, blk->col);
     134        8792 :       assert(0 <= rank && rank < nranks);
     135        8792 :       send_count[rank] += 2 + block_size;
     136             :     }
     137             :   }
     138             : 
     139             :   // 1st communication: Exchange counts.
     140         144 :   int recv_count[nranks];
     141         144 :   dbm_mpi_alltoall_int(send_count, 1, recv_count, 1, comm);
     142             : 
     143             :   // Compute displacements and allocate data buffers.
     144         144 :   int send_displ[nranks + 1], recv_displ[nranks + 1];
     145         144 :   send_displ[0] = recv_displ[0] = 0;
     146         432 :   for (int irank = 1; irank < nranks + 1; irank++) {
     147         288 :     send_displ[irank] = send_displ[irank - 1] + send_count[irank - 1];
     148         288 :     recv_displ[irank] = recv_displ[irank - 1] + recv_count[irank - 1];
     149             :   }
     150         144 :   const int total_send_count = send_displ[nranks];
     151         144 :   const int total_recv_count = recv_displ[nranks];
     152         144 :   double *data_send = malloc(total_send_count * sizeof(double));
     153         144 :   double *data_recv = malloc(total_recv_count * sizeof(double));
     154             : 
     155             :   // 2nd pass: Fill send_data.
     156         144 :   int send_data_positions[nranks];
     157         144 :   memcpy(send_data_positions, send_displ, nranks * sizeof(int));
     158         288 :   for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
     159         144 :     dbm_shard_t *shard = &matrix->shards[ishard];
     160        8936 :     for (int iblock = 0; iblock < shard->nblocks; iblock++) {
     161        8792 :       const dbm_block_t *blk = &shard->blocks[iblock];
     162        8792 :       const double *blk_data = &shard->data[blk->offset];
     163        8792 :       const int row_size = matrix->row_sizes[blk->row];
     164        8792 :       const int col_size = matrix->col_sizes[blk->col];
     165        8792 :       const int block_size = row_size * col_size;
     166        8792 :       const int rank = dbm_get_stored_coordinates(redist, blk->row, blk->col);
     167        8792 :       const int pos = send_data_positions[rank];
     168        8792 :       data_send[pos + 0] = blk->row; // send integers as doubles
     169        8792 :       data_send[pos + 1] = blk->col;
     170        8792 :       memcpy(&data_send[pos + 2], blk_data, block_size * sizeof(double));
     171        8792 :       send_data_positions[rank] += 2 + block_size;
     172             :     }
     173             :   }
     174         432 :   for (int irank = 0; irank < nranks; irank++) {
     175         288 :     assert(send_data_positions[irank] == send_displ[irank + 1]);
     176             :   }
     177             : 
     178             :   // 2nd communication: Exchange data.
     179         144 :   dbm_mpi_alltoallv_double(data_send, send_count, send_displ, data_recv,
     180             :                            recv_count, recv_displ, comm);
     181         144 :   free(data_send);
     182             : 
     183             :   // 3rd pass: Unpack data.
     184         144 :   dbm_clear(redist);
     185         144 :   int recv_data_pos = 0;
     186        8936 :   while (recv_data_pos < total_recv_count) {
     187        8792 :     const int row = (int)data_recv[recv_data_pos + 0];
     188        8792 :     const int col = (int)data_recv[recv_data_pos + 1];
     189        8792 :     assert(data_recv[recv_data_pos + 0] - (double)row == 0.0);
     190        8792 :     assert(data_recv[recv_data_pos + 1] - (double)col == 0.0);
     191        8792 :     dbm_put_block(redist, row, col, false, &data_recv[recv_data_pos + 2]);
     192        8792 :     const int row_size = matrix->row_sizes[row];
     193        8792 :     const int col_size = matrix->col_sizes[col];
     194        8792 :     const int block_size = row_size * col_size;
     195        8792 :     recv_data_pos += 2 + block_size;
     196             :   }
     197         144 :   assert(recv_data_pos == total_recv_count);
     198         144 :   free(data_recv);
     199         144 : }
     200             : 
     201             : /*******************************************************************************
     202             :  * \brief Looks up a block from given matrics. This routine is thread-safe.
     203             :  *        If the block is not found then a null pointer is returned.
     204             :  * \author Ole Schuett
     205             :  ******************************************************************************/
     206    17482199 : void dbm_get_block_p(dbm_matrix_t *matrix, const int row, const int col,
     207             :                      double **block, int *row_size, int *col_size) {
     208    17482199 :   assert(0 <= row && row < matrix->nrows);
     209    17482199 :   assert(0 <= col && col < matrix->ncols);
     210    17482199 :   assert(dbm_get_stored_coordinates(matrix, row, col) == matrix->dist->my_rank);
     211    17482199 :   *row_size = matrix->row_sizes[row];
     212    17482199 :   *col_size = matrix->col_sizes[col];
     213    17482199 :   *block = NULL;
     214             : 
     215    17482199 :   const int ishard = dbm_get_shard_index(matrix, row, col);
     216    17482199 :   const dbm_shard_t *shard = &matrix->shards[ishard];
     217    17482199 :   dbm_block_t *blk = dbm_shard_lookup(shard, row, col);
     218    17482199 :   if (blk != NULL) {
     219    16053194 :     *block = &shard->data[blk->offset];
     220             :   }
     221    17482199 : }
     222             : 
     223             : /*******************************************************************************
     224             :  * \brief Adds a block to given matrix. This routine is thread-safe.
     225             :  *        If block already exist then it gets overwritten (or summed).
     226             :  * \author Ole Schuett
     227             :  ******************************************************************************/
     228    27865223 : void dbm_put_block(dbm_matrix_t *matrix, const int row, const int col,
     229             :                    const bool summation, const double *block) {
     230    27865223 :   assert(0 <= row && row < matrix->nrows);
     231    27865223 :   assert(0 <= col && col < matrix->ncols);
     232    27865223 :   assert(dbm_get_stored_coordinates(matrix, row, col) == matrix->dist->my_rank);
     233    27865223 :   const int row_size = matrix->row_sizes[row];
     234    27865223 :   const int col_size = matrix->col_sizes[col];
     235    27865223 :   const int block_size = row_size * col_size;
     236             : 
     237    27865223 :   const int ishard = dbm_get_shard_index(matrix, row, col);
     238    27865223 :   dbm_shard_t *shard = &matrix->shards[ishard];
     239    27865223 :   omp_set_lock(&shard->lock);
     240    27865223 :   dbm_block_t *blk =
     241    27865223 :       dbm_shard_get_or_allocate_block(shard, row, col, block_size);
     242    27865223 :   double *blk_data = &shard->data[blk->offset];
     243    27865223 :   if (summation) {
     244  1104768087 :     for (int i = 0; i < block_size; i++) {
     245  1101751802 :       blk_data[i] += block[i];
     246             :     }
     247             :   } else {
     248    24848938 :     memcpy(blk_data, block, block_size * sizeof(double));
     249             :   }
     250    27865223 :   omp_unset_lock(&shard->lock);
     251    27865223 : }
     252             : 
     253             : /*******************************************************************************
     254             :  * \brief Remove all blocks from matrix, but does not release underlying memory.
     255             :  * \author Ole Schuett
     256             :  ******************************************************************************/
     257     1395192 : void dbm_clear(dbm_matrix_t *matrix) {
     258     1395192 :   assert(omp_get_num_threads() == 1);
     259             : 
     260     1395192 : #pragma omp parallel for schedule(dynamic)
     261             :   for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
     262             :     dbm_shard_t *shard = &matrix->shards[ishard];
     263             :     shard->nblocks = 0;
     264             :     shard->data_size = 0;
     265             :     shard->data_promised = 0;
     266             :     // Does not deallocate memory, hence data_allocated remains unchanged.
     267             :     memset(shard->hashtable, 0, shard->hashtable_size * sizeof(int));
     268             :   }
     269     1395192 : }
     270             : 
     271             : /*******************************************************************************
     272             :  * \brief Removes all blocks from the matrix whose norm is below the threshold.
     273             :  *        Blocks of size zero are always kept.
     274             :  * \author Ole Schuett
     275             :  ******************************************************************************/
     276      425987 : void dbm_filter(dbm_matrix_t *matrix, const double eps) {
     277      425987 :   assert(omp_get_num_threads() == 1);
     278             : 
     279      425987 :   if (eps == 0.0) {
     280             :     return;
     281             :   }
     282             : 
     283      419181 : #pragma omp parallel for schedule(dynamic)
     284             :   for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
     285             :     dbm_shard_t *shard = &matrix->shards[ishard];
     286             :     const int old_nblocks = shard->nblocks;
     287             :     shard->nblocks = 0;
     288             :     shard->data_promised = 0;
     289             :     memset(shard->hashtable, 0, shard->hashtable_size * sizeof(int));
     290             : 
     291             :     for (int iblock = 0; iblock < old_nblocks; iblock++) {
     292             :       const dbm_block_t old_blk = shard->blocks[iblock];
     293             :       const double *old_blk_data = &shard->data[old_blk.offset];
     294             :       const int row_size = matrix->row_sizes[old_blk.row];
     295             :       const int col_size = matrix->col_sizes[old_blk.col];
     296             :       const int block_size = row_size * col_size;
     297             :       double norm = 0.0; // Compute norm as double, but compare as float.
     298             :       for (int i = 0; i < block_size; i++) {
     299             :         norm += old_blk_data[i] * old_blk_data[i];
     300             :       }
     301             :       // For historic reasons zero-sized blocks are never filtered.
     302             :       if (sqrt((float)norm) < eps && block_size > 0) {
     303             :         continue; // filter the block
     304             :       }
     305             :       // Re-create block.
     306             :       dbm_block_t *new_blk = dbm_shard_promise_new_block(
     307             :           shard, old_blk.row, old_blk.col, block_size);
     308             :       assert(new_blk->offset <= old_blk.offset);
     309             :       if (new_blk->offset != old_blk.offset) {
     310             :         // Using memmove instead of memcpy because it handles overlap correctly.
     311             :         double *new_blk_data = &shard->data[new_blk->offset];
     312             :         memmove(new_blk_data, old_blk_data, block_size * sizeof(double));
     313             :       }
     314             :     }
     315             :     shard->data_size = shard->data_promised;
     316             :     // TODO: Could call realloc to release excess memory.
     317             :   }
     318             : }
     319             : 
     320             : /*******************************************************************************
     321             :  * \brief Adds list of blocks efficiently. The blocks will be filled with zeros.
     322             :  *        This routine must always be called within an OpenMP parallel region.
     323             :  * \author Ole Schuett
     324             :  ******************************************************************************/
     325     1045521 : void dbm_reserve_blocks(dbm_matrix_t *matrix, const int nblocks,
     326             :                         const int rows[], const int cols[]) {
     327     1045521 :   assert(omp_get_num_threads() == omp_get_max_threads() &&
     328             :          "Please call dbm_reserve_blocks within an OpenMP parallel region.");
     329     1045521 :   const int my_rank = matrix->dist->my_rank;
     330    24327930 :   for (int i = 0; i < nblocks; i++) {
     331    23282409 :     const int ishard = dbm_get_shard_index(matrix, rows[i], cols[i]);
     332    23282409 :     dbm_shard_t *shard = &matrix->shards[ishard];
     333    23282409 :     omp_set_lock(&shard->lock);
     334    23282409 :     assert(0 <= rows[i] && rows[i] < matrix->nrows);
     335    23282409 :     assert(0 <= cols[i] && cols[i] < matrix->ncols);
     336    23282409 :     assert(dbm_get_stored_coordinates(matrix, rows[i], cols[i]) == my_rank);
     337    23282409 :     const int row_size = matrix->row_sizes[rows[i]];
     338    23282409 :     const int col_size = matrix->col_sizes[cols[i]];
     339    23282409 :     const int block_size = row_size * col_size;
     340    23282409 :     dbm_shard_get_or_promise_block(shard, rows[i], cols[i], block_size);
     341    23282409 :     omp_unset_lock(&shard->lock);
     342             :   }
     343     1045521 : #pragma omp barrier
     344             : 
     345             : #pragma omp for schedule(dynamic)
     346     1045521 :   for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
     347     1045521 :     dbm_shard_t *shard = &matrix->shards[ishard];
     348     1045521 :     dbm_shard_allocate_promised_blocks(shard);
     349             :   }
     350     1045521 : }
     351             : 
     352             : /*******************************************************************************
     353             :  * \brief Multiplies all entries in the given matrix by the given factor alpha.
     354             :  * \author Ole Schuett
     355             :  ******************************************************************************/
     356      261613 : void dbm_scale(dbm_matrix_t *matrix, const double alpha) {
     357      261613 :   assert(omp_get_num_threads() == 1);
     358      261613 :   if (alpha == 1.0) {
     359             :     return;
     360             :   }
     361      185159 :   if (alpha == 0.0) {
     362      176263 :     dbm_zero(matrix);
     363      176263 :     return;
     364             :   }
     365             : 
     366        8896 : #pragma omp parallel for schedule(dynamic)
     367             :   for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
     368             :     dbm_shard_t *shard = &matrix->shards[ishard];
     369             :     for (int i = 0; i < shard->data_size; i++) {
     370             :       shard->data[i] *= alpha;
     371             :     }
     372             :   }
     373             : }
     374             : 
     375             : /*******************************************************************************
     376             :  * \brief Sets all blocks in the given matrix to zero.
     377             :  * \author Ole Schuett
     378             :  ******************************************************************************/
     379      176263 : void dbm_zero(dbm_matrix_t *matrix) {
     380      176263 :   assert(omp_get_num_threads() == 1);
     381             : 
     382      176263 : #pragma omp parallel for schedule(dynamic)
     383             :   for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
     384             :     dbm_shard_t *shard = &matrix->shards[ishard];
     385             :     memset(shard->data, 0, shard->data_size * sizeof(double));
     386             :   }
     387      176263 : }
     388             : 
     389             : /*******************************************************************************
     390             :  * \brief Adds matrix_b to matrix_a.
     391             :  * \author Ole Schuett
     392             :  ******************************************************************************/
     393      125667 : void dbm_add(dbm_matrix_t *matrix_a, const dbm_matrix_t *matrix_b) {
     394      125667 :   assert(omp_get_num_threads() == 1);
     395      125667 :   assert(matrix_a->dist == matrix_b->dist);
     396             : 
     397      125667 : #pragma omp parallel for schedule(dynamic)
     398             :   for (int ishard = 0; ishard < dbm_get_num_shards(matrix_b); ishard++) {
     399             :     dbm_shard_t *shard_a = &matrix_a->shards[ishard];
     400             :     const dbm_shard_t *shard_b = &matrix_b->shards[ishard];
     401             :     for (int iblock = 0; iblock < shard_b->nblocks; iblock++) {
     402             :       const dbm_block_t blk_b = shard_b->blocks[iblock];
     403             : 
     404             :       const int row_size = matrix_b->row_sizes[blk_b.row];
     405             :       const int col_size = matrix_b->col_sizes[blk_b.col];
     406             :       assert(row_size == matrix_a->row_sizes[blk_b.row]);
     407             :       assert(col_size == matrix_a->col_sizes[blk_b.col]);
     408             :       const int block_size = row_size * col_size;
     409             :       dbm_block_t *blk_a = dbm_shard_get_or_allocate_block(
     410             :           shard_a, blk_b.row, blk_b.col, block_size);
     411             :       double *data_a = &shard_a->data[blk_a->offset];
     412             :       const double *data_b = &shard_b->data[blk_b.offset];
     413             :       for (int i = 0; i < block_size; i++) {
     414             :         data_a[i] += data_b[i];
     415             :       }
     416             :     }
     417             :   }
     418      125667 : }
     419             : 
     420             : /*******************************************************************************
     421             :  * \brief Creates an iterator for the blocks of the given matrix.
     422             :  *        The iteration order is not stable.
     423             :  *        This routine must always be called within an OpenMP parallel region.
     424             :  * \author Ole Schuett
     425             :  ******************************************************************************/
     426     2486837 : void dbm_iterator_start(dbm_iterator_t **iter_out, const dbm_matrix_t *matrix) {
     427     2486837 :   assert(omp_get_num_threads() == omp_get_max_threads() &&
     428             :          "Please call dbm_iterator_start within an OpenMP parallel region.");
     429     2486837 :   dbm_iterator_t *iter = malloc(sizeof(dbm_iterator_t));
     430     2486837 :   iter->matrix = matrix;
     431     2486837 :   iter->next_block = 0;
     432     2486837 :   iter->next_shard = omp_get_thread_num();
     433     2931255 :   while (iter->next_shard < dbm_get_num_shards(matrix) &&
     434     2486837 :          matrix->shards[iter->next_shard].nblocks == 0) {
     435      444418 :     iter->next_shard += omp_get_num_threads();
     436             :   }
     437     2486837 :   assert(*iter_out == NULL);
     438     2486837 :   *iter_out = iter;
     439     2486837 : }
     440             : 
     441             : /*******************************************************************************
     442             :  * \brief Returns number of blocks the iterator will provide to calling thread.
     443             :  * \author Ole Schuett
     444             :  ******************************************************************************/
     445      491958 : int dbm_iterator_num_blocks(const dbm_iterator_t *iter) {
     446      491958 :   int num_blocks = 0;
     447      491958 :   for (int ishard = omp_get_thread_num();
     448      983916 :        ishard < dbm_get_num_shards(iter->matrix);
     449      491958 :        ishard += omp_get_num_threads()) {
     450      491958 :     num_blocks += iter->matrix->shards[ishard].nblocks;
     451             :   }
     452      491958 :   return num_blocks;
     453             : }
     454             : 
     455             : /*******************************************************************************
     456             :  * \brief Tests whether the given iterator has any block left.
     457             :  * \author Ole Schuett
     458             :  ******************************************************************************/
     459    46564361 : bool dbm_iterator_blocks_left(const dbm_iterator_t *iter) {
     460    46564361 :   return iter->next_shard < dbm_get_num_shards(iter->matrix);
     461             : }
     462             : 
     463             : /*******************************************************************************
     464             :  * \brief Returns the next block from the given iterator.
     465             :  * \author Ole Schuett
     466             :  ******************************************************************************/
     467    52132688 : void dbm_iterator_next_block(dbm_iterator_t *iter, int *row, int *col,
     468             :                              double **block, int *row_size, int *col_size) {
     469    52132688 :   const dbm_matrix_t *matrix = iter->matrix;
     470    52132688 :   assert(iter->next_shard < dbm_get_num_shards(matrix));
     471    52132688 :   const dbm_shard_t *shard = &matrix->shards[iter->next_shard];
     472    52132688 :   assert(iter->next_block < shard->nblocks);
     473    52132688 :   dbm_block_t *blk = &shard->blocks[iter->next_block];
     474             : 
     475    52132688 :   *row = blk->row;
     476    52132688 :   *col = blk->col;
     477    52132688 :   *row_size = matrix->row_sizes[blk->row];
     478    52132688 :   *col_size = matrix->col_sizes[blk->col];
     479    52132688 :   *block = &shard->data[blk->offset];
     480             : 
     481    52132688 :   iter->next_block++;
     482    52132688 :   if (iter->next_block >= shard->nblocks) {
     483             :     // Advance to the next non-empty shard...
     484     2042419 :     iter->next_shard += omp_get_num_threads();
     485     2042419 :     while (iter->next_shard < dbm_get_num_shards(matrix) &&
     486           0 :            matrix->shards[iter->next_shard].nblocks == 0) {
     487           0 :       iter->next_shard += omp_get_num_threads();
     488             :     }
     489     2042419 :     iter->next_block = 0; // ...and continue with its first block.
     490             :   }
     491    52132688 : }
     492             : 
     493             : /*******************************************************************************
     494             :  * \brief Releases the given iterator.
     495             :  * \author Ole Schuett
     496             :  ******************************************************************************/
     497     2486837 : void dbm_iterator_stop(dbm_iterator_t *iter) { free(iter); }
     498             : 
     499             : /*******************************************************************************
     500             :  * \brief Computes a checksum of the given matrix.
     501             :  * \author Ole Schuett
     502             :  ******************************************************************************/
     503         190 : double dbm_checksum(const dbm_matrix_t *matrix) {
     504         190 :   double checksum = 0.0;
     505         380 :   for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
     506         190 :     const dbm_shard_t *shard = &matrix->shards[ishard];
     507    14383658 :     for (int i = 0; i < shard->data_size; i++) {
     508    14383468 :       checksum += shard->data[i] * shard->data[i];
     509             :     }
     510             :   }
     511         190 :   dbm_mpi_sum_double(&checksum, 1, matrix->dist->comm);
     512         190 :   return checksum;
     513             : }
     514             : 
     515             : /*******************************************************************************
     516             :  * \brief Returns the absolute value of the larges element of the entire matrix.
     517             :  * \author Ole Schuett
     518             :  ******************************************************************************/
     519          48 : double dbm_maxabs(const dbm_matrix_t *matrix) {
     520          48 :   double maxabs = 0.0;
     521          96 :   for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
     522          48 :     dbm_shard_t *shard = &matrix->shards[ishard];
     523     1599050 :     for (int i = 0; i < shard->data_size; i++) {
     524     1599002 :       maxabs = fmax(maxabs, fabs(shard->data[i]));
     525             :     }
     526             :   }
     527          48 :   dbm_mpi_max_double(&maxabs, 1, matrix->dist->comm);
     528          48 :   return maxabs;
     529             : }
     530             : 
     531             : /*******************************************************************************
     532             :  * \brief Returns the name of the matrix of the given matrix.
     533             :  * \author Ole Schuett
     534             :  ******************************************************************************/
     535     1086610 : const char *dbm_get_name(const dbm_matrix_t *matrix) { return matrix->name; }
     536             : 
     537             : /*******************************************************************************
     538             :  * \brief Returns the number of local Non-Zero Elements of the given matrix.
     539             :  * \author Ole Schuett
     540             :  ******************************************************************************/
     541     1322990 : int dbm_get_nze(const dbm_matrix_t *matrix) {
     542     1322990 :   int nze = 0;
     543     2645980 :   for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
     544     1322990 :     nze += matrix->shards[ishard].data_size;
     545             :   }
     546     1322990 :   return nze;
     547             : }
     548             : 
     549             : /*******************************************************************************
     550             :  * \brief Returns the number of local blocks of the given matrix.
     551             :  * \author Ole Schuett
     552             :  ******************************************************************************/
     553      825008 : int dbm_get_num_blocks(const dbm_matrix_t *matrix) {
     554      825008 :   int nblocks = 0;
     555     1650016 :   for (int ishard = 0; ishard < dbm_get_num_shards(matrix); ishard++) {
     556      825008 :     nblocks += matrix->shards[ishard].nblocks;
     557             :   }
     558      825008 :   return nblocks;
     559             : }
     560             : 
     561             : /*******************************************************************************
     562             :  * \brief Returns the row block sizes of the given matrix.
     563             :  * \author Ole Schuett
     564             :  ******************************************************************************/
     565     3688572 : void dbm_get_row_sizes(const dbm_matrix_t *matrix, int *nrows,
     566             :                        const int **row_sizes) {
     567     3688572 :   *nrows = matrix->nrows;
     568     3688572 :   *row_sizes = matrix->row_sizes;
     569     3688572 : }
     570             : 
     571             : /*******************************************************************************
     572             :  * \brief Returns the column block sizes of the given matrix.
     573             :  * \author Ole Schuett
     574             :  ******************************************************************************/
     575     2831352 : void dbm_get_col_sizes(const dbm_matrix_t *matrix, int *ncols,
     576             :                        const int **col_sizes) {
     577     2831352 :   *ncols = matrix->ncols;
     578     2831352 :   *col_sizes = matrix->col_sizes;
     579     2831352 : }
     580             : 
     581             : /*******************************************************************************
     582             :  * \brief Returns the local row block sizes of the given matrix.
     583             :  * \author Ole Schuett
     584             :  ******************************************************************************/
     585      163452 : void dbm_get_local_rows(const dbm_matrix_t *matrix, int *nlocal_rows,
     586             :                         const int **local_rows) {
     587      163452 :   *nlocal_rows = matrix->dist->rows.nlocals;
     588      163452 :   *local_rows = matrix->dist->rows.local_indicies;
     589      163452 : }
     590             : 
     591             : /*******************************************************************************
     592             :  * \brief Returns the local column block sizes of the given matrix.
     593             :  * \author Ole Schuett
     594             :  ******************************************************************************/
     595       73518 : void dbm_get_local_cols(const dbm_matrix_t *matrix, int *nlocal_cols,
     596             :                         const int **local_cols) {
     597       73518 :   *nlocal_cols = matrix->dist->cols.nlocals;
     598       73518 :   *local_cols = matrix->dist->cols.local_indicies;
     599       73518 : }
     600             : 
     601             : /*******************************************************************************
     602             :  * \brief Returns the MPI rank on which the given block should be stored.
     603             :  * \author Ole Schuett
     604             :  ******************************************************************************/
     605    74805979 : int dbm_get_stored_coordinates(const dbm_matrix_t *matrix, const int row,
     606             :                                const int col) {
     607    74805979 :   return dbm_distribution_stored_coords(matrix->dist, row, col);
     608             : }
     609             : 
     610             : /*******************************************************************************
     611             :  * \brief Returns the distribution of the given matrix.
     612             :  * \author Ole Schuett
     613             :  ******************************************************************************/
     614      868777 : const dbm_distribution_t *dbm_get_distribution(const dbm_matrix_t *matrix) {
     615      868777 :   return matrix->dist;
     616             : }
     617             : 
     618             : // EOF