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  SolverCblk *cblk;
63  bcsc_cblk_t *blockcol;
64  pastix_int_t cblknum, bcscnum, iter, idxcol, nodeidx, colsize;
65 
66  bcsc->cscfnbr = solvmtx->cblknbr - solvmtx->faninnbr - solvmtx->recvnbr;
67  MALLOC_INTERN( bcsc->cscftab, bcsc->cscfnbr, bcsc_cblk_t );
68 
69  idxcol = 0;
70  bcscnum = 0;
71  cblk = solvmtx->cblktab;
72  blockcol = bcsc->cscftab;
73  for ( cblknum = 0; cblknum < solvmtx->cblknbr; cblknum++, cblk++ )
74  {
75  if ( cblk->cblktype & (CBLK_FANIN|CBLK_RECV) ) {
76  continue;
77  }
78 
79  blockcol->cblknum = cblknum;
80  blockcol->colnbr = cblk_colnbr( cblk );
81  assert( cblk->bcscnum == bcscnum );
82  MALLOC_INTERN( blockcol->coltab, blockcol->colnbr + 1, pastix_int_t );
83 
84  /*
85  * Works only for DoF constant
86  * TODO : rework this part to work in variadic dof
87  */
88  assert( cblk->fcolnum % dof == 0 );
89 
90  blockcol->coltab[0] = idxcol;
91  for ( iter=0; iter < blockcol->colnbr; iter++ )
92  {
93  nodeidx = ( cblk->fcolnum + (iter-iter%dof) ) / dof;
94 
95  colsize = (newcoltab[nodeidx+1] - newcoltab[nodeidx]) * dof;
96  blockcol->coltab[iter+1] = blockcol->coltab[iter] + colsize;
97  }
98 
99  idxcol = blockcol->coltab[blockcol->colnbr];
100 
101  blockcol++;
102  bcscnum++;
103  }
104  assert( (blockcol - bcsc->cscftab) == bcsc->cscfnbr );
105  assert( bcscnum == bcsc->cscfnbr );
106 
107  if ( idxcol > 0 ) {
108  MALLOC_INTERN( bcsc->rowtab, idxcol, pastix_int_t);
109  MALLOC_INTERN( bcsc->Lvalues, idxcol * pastix_size_of( bcsc->flttype ), char );
110  }
111  else {
112  bcsc->rowtab = NULL;
113  bcsc->Lvalues = NULL;
114  }
115  bcsc->Uvalues = NULL;
116 
117  return idxcol;
118 }
119 
120 /**
121  *******************************************************************************
122  *
123  * @ingroup bcsc_internal
124  *
125  * @brief Restore the coltab array
126  *
127  * Function to restore the coltab array when it has been modified to initialize
128  * the row and values arrays.
129  *
130  *******************************************************************************
131  *
132  * @param[inout] bcsc
133  * On entry, the bcsc to restore.
134  * On exit, the coltab array of the bcsc is restored to the correct
135  * indexes.
136  *
137  *******************************************************************************/
138 void
140 {
141  bcsc_cblk_t *blockcol;
142  pastix_int_t index, iter, idxcol, idxcoltmp;
143 
144  idxcol = 0;
145  blockcol = bcsc->cscftab;
146  for ( index=0; index<bcsc->cscfnbr; index++, blockcol++ )
147  {
148  for ( iter=0; iter <= blockcol->colnbr; iter++ )
149  {
150  idxcoltmp = blockcol->coltab[iter];
151  blockcol->coltab[iter] = idxcol;
152  idxcol = idxcoltmp;
153  }
154  }
155  return;
156 }
157 
158 /**
159  *******************************************************************************
160  *
161  * @brief Initialize the coltab of a centralized block csc matrix.
162  *
163  *******************************************************************************
164  *
165  * @param[in] spm
166  * The initial sparse matrix in the spm format.
167  *
168  * @param[in] ord
169  * The ordering that needs to be applied on the spm to generate the
170  * block csc.
171  *
172  * @param[in] solvmtx
173  * The solver matrix structure that describe the data distribution.
174  *
175  * @param[inout] bcsc
176  * On entry, the pointer to an allocated bcsc.
177  * On exit, the bcsc stores the input spm with the permutation applied
178  * and grouped accordingly to the distribution described in solvmtx.
179  *
180  *******************************************************************************
181  *
182  * @return The number of non zero unknowns in the matrix.
183  *
184  *******************************************************************************/
185 pastix_int_t
186 bcsc_init_global_coltab( const spmatrix_t *spm,
187  const pastix_order_t *ord,
188  const SolverMatrix *solvmtx,
189  pastix_bcsc_t *bcsc )
190 {
191  pastix_int_t valuesize, baseval;
192  pastix_int_t *globcol = NULL;
193  pastix_int_t *colptr = spm->colptr;
194  pastix_int_t *rowptr = spm->rowptr;
195  pastix_int_t dof = spm->dof;
196  int sym = (spm->mtxtype == SpmSymmetric) || (spm->mtxtype == SpmHermitian);
197 
198  bcsc->mtxtype = spm->mtxtype;
199  baseval = spm->colptr[0];
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 permutation 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->n == spm->gN );
211 
212  {
213  pastix_int_t *loc2glob;
214  pastix_int_t frow, lrow;
215  pastix_int_t k, j, ig, jg, ip, jp;
216 
217  loc2glob = spm->loc2glob;
218  for ( j=0; j<spm->n; j++, colptr++, loc2glob++ )
219  {
220  jg = (spm->loc2glob == NULL) ? j : *loc2glob - baseval;
221  jp = ord->permtab[jg];
222  frow = colptr[0] - baseval;
223  lrow = colptr[1] - baseval;
224 
225  globcol[jp] += lrow - frow;
226  assert( (lrow - frow) >= 0 );
227  if (sym) {
228  for ( k=frow; k<lrow; k++ )
229  {
230  ig = rowptr[k] - baseval;
231  if ( ig != jg ) {
232  ip = ord->permtab[ig];
233  globcol[ip]++;
234  }
235  }
236  }
237  }
238 
239  /* Compute displacements to update the colptr array */
240  {
241  pastix_int_t tmp, idx;
242 
243  idx = 0;
244  for (j=0; j<=spm->gN; j++)
245  {
246  tmp = globcol[j];
247  globcol[j] = idx;
248  idx += tmp;
249  }
250  }
251  }
252 
253  valuesize = bcsc_init_coltab( solvmtx, globcol, dof, bcsc );
254  memFree_null( globcol );
255 
256  return valuesize;
257 }
258 
259 /**
260  *******************************************************************************
261  *
262  * @brief Initialize a block csc.
263  *
264  *******************************************************************************
265  *
266  * @param[in] spm
267  * The initial sparse matrix in the spm format.
268  *
269  * @param[in] ord
270  * The ordering that needs to be applied on the spm to generate the
271  * block csc.
272  *
273  * @param[in] solvmtx
274  * The solver matrix structure that describe the data distribution.
275  *
276  * @param[in] initAt
277  * A flag to enable/disable the initialization of A'
278  *
279  * @param[inout] bcsc
280  * On entry, the pointer to an allocated bcsc.
281  * On exit, the bcsc stores the input spm with the permutation applied
282  * and grouped accordingly to the distribution described in solvmtx.
283  *
284  *******************************************************************************/
285 static inline void
286 bcsc_init( const spmatrix_t *spm,
287  const pastix_order_t *ord,
288  const SolverMatrix *solvmtx,
289  pastix_int_t initAt,
290  pastix_bcsc_t *bcsc )
291 {
292  pastix_int_t *col2cblk = NULL;
293 
294  bcsc->mtxtype = spm->mtxtype;
295  bcsc->flttype = spm->flttype;
296  bcsc->gN = spm->gN;
297  bcsc->n = spm->n;
298 
299  /*
300  * Initialize the col2cblk array. col2cblk[i] contains the cblk index of the
301  * i-th column. col2cblk[i] = -1 if not local.
302  */
303  {
304  SolverCblk *cblk = solvmtx->cblktab;
305  pastix_int_t cblknbr = solvmtx->cblknbr;
306  pastix_int_t j, cblknum;
307 
308 
309  MALLOC_INTERN( col2cblk, spm->gNexp, pastix_int_t );
310  memset( col2cblk, 0xff, spm->gNexp * sizeof(pastix_int_t) );
311 
312  for ( cblknum=0; cblknum<cblknbr; cblknum++, cblk++ )
313  {
314  if ( cblk->cblktype & (CBLK_FANIN|CBLK_RECV) ){
315  continue;
316  }
317  for ( j=cblk->fcolnum; j<=cblk->lcolnum; j++ )
318  {
319  col2cblk[j] = cblknum;
320  }
321  }
322  }
323 
324  /*
325  * Fill in the lower triangular part of the blocked csc with values and
326  * rows. The upper triangular part is done later if required through LU
327  * factorization.
328  */
329  switch( spm->flttype ) {
330  case SpmFloat:
331  bcsc_sinit( spm, ord, solvmtx, col2cblk, initAt, bcsc );
332  break;
333  case SpmDouble:
334  bcsc_dinit( spm, ord, solvmtx, col2cblk, initAt, bcsc );
335  break;
336  case SpmComplex32:
337  bcsc_cinit( spm, ord, solvmtx, col2cblk, initAt, bcsc );
338  break;
339  case SpmComplex64:
340  bcsc_zinit( spm, ord, solvmtx, col2cblk, initAt, bcsc );
341  break;
342  case SpmPattern:
343  default:
344  fprintf(stderr, "bcsc_init: Error unknown floating type for input spm\n");
345  }
346 
347  memFree_null(col2cblk);
348 }
349 
350 /**
351  *******************************************************************************
352  *
353  * @brief Initialize the block csc matrix.
354  *
355  * The block csc matrix is used to initialize the factorized matrix, and to
356  * perform the matvec operations in refinement.
357  *
358  *******************************************************************************
359  *
360  * @param[in] spm
361  * The initial sparse matrix in the spm format.
362  *
363  * @param[in] ord
364  * The ordering that needs to be applied on the spm to generate the
365  * block csc.
366  *
367  * @param[in] solvmtx
368  * The solver matrix structure that describe the data distribution.
369  *
370  * @param[in] initAt
371  * A flag to enable/disable the initialization of A'
372  *
373  * @param[inout] bcsc
374  * On entry, the pointer to an allocated bcsc.
375  * On exit, the bcsc stores the input spm with the permutation applied
376  * and grouped accordingly to the distribution described in solvmtx.
377  *
378  *******************************************************************************
379  *
380  * @return The time spent to initialize the bcsc structure.
381  *
382  *******************************************************************************/
383 double
384 bcscInit( const spmatrix_t *spm,
385  const pastix_order_t *ord,
386  const SolverMatrix *solvmtx,
387  pastix_int_t initAt,
388  pastix_bcsc_t *bcsc )
389 {
390  double time = 0.;
391 
392  assert( ord->baseval == 0 );
393  assert( ord->vertnbr == spm->n );
394 
395  clockStart(time);
396  bcsc_init( spm, ord, solvmtx, initAt, bcsc );
397  clockStop(time);
398 
399  return time;
400 }
401 
402 /**
403  *******************************************************************************
404  *
405  * @brief Cleanup the block csc structure but do not free the bcsc pointer.
406  *
407  *******************************************************************************
408  *
409  * @param[inout] bcsc
410  * The block csc matrix to free.
411  *
412  *******************************************************************************/
413 void
415 {
416  bcsc_cblk_t *cblk;
417  pastix_int_t i;
418 
419  if ( bcsc->cscftab == NULL ) {
420  return;
421  }
422 
423  for ( i=0, cblk=bcsc->cscftab; i < bcsc->cscfnbr; i++, cblk++ ) {
424  memFree_null( cblk->coltab );
425  }
426 
427  memFree_null( bcsc->cscftab );
428  memFree_null( bcsc->rowtab );
429 
430  if ( (bcsc->Uvalues != NULL) &&
431  (bcsc->Uvalues != bcsc->Lvalues) ) {
432  memFree_null( bcsc->Uvalues );
433  }
434 
435  memFree_null( bcsc->Lvalues );
436 }
solver.h
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_dinit
void bcsc_dinit(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:347
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:414
pastix_order_s::permtab
pastix_int_t * permtab
Definition: order.h:49
bcsc_cblk_s::coltab
pastix_int_t * coltab
Definition: bcsc.h:31
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
bcsc_init
static void bcsc_init(const spmatrix_t *spm, const pastix_order_t *ord, const SolverMatrix *solvmtx, pastix_int_t initAt, pastix_bcsc_t *bcsc)
Initialize a block csc.
Definition: bcsc.c:286
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
bcsc_zinit
void bcsc_zinit(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:347
pastix_order_s::vertnbr
pastix_int_t vertnbr
Definition: order.h:47
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_cblk_s
Compressed colptr format for the bcsc.
Definition: bcsc.h:28
bcsc_cinit
void bcsc_cinit(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:347
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:139
pastix_order_s::baseval
pastix_int_t baseval
Definition: order.h:46
bcsc_init_global_coltab
pastix_int_t bcsc_init_global_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:186
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_sinit
void bcsc_sinit(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:347
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:384
pastix_bcsc_s::gN
int gN
Definition: bcsc.h:38
pastix_bcsc_s::Uvalues
void * Uvalues
Definition: bcsc.h:46