PaStiX Handbook  6.2.1
symbol.c
Go to the documentation of this file.
1 /**
2  *
3  * @file symbol.c
4  *
5  * PaStiX symbol structure routines
6  *
7  * @copyright 2004-2021 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
8  * Univ. Bordeaux. All rights reserved.
9  *
10  * @version 6.2.0
11  * @author David Goudin
12  * @author Francois Pellegrini
13  * @author Mathieu Faverge
14  * @author Pascal Henon
15  * @author Pierre Ramet
16  * @author Tony Delarue
17  * @date 2021-01-03
18  *
19  **/
20 #include "common.h"
21 #include "graph/graph.h"
22 #include "pastix/order.h"
23 #include "symbol/symbol.h"
24 
25 /**
26  *******************************************************************************
27  *
28  * @ingroup symbol_dev
29  *
30  * @brief Add a dof array to the symbol matrix if any.
31  *
32  * If the dof parameter is variadic, then a permuted version of the initial dof
33  * array is constructed to match the symbol matrix that is working the permuted
34  * matrix A.
35  *
36  *******************************************************************************
37  *
38  * @param[in] graph
39  * The original graph of the matrix
40  *
41  * @param[in] order
42  * The ordering structure describing the permutation of the unknowns in
43  * the compressed form.
44  *
45  * @param[inout] symbptr
46  * The symbol structure to which dof array must be added.
47  *
48  *******************************************************************************/
49 static inline void
50 symbol_init_adddofs( const pastix_graph_t *graph,
51  const pastix_order_t *order,
52  symbol_matrix_t *symbptr )
53 {
54  symbptr->dof = graph->dof;
55  symbptr->dofs = NULL;
56 
57  if ( symbptr->dof < 1 ) {
58  pastix_int_t symbbase = symbptr->baseval;
59  pastix_int_t ordebase = order->baseval;
60  pastix_int_t i, ip, n, d, *dofs;
61 
62  n = graph->gN;
63 
64  MALLOC_INTERN( symbptr->dofs, n+1, pastix_int_t );
65 
66  dofs = symbptr->dofs;
67  dofs[0] = symbbase;
68 
69  for(ip=0; ip<n; ip++, dofs++) {
70  i = order->peritab[ip] - ordebase;
71 
72  assert( i < n );
73  d = graph->dofs[i+1] - graph->dofs[i];
74 
75  dofs[1] = dofs[0] + d;
76  }
77  assert( (symbptr->dofs[n] - symbbase) == (graph->dofs[n] - graph->dofs[0]) );
78  }
79 
80  return;
81 }
82 
83 /**
84  * @addtogroup pastix_symbol
85  * @{
86  *
87  *******************************************************************************
88  *
89  * @brief Initialize the symbol structure.
90  *
91  * Initialized the permuted dof array if graph and order are provided. The
92  * symbol is considered as dof = 1 otherwise.
93  *
94  *******************************************************************************
95  *
96  * @param[in] graph
97  * The original graph of the matrix
98  *
99  * @param[in] order
100  * The ordering structure describing the permutation of the unknowns in
101  * the compressed form.
102  *
103  * @param[inout] symbptr
104  * The symbol structure to initialize.
105  *
106  *******************************************************************************/
107 void
108 pastixSymbolInit ( const pastix_graph_t *graph,
109  const pastix_order_t *order,
110  symbol_matrix_t *symbptr )
111 {
112  memset (symbptr, 0, sizeof (symbol_matrix_t));
113  symbptr->dof = 1;
114  symbptr->schurfcol = -1;
115 
116  if ( (graph != NULL) && (order != NULL) ) {
117  symbol_init_adddofs( graph, order, symbptr );
118  }
119 
120  return;
121 }
122 
123 /**
124  *******************************************************************************
125  *
126  * @brief Free the content of symbolic matrix.
127  *
128  * All the arrays from the structure are freed and the structure is memset to 0
129  * at exit, but the symbol itself is not freed. It will require a new call to
130  * pastixSymbolInit() if the memory space area needs to be reused for a new
131  * symbol matrix.
132  *
133  *******************************************************************************
134  *
135  * @param[inout] symbptr
136  * The pointer to the structure to free.
137  *
138  *******************************************************************************/
139 void
141 {
142  if (symbptr->dofs != NULL) {
143  memFree_null( symbptr->dofs );
144  }
145  if (symbptr->cblktab != NULL) {
146  memFree_null( symbptr->cblktab );
147  }
148  if (symbptr->bloktab != NULL) {
149  memFree_null( symbptr->bloktab );
150  }
151  if (symbptr->browtab != NULL) {
152  memFree_null( symbptr->browtab );
153  }
154  pastixSymbolInit( NULL, NULL, symbptr );
155 }
156 
157 /**
158  *******************************************************************************
159  *
160  * @brief Reallocate the data structure to optimize the memory alignment.
161  *
162  * This function is used when the symbol need to be shrinked to a smaller or
163  * larger set of blocks and column blocks. The original data is copied in the
164  * new arrays.
165  *
166  *******************************************************************************
167  *
168  * @param[inout] symbptr
169  * The pointer to the structure that needs to be reallocated.
170  *
171  *******************************************************************************/
172 void
174 {
175  symbol_cblk_t *cblktab = NULL;
176  symbol_blok_t *bloktab = NULL;
177 
178  /* Move column block array */
179  MALLOC_INTERN( cblktab, symbptr->cblknbr+1, symbol_cblk_t );
180  memcpy(cblktab, symbptr->cblktab, (symbptr->cblknbr + 1) * sizeof (symbol_cblk_t));
181  memFree(symbptr->cblktab);
182  symbptr->cblktab = cblktab;
183 
184  /* Move block array */
185  MALLOC_INTERN( bloktab, symbptr->bloknbr, symbol_blok_t );
186  memcpy(bloktab, symbptr->bloktab, (symbptr->bloknbr) * sizeof (symbol_blok_t));
187  memFree(symbptr->bloktab);
188  symbptr->bloktab = bloktab;
189 }
190 
191 /**
192  *******************************************************************************
193  *
194  * @brief Search the targeted block C for a couple of blocks A and B.
195  *
196  * When executing the simulation run to map the data on the cores, it requires
197  * to compute dependencies of each block. In that case for each couple of blocks
198  * A (defined by bloknum), and B (defined by bloksrc), we need to find the block
199  * that will receive the contribution for the matrix-matrix product. To speedup
200  * the search, the startsearch parameter can be given to specify that the index
201  * of the block searched is larger than this parameter. It returns the index of
202  * the C block when found, or -1 if no block is facing the update. This happens,
203  * only if ILU(k) factorization is applied and ricar is set to true. Otherwise,
204  * returning -1 is an error.
205  *
206  *******************************************************************************
207  *
208  * @param[in] symbptr
209  * The pointer to the structure that needs to be reallocated.
210  *
211  * @param[in] bloksrc
212  * Index of the block used as A in teh matrix-matrix product.
213  *
214  * @param[in] bloknum
215  * Index of the block used as B in teh matrix-matrix product.
216  *
217  * @param[in] startsearch
218  * Index of a block belonging to the facing cblk of B used as a
219  * starting point to find the index of C.
220  *
221  * @param[in] ricar
222  * Booleen to enable ILU(k) factorization or not
223  *
224  *******************************************************************************
225  *
226  * @retval -1 No block where found. This is an error if ricar is disabled,
227  * otherwise it means that the path is longer than the k parameter in
228  * the ILU(k) factorization.
229  *
230  * @retval i The index of the block C that will receive the contribution from
231  * A * B^t
232  *
233  *******************************************************************************/
234 pastix_int_t
236  pastix_int_t bloksrc,
237  pastix_int_t bloknum,
238  pastix_int_t startsearch,
239  int ricar )
240 {
241  symbol_blok_t *bsrc;
242  symbol_blok_t *bdst;
243  pastix_int_t i, fcblknum, fbloknum, lbloknum;
244 
245  fcblknum = symbptr->bloktab[bloksrc].fcblknm;
246  fbloknum = symbptr->cblktab[fcblknum].bloknum;
247  lbloknum = symbptr->cblktab[fcblknum+1].bloknum;
248 
249  if(startsearch < fbloknum )
250  startsearch = fbloknum;
251 
252  assert( startsearch < lbloknum );
253 
254  /* Block in original column block */
255  bsrc = (symbptr->bloktab) + bloknum;
256 
257  /* Search for the facing block in the facing column block */
258  bdst = (symbptr->bloktab) + startsearch;
259 
260  if(ricar == 0)
261  {
262  for(i=startsearch; i<lbloknum; i++, bdst++ )
263  if( bdst->lrownum >= bsrc->frownum)
264  break;
265 
266  /* We should always exit the loop in non ilu(k) mode */
267  assert( (bdst->frownum <= bsrc->frownum) &&
268  (bdst->lrownum >= bsrc->lrownum) );
269 
270  return i;
271  }
272  else
273  {
274  for(i=startsearch; i<lbloknum; i++, bdst++)
275  {
276  if( ((bsrc->frownum >= bdst->frownum) && (bsrc->frownum <= bdst->lrownum)) ||
277  ((bsrc->lrownum >= bdst->frownum) && (bsrc->lrownum <= bdst->lrownum)) ||
278  ((bsrc->frownum <= bdst->frownum) && (bsrc->lrownum >= bdst->lrownum)) )
279  return i; /** We found the first block that matches **/
280 
281  if(bsrc->lrownum < bdst->frownum)
282  {
283  return -1;
284  }
285  }
286  }
287  return -1;
288 }
289 
290 /**
291  *******************************************************************************
292  *
293  * @brief Construct the browtab array that stores the blocks in a CSR way.
294  *
295  * The browtab is an equivalent of the columns array in a CSR for the symbolic
296  * structure in terms of block indexes.
297  *
298  *******************************************************************************
299  *
300  * @param[inout] symbptr
301  * The pointer to the symbolic structure to update.
302  *
303  *******************************************************************************/
304 void
306 {
307  symbol_cblk_t *cblk;
308  symbol_blok_t *blok;
309  pastix_int_t *innbr, *intmp, *browtab;
310  pastix_int_t itercblk;
311  pastix_int_t cblknbr;
312  pastix_int_t edgenbr = symbptr->bloknbr - symbptr->cblknbr;
313 
314  symbptr->browmax = 0;
315  cblknbr = symbptr->cblknbr;
316 
317  MALLOC_INTERN(innbr, cblknbr, pastix_int_t );
318  memset( innbr, 0, cblknbr * sizeof(pastix_int_t) );
319 
320  /* Count the number of input edge per cblk */
321  cblk = symbptr->cblktab;
322  blok = symbptr->bloktab;
323  for(itercblk=0; itercblk<cblknbr; itercblk++, cblk++)
324  {
325  pastix_int_t iterblok = cblk[0].bloknum + 1;
326  pastix_int_t lbloknum = cblk[1].bloknum;
327 
328  /* Skip diagonal block */
329  blok++;
330 
331  /* Off-diagonal blocks */
332  for( ; iterblok < lbloknum; iterblok++, blok++)
333  {
334  innbr[ blok->fcblknm ]++;
335  }
336  }
337 
338  /* Initialize the brownum fields */
339  cblk = symbptr->cblktab;
340  cblk->brownum = 0;
341  for(itercblk=0; itercblk<cblknbr; itercblk++, cblk++)
342  {
343  symbptr->browmax = pastix_imax( symbptr->browmax, innbr[ itercblk ] );
344  cblk[1].brownum = cblk[0].brownum + innbr[ itercblk ];
345  innbr[itercblk] = cblk[0].brownum;
346  }
347  assert( cblk[0].brownum == edgenbr );
348 
349  /* Initialize the browtab */
350  MALLOC_INTERN( browtab, edgenbr, pastix_int_t );
351 
352  cblk = symbptr->cblktab;
353  blok = symbptr->bloktab;
354  for(itercblk=0; itercblk<cblknbr; itercblk++, cblk++)
355  {
356  pastix_int_t iterblok = cblk[0].bloknum + 1;
357  pastix_int_t lbloknum = cblk[1].bloknum;
358 
359  /* Skip diagonal block */
360  blok++;
361 
362  /* Off-diagonal blocks */
363  for( ; iterblok < lbloknum; iterblok++, blok++)
364  {
365  intmp = innbr + blok->fcblknm;
366  browtab[ *intmp ] = iterblok;
367  (*intmp)++;
368  }
369  }
370 
371  if (symbptr->browtab == NULL) {
372  memFree(symbptr->browtab);
373  }
374  symbptr->browtab = browtab;
375 
376  memFree( innbr );
377  return;
378 }
379 
380 /**
381  *******************************************************************************
382  *
383  * @brief Print statistical information about the symbolic matrix structure
384  *
385  *******************************************************************************
386  *
387  * @param[in] symbptr
388  * The pointer to the symbolic structure.
389  *
390  *******************************************************************************/
391 void
393 {
394  symbol_cblk_t *cblk;
395  symbol_blok_t *blok;
396  pastix_int_t itercblk, dof;
397  pastix_int_t cblknbr, bloknbr, cblksel;
398  pastix_int_t cblkmin, cblkmax;
399  pastix_int_t blokmin, blokmax;
400  double cblkavg1, blokavg1;
401  double cblkavg2, blokavg2;
402  size_t mem = 0;
403 
404  cblknbr = symbptr->cblknbr;
405  cblksel = 0;
406  bloknbr = symbptr->bloknbr - cblknbr;
407  cblkmin = PASTIX_INT_MAX;
408  cblkmax = 0;
409  cblkavg1 = 0;
410  cblkavg2 = 0;
411  blokmin = PASTIX_INT_MAX;
412  blokmax = 0;
413  blokavg1 = 0;
414  blokavg2 = 0;
415 
416  cblk = symbptr->cblktab;
417  blok = symbptr->bloktab;
418 
419  for(itercblk=0; itercblk<cblknbr; itercblk++, cblk++)
420  {
421  pastix_int_t iterblok = cblk[0].bloknum + 1;
422  pastix_int_t lbloknum = cblk[1].bloknum;
423 
424  pastix_int_t colnbr = cblk->lcolnum - cblk->fcolnum + 1;
425 
426  cblksel += cblk->selevtx;
427  cblkmin = pastix_imin( cblkmin, colnbr );
428  cblkmax = pastix_imax( cblkmax, colnbr );
429  cblkavg1 += colnbr;
430  cblkavg2 += colnbr * colnbr;
431  blok++;
432 
433  /* Only extra diagonal */
434  for( ; iterblok < lbloknum; iterblok++, blok++)
435  {
436  pastix_int_t rownbr = blok->lrownum - blok->frownum + 1;
437 
438  blokmin = pastix_imin( blokmin, rownbr );
439  blokmax = pastix_imax( blokmax, rownbr );
440  blokavg1 += rownbr;
441  blokavg2 += rownbr * rownbr;
442  }
443  }
444 
445  dof = symbptr->dof;
446  cblkmin *= dof;
447  cblkmax *= dof;
448  cblkavg1 = (cblkavg1 * (double)dof ) / (double)cblknbr;
449  cblkavg2 = sqrt( ((cblkavg2 * (double)dof * (double)dof) / (double)cblknbr) - cblkavg1 * cblkavg1 );
450 
451  /* Check if we have at least one off-diagonal block */
452  if ( bloknbr > 0 ) {
453  blokmin *= dof;
454  blokmax *= dof;
455  blokavg1 = (blokavg1 * (double)dof ) / (double)bloknbr;
456  blokavg2 = sqrt( ((blokavg2 * (double)dof * (double)dof) / (double)bloknbr) - blokavg1 * blokavg1 );
457  }
458  else {
459  blokmin = 0;
460  blokmax = 0;
461  blokavg1 = 0.;
462  blokavg2 = 0.;
463  }
464 
465  /* Compute symbol matrix space */
466  mem = sizeof( symbol_matrix_t );
467  mem += sizeof( symbol_cblk_t ) * (cblknbr + 1);
468  mem += sizeof( symbol_blok_t ) * symbptr->bloknbr;
469  mem += sizeof( pastix_int_t ) * bloknbr;
470 
471  fprintf(stdout,
472  " Symbol Matrix statistics:\n"
473  " Number of cblk %10ld\n"
474  " Number of blok %10ld\n"
475  " Cblk width min %10ld\n"
476  " Cblk width max %10ld\n"
477  " Cblk width avg %11.2lf\n"
478  " Cblk width stdev %11.2lf\n"
479  " Blok height min %10ld\n"
480  " Blok height max %10ld\n"
481  " Blok height avg %11.2lf\n"
482  " Blok height stdev %11.2lf\n"
483  " Memory space %11.2lf %co\n",
484  (long)cblknbr, (long)bloknbr,
485  (long)cblkmin, (long)cblkmax, cblkavg1, cblkavg2,
486  (long)blokmin, (long)blokmax, blokavg1, blokavg2,
487  pastix_print_value( mem ),
488  pastix_print_unit( mem ) );
489 
490  if ( cblksel > 0 ) {
491  fprintf( stdout,
492  " Number of selected cblk %10ld\n",
493  (long)cblksel );
494  }
495 }
496 
497 /**
498  * @}
499  */
symbol_matrix_s
Symbol matrix structure.
Definition: symbol.h:75
symbol_matrix_s::schurfcol
pastix_int_t schurfcol
Definition: symbol.h:80
symbol_cblk_t
struct symbol_cblk_s symbol_cblk_t
Symbol column block structure.
symbol_blok_t
struct symbol_blok_s symbol_blok_t
Symbol block structure.
symbol_matrix_s::bloktab
symbol_blok_t * bloktab
Definition: symbol.h:82
pastix_order_s::peritab
pastix_int_t * peritab
Definition: order.h:50
pastixSymbolExit
void pastixSymbolExit(symbol_matrix_t *symbptr)
Free the content of symbolic matrix.
Definition: symbol.c:140
symbol_matrix_s::dofs
pastix_int_t * dofs
Definition: symbol.h:87
symbol_cblk_s::fcolnum
pastix_int_t fcolnum
Definition: symbol.h:44
pastix_order_s
Order structure.
Definition: order.h:45
symbol_matrix_s::cblknbr
pastix_int_t cblknbr
Definition: symbol.h:77
symbol_matrix_s::baseval
pastix_int_t baseval
Definition: symbol.h:76
symbol_init_adddofs
static void symbol_init_adddofs(const pastix_graph_t *graph, const pastix_order_t *order, symbol_matrix_t *symbptr)
Add a dof array to the symbol matrix if any.
Definition: symbol.c:50
symbol_matrix_s::browtab
pastix_int_t * browtab
Definition: symbol.h:83
symbol_matrix_s::dof
pastix_int_t dof
Definition: symbol.h:85
symbol_matrix_s::bloknbr
pastix_int_t bloknbr
Definition: symbol.h:78
symbol_blok_s::fcblknm
pastix_int_t fcblknm
Definition: symbol.h:61
symbol_cblk_s::bloknum
pastix_int_t bloknum
Definition: symbol.h:46
pastixSymbolRealloc
void pastixSymbolRealloc(symbol_matrix_t *symbptr)
Reallocate the data structure to optimize the memory alignment.
Definition: symbol.c:173
symbol_matrix_t
struct symbol_matrix_s symbol_matrix_t
Symbol matrix structure.
symbol_cblk_s::brownum
pastix_int_t brownum
Definition: symbol.h:47
symbol.h
symbol_cblk_s::lcolnum
pastix_int_t lcolnum
Definition: symbol.h:45
pastix_order_s::baseval
pastix_int_t baseval
Definition: order.h:46
symbol_blok_s::frownum
pastix_int_t frownum
Definition: symbol.h:58
graph.h
symbol_matrix_s::cblktab
symbol_cblk_t * cblktab
Definition: symbol.h:81
symbol_blok_s
Symbol block structure.
Definition: symbol.h:57
pastixSymbolBuildRowtab
void pastixSymbolBuildRowtab(symbol_matrix_t *symbptr)
Construct the browtab array that stores the blocks in a CSR way.
Definition: symbol.c:305
order.h
symbol_cblk_s
Symbol column block structure.
Definition: symbol.h:43
symbol_matrix_s::browmax
pastix_int_t browmax
Definition: symbol.h:84
pastixSymbolPrintStats
void pastixSymbolPrintStats(const symbol_matrix_t *symbptr)
Print statistical information about the symbolic matrix structure.
Definition: symbol.c:392
symbol_blok_s::lrownum
pastix_int_t lrownum
Definition: symbol.h:59
pastixSymbolGetFacingBloknum
pastix_int_t pastixSymbolGetFacingBloknum(const symbol_matrix_t *symbptr, pastix_int_t bloksrc, pastix_int_t bloknum, pastix_int_t startsearch, int ricar)
Search the targeted block C for a couple of blocks A and B.
Definition: symbol.c:235