PaStiX Handbook  6.2.1
pastix_task_refine.c
Go to the documentation of this file.
1 /**
2  *
3  * @file pastix_task_refine.c
4  *
5  * PaStiX refinement functions implementations.
6  *
7  * @copyright 2015-2021 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
8  * Univ. Bordeaux. All rights reserved.
9  *
10  * @version 6.2.0
11  * @author Mathieu Faverge
12  * @author Pierre Ramet
13  * @author Xavier Lacoste
14  * @author Gregoire Pichon
15  * @author Tony Delarue
16  * @date 2021-03-30
17  *
18  **/
19 #include "common.h"
20 #include "bcsc/bcsc.h"
25 #include "pastix/order.h"
26 
27 /**
28  *******************************************************************************
29  *
30  * @ingroup pastix_dev_refine
31  *
32  * @brief Select the refinement function to call depending on the matrix type
33  * and the precision
34  *
35  *******************************************************************************/
36 static pastix_int_t (*sopalinRefine[4][4])(pastix_data_t *pastix_data, void *x, void *b) =
37 {
38  // PastixRefineGMRES
39  {
44  },
45  // PastixRefineCG
46  {
47  s_grad_smp,
48  d_grad_smp,
49  c_grad_smp,
51  },
52  // PastixRefineSR
53  {
58  },
59  // PastixRefineBiCGSTAB
60  {
65  }
66 };
67 
68 /**
69  *******************************************************************************
70  *
71  * @ingroup pastix_refine
72  *
73  * @brief Perform the iterative refinement without apply the permutations.
74  *
75  * This routine is affected by the following parameters:
76  * IPARM_REFINEMENT, DPARM_EPSILON_REFINEMENT
77  *
78  *******************************************************************************
79  *
80  * @param[in] pastix_data
81  * The PaStiX data structure that describes the solver instance.
82  *
83  * @param[in] n
84  * The size of system to solve, and the number of rows of both
85  * matrices b and x.
86  *
87  * @param[in] nrhs
88  * The number of right hand side members, and the number of columns of
89  * b and x.
90  *
91  * @param[inout] b
92  * The right hand side matrix of size ldb-by-nrhs.
93  * B is noted as inout, as permutation might be performed on the
94  * matrix. On exit, the matrix is restored as it was on entry.
95  *
96  * @param[in] ldb
97  * The leading dimension of the matrix b. ldb >= n.
98  *
99  * @param[inout] x
100  * The matrix x of size ldx-by-nrhs.
101  * On entry, the initial guess x0 for the refinement step, that may be
102  * the solution returned by the solve step or any other initial guess.
103  * On exit, contains the final solution after the iterative refinement.
104  *
105  * @param[in] ldx
106  * The leading dimension of the matrix x. ldx >= n.
107  *
108  *******************************************************************************
109  *
110  * @retval PASTIX_SUCCESS on successful exit,
111  * @retval PASTIX_ERR_BADPARAMETER if one parameter is incorrect,
112  *
113  *******************************************************************************/
114 int
115 pastix_subtask_refine( pastix_data_t *pastix_data,
116  pastix_int_t n, pastix_int_t nrhs,
117  const void *b, pastix_int_t ldb,
118  void *x, pastix_int_t ldx )
119 {
120  pastix_int_t *iparm = pastix_data->iparm;
121  pastix_bcsc_t *bcsc = pastix_data->bcsc;
122  double timer;
123 
124  if (nrhs > 1)
125  {
126  errorPrintW("Refinement works only with 1 rhs, We will iterate on each RHS one by one\n");
127  }
128 
129  if ( (pastix_data->schur_n > 0) && (iparm[IPARM_SCHUR_SOLV_MODE] != PastixSolvModeLocal))
130  {
131  fprintf(stderr, "Refinement is not available with Schur complement when non local solve is required\n");
133  }
134 
135  /* Prepare the refinement threshold, if not set by the user */
136  if ( pastix_data->dparm[DPARM_EPSILON_REFINEMENT] < 0. ) {
137  if ( (bcsc->flttype == PastixFloat) ||
138  (bcsc->flttype == PastixComplex32) ) {
139  pastix_data->dparm[DPARM_EPSILON_REFINEMENT] = 1e-6;
140  }
141  else {
142  pastix_data->dparm[DPARM_EPSILON_REFINEMENT] = 1e-12;
143  }
144  }
145 
146  clockSyncStart( timer, pastix_data->inter_node_comm );
147  {
148  pastix_int_t (*refinefct)(pastix_data_t *, void *, void *) = sopalinRefine[iparm[IPARM_REFINEMENT]][pastix_data->bcsc->flttype -2];
149  char *xptr = (char *)x;
150  char *bptr = (char *)b;
151  size_t shiftx, shiftb;
152  int i;
153 
154  shiftx = ldx * pastix_size_of( pastix_data->bcsc->flttype );
155  shiftb = ldb * pastix_size_of( pastix_data->bcsc->flttype );
156 
157  for(i=0; i<nrhs; i++, xptr += shiftx, bptr += shiftb ) {
158  pastix_int_t it;
159  it = refinefct( pastix_data, xptr, bptr );
160  pastix_data->iparm[IPARM_NBITER] = pastix_imax( it, pastix_data->iparm[IPARM_NBITER] );
161  }
162  }
163  clockSyncStop( timer, pastix_data->inter_node_comm );
164 
165  pastix_data->dparm[DPARM_REFINE_TIME] = clockVal(timer);
166  if ( iparm[IPARM_VERBOSE] > PastixVerboseNot ) {
167  pastix_print( pastix_data->inter_node_procnum,
168  0, OUT_TIME_REFINE,
169  pastix_data->dparm[DPARM_REFINE_TIME] );
170  }
171 
172  (void)n;
173  return PASTIX_SUCCESS;
174 }
175 
176 /**
177  *******************************************************************************
178  *
179  * @ingroup pastix_users
180  *
181  * @brief Perform iterative refinement.
182  *
183  * This routine performs the permutation of x, and b before and after the
184  * iterative refinement solution. To prevent extra permuation to happen, see
185  * pastix_subtask_refine().
186  * This routine is affected by the following parameters:
187  * IPARM_REFINEMENT, DPARM_EPSILON_REFINEMENT
188  *
189  *******************************************************************************
190  *
191  * @param[in] pastix_data
192  * The PaStiX data structure that describes the solver instance.
193  *
194  * @param[in] n
195  * The size of system to solve, and the number of rows of both
196  * matrices b and x.
197  *
198  * @param[in] nrhs
199  * The number of right hand side members, and the number of columns of
200  * b and x.
201  *
202  * @param[inout] b
203  * The right hand side matrix of size ldb-by-nrhs.
204  * B is noted as inout, as permutation might be performed on the
205  * matrix. On exit, the matrix is restored as it was on entry.
206  *
207  * @param[in] ldb
208  * The leading dimension of the matrix b. ldb >= n.
209  *
210  * @param[inout] x
211  * The matrix x of size ldx-by-nrhs.
212  * On entry, the initial guess x0 for the refinement step, that may be
213  * the solution returned by the solve step or any other initial guess.
214  * On exit, contains the final solution after the iterative refinement.
215  *
216  * @param[in] ldx
217  * The leading dimension of the matrix x. ldx >= n.
218  *
219  *******************************************************************************
220  *
221  * @retval PASTIX_SUCCESS on successful exit,
222  * @retval PASTIX_ERR_BADPARAMETER if one parameter is incorrect,
223  *
224  *******************************************************************************/
225 int
226 pastix_task_refine( pastix_data_t *pastix_data,
227  pastix_int_t n, pastix_int_t nrhs,
228  void *b, pastix_int_t ldb,
229  void *x, pastix_int_t ldx )
230 {
231  pastix_int_t *iparm = pastix_data->iparm;
232  pastix_bcsc_t *bcsc = pastix_data->bcsc;
233  int rc;
234  void *bglob, *btmp = NULL;
235  void *xglob, *xtmp = NULL;
236 
237  if ( (pastix_data->schur_n > 0) && (iparm[IPARM_SCHUR_SOLV_MODE] != PastixSolvModeLocal))
238  {
239  fprintf(stderr, "Refinement is not available with Schur complement when non local solve is required\n");
241  }
242 
243  /* Prepare the refinement threshold, if not set by the user */
244  if ( pastix_data->dparm[DPARM_EPSILON_REFINEMENT] < 0. ) {
245  if ( (bcsc->flttype == PastixFloat) ||
246  (bcsc->flttype == PastixComplex32) ) {
247  pastix_data->dparm[DPARM_EPSILON_REFINEMENT] = 1e-6;
248  }
249  else {
250  pastix_data->dparm[DPARM_EPSILON_REFINEMENT] = 1e-12;
251  }
252  }
253 
254  /* The spm is distributed, we have to gather it for the moment */
255  if ( pastix_data->csc->loc2glob != NULL ) {
256  if( iparm[IPARM_VERBOSE] > PastixVerboseNo ) {
257  pastix_print( pastix_data->procnum, 0, "pastix_task_refine: the RHS has to be centralized for the moment\n" );
258  }
259  btmp = b;
260  xtmp = x;
261 
262  spmGatherRHS( nrhs, pastix_data->csc, b, ldb, &bglob,-1 );
263  spmGatherRHS( nrhs, pastix_data->csc, x, ldx, &xglob,-1 );
264 
265  b = bglob;ldb = pastix_data->csc->gNexp;
266  x = xglob;ldx = pastix_data->csc->gNexp;
267  n = pastix_data->csc->gNexp;
268  }
269 
270  /* Compute P * b */
271  rc = pastix_subtask_applyorder( pastix_data, bcsc->flttype,
272  PastixDirForward, bcsc->gN, nrhs, b, ldb );
273  if( rc != PASTIX_SUCCESS ) {
274  return rc;
275  }
276 
277  /* Compute P * x */
278  rc = pastix_subtask_applyorder( pastix_data, bcsc->flttype,
279  PastixDirForward, bcsc->gN, nrhs, x, ldx );
280  if( rc != PASTIX_SUCCESS ) {
281  return rc;
282  }
283 
284  /* Performe the iterative refinement */
285  rc = pastix_subtask_refine( pastix_data, n, nrhs, b, ldb, x, ldx );
286  if( rc != PASTIX_SUCCESS ) {
287  return rc;
288  }
289 
290  /* Compute P * b */
291  rc = pastix_subtask_applyorder( pastix_data, bcsc->flttype,
292  PastixDirBackward, bcsc->gN, nrhs, b, ldb );
293  if( rc != PASTIX_SUCCESS ) {
294  return rc;
295  }
296 
297  /* Compute P * x */
298  rc = pastix_subtask_applyorder( pastix_data, bcsc->flttype,
299  PastixDirBackward, bcsc->gN, nrhs, x, ldx );
300  if( rc != PASTIX_SUCCESS ) {
301  return rc;
302  }
303 
304  if ( btmp != NULL) {
305  pastix_int_t ldbglob = ldb;
306  ldb = pastix_data->csc->nexp;
307  b = btmp;
308 
309  spmExtractLocalRHS( nrhs, pastix_data->csc, bglob, ldbglob, b, ldb );
310  memFree_null( bglob );
311  }
312  if ( xtmp != NULL ) {
313  pastix_int_t ldxglob = ldx;
314  ldx = pastix_data->csc->nexp;
315  x = xtmp;
316 
317  spmExtractLocalRHS( nrhs, pastix_data->csc, xglob, ldxglob, x, ldx );
318  memFree_null( xglob );
319  }
320 
321  (void)n;
322  return PASTIX_SUCCESS;
323 }
z_gmres_smp
pastix_int_t z_gmres_smp(pastix_data_t *pastix_data, void *x, void *b)
Definition: z_refine_gmres.c:48
s_refine_functions.h
IPARM_NBITER
@ IPARM_NBITER
Definition: api.h:111
pastix_subtask_applyorder
int pastix_subtask_applyorder(pastix_data_t *pastix_data, pastix_coeftype_t flttype, pastix_dir_t dir, pastix_int_t m, pastix_int_t n, void *b, pastix_int_t ldb)
Apply a permutation on the right-and-side vector before the solve step.
Definition: pastix_task_solve.c:122
IPARM_SCHUR_SOLV_MODE
@ IPARM_SCHUR_SOLV_MODE
Definition: api.h:106
pastix_bcsc_s
Internal column block distributed CSC matrix.
Definition: bcsc.h:37
z_pivot_smp
pastix_int_t z_pivot_smp(pastix_data_t *pastix_data, void *x, void *b)
Definition: z_refine_pivot.c:48
DPARM_EPSILON_REFINEMENT
@ DPARM_EPSILON_REFINEMENT
Definition: api.h:154
c_refine_functions.h
DPARM_REFINE_TIME
@ DPARM_REFINE_TIME
Definition: api.h:171
PastixDirForward
@ PastixDirForward
Definition: api.h:489
z_bicgstab_smp
pastix_int_t z_bicgstab_smp(pastix_data_t *pastix_data, void *x, void *b)
Definition: z_refine_bicgstab.c:48
d_gmres_smp
pastix_int_t d_gmres_smp(pastix_data_t *pastix_data, void *x, void *b)
Definition: d_refine_gmres.c:48
sopalinRefine
static pastix_int_t(* sopalinRefine[4][4])(pastix_data_t *pastix_data, void *x, void *b)
Select the refinement function to call depending on the matrix type and the precision.
Definition: pastix_task_refine.c:36
pastix_bcsc_s::flttype
int flttype
Definition: bcsc.h:41
s_gmres_smp
pastix_int_t s_gmres_smp(pastix_data_t *pastix_data, void *x, void *b)
Definition: s_refine_gmres.c:48
bcsc.h
d_pivot_smp
pastix_int_t d_pivot_smp(pastix_data_t *pastix_data, void *x, void *b)
Definition: d_refine_pivot.c:48
c_bicgstab_smp
pastix_int_t c_bicgstab_smp(pastix_data_t *pastix_data, void *x, void *b)
Definition: c_refine_bicgstab.c:48
c_pivot_smp
pastix_int_t c_pivot_smp(pastix_data_t *pastix_data, void *x, void *b)
Definition: c_refine_pivot.c:48
PastixVerboseNot
@ PastixVerboseNot
Definition: api.h:207
s_grad_smp
pastix_int_t s_grad_smp(pastix_data_t *pastix_data, void *x, void *b)
Definition: s_refine_grad.c:48
PASTIX_SUCCESS
@ PASTIX_SUCCESS
Definition: api.h:344
PastixVerboseNo
@ PastixVerboseNo
Definition: api.h:208
PastixDirBackward
@ PastixDirBackward
Definition: api.h:490
c_gmres_smp
pastix_int_t c_gmres_smp(pastix_data_t *pastix_data, void *x, void *b)
Definition: c_refine_gmres.c:48
d_bicgstab_smp
pastix_int_t d_bicgstab_smp(pastix_data_t *pastix_data, void *x, void *b)
Definition: d_refine_bicgstab.c:48
d_grad_smp
pastix_int_t d_grad_smp(pastix_data_t *pastix_data, void *x, void *b)
Definition: d_refine_grad.c:48
z_refine_functions.h
IPARM_VERBOSE
@ IPARM_VERBOSE
Definition: api.h:36
order.h
s_pivot_smp
pastix_int_t s_pivot_smp(pastix_data_t *pastix_data, void *x, void *b)
Definition: s_refine_pivot.c:48
pastix_task_refine
int pastix_task_refine(pastix_data_t *pastix_data, pastix_int_t n, pastix_int_t nrhs, void *b, pastix_int_t ldb, void *x, pastix_int_t ldx)
Perform iterative refinement.
Definition: pastix_task_refine.c:226
z_grad_smp
pastix_int_t z_grad_smp(pastix_data_t *pastix_data, void *x, void *b)
Definition: z_refine_grad.c:48
pastix_bcsc_s::gN
int gN
Definition: bcsc.h:38
PASTIX_ERR_BADPARAMETER
@ PASTIX_ERR_BADPARAMETER
Definition: api.h:351
d_refine_functions.h
IPARM_REFINEMENT
@ IPARM_REFINEMENT
Definition: api.h:110
s_bicgstab_smp
pastix_int_t s_bicgstab_smp(pastix_data_t *pastix_data, void *x, void *b)
Definition: s_refine_bicgstab.c:48
c_grad_smp
pastix_int_t c_grad_smp(pastix_data_t *pastix_data, void *x, void *b)
Definition: c_refine_grad.c:48
pastix_subtask_refine
int pastix_subtask_refine(pastix_data_t *pastix_data, pastix_int_t n, pastix_int_t nrhs, const void *b, pastix_int_t ldb, void *x, pastix_int_t ldx)
Perform the iterative refinement without apply the permutations.
Definition: pastix_task_refine.c:115