PaStiX Handbook  6.3.2
bvec.c
Go to the documentation of this file.
1 /**
2  *
3  * @file bvec.c
4  *
5  * @copyright 2004-2023 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
6  * Univ. Bordeaux. All rights reserved.
7  *
8  * @version 6.3.2
9  * @author Mathieu Faverge
10  * @author Pierre Ramet
11  * @author Vincent Bridonneau
12  * @author Alycia Lisito
13  * @date 2023-07-21
14  *
15  **/
16 #include "common.h"
17 #include <math.h>
18 #include "lapacke.h"
19 #include "bcsc/bcsc.h"
20 #include "bcsc/bvec.h"
21 #include "order/order_internal.h"
22 #include "cblas.h"
23 #include "blend/solver.h"
24 
25 /**
26  *******************************************************************************
27  *
28  * @ingroup bcsc
29  *
30  * @brief Allocate a vector
31  *
32  *******************************************************************************
33  *
34  * @param[in] size
35  * The size of the vector
36  *
37  *******************************************************************************
38  *
39  * @return The allocated vector
40  *
41  *******************************************************************************/
42 void *bvec_malloc( size_t size )
43 {
44  void *x = NULL;
45  MALLOC_INTERN(x, size, char);
46  return x;
47 }
48 
49 /**
50  *******************************************************************************
51  *
52  * @ingroup bcsc
53  *
54  * @brief Free a vector
55  *
56  *******************************************************************************
57  *
58  * @param[inout] x
59  * The vector to be free
60  *
61  *******************************************************************************/
62 void bvec_free( void *x )
63 {
64  memFree_null(x);
65 }
66 
67 #if defined( PASTIX_WITH_MPI )
68 /**
69  *******************************************************************************
70  *
71  * @ingroup bcsc_internal
72  *
73  * @brief Initializes the rhs_comm field of an RHS data structure.
74  *
75  *******************************************************************************
76  *
77  * @param[in] pastix_data
78  * The pastix_data structure.
79  *
80  * @param[inout] Pb
81  * On entry, the initialized pastix_rhs_t data structure but the rhs_comm
82  * field is NULL.
83  * On exit, rhs_comm is initialized.
84  *
85  *******************************************************************************
86  *
87  * @retval PASTIX_SUCCESS on successful exit.
88  *
89  *******************************************************************************/
90 int
91 bvec_handle_comm_init( const pastix_data_t *pastix_data,
92  pastix_rhs_t Pb )
93 {
94  const SolverMatrix *solvmtx = pastix_data->solvmatr;
95  bvec_handle_comm_t *rhs_comm = NULL;
96 
97  /* Initializes bvec_comm. */
98  if ( Pb->rhs_comm == NULL ) {
99  size_t size = sizeof(bvec_handle_comm_t) + (solvmtx->clustnbr-1) * sizeof(bvec_proc_comm_t);
100 
101  Pb->rhs_comm = (bvec_handle_comm_t *)malloc( size );
102  }
103 
104  rhs_comm = Pb->rhs_comm;
105  rhs_comm->clustnbr = solvmtx->clustnbr;
106  rhs_comm->clustnum = solvmtx->clustnum;
107  rhs_comm->comm = solvmtx->solv_comm;
108  rhs_comm->flttype = Pb->flttype;
109  rhs_comm->max_idx = 0;
110  rhs_comm->max_val = 0;
111 
112  memset( rhs_comm->data_comm, 0, rhs_comm->clustnbr * sizeof(bvec_proc_comm_t) );
113 
114  return PASTIX_SUCCESS;
115 }
116 
117 /**
118  *******************************************************************************
119  *
120  * @ingroup bcsc_internal
121  *
122  * @brief Cleans-up a bvec_handle_comm_t structure.
123  *
124  *******************************************************************************
125 
126  * @param[inout] rhs_comm
127  * On entry, the initialized bvec_handle_comm_t data structure.
128  * On exit, the structure is destroyed and should no longer be used.
129  *
130  *******************************************************************************
131  *
132  * @retval PASTIX_SUCCESS on successful exit.
133  *
134  *******************************************************************************/
135 int
136 bvec_handle_comm_exit( bvec_handle_comm_t *rhs_comm )
137 {
138  int c;
139  int clustnbr = rhs_comm->clustnbr;
140  bvec_proc_comm_t *data;
141 
142  for ( c = 0; c < clustnbr; c++ ) {
143  data = rhs_comm->data_comm + c;
144 
145  if ( data->send_idxbuf != NULL ) {
146  memFree_null( data->send_idxbuf );
147  }
148  if ( data->send_valbuf != NULL ) {
149  memFree_null( data->send_valbuf );
150  }
151  }
152  rhs_comm->max_idx = 0;
153  rhs_comm->max_val = 0;
154 
155  return PASTIX_SUCCESS;
156 }
157 
158 /**
159  *******************************************************************************
160  *
161  * @ingroup bcsc_internal
162  *
163  * @brief Converts the non permuted global index into the permuted local one.
164  *
165  *******************************************************************************
166  *
167  * @param[in] pastix_data
168  * The pastix_data structure.
169  *
170  * @param[in] ig
171  * Index global not permuted.
172  *
173  *******************************************************************************
174  *
175  * @retval Returns the local permuted index or c = -(p+1) with p the process to
176  * which the local permuted column/row belongs to.
177  *
178  *******************************************************************************/
180 bvec_glob2Ploc( const pastix_data_t *pastix_data,
181  pastix_int_t ig )
182 {
183  const spmatrix_t *spm = pastix_data->csc;
184  const pastix_bcsc_t *bcsc = pastix_data->bcsc;
185  const pastix_order_t *ord = pastix_data->ordemesh;
186  const SolverMatrix *solvmatr = pastix_data->solvmatr;
187  const pastix_int_t *col2cblk = bcsc->col2cblk;
188  pastix_int_t basespm = spm->baseval;
189  pastix_int_t dof = spm->dof;
190  const pastix_int_t *dofs = spm->dofs;
191  const SolverCblk *cblk;
192  pastix_int_t igp, igpe, ilpe, cblknum;
193 
194  assert( ord->baseval == 0 );
195 
196  igp = ord->permtab[ ig ];
197  igpe = ( dof > 0 ) ? igp * dof : dofs[ ig ] - basespm;/* vdof incorrect */
198  cblknum = col2cblk[ igpe ];
199  if ( cblknum >= 0 ) {
200  cblk = solvmatr->cblktab + cblknum;
201  ilpe = cblk->lcolidx + igpe - cblk->fcolnum;
202  }
203  else {
204  ilpe = cblknum;
205  }
206 
207  return ilpe;
208 }
209 
210 /**
211  *******************************************************************************
212  *
213  * @ingroup bcsc_internal
214  *
215  * @brief Converts the permuted global index into the non permuted local one.
216  *
217  *******************************************************************************
218  *
219  * @param[in] pastix_data
220  * The pastix_data structure.
221  *
222  * @param[in] glob2loc
223  * The glob2loc array associated to the user spm stored in pastix_data.
224  *
225  * @param[in] igp
226  * Index global permuted.
227  *
228  *******************************************************************************
229  *
230  * @retval Returns the local non permuted index or c = -(p+1) with p the process
231  * to which the local permuted column/row belongs to.
232  *
233  *******************************************************************************/
235 bvec_Pglob2loc( const pastix_data_t *pastix_data,
236  const pastix_int_t *glob2loc,
237  pastix_int_t igp )
238 {
239  const spmatrix_t *spm = pastix_data->csc;
240  const pastix_order_t *ord = pastix_data->ordemesh;
241  pastix_int_t basespm = spm->baseval;
242  pastix_int_t dof = spm->dof;
243  const pastix_int_t *dofs = spm->dofs;
244  pastix_int_t ig, il, ile;
245 
246  assert( ord->baseval == 0 );
247 
248  ig = ord->peritab[ igp ];
249  il = glob2loc[ig];
250  if ( il >= 0 ) {
251  ile = ( dof > 0 ) ? il * dof : dofs[ ig ] - basespm;/* vdof incorrect */
252  }
253  else {
254  ile = il;
255  }
256 
257  return ile;
258 }
259 
260 /**
261  *******************************************************************************
262  *
263  * @ingroup bcsc_internal
264  *
265  * @brief Computes the Ploc2Pglob array.
266  *
267  *******************************************************************************
268  *
269  * @param[in] pastix_data
270  * The pastix_data structure.
271  *
272  * @param[inout] Pb
273  * On entry, the initialized pastix_rhs_t data structure but the
274  * Ploc2Pglob field is NULL.
275  * On exit, Ploc2Pglob is computed.
276  *
277  *******************************************************************************
278  *
279  * @retval PASTIX_SUCCESS
280  *
281  *******************************************************************************/
282 int
283 bvec_compute_Ploc2Pglob( pastix_data_t *pastix_data,
284  pastix_rhs_t Pb )
285 {
286  const spmatrix_t *spm = pastix_data->csc;
287  pastix_int_t *col2cblk = pastix_data->bcsc->col2cblk;
288  pastix_int_t *Ploc2Pglob = NULL;
289  pastix_int_t dof = spm->dof;
290  pastix_int_t bcsc_n = ( dof > 0 ) ? Pb->m / dof : Pb->m ; // ToDo : mvdof
291  pastix_int_t kgpe, kgp, klp;
292 
293  if ( Pb->Ploc2Pglob != NULL ) {
294  memFree_null( Pb->Ploc2Pglob );
295  }
296 
297  /*
298  * Creates Ploc2Pglob with col2cblk: gives the local index corresponding to the
299  * global one (Ploc2Pglob[ig] = il).
300  */
301  if ( bcsc_n > 0 ) {
302  MALLOC_INTERN( Ploc2Pglob, bcsc_n, pastix_int_t );
303  klp = 0;
304  kgp = 0;
305  for ( kgpe = 0; kgpe < spm->gNexp; kgpe += dof, kgp ++ ) {
306  if ( col2cblk[ kgpe ] >= 0 ) {
307  Ploc2Pglob[ klp ] = kgp;
308  klp ++;
309  }
310  }
311  assert( klp == bcsc_n );
312  }
313  Pb->Ploc2Pglob = Ploc2Pglob;
314 
315  return PASTIX_SUCCESS;
316 }
317 
318 /**
319  *******************************************************************************
320  *
321  * @ingroup bcsc_internal
322  *
323  * @brief Exchanges the amount of data in the replicated case.
324  *
325  *******************************************************************************
326  *
327  * @param[inout] rhs_comm
328  * On entry the rhs_comm of the permuted right hand side initialised.
329  * At exit rhs_comm is filled with the amount of data exchanged.
330  *
331  *******************************************************************************
332  *
333  * @retval PASTIX_SUCCESS
334  *
335  *******************************************************************************/
336 int
337 bvec_exchange_amount_rep( bvec_handle_comm_t *rhs_comm )
338 {
339  bvec_proc_comm_t *data_comm = rhs_comm->data_comm;
340  pastix_int_t clustnbr = rhs_comm->clustnbr;
341  pastix_int_t clustnum = rhs_comm->clustnum;
342  pastix_int_t max_idx = 0;
343  pastix_int_t max_val = 0;
344  pastix_int_t idx_cnt, val_cnt;
345  pastix_int_t c;
346 
347  /* Receives the amount of indexes and values. */
348  for ( c = 0; c < clustnbr; c++ ) {
349  data_comm = rhs_comm->data_comm + c;
350  if ( c == clustnum ) {
351  MPI_Bcast( &(data_comm->nsends), 2, PASTIX_MPI_INT, c, rhs_comm->comm );
352  continue;
353  }
354 
355  MPI_Bcast( &(data_comm->nrecvs), 2, PASTIX_MPI_INT, c, rhs_comm->comm );
356 
357  idx_cnt = data_comm->nrecvs.idxcnt;
358  val_cnt = data_comm->nrecvs.valcnt;
359 
360  assert( idx_cnt <= val_cnt );
361  if ( idx_cnt == 0 ) {
362  assert( val_cnt == 0 );
363  }
364 
365  if ( max_idx < idx_cnt ) {
366  max_idx = idx_cnt;
367  }
368  if ( max_val < val_cnt ) {
369  max_val = val_cnt;
370  }
371  }
372 
373  assert( max_idx <= max_val );
374 
375  rhs_comm->max_idx = max_idx;
376  rhs_comm->max_val = max_val;
377 
378  return PASTIX_SUCCESS;
379 }
380 
381 /**
382  *******************************************************************************
383  *
384  * @ingroup bcsc_internal
385  *
386  * @brief Computes the amount of sending data in the distributed case.
387  *
388  *******************************************************************************
389  *
390  * @param[in] pastix_data
391  * The pastix_data structure.
392  *
393  * @param[in] m
394  * The number of rows of the right hand side b.
395  *
396  * @param[in] nrhs
397  * The number of columns in the right hand side b.
398  *
399  * @param[inout] Pb
400  * On entry the rhs structure initialised.
401  * At exit the field rhs_comm is filled with the amount of sending
402  * data.
403  *
404  *******************************************************************************
405  *
406  * @retval PASTIX_SUCCESS
407  *
408  *******************************************************************************/
409 static inline int
410 bvec_compute_amount_dst( const pastix_data_t *pastix_data,
411  pastix_int_t m,
412  pastix_int_t nrhs,
413  pastix_rhs_t Pb )
414 {
415  const spmatrix_t *spm = pastix_data->csc;
416  pastix_int_t baseval_spm = spm->baseval;
417  pastix_int_t dof = spm->dof;
418  pastix_int_t *dofs = spm->dofs;
419  pastix_int_t *loc2glob = spm->loc2glob;
420  bvec_handle_comm_t *comm_rhs = NULL;
421  bvec_proc_comm_t *data = NULL;
422  pastix_int_t il, ile, ilpe, ig, dofi, c;
423 
424  bvec_handle_comm_init( pastix_data, Pb );
425  comm_rhs = Pb->rhs_comm;
426 
427  /*
428  * Goes through b to fill the data_comm with the data to send and
429  * fills pb with the local data.
430  */
431  for ( il = 0, ile = 0; ile < m; ile += dofi, il++ ) {
432  ig = loc2glob[ il ] - baseval_spm;
433  ilpe = bvec_glob2Ploc( pastix_data, ig );
434  dofi = ( dof > 0 ) ? dof : dofs[ ig+1 ] - dofs[ ig ];
435 
436  if ( ilpe < 0 ) {
437  c = - ( ilpe + 1 );
438  data = comm_rhs->data_comm + c;
439 
440  data->nsends.idxcnt += 1;
441  data->nsends.valcnt += dofi * nrhs;
442  }
443  }
444 
445  return PASTIX_SUCCESS;
446 }
447 
448 /**
449  *******************************************************************************
450  *
451  * @ingroup bcsc
452  *
453  * @brief Computes the maximum size of the sending indexes and values buffers.
454  *
455  *******************************************************************************
456  *
457  * @param[inout] rhs_comm
458  * On entry the rhs_comm of the permuted right hand side initialized.
459  * At exit the fields max_idx and max_val of rhs_comm are updated.
460  *
461  *******************************************************************************
462  *
463  * @retval PASTIX_SUCCESS
464  *
465  *******************************************************************************/
466 static inline int
467 bvec_compute_max( bvec_handle_comm_t *rhs_comm )
468 {
469  bvec_proc_comm_t *data = NULL;
470  pastix_int_t clustnbr = rhs_comm->clustnbr;
471  pastix_int_t clustnum = rhs_comm->clustnum;
472  pastix_int_t max_idx = 0;
473  pastix_int_t max_val = 0;
474  pastix_int_t idx_cnt, val_cnt, c;
475 
476  /* Receives the amount of indexes and values. */
477  for ( c = 0; c < clustnbr; c++ ) {
478  data = rhs_comm->data_comm + c;
479  if ( c == clustnum ) {
480  continue;
481  }
482 
483  idx_cnt = data->nrecvs.idxcnt;
484  val_cnt = data->nrecvs.valcnt;
485 
486  if ( max_idx < idx_cnt ) {
487  max_idx = idx_cnt;
488  }
489  if ( max_val < val_cnt ) {
490  max_val = val_cnt;
491  }
492  }
493 
494  assert( max_idx <= max_val );
495 
496  rhs_comm->max_idx = max_idx;
497  rhs_comm->max_val = max_val;
498 
499  return PASTIX_SUCCESS;
500 }
501 
502 /**
503  *******************************************************************************
504  *
505  * @ingroup bcsc
506  *
507  * @brief Switches the amount of the sending and receiving data in the
508  * distributed case.
509  *
510  *******************************************************************************
511  *
512  * @param[inout] rhs_comm
513  * On entry the rhs_comm of the permuted right hand side initialized.
514  * At exit rhs_comm is updated with the right amount of sending and
515  * receiving data.
516  *
517  *******************************************************************************
518  *
519  * @retval PASTIX_SUCCESS
520  *
521  *******************************************************************************/
522 static inline int
523 bvec_switch_amount_dst( bvec_handle_comm_t *rhs_comm )
524 {
525  bvec_proc_comm_t *data = NULL;
526  pastix_int_t clustnbr = rhs_comm->clustnbr;
527  pastix_int_t clustnum = rhs_comm->clustnum;
528  pastix_int_t c, idx_tmp, val_tmp;
529 
530  /* Sends the same amout of data to all process. */
531  for ( c = 0; c < clustnbr; c ++ ) {
532 
533  data = rhs_comm->data_comm + c;
534 
535  if ( c == clustnum ) {
536  continue;
537  }
538 
539  idx_tmp = data->nsends.idxcnt;
540  data->nsends.idxcnt = data->nrecvs.idxcnt;
541  data->nrecvs.idxcnt = idx_tmp;
542 
543  val_tmp = data->nsends.valcnt;
544  data->nsends.valcnt = data->nrecvs.valcnt;
545  data->nrecvs.valcnt = val_tmp;
546 
547  }
548 
549  return PASTIX_SUCCESS;
550 }
551 
552 /**
553  *******************************************************************************
554  *
555  * @ingroup bcsc_internal
556  *
557  * @brief Exchanges the amount of data in the distributed case.
558  *
559  *******************************************************************************
560  *
561  * @param[in] pastix_data
562  * The pastix_data structure.
563  *
564  * @param[in] dir
565  * The direction of the permutation.
566  * If PastixDirForward, b is permuted into Pb.
567  * If PastixDirBackward, Pb is permuted into b.
568  *
569  * @param[in] m
570  * If dir == PastixDirForward:
571  * m is the number of rows of the right hand side b.
572  * If dir == PastixDirBackward:
573  * m is not used.
574  *
575  * @param[in] nrhs
576  * The number of columns in the right hand side b.
577  *
578  * @param[inout] Pb
579  * On entry the rhs structure initialised.
580  * At exit the field rhs_comm is filled with the amount of data
581  * exchanged.
582  *
583  *******************************************************************************
584  *
585  * @retval PASTIX_SUCCESS
586  *
587  *******************************************************************************/
588 int
589 bvec_exchange_amount_dst( pastix_data_t *pastix_data,
590  pastix_dir_t dir,
591  pastix_int_t m,
592  pastix_int_t nrhs,
593  pastix_rhs_t Pb )
594 {
595  bvec_handle_comm_t *rhs_comm = Pb->rhs_comm;
596  pastix_int_t clustnbr = rhs_comm->clustnbr;
597  pastix_int_t clustnum = rhs_comm->clustnum;
598  bvec_proc_comm_t *data_comm = NULL;
599  bvec_proc_comm_t *data_send = NULL;
600  bvec_proc_comm_t *data_recv = NULL;
601  pastix_int_t counter_req = 0;
602  MPI_Status statuses[(clustnbr-1)*2];
603  MPI_Request requests[(clustnbr-1)*2];
604  bvec_data_amount_t *sends, *recvs;
605  pastix_int_t c_send, c_recv, k;
606 
607  if ( dir == PastixDirBackward ) {
608  bvec_switch_amount_dst( Pb->rhs_comm );
609  bvec_compute_max( Pb->rhs_comm );
610  return PASTIX_SUCCESS;
611  }
612 
613  bvec_compute_amount_dst( pastix_data, m, nrhs, Pb );
614  rhs_comm = Pb->rhs_comm;
615  data_comm = rhs_comm->data_comm;
616 
617  /* Receives the amount of indexes and values. */
618  c_send = (clustnum+1) % clustnbr;
619  c_recv = (clustnum-1+clustnbr) % clustnbr;
620  for ( k = 0; k < clustnbr-1; k++ ) {
621  data_send = data_comm + c_send;
622  data_recv = data_comm + c_recv;
623  sends = &( data_send->nsends );
624  recvs = &( data_recv->nrecvs );
625  if ( c_send == clustnum ) {
626  continue;
627  }
628 
629  MPI_Irecv( recvs, 2, PASTIX_MPI_INT, c_recv, PastixTagAmount,
630  rhs_comm->comm, &requests[counter_req++] );
631 
632  MPI_Isend( sends, 2, PASTIX_MPI_INT, c_send, PastixTagAmount,
633  rhs_comm->comm, &requests[counter_req++] );
634 
635  c_send = (c_send+1) % clustnbr;
636  c_recv = (c_recv-1+clustnbr) % clustnbr;
637  }
638 
639  MPI_Waitall( counter_req, requests, statuses );
640 
641  bvec_compute_max( rhs_comm );
642 
643  return PASTIX_SUCCESS;
644 }
645 #endif
BEGIN_C_DECLS typedef int pastix_int_t
Definition: datatypes.h:51
pastix_int_t max_idx
Definition: bvec.h:71
pastix_int_t valcnt
Definition: bvec.h:42
pastix_int_t max_val
Definition: bvec.h:72
pastix_int_t idxcnt
Definition: bvec.h:41
pastix_int_t * send_idxbuf
Definition: bvec.h:52
bvec_data_amount_t nrecvs
Definition: bvec.h:51
bvec_data_amount_t nsends
Definition: bvec.h:50
pastix_int_t clustnbr
Definition: bvec.h:67
void * send_valbuf
Definition: bvec.h:53
pastix_int_t clustnum
Definition: bvec.h:68
PASTIX_Comm comm
Definition: bvec.h:69
pastix_coeftype_t flttype
Definition: bvec.h:70
bvec_proc_comm_t data_comm[1]
Definition: bvec.h:73
void bvec_free(void *x)
Free a vector.
Definition: bvec.c:62
struct bvec_proc_comm_s bvec_proc_comm_t
Informations of the data exchanged with other processes.
void * bvec_malloc(size_t size)
Allocate a vector.
Definition: bvec.c:42
struct bvec_handle_comm_s bvec_handle_comm_t
Structure to manage communications with distributed rhs.
Information about the amount of data exchanged to permute the pivots.
Definition: bvec.h:38
Structure to manage communications with distributed rhs.
Definition: bvec.h:66
Informations of the data exchanged with other processes.
Definition: bvec.h:49
enum pastix_dir_e pastix_dir_t
Direction.
@ PastixDirBackward
Definition: api.h:514
@ PASTIX_SUCCESS
Definition: api.h:367
pastix_int_t baseval
Definition: order.h:48
pastix_int_t * permtab
Definition: order.h:51
pastix_int_t * peritab
Definition: order.h:52
Order structure.
Definition: order.h:47
SolverMatrix * solvmatr
Definition: pastixdata.h:102
pastix_int_t * Ploc2Pglob
Definition: pastixdata.h:159
pastix_order_t * ordemesh
Definition: pastixdata.h:97
bvec_handle_comm_t * rhs_comm
Definition: pastixdata.h:158
const spmatrix_t * csc
Definition: pastixdata.h:89
pastix_coeftype_t flttype
Definition: pastixdata.h:152
pastix_int_t m
Definition: pastixdata.h:153
pastix_bcsc_t * bcsc
Definition: pastixdata.h:101
Main PaStiX data structure.
Definition: pastixdata.h:67
Main PaStiX RHS structure.
Definition: pastixdata.h:150
pastix_int_t lcolidx
Definition: solver.h:165
SolverCblk *restrict cblktab
Definition: solver.h:222
pastix_int_t fcolnum
Definition: solver.h:161
Solver column block structure.
Definition: solver.h:156
Solver column block structure.
Definition: solver.h:200