PaStiX Handbook  6.2.1
fax_csr_iluk.c
Go to the documentation of this file.
1 /**
2  *
3  * @file fax_csr_iluk.c
4  *
5  * This file contains croutines to create the graph of the incomplete
6  * factorization of a graph A.
7  *
8  * @copyright 2004-2021 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
9  * Univ. Bordeaux. All rights reserved.
10  *
11  * @version 6.0.3
12  * @author Pascal Henon
13  * @author Mathieu Faverge
14  * @date 2019-11-12
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 levelized
27  * incomplete factor for a sparse lower triangular matrix in CSC format. This
28  * pattern is exact iff the matrix has a SYMMETRIC non zero structure.
29  *
30  * @remark This algorithm has been implemented according to the paper of David
31  * Hysom and Alex Pothen : Level-based Incomplete LU factorization: Graph
32  * Model and Algorithm
33  *
34  *******************************************************************************
35  *
36  * @param[in] graphA
37  * The graph structure of the original matrix A to factorize.
38  *
39  * @param[in] order
40  * The order structure holding the number of cblk and the partition
41  *
42  * @param[in] level
43  * It is the level desired for the ilu(k) factorization. (level >= 0)
44  *
45  * @param[inout] graphL
46  * The graph the structure of the non zero pattern of the factorized matrix.
47  * On entry, a pointer to a graph structure. No need for initialization.
48  * On exit, the structure contains the computed graph.
49  *
50  *******************************************************************************
51  *
52  * @retval >=0, the number of non zero entries in the generated graph.
53  * @retval -i, if the i^th parameter is incorrect
54  *
55  *******************************************************************************/
56 pastix_int_t
57 faxCSRFactILUk( const fax_csr_t *graphA,
58  const pastix_order_t *order,
59  pastix_int_t level,
60  fax_csr_t *graphL )
61 {
62  pastix_int_t *visited = NULL;
63  pastix_int_t *length = NULL;
64  pastix_int_t *stack = NULL;
65  pastix_int_t *adj = NULL;
66  pastix_int_t *ja = NULL;
67  pastix_int_t used;
68  pastix_int_t h, i, j, k, t;
69  long nnz;
70  fax_csr_t tmpgrah;
71 
72  /* Check parameters */
73  if ( graphA == NULL ) {
74  return -1;
75  }
76  if ( level < 0 ) {
77  return -2;
78  }
79  if ( graphL == NULL ) {
80  return -3;
81  }
82 
83  /* Quick return */
84  if ( graphA->n == 0 ) {
85  return 0;
86  }
87 
88  /** Allocated the working array **/
89  MALLOC_INTERN( visited, graphA->n, pastix_int_t );
90  MALLOC_INTERN( length, graphA->n, pastix_int_t );
91  MALLOC_INTERN( stack, graphA->n, pastix_int_t );
92  MALLOC_INTERN( ja, graphA->n, pastix_int_t );
93  nnz = 0;
94 
95  /** Initialized visited ***/
96  for ( j = 0; j < graphA->n; j++ ) {
97  visited[j] = -1;
98  length[j] = 0;
99  }
100 
101  /** Apply GS_Urow for each row **/
102  faxCSRInit( graphA->n, graphL );
103  for ( i = 0; i < graphA->n; i++ ) {
104  /** Reset the stack number of elements **/
105  stack[0] = i;
106  used = 1;
107 
108  length[i] = 0;
109  visited[i] = i;
110 
111  ja[0] = i; /** Put the diagonal term **/
112  k = 1;
113 
114  /** BFS phase **/
115  while ( used > 0 ) {
116  used--;
117  h = stack[used];
118  adj = graphA->rows[h];
119  for ( j = 0; j < graphA->nnz[h]; j++ ) {
120  t = adj[j];
121  if ( visited[t] != i ) {
122  visited[t] = i;
123  if ( ( t < i ) && ( length[h] < level ) ) {
124  stack[used] = t;
125  used++;
126  length[t] = length[h] + 1;
127  }
128  if ( t > i ) {
129  ja[k++] = t;
130  }
131  }
132  }
133  }
134 
135  assert( k > 0 );
136 
137  graphL->nnz[i] = k;
138  MALLOC_INTERN( graphL->rows[i], k, pastix_int_t );
139  memcpy( graphL->rows[i], ja, k * sizeof( pastix_int_t ) );
140 
141  intSort1asc1( graphL->rows[i], graphL->nnz[i] );
142 
143  nnz += k;
144  }
145 
146  /* Compress with the elementary partition */
147  tmpgrah = *graphL;
148  faxCSRCblkCompress( &tmpgrah, order, graphL, visited );
149  faxCSRClean( &tmpgrah );
150 
151  memFree_null( ja );
152  memFree_null( visited );
153  memFree_null( length );
154  memFree_null( stack );
155 
156  graphL->total_nnz = nnz;
157  return nnz;
158 }
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
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
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
fax_csr_s
Fax blocked csr structure.
Definition: fax_csr.h:25
fax_csr.h
order.h