PaStiX Handbook  6.2.1
fax_csr.c
Go to the documentation of this file.
1 /**
2  *
3  * @file fax_csr.c
4  *
5  * PaStiX Fax csr routines to handle the graphs of PA and L
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  * @author Tony Delarue
14  * @date 2021-01-25
15  *
16  **/
17 #include "common.h"
18 #include "graph/graph.h"
19 #include "pastix/order.h"
20 #include "fax_csr.h"
21 
22 /**
23  *******************************************************************************
24  *
25  * @ingroup symbol_dev_csr
26  *
27  * @brief Initialize the data structure by doing the first allocations within
28  * the structure and initializing the fields.
29  *
30  *******************************************************************************
31  *
32  * @param[in] n
33  * The size of the graph that needs to be initialized.
34  *
35  * @param[out] csr
36  * The graph to initialize.
37  *
38  *******************************************************************************/
39 void
40 faxCSRInit( pastix_int_t n, fax_csr_t *csr )
41 {
42  csr->n = n;
43  csr->total_nnz = 0;
44  MALLOC_INTERN( csr->nnz, n, pastix_int_t );
45  MALLOC_INTERN( csr->rows, n, pastix_int_t * );
46 
47  memset( csr->nnz, 0, n * sizeof( pastix_int_t ) );
48  memset( csr->rows, 0, n * sizeof( pastix_int_t * ) );
49 }
50 
51 /**
52  *******************************************************************************
53  *
54  * @ingroup symbol_dev_csr
55  *
56  * @brief Free the data store in the structure.
57  *
58  *******************************************************************************
59  *
60  * @param[inout] csr
61  * The graph to clean.
62  *
63  *******************************************************************************/
64 void
66 {
67  pastix_int_t i;
68  for ( i = 0; i < csr->n; i++ ) {
69  if ( csr->nnz[i] != 0 ) {
70  memFree_null( csr->rows[i] );
71  }
72  }
73  memFree_null( csr->rows );
74  memFree_null( csr->nnz );
75 }
76 
77 /**
78  *******************************************************************************
79  *
80  * @ingroup symbol_dev_csr
81  *
82  * @brief Computes the number of non zero entries in the graph.
83  *
84  * It is using the following formula: nnz = sum( i=0..n, nnz[n] )
85  * The formula must be post computed to adapt to presence of diagonal elements
86  * or not, and to the symmetry of the graph.
87  *
88  *******************************************************************************
89  *
90  * @param[in] csr
91  * The graph on which the number of non zero entries is computed.
92  *
93  *******************************************************************************
94  *
95  * @retval The number of non zero entries.
96  *
97  *******************************************************************************/
98 pastix_int_t
99 faxCSRGetNNZ( const fax_csr_t *csr )
100 {
101  pastix_int_t i, nnz;
102  nnz = 0;
103  for ( i = 0; i < csr->n; i++ ) {
104  nnz += csr->nnz[i];
105  }
106  return nnz;
107 }
108 
109 /**
110  *******************************************************************************
111  *
112  * @ingroup symbol_dev_csr
113  *
114  * @brief Compact a compressed graph.
115  *
116  * All nodes with no non zero entries are removed from the graph, the allocated
117  * space is not adjusted.
118  *
119  *******************************************************************************
120  *
121  * @param[inout] csr
122  * The graph to compact.
123  * On entry, graph which might contain nodes with no non zero entries.
124  * On exit, all those nodes are suppressed and the compressed graph is
125  * returned.
126  *
127  *******************************************************************************/
128 void
130 {
131  pastix_int_t n = csr->n;
132  pastix_int_t i, j;
133 
134  /* Look for first deleted node */
135  for ( i = 0, j = 0; i < n; i++, j++ ) {
136  if ( csr->nnz[i] == 0 )
137  break;
138  }
139 
140  /* Compact the nodes */
141  for ( ; i < n; i++ ) {
142  if ( csr->nnz[i] > 0 ) {
143  assert( j < i );
144  csr->nnz[j] = csr->nnz[i];
145  csr->rows[j] = csr->rows[i];
146  csr->nnz[i] = 0;
147  csr->rows[i] = NULL;
148  j++;
149  }
150  }
151 
152  csr->n = j;
153 }
154 
155 /**
156  *******************************************************************************
157  *
158  * @ingroup symbol_dev_csr
159  *
160  * @brief Generate the graph of P*A from the graph of A and the
161  * permutation vector.
162  *
163  *******************************************************************************
164  *
165  * @param[in] graphA
166  * The original graph Aon which the permutation will be applied.
167  *
168  * @param[in] perm
169  * Integer array of size graphA->n. Contains the permutation to apply to A.
170  *
171  * @param[inout] graphPA
172  * On entry, the initialized graph with size graphA->n.
173  * On exit, contains the graph of P A.
174  *
175  *******************************************************************************
176  *
177  * @retval PASTIX_SUCCESS on success.
178  * @retval PASTIX_ERR_ALLOC if allocation went wrong.
179  * @retval PASTIX_ERR_BADPARAMETER if incorrect parameters are given.
180  *
181  *******************************************************************************/
182 int
183 faxCSRGenPA( const pastix_graph_t *graphA, const pastix_int_t *perm, fax_csr_t *graphPA )
184 {
185  pastix_int_t *rowsPA, *rowsA;
186  pastix_int_t i, j, ip;
187  pastix_int_t baseval;
188  pastix_int_t n = graphPA->n = graphA->n;
189 
190  baseval = graphA->colptr[0];
191 
192  faxCSRInit( graphA->n, graphPA );
193 
194  /* Compute the number of nnz per vertex */
195  for ( i = 0; i < n; i++ ) {
196  ip = perm[i];
197  /* Add diagonal (could be removed for direct) */
198  graphPA->nnz[ip] = graphA->colptr[i + 1] - graphA->colptr[i] + 1;
199  }
200 
201  /* Create the row vector */
202  for ( i = 0; i < n; i++ ) {
203  ip = perm[i] - baseval;
204 
205  MALLOC_INTERN( graphPA->rows[ip], graphPA->nnz[ip], pastix_int_t );
206 
207  rowsPA = graphPA->rows[ip];
208  rowsA = graphA->rowptr + graphA->colptr[i] - baseval;
209 
210  /* Add diagonal */
211  *rowsPA = ip;
212  rowsPA++;
213 
214  for ( j = 1; j < graphPA->nnz[ip]; j++, rowsPA++ ) {
215  *rowsPA = perm[*rowsA];
216  rowsA++;
217  }
218 
219  intSort1asc1( graphPA->rows[ip], graphPA->nnz[ip] );
220  }
221  return PASTIX_SUCCESS;
222 }
223 
224 /**
225  *******************************************************************************
226  *
227  * @ingroup symbol_dev_csr
228  *
229  * @brief Compact a element wise graph of a matrix A, according to the
230  * associated partition.
231  *
232  *******************************************************************************
233  *
234  * @param[in] graphA
235  * The original graph of A element wise to compress.
236  *
237  * @param[in] order
238  * The ordering structure that holds the partition associated to graphAo.
239  *
240  * @param[out] graphL
241  * On entry, the block wise initialized graph of A with size order->cblknbr.
242  * On exit, contains the compressed graph of A.
243  *
244  * @param[inout] work
245  * Workspace of size >= max( degree(L_i) ), so >= grapA->n.
246  *
247  *******************************************************************************/
248 void
250  const pastix_order_t *order,
251  fax_csr_t *graphL,
252  pastix_int_t *work )
253 {
254  pastix_int_t i, j, k, nnznbr;
255  const pastix_int_t cblknbr = order->cblknbr;
256  const pastix_int_t *rangtab = order->rangtab;
257  pastix_int_t *work2;
258  pastix_int_t *tmp, *tmp1, *tmp2;
259 
260  MALLOC_INTERN( work2, graphA->n, pastix_int_t );
261  tmp1 = work;
262  tmp2 = work2;
263 
264  assert( order->baseval == 0 );
265  faxCSRInit( cblknbr, graphL );
266 
267  /*
268  * Let's accumulate the row presents in the column blok k in tmp1
269  * Then use tmp2 to merge the elements of the next column, and tmp to switch
270  * pointers.
271  */
272  for ( k = 0; k < cblknbr; k++ ) {
273  /* Put the diagonal elements (In case A does not contains them) */
274  j = 0;
275  for ( i = rangtab[k]; i < rangtab[k + 1]; i++ ) {
276  tmp1[j++] = i;
277  }
278  nnznbr = j;
279 
280  for ( i = rangtab[k]; i < rangtab[k + 1]; i++ ) {
281  j = 0;
282 
283  /* We take only the elements greater than i */
284  while ( ( j < graphA->nnz[i] ) && ( graphA->rows[i][j] <= i ) ) {
285  j++;
286  }
287 
288  /* Merge the actual list with the edges of the ith vertex */
289  nnznbr =
290  pastix_intset_union( nnznbr, tmp1, graphA->nnz[i] - j, graphA->rows[i] + j, tmp2 );
291 
292  /* Swap tmp1 and the merged set tmp */
293  tmp = tmp1;
294  tmp1 = tmp2;
295  tmp2 = tmp;
296  }
297 
298 #if !defined( NDEBUG ) && defined( PASTIX_DEBUG_SYMBOL )
299  /* Check that the first elements are the diagonal ones */
300  {
301  pastix_int_t ind;
302  ind = 0;
303  assert( nnznbr >= ( rangtab[k + 1] - rangtab[k] ) );
304  for ( j = rangtab[k]; j < rangtab[k + 1]; j++ ) {
305  assert( tmp1[ind] == j );
306  ind++;
307  }
308  assert( nnznbr > 0 );
309  }
310 #endif
311 
312  /* Update graphL */
313  graphL->nnz[k] = nnznbr;
314  MALLOC_INTERN( graphL->rows[k], nnznbr, pastix_int_t );
315  memcpy( graphL->rows[k], tmp1, sizeof( pastix_int_t ) * nnznbr );
316  }
317 
318  memFree( work2 );
319 }
pastix_order_s::cblknbr
pastix_int_t cblknbr
Definition: order.h:48
faxCSRClean
void faxCSRClean(fax_csr_t *csr)
Free the data store in the structure.
Definition: fax_csr.c:65
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
faxCSRCompact
void faxCSRCompact(fax_csr_t *csr)
Compact a compressed graph.
Definition: fax_csr.c:129
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
PASTIX_SUCCESS
@ PASTIX_SUCCESS
Definition: api.h:346
pastix_order_s::baseval
pastix_int_t baseval
Definition: order.h:46
graph.h
faxCSRCblkCompress
void faxCSRCblkCompress(const fax_csr_t *graphA, const pastix_order_t *order, fax_csr_t *graphL, pastix_int_t *work)
Compact a element wise graph of a matrix A, according to the associated partition.
Definition: fax_csr.c:249
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
pastix_order_s::rangtab
pastix_int_t * rangtab
Definition: order.h:51
order.h