PaStiX Handbook  6.4.0
bvec_zlapmr.c
Go to the documentation of this file.
1 /**
2  *
3  * @file bvec_zlapmr.c
4  *
5  * Functions to compute the permutation of the right hand side.
6  *
7  * @copyright 2004-2024 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
8  * Univ. Bordeaux. All rights reserved.
9  *
10  * @version 6.4.0
11  * @author Mathieu Faverge
12  * @author Pierre Ramet
13  * @author Xavier Lacoste
14  * @author Gregoire Pichon
15  * @author Theophile Terraz
16  * @author Tony Delarue
17  * @author Vincent Bridonneau
18  * @date 2024-07-05
19  * @generated from /builds/solverstack/pastix/bcsc/bvec_zlapmr.c, normal z -> z, Tue Oct 8 14:17:55 2024
20  *
21  * This file implements the function bvec_zlapmr with the following hierarchy:
22  *
23  * bvec_zlapmr():
24  * - bvec_zlapmr_shm() for shared memory case
25  * - bvec_zlapmr_rep() for replicated rhs case
26  * - bvec_zlapmr_rep_{vec2bvec,bvec2vec}()
27  * - bvec_zlapmr_dst() for distributed rhs case
28  * - bvec_zlapmr_dst_{vec2bvec,bvec2vec}()
29  *
30  **/
31 #include "common.h"
32 #include <math.h>
33 #include "lapacke.h"
34 #include "bcsc/bcsc.h"
35 #include "bcsc/bvec.h"
36 #include "bcsc_z.h"
37 #include "order/order_internal.h"
38 #include "cblas.h"
39 #include "blend/solver.h"
40 #include "spm/src/common.h"
41 
42 #if defined( PASTIX_WITH_MPI )
43 /**
44  *******************************************************************************
45  *
46  * @ingroup bcsc_internal
47  *
48  * @brief Applies a row permutation (permtab) to the right hand side b and stores it
49  * in Pb in the distributed case. It also sends and receives the part of pb
50  * according to the bcsc partition.
51  *
52  *******************************************************************************
53  *
54  * @param[in] pastix_data
55  * The pastix_data structure.
56  *
57  * @param[in] m
58  * The number of rows in the right hand side b.
59  *
60  * @param[in] nrhs
61  * The number of columns in the right hand side b.
62  *
63  * @param[in] b
64  * A right hand side of size ldb-by-n.
65  *
66  * @param[in] ldb
67  * The leading dimension of b >= m.
68  *
69  * @param[inout] Pb
70  * The structure of the permuted right hand side b.
71  * On entry, the structure is initialized. On exit, contains the
72  * permuted right hand side b.
73  *
74  *******************************************************************************
75  *
76  * @retval PASTIX_SUCCESS
77  *
78  *******************************************************************************/
79 static inline int
80 bvec_zlapmr_dst_vec2bvec( pastix_data_t *pastix_data,
81  pastix_int_t m,
82  pastix_int_t nrhs,
83  pastix_complex64_t *b,
84  pastix_int_t ldb,
85  pastix_rhs_t Pb )
86 {
87  pastix_complex64_t *pb = Pb->b;
88  const spmatrix_t *spm = pastix_data->csc;
89  pastix_bcsc_t *bcsc = pastix_data->bcsc;
90  const SolverMatrix *solvmatr = pastix_data->solvmatr;
91  pastix_int_t baseval_spm = spm->baseval;
92  pastix_int_t dof = spm->dof;
93  pastix_int_t *dofs = spm->dofs;
94  pastix_int_t *loc2glob = spm->loc2glob;
95  pastix_int_t clustnbr = solvmatr->clustnbr;
96  pastix_int_t bcsc_n = bcsc->n;
97  bvec_handle_comm_t *comm_rhs = NULL;
98  bvec_proc_comm_t *data_comm = NULL;
99  pastix_int_t il, ile, ilpe, ig, dofi;
100  pastix_int_t j, c;
101  pastix_int_t *idx_cnt, *val_cnt;
102  pastix_complex64_t *values_c;
103 
104  assert( m == spm->nexp );
105  assert( Pb->m == bcsc_n );
106  assert( Pb->allocated == 1 );
107  assert( Pb->ld == bcsc_n );
108  assert( b != pb );
109 
110  bvec_handle_comm_init( pastix_data, Pb );
111 
112  bvec_exchange_amount_dst( pastix_data, PastixDirForward, m, nrhs, Pb );
113  comm_rhs = Pb->rhs_comm;
114 
115  /* Allocates the indexes and values buffers. */
116  bvec_zallocate_buf_dst( comm_rhs );
117  data_comm = comm_rhs->data_comm;
118 
119  /*
120  * Allocates and initialises the counters used to fill rhs_comm->values
121  * and rhs_comm->indexes.
122  */
123  MALLOC_INTERN( idx_cnt, clustnbr, pastix_int_t );
124  MALLOC_INTERN( val_cnt, clustnbr, pastix_int_t );
125  memset( idx_cnt, 0, clustnbr * sizeof(pastix_int_t) );
126  memset( val_cnt, 0, clustnbr * sizeof(pastix_int_t) );
127 
128  /*
129  * Goes through b to fill the data_comm with the data to send and
130  * fills pb with the local data.
131  */
132  for ( il = 0, ile = 0; ile < m; ile += dofi, il++ ) {
133  ig = loc2glob[ il ] - baseval_spm;
134  ilpe = bvec_glob2Ploc( pastix_data, ig );
135  dofi = ( dof > 0 ) ? dof : dofs[ ig+1 ] - dofs[ ig ];
136 
137  if ( ilpe < 0 ) {
138  c = - ( ilpe + 1 );
139  data_comm = comm_rhs->data_comm + c;
140 
141  /* Stores the indexes to send to c: (ipe, j). */
142  data_comm->send_idxbuf[ idx_cnt[ c ] ] = ig;
143  idx_cnt[ c ] ++;
144 
145  /* Stores the value to send to c. */
146  for ( j = 0; j < nrhs; j++ ) {
147  values_c = ((pastix_complex64_t*)(data_comm->send_valbuf)) + val_cnt[ c ];
148  memcpy( values_c, b + ile + j * ldb, dofi * sizeof(pastix_complex64_t) );
149  val_cnt[ c ] += dofi;
150  }
151  }
152  else {
153  for ( j = 0; j < nrhs; j++ ) {
154  memcpy( pb + ilpe + j * bcsc_n, b + ile + j * ldb, dofi * sizeof(pastix_complex64_t) );
155  }
156  }
157  }
158 
159  bvec_zexchange_data_dst( pastix_data, PastixDirForward, nrhs, b, ldb, Pb, NULL );
160 
161  bvec_handle_comm_exit( Pb->rhs_comm );
162 
163  memFree_null( idx_cnt );
164  memFree_null( val_cnt );
165  return PASTIX_SUCCESS;
166 }
167 
168 /**
169  *******************************************************************************
170  *
171  * @ingroup bcsc
172  *
173  * @brief Applies a row permutation (peritab) to the right hand side pb and stores it in
174  * b in the distributed case. It also sends and receives the part of b according
175  * to the spm repartition.
176  *
177  *******************************************************************************
178  *
179  * @param[in] pastix_data
180  * The pastix_data structure.
181  *
182  * @param[in] m
183  * The number of rows in the right hand side b.
184  *
185  * @param[in] nrhs
186  * The number of columns in the right hand side b.
187  *
188  * @param[inout] b
189  * A right hand side of size ldb-by-n.
190  * On entry, the allocated right hand side.
191  * At exit, contains the reverse permutation of Pb.
192  *
193  * @param[in] ldb
194  * The leading dimension of b >= m.
195  *
196  * @param[inout] Pb
197  * The structure of the permuted right hand side b.
198  *
199  *******************************************************************************
200  *
201  * @retval PASTIX_SUCCESS
202  *
203  *******************************************************************************/
204 static inline int
205 bvec_zlapmr_dst_bvec2vec( pastix_data_t *pastix_data,
206  pastix_int_t m,
207  pastix_int_t nrhs,
208  pastix_complex64_t *b,
209  pastix_int_t ldb,
210  pastix_rhs_t Pb )
211 {
212  pastix_complex64_t *bp = Pb->b;
213  const spmatrix_t *spm = pastix_data->csc;
214  pastix_int_t dof = spm->dof;
215  bvec_handle_comm_t *comm_rhs = NULL;
216  pastix_int_t clustnbr = pastix_data->solvmatr->clustnbr;
217  pastix_int_t bcsc_n = Pb->m;
218  pastix_int_t *Ploc2Pglob = Pb->Ploc2Pglob;
219  pastix_int_t *glob2loc = NULL;
220  pastix_int_t ilp, ilpe, igp, ile;
221  pastix_int_t j, c, dofi;
222  bvec_proc_comm_t *data_comm;
223  pastix_int_t *idx_cnt, *val_cnt;
224  pastix_complex64_t *values_c;
225 
226  assert( Pb->m == pastix_data->bcsc->n );
227  assert( m == spm->nexp );
228  assert( b != NULL );
229  memset( b, 0, ldb * nrhs * sizeof( pastix_complex64_t) );
230 
231  bvec_compute_Ploc2Pglob( pastix_data, Pb );
232  Ploc2Pglob = Pb->Ploc2Pglob;
233 
234  if ( spm->glob2loc == NULL ) {
235  spmatrix_t spmcpy;
236 
237  memcpy( &spmcpy, spm, sizeof( spmatrix_t ) );
238  glob2loc = spm_get_glob2loc( &spmcpy );
239  }
240  else {
241  glob2loc = spm->glob2loc;
242  }
243  bvec_exchange_amount_dst( pastix_data, PastixDirBackward, m, nrhs, Pb );
244  comm_rhs = Pb->rhs_comm;
245 
246  /* Allocates the indexes and values buffers. */
247  bvec_zallocate_buf_dst( comm_rhs );
248  data_comm = comm_rhs->data_comm;
249 
250  /*
251  * Allocates and initialises the counters used to fill bcsc_comm->values
252  * and bcsc_comm->indexes.
253  */
254  MALLOC_INTERN( idx_cnt, clustnbr, pastix_int_t );
255  MALLOC_INTERN( val_cnt, clustnbr, pastix_int_t );
256  memset( idx_cnt, 0, clustnbr * sizeof(pastix_int_t) );
257  memset( val_cnt, 0, clustnbr * sizeof(pastix_int_t) );
258 
259  /*
260  * Goes through b to fill the data_comm with the data to send and
261  * fills pb with the local data.
262  */
263  ilp = 0;
264  dofi = dof; /* vdof incorrect */
265  for ( ilpe = 0; ilpe < bcsc_n; ilpe += dofi, ilp ++ ) {
266  igp = Ploc2Pglob[ ilp ];
267  ile = bvec_Pglob2loc( pastix_data, glob2loc, igp );
268 
269  if ( ile < 0 ) {
270  c = - ( ile + 1 );
271  data_comm = comm_rhs->data_comm + c;
272 
273  /* Stores the indexes to send to c: (ipe, j). */
274  data_comm->send_idxbuf[ idx_cnt[ c ] ] = igp;
275  idx_cnt[ c ] ++;
276  /* Stores the value to send to c. */
277  for ( j = 0; j < nrhs; j++ ) {
278  values_c = ((pastix_complex64_t*)(data_comm->send_valbuf)) + val_cnt[ c ];
279  memcpy( values_c, bp + ilpe + j * Pb->ld, dofi * sizeof(pastix_complex64_t) );
280  val_cnt[ c ] += dofi;
281  }
282  }
283  else {
284  for ( j = 0; j < nrhs; j++ ) {
285  memcpy( b + ile + j * ldb, bp + ilpe + j * Pb->ld, dofi * sizeof(pastix_complex64_t) );
286  }
287  }
288  }
289 
290  bvec_zexchange_data_dst( pastix_data, PastixDirBackward, nrhs, b, ldb, Pb, glob2loc );
291 
292  bvec_handle_comm_exit( Pb->rhs_comm );
293 
294  if ( spm->glob2loc == NULL ) {
295  free( glob2loc );
296  }
297  memFree_null( idx_cnt );
298  memFree_null( val_cnt );
299  return PASTIX_SUCCESS;
300 }
301 
302 /**
303  *
304  *******************************************************************************
305  *
306  * @ingroup bcsc_internal
307  *
308  * @brief Apply a row permutation to a right hand side A (LAPACK xlatmr) in the
309  * distributed case.
310  *
311  *******************************************************************************
312  *
313  * @param[in] pastix_data
314  * The pastix_data structure.
315  *
316  * @param[in] dir
317  * The direction of the permutation.
318  * If PastixDirForward, A is permuted into PA.
319  * If PastixDirBackward, PA is permuted into A.
320  *
321  * @param[in] m
322  * The number of rows in the right hand side A, and the number of elements in
323  * perm.
324  *
325  * @param[in] n
326  * The number of columns in the right hand side A.
327  *
328  * @param[inout] A
329  * A right hand side of size lda-by-n.
330  * Referenced as input if dir is PastixDirForward, as output otherwise.
331  *
332  * @param[in] lda
333  * The leading dimension of A.
334  *
335  * @param[inout] PA
336  * The structure of the permuted right hand side A.
337  * Referenced as inout if dir is PastixDirForward, as input otherwise.
338  *
339  *******************************************************************************
340  *
341  * @retval PASTIX_SUCCESS
342  *
343  *******************************************************************************/
344 static inline int
345 bvec_zlapmr_dst( pastix_data_t *pastix_data,
346  pastix_dir_t dir,
347  pastix_int_t m,
348  pastix_int_t n,
349  pastix_complex64_t *A,
350  pastix_int_t lda,
351  pastix_rhs_t PA )
352 {
353  if ( dir == PastixDirForward ) {
354  return bvec_zlapmr_dst_vec2bvec( pastix_data, m, n, A, lda, PA );
355  }
356  else {
357  return bvec_zlapmr_dst_bvec2vec( pastix_data, m, n, A, lda, PA );
358  }
359 }
360 
361 /**
362  *******************************************************************************
363  *
364  * @ingroup bcsc_internal
365  *
366  * @brief Applies a row permutation (permtab) to the right hand side b and stores it
367  * in Pb in the replicated case.
368  *
369  *******************************************************************************
370  *
371  * @param[in] pastix_data
372  * The pastix_data structure.
373  *
374  * @param[in] m
375  * The number of rows in the right hand side b.
376  *
377  * @param[in] nrhs
378  * The number of columns in the right hand side b.
379  *
380  * @param[in] b
381  * A right hand side of size ldb-by-n.
382  *
383  * @param[in] ldb
384  * The leading dimension of b >= m.
385  *
386  * @param[inout] Pb
387  * The structure of the permuted right hand side b.
388  * On entry, the structure is initialized. On exit, contains the
389  * permuted right hand side b.
390  *
391  *******************************************************************************
392  *
393  * @retval PASTIX_SUCCESS
394  *
395  *******************************************************************************/
396 static inline int
397 bvec_zlapmr_rep_vec2bvec( const pastix_data_t *pastix_data,
398  pastix_int_t m,
399  pastix_int_t nrhs,
400  const pastix_complex64_t *b,
401  pastix_int_t ldb,
402  pastix_rhs_t Pb )
403 {
404  pastix_complex64_t *pb;
405  const spmatrix_t *spm = pastix_data->csc;
406  pastix_int_t dof = spm->dof;
407  const pastix_int_t *dofs = spm->dofs;
408  pastix_int_t ldpb = Pb->ld;
409  pastix_int_t idx_cnt = 0;
410  pastix_int_t val_cnt = 0;
411  pastix_int_t j, ig, ige, ilpe, dofi;
412  bvec_data_amount_t *data;
413  bvec_handle_comm_t *comm_rhs;
414 
415  /* Check on b */
416  assert( m == spm->gNexp );
417  assert( m <= ldb );
418  assert( m == spm->nexp );
419  assert( pastix_data->bcsc->n == Pb->m );
420  assert( pastix_data->bcsc->n == Pb->ld );
421 
422  /*
423  * Goes through b to fill the data_comm with the data to send and
424  * fills pb with the local data.
425  */
426  pb = Pb->b;
427  ige = 0;
428  for ( ig = 0; ig < spm->gN; ig++, ige += dofi ) {
429  dofi = ( dof > 0 ) ? dof : dofs[ ig+1 ] - dofs[ ig ];
430  ilpe = bvec_glob2Ploc( pastix_data, ig );
431 
432  if ( ilpe < 0 ) {
433  continue;
434  }
435 
436  idx_cnt += 1;
437  val_cnt += dofi * nrhs;
438  for ( j = 0; j < nrhs; j++ ) {
439  memcpy( pb + ilpe + j * ldpb,
440  b + ige + j * ldb,
441  dofi * sizeof(pastix_complex64_t) );
442  }
443  }
444 
445  assert( idx_cnt <= val_cnt );
446 
447  bvec_handle_comm_init( pastix_data, Pb );
448 
449  comm_rhs = Pb->rhs_comm;
450  data = &(comm_rhs->data_comm[comm_rhs->clustnum].nsends);
451  data->idxcnt = idx_cnt;
452  data->valcnt = val_cnt;
453 
454  bvec_exchange_amount_rep( Pb->rhs_comm );
455 
456  (void)m;
457  return PASTIX_SUCCESS;
458 }
459 
460 /**
461  *******************************************************************************
462  *
463  * @ingroup bcsc_internal
464  *
465  * @brief Applies a row permutation (permtab) to the right hand side b. and stores it
466  * in Pb.
467  *
468  *******************************************************************************
469  *
470  * @param[in] pastix_data
471  * The pastix_data structure.
472  *
473  * @param[in] m
474  * The number of rows in the right hand side b.
475  *
476  * @param[in] nrhs
477  * The number of columns in the right hand side b.
478  *
479  * @param[inout] b
480  * A right hand side of size ldb-by-n.
481  * On entry, the allocated right hand side.
482  * On exit, contains the revers permutation of Pb.
483  *
484  * @param[in] ldb
485  * The leading dimension of b >= m.
486  *
487  * @param[inout] Pb
488  * The structure of the permuted right hand side b.
489  *
490  *******************************************************************************/
491 static inline int
492 bvec_zlapmr_rep_bvec2vec( pastix_data_t *pastix_data,
493  pastix_int_t m,
494  pastix_int_t nrhs,
495  pastix_complex64_t *b,
496  pastix_int_t ldb,
497  pastix_rhs_t Pb )
498 {
499  pastix_complex64_t *bp = Pb->b;
500  bvec_handle_comm_t *comm_rhs = Pb->rhs_comm;
501  const pastix_order_t *ord = pastix_data->ordemesh;
502  const spmatrix_t *spm = pastix_data->csc;
503  pastix_int_t dof = spm->dof;
504  const pastix_int_t *dofs = spm->dofs;
505  pastix_int_t *Ploc2Pglob = Pb->Ploc2Pglob;
506  pastix_int_t clustnum = pastix_data->solvmatr->clustnum;
507  pastix_int_t bcsc_n = Pb->m;
508  pastix_int_t ilpe, ilp, igp, ig, ige;
509  pastix_int_t j, dofi;
510  bvec_proc_comm_t *data_send;
511  bvec_data_amount_t *sends;
512  pastix_int_t *idxptr;
513 
514  assert( Pb->m == pastix_data->bcsc->n );
515  assert( m == pastix_data->csc->nexp );
516  assert( b != NULL );
517 
518  bvec_compute_Ploc2Pglob( pastix_data, Pb );
519 
520  Ploc2Pglob = Pb->Ploc2Pglob;
521  data_send = comm_rhs->data_comm + clustnum;
522  sends = &( data_send->nsends );
523 
524  /*
525  * Allocates the sending indexes buffer.
526  */
527  if ( ( sends->idxcnt != 0 ) && ( data_send->send_idxbuf == NULL ) ) {
528  MALLOC_INTERN( data_send->send_idxbuf, sends->idxcnt, pastix_int_t );
529  }
530 
531  /*
532  * Sets the counters used to fill indexes and values buffers.
533  */
534  idxptr = data_send->send_idxbuf;
535 
536  /*
537  * Goes through b to fill the data_comm with the data to send and
538  * fills pb with the local data.
539  */
540  ilp = 0;
541  dofi = ( dof > 0 ) ? dof : dofs[ 1 ] - dofs[ 0 ];
542  for ( ilpe = 0; ilpe < bcsc_n; ilpe += dofi, ilp ++ ) {
543  igp = Ploc2Pglob[ ilp ];
544  ig = ord->peritab[ igp ];
545  ige = ( dof > 0 ) ? ig * dof : dofs[ ig ];
546  dofi = ( dof > 0 ) ? dof : dofs[ ig+1 ] - dofs[ ig ];
547 
548  /* Stores the indexes to send to c: (ig, j). */
549  *idxptr = ig;
550  idxptr++;
551 
552  /* Stores the values to send to c. */
553  for ( j = 0; j < nrhs; j++ ) {
554  memcpy( b + ige + j * ldb, bp + ilpe + j * Pb->ld, dofi * sizeof(pastix_complex64_t) );
555  }
556  }
557  assert( (idxptr - data_send->send_idxbuf) == sends->idxcnt );
558 
559  bvec_zexchange_data_rep( pastix_data, nrhs, b, ldb, Pb );
560 
561  bvec_handle_comm_exit( Pb->rhs_comm );
562 
563  (void)m;
564  return PASTIX_SUCCESS;
565 }
566 
567 /**
568  *******************************************************************************
569  *
570  * @ingroup bcsc_internal
571  *
572  * @brief Apply a row permutation to a right hand side A (LAPACK xlatmr) in the
573  * replicated case.
574  *
575  *******************************************************************************
576  *
577  * @param[in] pastix_data
578  * The pastix_data structure.
579  *
580  * @param[in] dir
581  * The direction of the permutation.
582  * If PastixDirForward, A is permuted into PA.
583  * If PastixDirBackward, PA is permuted into A.
584  *
585  * @param[in] m
586  * The number of rows in the right hand side A, and the number of elements in
587  * perm.
588  *
589  * @param[in] n
590  * The number of columns in the right hand side A.
591  *
592  * @param[inout] A
593  * A right hand side of size lda-by-n.
594  * Referenced as input if dir is PastixDirForward, as output otherwise.
595  *
596  * @param[in] lda
597  * The leading dimension of A.
598  *
599  * @param[inout] PA
600  * The structure of the permuted right hand side A.
601  * Referenced as inout if dir is PastixDirForward, as input otherwise.
602  *
603  *******************************************************************************
604  *
605  * @retval PASTIX_SUCCESS
606  *
607  *******************************************************************************/
608 static inline int
609 bvec_zlapmr_rep( pastix_data_t *pastix_data,
610  pastix_dir_t dir,
611  pastix_int_t m,
612  pastix_int_t n,
613  pastix_complex64_t *A,
614  pastix_int_t lda,
615  pastix_rhs_t PA )
616 {
617  if ( dir == PastixDirForward ) {
618  return bvec_zlapmr_rep_vec2bvec( pastix_data, m, n, A, lda, PA );
619  }
620  else {
621  return bvec_zlapmr_rep_bvec2vec( pastix_data, m, n, A, lda, PA );
622  }
623 }
624 #endif
625 
626 /**
627  *******************************************************************************
628  *
629  * @ingroup bcsc_internal
630  *
631  * @brief Apply a row permutation to a right hand side A (LAPACK xlatmr) in the shared
632  * memory case.
633  *
634  *******************************************************************************
635  *
636  * @param[in] pastix_data
637  * The pastix_data structure.
638  *
639  * @param[in] dir
640  * The direction of the permutation.
641  * If PastixDirForward, A is permuted into PA.
642  * If PastixDirBackward, PA is permuted into A.
643  *
644  * @param[in] m
645  * The number of rows in the right hand side A, and the number of elements in
646  * perm.
647  *
648  * @param[in] n
649  * The number of columns in the right hand side A.
650  *
651  * @param[inout] A
652  * A right hand side of size lda-by-n.
653  * Referenced as input if dir is PastixDirForward, as output otherwise.
654  *
655  * @param[in] lda
656  * The leading dimension of A.
657  *
658  * @param[inout] PA
659  * The structure of the permuted right hand side A.
660  * Referenced as inout if dir is PastixDirForward, as input otherwise.
661  *
662  *******************************************************************************
663  *
664  * @retval PASTIX_SUCCESS
665  *
666  *******************************************************************************/
667 static inline int
669  pastix_dir_t dir,
670  pastix_int_t m,
671  pastix_int_t n,
672  pastix_complex64_t *A,
673  pastix_int_t lda,
674  pastix_rhs_t PA )
675 {
676  pastix_complex64_t tmp;
677  pastix_int_t i, j, k, jj;
678  pastix_int_t *perm, *perm_cpy;
679  int thread_safe = pastix_data->iparm[IPARM_APPLYPERM_WS];
680 
681  if ( PA->b != A ) {
682  pastix_print_error( "Incorrect definition of the right hand side for in place permutation\n" );
684  }
685  assert( PA->allocated == 0 );
686  assert( PA->m == m );
687  assert( PA->n == n );
688  assert( PA->ld == lda );
689 
690  perm = orderGetExpandedPeritab( pastix_data->ordemesh, pastix_data->csc );
691  assert( perm != NULL );
692 
693  if ( thread_safe ) {
694  perm_cpy = malloc( m * sizeof(pastix_int_t) );
695  memcpy( perm_cpy, perm, m * sizeof(pastix_int_t) );
696  }
697  else {
698  perm_cpy = perm;
699  }
700 
701  if ( dir == PastixDirBackward ) {
702  for( k = 0; k < m; k++ ) {
703  i = k;
704  j = perm_cpy[i];
705 
706  /* Cycle already seen */
707  if ( j < 0 ) {
708  continue;
709  }
710 
711  /* Mark the i^th element as being seen */
712  perm_cpy[i] = -j-1;
713 
714  while( j != k ) {
715 
716  for( jj = 0; jj < n; jj++ ) {
717  tmp = A[j + jj * lda];
718  A[j + jj * lda] = A[k + jj * lda];
719  A[k + jj * lda] = tmp;
720  }
721 
722  i = j;
723  j = perm_cpy[i];
724  perm_cpy[i] = -j-1;
725 
726  assert( (j != i) && (j >= 0) );
727  }
728  }
729  }
730  else {
731  for( k = 0; k < m; k++ ) {
732  i = k;
733  j = perm_cpy[i];
734  perm_cpy[i] = -j-1;
735 
736  /* Cycle already seen */
737  if ( j < 0 ) {
738  continue;
739  }
740 
741  i = perm_cpy[j];
742 
743  /* Mark the i^th element as being seen */
744  while( i >= 0 ) {
745 
746  for( jj = 0; jj < n; jj++ ) {
747  tmp = A[j + jj * lda];
748  A[j + jj * lda] = A[i + jj * lda];
749  A[i + jj * lda] = tmp;
750  }
751 
752  perm_cpy[j] = -i-1;
753  j = i;
754  i = perm_cpy[j];
755 
756  assert( j != i );
757  }
758  }
759  }
760 
761  if ( thread_safe ) {
762  free( perm_cpy );
763  }
764  else {
765  /* Restore perm array */
766  for( k = 0; k < m; k++ ) {
767  assert(perm[k] < 0);
768  perm[k] = - perm[k] - 1;
769  }
770  }
771 
772  return PASTIX_SUCCESS;
773 }
774 
775 /**
776  *******************************************************************************
777  *
778  * @ingroup bcsc_internal
779  *
780  * @brief Apply a row permutation to a right hand side A (LAPACK xlatmr)
781  *
782  *******************************************************************************
783  *
784  * @param[in] pastix_data
785  * The pastix_data structure.
786  *
787  * @param[in] dir
788  * The direction of the permutation.
789  * If PastixDirForward, A is permuted into PA.
790  * If PastixDirBackward, PA is permuted into A.
791  *
792  * @param[in] m
793  * The number of rows in the right hand side A, and the number of elements in
794  * perm.
795  *
796  * @param[in] n
797  * The number of columns in the right hand side A.
798  *
799  * @param[inout] A
800  * A right hand side of size lda-by-n.
801  * Referenced as input if dir is PastixDirForward, as output otherwise.
802  *
803  * @param[in] lda
804  * The leading dimension of A.
805  *
806  * @param[inout] PA
807  * The structure of the permuted right hand side A.
808  * Referenced as inout if dir is PastixDirForward, as input otherwise.
809  *
810  *******************************************************************************
811  *
812  * @retval PASTIX_SUCCESS
813  *
814  *******************************************************************************/
815 int
817  pastix_dir_t dir,
818  pastix_int_t m,
819  pastix_int_t n,
820  pastix_complex64_t *A,
821  pastix_int_t lda,
822  pastix_rhs_t PA )
823 {
824  const spmatrix_t *spm = pastix_data->csc;
825  const pastix_bcsc_t *bcsc = pastix_data->bcsc;
826  const SolverMatrix *solvmatr = pastix_data->solvmatr;
827  int rc;
828 
829  assert( lda >= m );
830  if ( dir == PastixDirForward ) {
831  PA->flttype = PastixComplex64;
832  PA->m = bcsc->n;
833  PA->n = n;
834 
835  if ( solvmatr->clustnbr > 1 ) {
836  PA->allocated = 1;
837  PA->ld = PA->m;
838  PA->b = malloc( PA->ld * PA->n * pastix_size_of( PA->flttype ) );
839  }
840  else {
841  assert( m == PA->m );
842  PA->allocated = 0;
843  PA->ld = lda;
844  PA->b = A;
845  }
846  }
847 #if !defined(NDEBUG)
848  else {
849  assert( PA->allocated >= 0 );
850  assert( PA->flttype == PastixComplex64 );
851  assert( PA->m == bcsc->n );
852  assert( PA->n == n );
853 
854  if ( PA->allocated == 0 )
855  {
856  assert( PA->b == A );
857  assert( PA->ld == lda );
858  }
859  else {
860  assert( PA->b != A );
861  assert( PA->ld == PA->m );
862  }
863  }
864 #endif
865 
866 #if defined(PASTIX_WITH_MPI)
867  if ( spm->clustnbr > 1 ) {
868  if ( spm->loc2glob != NULL ) {
869  /*
870  * The input right hand side is distributed, we redispatch it following the
871  * ordering.
872  */
873  rc = bvec_zlapmr_dst( pastix_data, dir, m, n, A, lda, PA );
874  }
875  else {
876  /*
877  * The input right hand side is replicated, we extract or collect the data
878  * following the ordering.
879  */
880  rc = bvec_zlapmr_rep( pastix_data, dir, m, n, A, lda, PA );
881  }
882  }
883  else
884 #endif
885  {
886  /*
887  * We are in the shared memory case, the permutation is applied in place.
888  */
889  rc = bvec_zlapmr_shm( pastix_data, dir, m, n, A, lda, PA );
890  (void)spm;
891  }
892 
893  if ( dir == PastixDirBackward ) {
894  if ( PA->allocated > 0 ) {
895  memFree_null( PA->b );
896  }
897 
898  PA->allocated = -1;
899  PA->flttype = PastixPattern;
900  PA->m = -1;
901  PA->n = -1;
902  PA->ld = -1;
903  PA->b = NULL;
904 
905  if ( PA->rhs_comm != NULL ) {
906  memFree_null( PA->rhs_comm );
907  }
908  }
909 
910  return rc;
911 }
BEGIN_C_DECLS typedef int pastix_int_t
Definition: datatypes.h:51
static int bvec_zlapmr_shm(pastix_data_t *pastix_data, pastix_dir_t dir, pastix_int_t m, pastix_int_t n, pastix_complex64_t *A, pastix_int_t lda, pastix_rhs_t PA)
Apply a row permutation to a right hand side A (LAPACK xlatmr) in the shared memory case.
Definition: bvec_zlapmr.c:668
int bvec_zlapmr(pastix_data_t *pastix_data, pastix_dir_t dir, pastix_int_t m, pastix_int_t n, pastix_complex64_t *A, pastix_int_t lda, pastix_rhs_t PA)
Apply a row permutation to a right hand side A (LAPACK xlatmr)
Definition: bvec_zlapmr.c:816
pastix_int_t valcnt
Definition: bvec.h:42
pastix_int_t idxcnt
Definition: bvec.h:41
pastix_int_t * send_idxbuf
Definition: bvec.h:52
bvec_data_amount_t nsends
Definition: bvec.h:50
void * send_valbuf
Definition: bvec.h:53
pastix_int_t clustnum
Definition: bvec.h:68
bvec_proc_comm_t data_comm[1]
Definition: bvec.h:73
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.
@ PastixDirForward
Definition: api.h:513
@ PastixDirBackward
Definition: api.h:514
@ IPARM_APPLYPERM_WS
Definition: api.h:108
@ PASTIX_SUCCESS
Definition: api.h:367
@ PASTIX_ERR_BADPARAMETER
Definition: api.h:374
pastix_int_t * peritab
Definition: order.h:52
pastix_int_t * orderGetExpandedPeritab(pastix_order_t *ordeptr, const spmatrix_t *spm)
This routine expand the peritab array for multi-dof matrices.
Definition: order.c:494
Order structure.
Definition: order.h:47
SolverMatrix * solvmatr
Definition: pastixdata.h:103
pastix_int_t * Ploc2Pglob
Definition: pastixdata.h:164
pastix_order_t * ordemesh
Definition: pastixdata.h:98
pastix_int_t * iparm
Definition: pastixdata.h:70
pastix_int_t ld
Definition: pastixdata.h:160
bvec_handle_comm_t * rhs_comm
Definition: pastixdata.h:163
const spmatrix_t * csc
Definition: pastixdata.h:90
pastix_coeftype_t flttype
Definition: pastixdata.h:157
pastix_int_t m
Definition: pastixdata.h:158
pastix_bcsc_t * bcsc
Definition: pastixdata.h:102
pastix_int_t n
Definition: pastixdata.h:159
int8_t allocated
Definition: pastixdata.h:156
Main PaStiX data structure.
Definition: pastixdata.h:68
Main PaStiX RHS structure.
Definition: pastixdata.h:155
Solver column block structure.
Definition: solver.h:203