PaStiX Handbook  6.2.1
symbol_fax_iluk.c
Go to the documentation of this file.
1 /**
2  *
3  * @file symbol_fax_iluk.c
4  *
5  * PaStiX symbolic factorization routines for ILU(k)
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 Pascal Henon
12  * @author Mathieu Faverge
13  * @date 2021-01-03
14  *
15  **/
16 #include "common.h"
17 #include "graph/graph.h"
18 #include "pastix/order.h"
19 #include "symbol/symbol.h"
20 #include "fax_csr.h"
21 
22 /**
23  *******************************************************************************
24  *
25  * @ingroup symbol_dev_fax
26  *
27  * @brief Create the symbol matrix from the graph of the non zero pattern of the
28  * factorized matrix and the supernode partition.
29  *
30  *******************************************************************************
31  *
32  * @param[inout] P
33  * The non zero pattern of the factorized matrix. WARNING: on exit, the
34  * graph is destroyed.
35  *
36  * @param[in] cblknbr
37  * The number of supernode. Must be equal to P->n.
38  *
39  * @param[in] rangtab
40  * Integer array of size cblknbr+1.
41  * Contains the supernode partition of the graph.
42  *
43  * @param[out] symbmtx
44  * On entry, an initialized structure of symbol matrix (see pastixSymbolInit()).
45  * On exit, contains the symbol matrix associated to the graph P and
46  * the supernode partition given.
47  *
48  *******************************************************************************
49  *
50  * @retval 0 on success.
51  * @retval !0 on failure.
52  *
53  *******************************************************************************/
54 int
56  pastix_int_t levelk,
57  const pastix_graph_t *graphA,
58  const pastix_order_t *ordeptr )
59 {
60  pastix_int_t i, j, k, l;
61  pastix_int_t cblknum;
62  pastix_int_t ind;
63  pastix_int_t *tmp = NULL;
64  pastix_int_t *node2cblk = NULL;
65  pastix_int_t *ja = NULL;
66  pastix_int_t n, bloknbr, tmpsize;
67  pastix_int_t cblknbr;
68  pastix_int_t *rangtab;
69  fax_csr_t graphPA;
70  fax_csr_t graphL;
71 
72  assert( ordeptr->baseval == 0 );
73  cblknbr = ordeptr->cblknbr;
74  rangtab = ordeptr->rangtab;
75 
76  n = rangtab[cblknbr];
77  tmpsize = pastix_iceil( n, 2 );
78  MALLOC_INTERN( tmp, tmpsize, pastix_int_t );
79  MALLOC_INTERN( node2cblk, n, pastix_int_t );
80 
81  /*
82  * Let's regenerate the compresse graph of L:
83  * 1) Apply the permutation
84  * 2) Apply the symbolic ILUk
85  * 3) Compress the obtained L, with the partition previsouly computed at order stage
86  */
87  /* Create the graph of P A */
88  faxCSRGenPA( graphA, ordeptr->permtab, &graphPA );
89 
90  /* Create the graph of L */
91  memset( &graphL, 0, sizeof(fax_csr_t) );
92  faxCSRFactILUk( &graphPA, ordeptr, levelk, &graphL );
93  faxCSRClean( &graphPA );
94 
95  assert( cblknbr == graphL.n );
96 
97  /**** First we transform the P matrix to find the block ****/
98  for ( k = 0; k < cblknbr; k++ ) {
99  for ( i = rangtab[k]; i < rangtab[k + 1]; i++ ) {
100  node2cblk[i] = k;
101  }
102  }
103 
104  /* Let's update P to store all the couples (frownum,lrownum) in rows arrays */
105  bloknbr = 0;
106  for ( k = 0; k < cblknbr; k++ ) {
107  assert( graphL.nnz[k] >= ( rangtab[k + 1] - rangtab[k] ) );
108 
109 #if defined( PASTIX_DEBUG_SYMBOL )
110  for ( l = 0; l < rangtab[k + 1] - rangtab[k]; l++ ) {
111  assert( graphL.rows[k][l] == rangtab[k] + l );
112  assert( node2cblk[graphL.rows[k][l]] == k );
113  }
114 #endif
115 
116  ja = graphL.rows[k];
117  j = 0;
118  ind = 0;
119  while ( j < graphL.nnz[k] ) {
120  cblknum = node2cblk[ja[j]];
121  l = j + 1;
122  while ( ( l < graphL.nnz[k] ) && ( ja[l] == ja[l - 1] + 1 ) &&
123  ( node2cblk[ja[l]] == cblknum ) ) {
124  l++;
125  }
126 
127  assert( ind < tmpsize );
128  tmp[ind++] = ja[j];
129  tmp[ind++] = ja[l - 1];
130  assert( ( ja[l - 1] - ja[j] + 1 ) == ( l - j ) );
131 
132  j = l;
133  }
134 
135  graphL.nnz[k] = ind;
136  bloknbr += ind / 2;
137 
138  memFree( graphL.rows[k] );
139  MALLOC_INTERN( graphL.rows[k], ind, pastix_int_t );
140  memcpy( graphL.rows[k], tmp, ind * sizeof( pastix_int_t ) );
141  }
142 
143  memFree( tmp );
144 
145  assert( bloknbr == faxCSRGetNNZ( &graphL ) / 2 );
146 
147 #if defined( PASTIX_DEBUG_SYMBOL )
148  for ( k = 0; k < cblknbr; k++ ) {
149  assert( graphL.nnz[k] > 0 );
150  assert( graphL.rows[k][0] == rangtab[k] );
151  assert( ( graphL.rows[k][1] - graphL.rows[k][0] + 1 ) == ( rangtab[k + 1] - rangtab[k] ) );
152  }
153 #endif
154 
155  /**********************************/
156  /*** Compute the symbol matrix ****/
157  /**********************************/
158  symbptr->baseval = 0;
159  symbptr->cblknbr = cblknbr;
160  symbptr->bloknbr = bloknbr;
161  symbptr->nodenbr = n;
162  symbptr->browtab = NULL;
163 
164  MALLOC_INTERN( symbptr->cblktab, cblknbr + 1, symbol_cblk_t );
165  MALLOC_INTERN( symbptr->bloktab, symbptr->bloknbr, symbol_blok_t );
166 
167  ind = 0;
168  for ( k = 0; k < cblknbr; k++ ) {
169  symbptr->cblktab[k].fcolnum = rangtab[k];
170  symbptr->cblktab[k].lcolnum = rangtab[k + 1] - 1;
171  symbptr->cblktab[k].bloknum = ind;
172  symbptr->cblktab[k].brownum = -1;
173 
174  for ( i = 0; i < graphL.nnz[k]; i += 2 ) {
175  j = graphL.rows[k][i];
176  symbptr->bloktab[ind].frownum = j;
177  symbptr->bloktab[ind].lrownum = graphL.rows[k][i + 1];
178  symbptr->bloktab[ind].lcblknm = k;
179  symbptr->bloktab[ind].fcblknm = node2cblk[j];
180  ind++;
181 
182  assert( node2cblk[j] == node2cblk[graphL.rows[k][i + 1]] );
183  }
184 
185 #if defined( PASTIX_DEBUG_SYMBOL )
186  assert( symbptr->bloktab[symbptr->cblktab[k].bloknum].frownum ==
187  symbptr->cblktab[k].fcolnum );
188  assert( symbptr->bloktab[symbptr->cblktab[k].bloknum].lrownum ==
189  symbptr->cblktab[k].lcolnum );
190  assert( symbptr->bloktab[symbptr->cblktab[k].bloknum].fcblknm == k );
191 #endif
192  }
193 
194  /* virtual cblk to avoid side effect in the loops on cblk bloks */
195  symbptr->cblktab[cblknbr].fcolnum = symbptr->cblktab[cblknbr - 1].lcolnum + 1;
196  symbptr->cblktab[cblknbr].lcolnum = symbptr->cblktab[cblknbr - 1].lcolnum + 1;
197  symbptr->cblktab[cblknbr].bloknum = ind;
198  symbptr->cblktab[cblknbr].brownum = -1;
199 
200  assert( ind == symbptr->bloknbr );
201  memFree( node2cblk );
202 
203  return 0;
204 }
205 
206 /**
207  *******************************************************************************
208  *
209  * @ingroup symbol_dev_fax
210  *
211  * @brief Patch the symbol matrix to add blocks in order to get a
212  * real elimination tree.
213  *
214  * This function is called when ILU(k) factorization is
215  * performed and the faxCSRBuildSymbol() function might have returned a symbol
216  * matrix that doesn't provide a real elimination tree.
217  *
218  *******************************************************************************
219  *
220  * @param[inout] symbmtx
221  * On entry, a generated symbol matrix with faxCSRBuildSymbol() for example.
222  * On exit, the patched symbol matrix with extra blocks to have a real
223  * elimination tree.
224  *
225  *******************************************************************************/
226 void
228 {
229  pastix_int_t i, j, k;
230  pastix_int_t *father = NULL; /** For the cblk of the symbol matrix **/
231  symbol_blok_t *newbloktab = NULL;
232  symbol_cblk_t *cblktab = NULL;
233  symbol_blok_t *bloktab = NULL;
234  fax_csr_t Q;
235 
236  cblktab = symbmtx->cblktab;
237  bloktab = symbmtx->bloktab;
238 
239  MALLOC_INTERN( father, symbmtx->cblknbr, pastix_int_t );
240  MALLOC_INTERN( newbloktab, symbmtx->cblknbr + symbmtx->bloknbr, symbol_blok_t );
241 
242  faxCSRInit( symbmtx->cblknbr, &Q );
243 
244  /*
245  * Count how many extra-diagonal bloks are facing each diagonal blok
246  * Double-loop to avoid diagonal blocks.
247  */
248  for ( i = 0; i < symbmtx->cblknbr; i++ ) {
249  for ( j = cblktab[i].bloknum + 1; j < cblktab[i + 1].bloknum; j++ ) {
250  Q.nnz[bloktab[j].fcblknm]++;
251  }
252  }
253 
254  /* Allocate nFacingBlok integer for each diagonal blok */
255  for ( i = 0; i < symbmtx->cblknbr; i++ ) {
256  if ( Q.nnz[i] > 0 ) {
257  MALLOC_INTERN( Q.rows[i], Q.nnz[i], pastix_int_t );
258  }
259  else {
260  Q.rows[i] = NULL;
261  }
262  }
263 
264  memset( Q.nnz, 0, symbmtx->cblknbr * sizeof( pastix_int_t ) );
265 
266  /*
267  * Q.rows[k] will contain, for each extra-diagonal facing blok
268  * of the column blok k, its column blok.
269  */
270  for ( i = 0; i < symbmtx->cblknbr; i++ ) {
271  for ( j = cblktab[i].bloknum + 1; j < cblktab[i + 1].bloknum; j++ ) {
272  k = bloktab[j].fcblknm;
273  Q.rows[k][Q.nnz[k]++] = i;
274  }
275  }
276 
277  for ( i = 0; i < Q.n; i++ ) {
278  father[i] = -1;
279  }
280 
281  for ( i = 0; i < Q.n; i++ ) {
282  /*
283  * For each blok facing diagonal blok i, belonging to column blok k.
284  */
285  for ( j = 0; j < Q.nnz[i]; j++ ) {
286  k = Q.rows[i][j];
287  assert( k < i );
288 
289  while ( ( father[k] != -1 ) && ( father[k] != i ) ) {
290  k = father[k];
291  }
292  father[k] = i;
293  }
294  }
295 
296  for ( i = 0; i < Q.n; i++ ) {
297  if ( father[i] == -1 ) {
298  father[i] = i + 1;
299  }
300  }
301 
302  faxCSRClean( &Q );
303 
304  k = 0;
305  for ( i = 0; i < symbmtx->cblknbr - 1; i++ ) {
306  pastix_int_t odb, fbloknum;
307 
308  fbloknum = cblktab[i].bloknum;
309  memcpy( newbloktab + k, bloktab + fbloknum, sizeof( symbol_blok_t ) );
310  cblktab[i].bloknum = k;
311  k++;
312  odb = cblktab[i + 1].bloknum - fbloknum;
313  if ( odb <= 1 || bloktab[fbloknum + 1].fcblknm != father[i] ) {
314  /** Add a blok toward the father **/
315  newbloktab[k].frownum = cblktab[father[i]].fcolnum;
316  newbloktab[k].lrownum = cblktab[father[i]].fcolnum; /** OIMBE try lcolnum **/
317  newbloktab[k].lcblknm = i;
318  newbloktab[k].fcblknm = father[i];
319 #if defined( PASTIX_DEBUG_SYMBOL )
320  if ( father[i] != i ) {
321  assert( cblktab[father[i]].fcolnum > cblktab[i].lcolnum );
322  }
323 #endif
324  k++;
325  }
326 
327  if ( odb > 1 ) {
328  memcpy( newbloktab + k, bloktab + fbloknum + 1, sizeof( symbol_blok_t ) * ( odb - 1 ) );
329  k += odb - 1;
330  }
331  }
332 
333  /** Copy the last one **/
334  memcpy( newbloktab + k,
335  bloktab + symbmtx->cblktab[symbmtx->cblknbr - 1].bloknum,
336  sizeof( symbol_blok_t ) );
337  cblktab[symbmtx->cblknbr - 1].bloknum = k;
338  k++;
339  /** Virtual cblk **/
340  symbmtx->cblktab[symbmtx->cblknbr].bloknum = k;
341 
342 #if defined( PASTIX_DEBUG_SYMBOL )
343  assert( k >= symbmtx->bloknbr );
344  assert( k < symbmtx->cblknbr + symbmtx->bloknbr );
345 #endif
346  symbmtx->bloknbr = k;
347  memFree( symbmtx->bloktab );
348  MALLOC_INTERN( symbmtx->bloktab, k, symbol_blok_t );
349  memcpy( symbmtx->bloktab, newbloktab, sizeof( symbol_blok_t ) * symbmtx->bloknbr );
350  /* virtual cblk to avoid side effect in the loops on cblk bloks */
351  cblktab[symbmtx->cblknbr].bloknum = k;
352 
353  memFree( father );
354  memFree( newbloktab );
355 }
symbol_matrix_s
Symbol matrix structure.
Definition: symbol.h:75
pastix_order_s::cblknbr
pastix_int_t cblknbr
Definition: order.h:48
symbol_blok_s::lcblknm
pastix_int_t lcblknm
Definition: symbol.h:60
symbol_matrix_s::bloktab
symbol_blok_t * bloktab
Definition: symbol.h:82
pastix_order_s::permtab
pastix_int_t * permtab
Definition: order.h:49
faxCSRClean
void faxCSRClean(fax_csr_t *csr)
Free the data store in the structure.
Definition: fax_csr.c:65
symbol_cblk_s::fcolnum
pastix_int_t fcolnum
Definition: symbol.h:44
faxCSRInit
void faxCSRInit(pastix_int_t n, fax_csr_t *csr)
Initialize the data structure by doing the first allocations within the structure and initializing th...
Definition: fax_csr.c:40
pastix_order_s
Order structure.
Definition: order.h:45
symbol_matrix_s::nodenbr
pastix_int_t nodenbr
Definition: symbol.h:79
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_matrix_s::browtab
pastix_int_t * browtab
Definition: symbol.h:83
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
faxCSRGetNNZ
pastix_int_t faxCSRGetNNZ(const fax_csr_t *csr)
Computes the number of non zero entries in the graph.
Definition: fax_csr.c:99
symbol_cblk_s::brownum
pastix_int_t brownum
Definition: symbol.h:47
pastixSymbolFaxILUk
int pastixSymbolFaxILUk(symbol_matrix_t *symbptr, pastix_int_t levelk, const pastix_graph_t *graphA, const pastix_order_t *ordeptr)
Create the symbol matrix from the graph of the non zero pattern of the factorized matrix and the supe...
Definition: symbol_fax_iluk.c:55
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
faxCSRFactILUk
pastix_int_t faxCSRFactILUk(const fax_csr_t *graphA, const pastix_order_t *order, pastix_int_t level, fax_csr_t *graphL)
Compute the non zero pattern of the levelized incomplete factor for a sparse lower triangular matrix ...
Definition: fax_csr_iluk.c:57
graph.h
symbol_matrix_s::cblktab
symbol_cblk_t * cblktab
Definition: symbol.h:81
faxCSRPatchSymbol
void faxCSRPatchSymbol(symbol_matrix_t *symbmtx)
Patch the symbol matrix to add blocks in order to get a real elimination tree.
Definition: symbol_fax_iluk.c:227
faxCSRGenPA
int faxCSRGenPA(const pastix_graph_t *graphA, const pastix_int_t *perm, fax_csr_t *graphPA)
Generate the graph of P*A from the graph of A and the permutation vector.
Definition: fax_csr.c:183
fax_csr_s
Fax blocked csr structure.
Definition: fax_csr.h:25
fax_csr.h
symbol_blok_s
Symbol block structure.
Definition: symbol.h:57
pastix_order_s::rangtab
pastix_int_t * rangtab
Definition: order.h:51
order.h
symbol_cblk_s
Symbol column block structure.
Definition: symbol.h:43
symbol_blok_s::lrownum
pastix_int_t lrownum
Definition: symbol.h:59