PaStiX Handbook  6.4.0
fax_csr_direct.c
Go to the documentation of this file.
1 /**
2  *
3  * @file fax_csr_direct.c
4  *
5  * This file contains routines to create the graph of the direct
6  * factorization of a graph A.
7  *
8  * @copyright 2004-2024 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
9  * Univ. Bordeaux. All rights reserved.
10  *
11  * @version 6.4.0
12  * @author Pascal Henon
13  * @author Mathieu Faverge
14  * @date 2024-07-05
15  *
16  **/
17 #include "common.h"
18 #include "pastix/order.h"
19 #include "fax_csr.h"
20 
21 /**
22  *******************************************************************************
23  *
24  * @ingroup symbol_dev_csr
25  *
26  * @brief Compute the non zero pattern of the direct factorization of a matrix
27  * A, given the supernode partition associated.
28  *
29  *******************************************************************************
30  *
31  * @param[in] graphA
32  * The graph structure of the original matrix A to factorize.
33  *
34  * @param[in] order
35  * The order structure holding the number of cblk and the partition
36  *
37  * @param[out] graphL
38  * The graph the structure of the non zero pattern of the factorized
39  * matrix. On entry, a pointer to a graph structure. No need for
40  * initialization. On exit, the structure contains the computed graph.
41  *
42  *******************************************************************************
43  *
44  * @retval >=0, the number of non zero entries in the generated graph.
45  * @retval -i, if the i^th parameter is incorrect
46  *
47  *******************************************************************************/
49 faxCSRFactDirect( const fax_csr_t *graphA, const pastix_order_t *order, fax_csr_t *graphL )
50 {
51  pastix_int_t i, k, nnz;
52  pastix_int_t nnznbr, father;
53  pastix_int_t *tmp = NULL;
54  pastix_int_t *ja = NULL;
55  pastix_int_t cblknbr;
56  const pastix_int_t *rangtab;
57  const pastix_int_t *treetab;
58 
59  /* Check parameters */
60  if ( graphA == NULL ) {
61  return -1;
62  }
63  if ( order == NULL ) {
64  return -2;
65  }
66  if ( graphL == NULL ) {
67  return -3;
68  }
69 
70  cblknbr = order->cblknbr;
71  rangtab = order->rangtab;
72  treetab = order->treetab;
73 
74  if ( ( cblknbr < 0 ) || ( cblknbr > graphA->n ) ) {
75  return -4;
76  }
77 
78  /* Quick return */
79  if ( graphA->n == 0 ) {
80  return 0;
81  }
82 
83  MALLOC_INTERN( tmp, graphA->n, pastix_int_t );
84 
85  /* Compute the nnz structure of each supernode in A */
86  faxCSRCblkCompress( graphA, order, graphL, tmp );
87 
88  /* Compute the symbolic factorization */
89  for ( k = 0; k < cblknbr; k++, treetab++ ) {
90  father = *treetab;
91 
92  /* Merge son's nodes into father's list */
93  if ( ( father != k ) && ( father > 0 ) ) {
94  i = 0;
95  ja = graphL->rows[k];
96 
97  /* Take only the trows outside the cblk */
98  while ( ( i < graphL->nnz[k] ) && ( ja[i] < rangtab[k+1] ) ) {
99  i++;
100  }
101 
102  nnznbr = pastix_intset_union( graphL->nnz[k] - i,
103  graphL->rows[k] + i,
104  graphL->nnz[father],
105  graphL->rows[father],
106  tmp );
107 
108  memFree( graphL->rows[father] );
109  MALLOC_INTERN( graphL->rows[father], nnznbr, pastix_int_t );
110  memcpy( graphL->rows[father], tmp, sizeof( pastix_int_t ) * nnznbr );
111  graphL->nnz[father] = nnznbr;
112  }
113  }
114 
115 #if defined( PASTIX_DEBUG_SYMBOL )
116  /* Check that all terms of A are in the pattern */
117  {
118  pastix_int_t ind;
119  pastix_int_t j;
120  for ( k = 0; k < cblknbr; k++ ) {
121  /* Put the diagonal elements (A does not contains them) */
122  for ( i = rangtab[k]; i < rangtab[k + 1]; i++ ) {
123  j = 0;
124  while ( ( j < graphA->nnz[i] ) && ( graphA->rows[i][j] < i ) ) {
125  j++;
126  }
127 
128  for ( ind = j; ind < graphA->nnz[i]; ind++ ) {
129  assert( graphA->rows[i][ind] >= i );
130  }
131  for ( ind = j + 1; ind < graphA->nnz[i]; ind++ ) {
132  assert( graphA->rows[i][ind] > graphA->rows[i][ind - 1] );
133  }
134 
135  ind = pastix_intset_union(
136  graphL->nnz[k], graphL->rows[k], graphA->nnz[i] - j, graphA->rows[i] + j, tmp );
137 
138  assert( ind <= graphL->nnz[k] );
139  }
140  }
141  }
142 #endif
143 
144  memFree( tmp );
145 
146  /*
147  * Computes the nnz in graphL
148  */
149  nnz = 0;
150  for ( i = 0; i < cblknbr; i++ ) {
151  pastix_int_t ncol, nrow;
152  ncol = rangtab[i + 1] - rangtab[i];
153  nrow = graphL->nnz[i];
154 
155  assert( nrow >= ncol );
156  assert( nrow <= graphA->n );
157 
158  nnz += ( ncol * ( ncol + 1 ) ) / 2;
159  nnz += ( ncol * ( nrow - ncol ) );
160  }
161 
162  graphL->total_nnz = nnz;
163  return nnz;
164 }
BEGIN_C_DECLS typedef int pastix_int_t
Definition: datatypes.h:51
pastix_int_t * treetab
Definition: order.h:54
pastix_int_t cblknbr
Definition: order.h:50
pastix_int_t * rangtab
Definition: order.h:53
Order structure.
Definition: order.h:47
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
pastix_int_t faxCSRFactDirect(const fax_csr_t *graphA, const pastix_order_t *order, fax_csr_t *graphL)
Compute the non zero pattern of the direct factorization of a matrix A, given the supernode partition...
Fax blocked csr structure.
Definition: fax_csr.h:25