PaStiX Handbook  6.2.1
bcsc.c
Go to the documentation of this file.
1 /**
2  *
3  * @file bcsc.c
4  *
5  * @copyright 2004-2021 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
6  * Univ. Bordeaux. All rights reserved.
7  *
8  * @version 6.2.0
9  * @author Mathieu Faverge
10  * @author Pierre Ramet
11  * @author Xavier Lacoste
12  * @author Theophile Terraz
13  * @author Tony Delarue
14  * @date 2021-03-30
15  *
16  **/
17 #include "common.h"
18 #include "pastix/order.h"
19 #include <spm.h>
20 #include "blend/solver.h"
21 #include "bcsc/bcsc.h"
22 
23 #include "bcsc/bcsc_z.h"
24 #include "bcsc/bcsc_c.h"
25 #include "bcsc/bcsc_d.h"
26 #include "bcsc/bcsc_s.h"
27 
28 /**
29  *******************************************************************************
30  *
31  * @ingroup bcsc_internal
32  *
33  * @brief Initialize the coltab of a block csc matrix.
34  *
35  *******************************************************************************
36  *
37  * @param[in] solvmtx
38  * The solver matrix structure that describe the data distribution.
39  *
40  * @param[in] newcoltab
41  * Array of size spm->gN+1. This array is global coltab with -1 for non
42  * local indexes.
43  *
44  * @param[in] dof
45  * The degree of freedom of each unknown.
46  *
47  * @param[inout] bcsc
48  * On entry, the pointer to an allocated bcsc.
49  * On exit, the bcsc stores the ininitialized coltab split per block.
50  *
51  *******************************************************************************
52  *
53  * @return The number of non zero unknowns in the matrix.
54  *
55  *******************************************************************************/
56 static inline pastix_int_t
57 bcsc_init_coltab( const SolverMatrix *solvmtx,
58  const pastix_int_t *newcoltab,
59  pastix_int_t dof,
60  pastix_bcsc_t *bcsc )
61 {
62  bcsc_cblk_t *blockcol;
63  pastix_int_t cblknum, bcscnum, iter, idxcol, nodeidx, colsize;
64 
65  SolverCblk *cblk = solvmtx->cblktab;
66  bcsc->cscfnbr = solvmtx->cblknbr - solvmtx->faninnbr - solvmtx->recvnbr;
67  MALLOC_INTERN( bcsc->cscftab, bcsc->cscfnbr, bcsc_cblk_t );
68 
69  idxcol = 0;
70  cblk = solvmtx->cblktab;
71  blockcol = bcsc->cscftab;
72  bcscnum = 0;
73  for (cblknum = 0; cblknum < solvmtx->cblknbr; cblknum++, cblk++)
74  {
75  pastix_int_t fcolnum = cblk->fcolnum;
76 
77  if ( cblk->cblktype & (CBLK_FANIN|CBLK_RECV) ) {
78  continue;
79  }
80 
81  blockcol->cblknum = cblknum;
82  blockcol->colnbr = cblk_colnbr( cblk );
83  assert( cblk->bcscnum == bcscnum );
84  MALLOC_INTERN( blockcol->coltab, blockcol->colnbr + 1, pastix_int_t );
85 
86  /* Works only for DoF constant */
87  assert( fcolnum % dof == 0 );
88 
89  blockcol->coltab[0] = idxcol;
90  for (iter=0; iter < blockcol->colnbr; iter++)
91  {
92  nodeidx = ( fcolnum + (iter-iter%dof) ) / dof;
93 
94  colsize = (newcoltab[nodeidx+1] - newcoltab[nodeidx]) * dof;
95  blockcol->coltab[iter+1] = blockcol->coltab[iter] + colsize;
96  }
97 
98  idxcol = blockcol->coltab[blockcol->colnbr];
99 
100  blockcol++;
101  bcscnum++;
102  }
103  assert( (blockcol - bcsc->cscftab) == bcsc->cscfnbr );
104  assert( bcscnum == bcsc->cscfnbr );
105 
106  if ( idxcol > 0 ) {
107  MALLOC_INTERN( bcsc->rowtab, idxcol, pastix_int_t);
108  MALLOC_INTERN( bcsc->Lvalues, idxcol * pastix_size_of( bcsc->flttype ), char );
109  }
110  else {
111  bcsc->rowtab = NULL;
112  bcsc->Lvalues = NULL;
113  }
114  bcsc->Uvalues = NULL;
115 
116  return idxcol;
117 }
118 
119 /**
120  *******************************************************************************
121  *
122  * @ingroup bcsc_internal
123  *
124  * @brief Restore the coltab array
125  *
126  * Function to restore the coltab array when it has been modified to initialize
127  * the row and values arrays.
128  *
129  *******************************************************************************
130  *
131  * @param[inout] bcsc
132  * On entry, the bcsc to restore.
133  * On exit, the coltab array of the bcsc is restored to the correct
134  * indexes.
135  *
136  *******************************************************************************/
137 void
139 {
140  bcsc_cblk_t *blockcol;
141  pastix_int_t index, iter, idxcol, idxcoltmp;
142 
143  idxcol = 0;
144  blockcol = bcsc->cscftab;
145  for (index=0; index<bcsc->cscfnbr; index++, blockcol++)
146  {
147  for (iter=0; iter <= blockcol->colnbr; iter++)
148  {
149  idxcoltmp = blockcol->coltab[iter];
150  blockcol->coltab[iter] = idxcol;
151  idxcol = idxcoltmp;
152  }
153  }
154  return;
155 }
156 
157 /**
158  *******************************************************************************
159  *
160  * @brief Initialize the coltab of a centralized block csc matrix.
161  *
162  *******************************************************************************
163  *
164  * @param[in] spm
165  * The initial sparse matrix in the spm format.
166  *
167  * @param[in] ord
168  * The ordering that needs to be applied on the spm to generate the
169  * block csc.
170  *
171  * @param[in] solvmtx
172  * The solver matrix structure that describe the data distribution.
173  *
174  * @param[inout] bcsc
175  * On entry, the pointer to an allocated bcsc.
176  * On exit, the bcsc stores the input spm with the permutation applied
177  * and grouped accordingly to the distribution described in solvmtx.
178  *
179  *******************************************************************************
180  *
181  * @return The number of non zero unknowns in the matrix.
182  *
183  *******************************************************************************/
184 pastix_int_t
185 bcsc_init_centralized_coltab( const spmatrix_t *spm,
186  const pastix_order_t *ord,
187  const SolverMatrix *solvmtx,
188  pastix_bcsc_t *bcsc )
189 {
190  pastix_int_t valuesize, baseval;
191  pastix_int_t *globcol = NULL;
192  pastix_int_t *colptr = spm->colptr;
193  pastix_int_t *rowptr = spm->rowptr;
194  int dof = spm->dof;
195  int sym = (spm->mtxtype == SpmSymmetric) || (spm->mtxtype == SpmHermitian);
196 
197  bcsc->mtxtype = spm->mtxtype;
198  baseval = spm->colptr[0];
199 
200  /*
201  * Allocate and initialize globcol that contains the number of elements in
202  * each column of the input matrix
203  * Globcol is equivalent to the classic colptr for the internal blocked
204  * csc. The blocked csc integrate the perumtation computed within order
205  * structure.
206  */
207  MALLOC_INTERN( globcol, spm->gN+1, pastix_int_t );
208  memset( globcol, 0, (spm->gN+1) * sizeof(pastix_int_t) );
209 
210  assert( spm->loc2glob == NULL );
211 
212  {
213  pastix_int_t itercol, newcol;
214 
215  for (itercol=0; itercol<spm->gN; itercol++)
216  {
217  pastix_int_t frow = colptr[itercol] - baseval;
218  pastix_int_t lrow = colptr[itercol+1] - baseval;
219  newcol = ord->permtab[itercol];
220  globcol[newcol] += lrow - frow;
221 
222  assert( (lrow - frow) >= 0 );
223  if (sym) {
224  pastix_int_t iterrow, newrow;
225 
226  for (iterrow=frow; iterrow<lrow; iterrow++)
227  {
228  pastix_int_t tmprow = rowptr[iterrow] - baseval;
229  if (tmprow != itercol) {
230  newrow = ord->permtab[tmprow];
231  globcol[newrow]++;
232  }
233  }
234  }
235  }
236 
237  /* Compute displacements to update the colptr array */
238  {
239  pastix_int_t tmp, idx;
240 
241  idx = 0;
242  for (itercol=0; itercol<=spm->gN; itercol++)
243  {
244  tmp = globcol[itercol];
245  globcol[itercol] = idx;
246  idx += tmp;
247  }
248  }
249  }
250 
251  valuesize = bcsc_init_coltab( solvmtx, globcol, dof, bcsc );
252  memFree_null( globcol );
253 
254  return valuesize;
255 }
256 
257 /**
258  *******************************************************************************
259  *
260  * @brief Initialize a centralized block csc when no MPI processes are involved.
261  *
262  *******************************************************************************
263  *
264  * @param[in] spm
265  * The initial sparse matrix in the spm format.
266  *
267  * @param[in] ord
268  * The ordering that needs to be applied on the spm to generate the
269  * block csc.
270  *
271  * @param[in] solvmtx
272  * The solver matrix structure that describe the data distribution.
273  *
274  * @param[in] initAt
275  * A flag to enable/disable the initialization of A'
276  *
277  * @param[inout] bcsc
278  * On entry, the pointer to an allocated bcsc.
279  * On exit, the bcsc stores the input spm with the permutation applied
280  * and grouped accordingly to the distribution described in solvmtx.
281  *
282  *******************************************************************************/
283 void
284 bcsc_init_centralized( const spmatrix_t *spm,
285  const pastix_order_t *ord,
286  const SolverMatrix *solvmtx,
287  pastix_int_t initAt,
288  pastix_bcsc_t *bcsc )
289 {
290  pastix_int_t itercol, itercblk;
291  pastix_int_t cblknbr = solvmtx->cblknbr;
292  pastix_int_t eltnbr = spm->gNexp;
293  pastix_int_t *col2cblk = NULL;
294 
295  bcsc->mtxtype = spm->mtxtype;
296  bcsc->flttype = spm->flttype;
297  bcsc->gN = spm->gN;
298  bcsc->n = spm->n;
299 
300  assert( spm->loc2glob == NULL );
301 
302  /*
303  * Initialize the col2cblk array. col2cblk[i] contains the cblk index of the
304  * i-th column. col2cblk[i] = -1 if not local.
305  */
306  {
307  SolverCblk *cblk = solvmtx->cblktab;
308 
309  MALLOC_INTERN( col2cblk, eltnbr, pastix_int_t );
310  for (itercol=0; itercol<eltnbr; itercol++)
311  {
312  col2cblk[itercol] = -1;
313  }
314 
315  for (itercblk=0; itercblk<cblknbr; itercblk++, cblk++)
316  {
317  if( cblk->cblktype & (CBLK_FANIN|CBLK_RECV) ){
318  continue;
319  }
320  for (itercol = cblk->fcolnum;
321  itercol <= cblk->lcolnum;
322  itercol++ )
323  {
324  col2cblk[itercol] = itercblk;
325  }
326  }
327  }
328 
329  /*
330  * Fill in the lower triangular part of the blocked csc with values and
331  * rows. The upper triangular part is done later if required through LU
332  * factorization.
333  */
334  switch( spm->flttype ) {
335  case SpmFloat:
336  bcsc_sinit_centralized( spm, ord, solvmtx, col2cblk, initAt, bcsc );
337  break;
338  case SpmDouble:
339  bcsc_dinit_centralized( spm, ord, solvmtx, col2cblk, initAt, bcsc );
340  break;
341  case SpmComplex32:
342  bcsc_cinit_centralized( spm, ord, solvmtx, col2cblk, initAt, bcsc );
343  break;
344  case SpmComplex64:
345  bcsc_zinit_centralized( spm, ord, solvmtx, col2cblk, initAt, bcsc );
346  break;
347  case SpmPattern:
348  default:
349  fprintf(stderr, "bcsc_init_centralized: Error unknown floating type for input spm\n");
350  }
351 
352  memFree_null(col2cblk);
353 }
354 
355 /**
356  *******************************************************************************
357  *
358  * @brief Initialize the block csc matrix.
359  *
360  * The block csc matrix is used to initialize the factorized matrix, and to
361  * perform the matvec operations in refinement.
362  *
363  *******************************************************************************
364  *
365  * @param[in] spm
366  * The initial sparse matrix in the spm format.
367  *
368  * @param[in] ord
369  * The ordering that needs to be applied on the spm to generate the
370  * block csc.
371  *
372  * @param[in] solvmtx
373  * The solver matrix structure that describe the data distribution.
374  *
375  * @param[in] initAt
376  * A flag to enable/disable the initialization of A'
377  *
378  * @param[inout] bcsc
379  * On entry, the pointer to an allocated bcsc.
380  * On exit, the bcsc stores the input spm with the permutation applied
381  * and grouped accordingly to the distribution described in solvmtx.
382  *
383  *******************************************************************************
384  *
385  * @return The time spent to initialize the bcsc structure.
386  *
387  *******************************************************************************/
388 double
389 bcscInit( const spmatrix_t *spm,
390  const pastix_order_t *ord,
391  const SolverMatrix *solvmtx,
392  pastix_int_t initAt,
393  pastix_bcsc_t *bcsc )
394 {
395  spmatrix_t *spmg;
396  double time = 0.;
397 
398  if ( spm->loc2glob == NULL ) {
399  spmg = (spmatrix_t *)spm;
400  }
401 #if defined(PASTIX_WITH_MPI)
402  else {
403  pastix_print(solvmtx->clustnum, 0, "bcscInit: the SPM has to be centralized for the moment\n");
404  spmg = spmGather( spm, -1 );
405  }
406 #endif
407 
408  assert( ord->baseval == 0 );
409  assert( ord->vertnbr == spmg->n );
410 
411  clockStart(time);
412 
413  bcsc_init_centralized( spmg, ord, solvmtx, initAt, bcsc );
414 
415  clockStop(time);
416 
417  if( spmg != spm ) {
418  spmExit(spmg);
419  memFree_null( spmg );
420  }
421 
422  return time;
423 }
424 
425 /**
426  *******************************************************************************
427  *
428  * @brief Cleanup the block csc structure but do not free the bcsc pointer.
429  *
430  *******************************************************************************
431  *
432  * @param[inout] bcsc
433  * The block csc matrix to free.
434  *
435  *******************************************************************************/
436 void
438 {
439  bcsc_cblk_t *cblk;
440  pastix_int_t i;
441 
442  if ( bcsc->cscftab == NULL ) {
443  return;
444  }
445 
446  for (i=0, cblk=bcsc->cscftab; i < bcsc->cscfnbr; i++, cblk++ ) {
447  memFree_null( cblk->coltab );
448  }
449 
450  memFree_null( bcsc->cscftab );
451  memFree_null( bcsc->rowtab );
452 
453  if ( (bcsc->Uvalues != NULL) &&
454  (bcsc->Uvalues != bcsc->Lvalues) ) {
455  memFree_null( bcsc->Uvalues );
456  }
457 
458  memFree_null( bcsc->Lvalues );
459 }
solver.h
bcsc_dinit_centralized
void bcsc_dinit_centralized(const spmatrix_t *spm, const pastix_order_t *ord, const SolverMatrix *solvmtx, const pastix_int_t *col2cblk, int initAt, pastix_bcsc_t *bcsc)
Initialize a centralize double block csc.
Definition: bcsc_dinit.c:488
cblk_colnbr
static pastix_int_t cblk_colnbr(const SolverCblk *cblk)
Compute the number of columns in a column block.
Definition: solver.h:247
bcsc_cinit_centralized
void bcsc_cinit_centralized(const spmatrix_t *spm, const pastix_order_t *ord, const SolverMatrix *solvmtx, const pastix_int_t *col2cblk, int initAt, pastix_bcsc_t *bcsc)
Initialize a centralize pastix_complex32_t block csc.
Definition: bcsc_cinit.c:488
bcsc_s.h
pastix_bcsc_s
Internal column block distributed CSC matrix.
Definition: bcsc.h:37
pastix_bcsc_s::cscftab
bcsc_cblk_t * cscftab
Definition: bcsc.h:43
bcscExit
void bcscExit(pastix_bcsc_t *bcsc)
Cleanup the block csc structure but do not free the bcsc pointer.
Definition: bcsc.c:437
pastix_order_s::permtab
pastix_int_t * permtab
Definition: order.h:49
bcsc_cblk_s::coltab
pastix_int_t * coltab
Definition: bcsc.h:31
bcsc_zinit_centralized
void bcsc_zinit_centralized(const spmatrix_t *spm, const pastix_order_t *ord, const SolverMatrix *solvmtx, const pastix_int_t *col2cblk, int initAt, pastix_bcsc_t *bcsc)
Initialize a centralize pastix_complex64_t block csc.
Definition: bcsc_zinit.c:488
solver_cblk_s
Solver column block structure.
Definition: solver.h:127
bcsc_init_coltab
static pastix_int_t bcsc_init_coltab(const SolverMatrix *solvmtx, const pastix_int_t *newcoltab, pastix_int_t dof, pastix_bcsc_t *bcsc)
Initialize the coltab of a block csc matrix.
Definition: bcsc.c:57
pastix_order_s
Order structure.
Definition: order.h:45
pastix_bcsc_s::Lvalues
void * Lvalues
Definition: bcsc.h:45
bcsc_cblk_s::colnbr
pastix_int_t colnbr
Definition: bcsc.h:29
solver_cblk_s::bcscnum
pastix_int_t bcscnum
Definition: solver.h:141
pastix_bcsc_s::n
int n
Definition: bcsc.h:39
pastix_order_s::vertnbr
pastix_int_t vertnbr
Definition: order.h:47
bcsc_sinit_centralized
void bcsc_sinit_centralized(const spmatrix_t *spm, const pastix_order_t *ord, const SolverMatrix *solvmtx, const pastix_int_t *col2cblk, int initAt, pastix_bcsc_t *bcsc)
Initialize a centralize float block csc.
Definition: bcsc_sinit.c:488
bcsc_cblk_s::cblknum
pastix_int_t cblknum
Definition: bcsc.h:30
pastix_bcsc_s::flttype
int flttype
Definition: bcsc.h:41
pastix_bcsc_s::mtxtype
int mtxtype
Definition: bcsc.h:40
pastix_bcsc_s::cscfnbr
pastix_int_t cscfnbr
Definition: bcsc.h:42
bcsc.h
bcsc_init_centralized_coltab
pastix_int_t bcsc_init_centralized_coltab(const spmatrix_t *spm, const pastix_order_t *ord, const SolverMatrix *solvmtx, pastix_bcsc_t *bcsc)
Initialize the coltab of a centralized block csc matrix.
Definition: bcsc.c:185
bcsc_cblk_s
Compressed colptr format for the bcsc.
Definition: bcsc.h:28
pastix_bcsc_s::rowtab
pastix_int_t * rowtab
Definition: bcsc.h:44
bcsc_d.h
bcsc_restore_coltab
void bcsc_restore_coltab(pastix_bcsc_t *bcsc)
Restore the coltab array.
Definition: bcsc.c:138
pastix_order_s::baseval
pastix_int_t baseval
Definition: order.h:46
bcsc_c.h
solver_cblk_s::cblktype
int8_t cblktype
Definition: solver.h:130
solver_cblk_s::fcolnum
pastix_int_t fcolnum
Definition: solver.h:132
order.h
bcsc_z.h
bcscInit
double bcscInit(const spmatrix_t *spm, const pastix_order_t *ord, const SolverMatrix *solvmtx, pastix_int_t initAt, pastix_bcsc_t *bcsc)
Initialize the block csc matrix.
Definition: bcsc.c:389
pastix_bcsc_s::gN
int gN
Definition: bcsc.h:38
bcsc_init_centralized
void bcsc_init_centralized(const spmatrix_t *spm, const pastix_order_t *ord, const SolverMatrix *solvmtx, pastix_int_t initAt, pastix_bcsc_t *bcsc)
Initialize a centralized block csc when no MPI processes are involved.
Definition: bcsc.c:284
pastix_bcsc_s::Uvalues
void * Uvalues
Definition: bcsc.h:46