PaStiX Handbook  6.2.1
pastix.c
Go to the documentation of this file.
1 /**
2  *
3  * @file pastix.c
4  *
5  * PaStiX main interface for compatibility with former releases
6  *
7  * @copyright 2011-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 Matias Hastaran
16  * @date 2021-01-25
17  *
18  **/
19 #include "common.h"
20 #include <spm.h>
21 #include "pastix/order.h"
22 
23 /**
24  *******************************************************************************
25  *
26  * @ingroup pastix_api
27  * @brief Main function for compatibility with former releases
28  *
29  *******************************************************************************
30  *
31  * @param[inout] pastix_data_ptr
32  * The pastix data structure of the solver to store the state of the
33  * solver at every call.
34  *
35  * @param[in] pastix_comm
36  * The MPI communicator to use for the distributed solver.
37  *
38  * @param[in] n
39  * The size of the sparse matrix system to solve.
40  *
41  * @param[inout] colptr
42  * The pointer to the column index of the sparse matrix in the CSC
43  * format.
44  * On exit, the base value of the array might have changed, and/or the
45  * pointer might have been freed if IPARM_FREE_CSCUSER is set, and the
46  * factorization step is performed.
47  *
48  * @param[inout] row
49  * The pointer to the row array of the sparse matrix in the CSC
50  * format.
51  * On exit, the base value of the array might have changed, and/or the
52  * pointer might have been freed if IPARM_FREE_CSCUSER is set, and the
53  * factorization step is performed.
54  *
55  * @param[inout] avals
56  * The pointer to the values array of the sparse matrix in the CSC
57  * format.
58  * On exit, the pointer might have been freed if IPARM_FREE_CSCUSER is
59  * set, and the factorization step is performed.
60  *
61  * @param[inout] perm
62  * The pointer to the permutation array.
63  * On entry: the pointer might be allocated to store the generated
64  * permutation on exit, or to provide the user permutation.
65  * On exit, the permutation used by the solver is returned if perm is
66  * not NULL.
67  *
68  * @param[inout] invp
69  * The pointer to the inverse permutation array.
70  * On entry: the pointer might be allocated to store the generated
71  * inverse permutation on exit, or to provide the user permutation.
72  * On exit, the inverse permutation used by the solver is returned if
73  * invp is not NULL.
74  *
75  * @param[inout] b
76  * Array of size n -by- nrhs
77  * On entry, contains the nrhs vectors of the problem.
78  * On exit, contains the nrhs solution vectors of the problem.
79  *
80  * @param[in] nrhs
81  * The number of right hand side in the problem.
82  *
83  * @param[inout] iparm
84  * Array of size IPARM_SIZE
85  * On entry, contains all the integer parameters of the solver.
86  * On exit, the aray is updated with integer outputs of the solver.
87  *
88  * @param[inout] dparm
89  * Array of size DPARM_SIZE
90  * On entry, contains all the double parameters of the solver.
91  * On exit, the aray is updated with double outputs of the solver.
92  *
93  *******************************************************************************
94  *
95  * @retval PASTIX_SUCCESS on succesful exit,
96  * @retval PASTIX_ERR_BADPARAMETER on incorrect input parameter,
97  * @retval PASTIX_ERR_NOTIMPLEMENTED on variadic dofs,
98  * @retval PASTIX_ERR_UNKNOWN on undefined behaviors.
99  *
100  *******************************************************************************
101  */
102 int
103 pastix( pastix_data_t **pastix_data_ptr,
104  PASTIX_Comm pastix_comm,
105  pastix_int_t n,
106  pastix_int_t *colptr,
107  pastix_int_t *row,
108  void *avals,
109  pastix_int_t *perm,
110  pastix_int_t *invp,
111  void *b,
112  pastix_int_t nrhs,
113  pastix_int_t *iparm,
114  double *dparm )
115 {
116  pastix_data_t *pastix_data;
117  spmatrix_t *spm = NULL;
118  int ret;
119  size_t size;
120 
121  /*
122  * Initialize iparm/dparm to default values and exit when called with
123  * IPARM_MODIFY_PARAMETER set to 0
124  */
125  if (!iparm[IPARM_MODIFY_PARAMETER])
126  {
127  pastixInitParam(iparm, dparm);
128  iparm[IPARM_MODIFY_PARAMETER] = 1;
129  return PASTIX_SUCCESS;
130  }
131 
132  /*
133  * Initialization step
134  * Create the pastix_data structure and initialize the runtimes
135  */
136  if (iparm[IPARM_END_TASK] < PastixTaskInit) {
137  return PASTIX_SUCCESS;
138  }
139 
140  if (iparm[IPARM_START_TASK] == PastixTaskInit) {
141  if (*pastix_data_ptr != NULL)
142  {
143  /*
144  * Let's consider the user want to restart pastix with different
145  * parameters
146  */
147  if (iparm[IPARM_VERBOSE] > PastixVerboseNo) {
148  pastix_print( (*pastix_data_ptr)->inter_node_procnum, 0, "WARNING: PaStiX schedulers restarted\n" );
149  }
150  pastixFinalize( pastix_data_ptr );
151  }
152 
153  /*
154  * Initialize pastix_data structure, and start scheduler(s)
155  */
156  pastixInit( pastix_data_ptr, pastix_comm, iparm, dparm );
157 
158  iparm[IPARM_START_TASK]++;
159  }
160 
161  /*
162  * Return now if only initialization is required
163  */
164  if (iparm[IPARM_END_TASK] < PastixTaskOrdering) {
165  return PASTIX_SUCCESS;
166  }
167 
168  pastix_data = *pastix_data_ptr;
169 
170  /*
171  * Initialize the internal spm structure.
172  * We perform if only if starting step is lower than numerical
173  * factorization, because further steps are using the internal bcsc for
174  * computations with A.
175  */
176  if (iparm[IPARM_START_TASK] <= PastixTaskNumfact) {
177  if ( (pastix_data->csc != NULL) &&
178  ((pastix_data->csc->n != n) ||
179  (pastix_data->csc->nnz != (colptr[n] - colptr[0])) ||
180  (pastix_data->csc->colptr != colptr) ||
181  (pastix_data->csc->rowptr != row)) )
182  {
183  /*
184  * This is a new csc, we need to delete the old one stored in pastix_data, and create a new one
185  * We do not use spmExit, because the user allocated the fields.
186  */
187  memFree_null( pastix_data->csc );
188  }
189 
190  if ( pastix_data->csc == NULL )
191  {
192  /*
193  * Check and set the matrix type
194  */
195  if (iparm[IPARM_FLOAT] == -1) {
196  printf("Pastix old interface: you have to set iparm[IPARM_FLOAT]\n");
198  }
199  if (iparm[IPARM_MTX_TYPE] == -1) {
200  printf("Pastix old interface: you have to set iparm[IPARM_MTX_TYPE]\n");
202  }
203  if (iparm[IPARM_DOF_NBR] < 1) {
204  fprintf(stderr,
205  "pastix: Variadic dofs are not supported in old pastix interface.\n"
206  " Please switch to the new interface to use this feature, \n"
207  " or set to a positive value\n");
209  }
210 
211  spm = malloc(sizeof( spmatrix_t ));
212  spmInit( spm );
213 
214  spm->mtxtype = iparm[IPARM_MTX_TYPE];
215  spm->flttype = iparm[IPARM_FLOAT];
216  spm->fmttype = SpmCSC;
217 
218  spm->n = n;
219  spm->nnz = colptr[n] - colptr[0];
220  spm->dof = iparm[IPARM_DOF_NBR];
221 
222  spm->colptr = colptr;
223  spm->rowptr = row;
224  spm->values = avals;
225 
226  spm->baseval = spmFindBase(spm);
227  spmUpdateComputedFields( spm );
228 
229  pastix_data->csc = spm;
230  }
231  else {
232  /* Cast to overwrite the const */
233  spm = (spmatrix_t*)(pastix_data->csc);
234  }
235 
236  /*
237  * Update value field if given only at numerical steps
238  */
239  if ( spm->values == NULL ) {
240  spm->values = avals;
241  }
242  }
243  else {
244  /* Cast to overwrite the const */
245  spm = (spmatrix_t*)(pastix_data->csc);
246  }
247 
248  /*
249  * Ordering
250  */
251  if (iparm[IPARM_START_TASK] == PastixTaskOrdering)
252  {
253  pastix_order_t *o = NULL;
254 
255  if ( (perm != NULL) || (invp != NULL) ) {
256  o = malloc( sizeof(pastix_order_t) );
257  ret = pastixOrderAlloc( o, 0, 0 );
258  if (PASTIX_SUCCESS != ret) {
259  free(o);
260  return ret;
261  }
262 
263  if ( perm != NULL ) {
264  MALLOC_INTERN(o->permtab, n, pastix_int_t);
265  memcpy( o->permtab, perm, n * sizeof(pastix_int_t) );
266  o->vertnbr = n;
267  }
268  if ( invp != NULL ) {
269  MALLOC_INTERN(o->peritab, n, pastix_int_t);
270  memcpy( o->peritab, invp, n * sizeof(pastix_int_t) );
271  o->vertnbr = n;
272  }
273  }
274 
275  ret = pastix_subtask_order( pastix_data, spm, o );
276  if (PASTIX_SUCCESS != ret) {
277  if ( o != NULL ) {
278  pastixOrderExit(o);
279  free(o);
280  }
281  return ret;
282  }
283 
284  if ( o != NULL ) {
285  if ( perm != NULL ) {
286  assert( o->permtab != NULL );
287  assert( o->vertnbr == n );
288  memcpy( perm, o->permtab, n * sizeof(pastix_int_t));
289  }
290  if ( invp != NULL ) {
291  assert( o->peritab != NULL );
292  assert( o->vertnbr == n );
293  memcpy( invp, o->peritab, n * sizeof(pastix_int_t));
294  }
295  pastixOrderExit(o);
296  free(o);
297  }
298  iparm[IPARM_START_TASK]++;
299  }
300 
301  /*
302  * Symbolic factorization
303  */
304  if (iparm[IPARM_END_TASK] < PastixTaskSymbfact) {
305  return PASTIX_SUCCESS;
306  }
307 
308  if (iparm[IPARM_START_TASK] == PastixTaskSymbfact)
309  {
310  ret = pastix_subtask_symbfact( pastix_data );
311  if (PASTIX_SUCCESS != ret)
312  {
313  return ret;
314  }
315  iparm[IPARM_START_TASK]++;
316  }
317 
318  /*
319  * Analyze step
320  */
321  if (iparm[IPARM_END_TASK] < PastixTaskAnalyze) {
322  return PASTIX_SUCCESS;
323  }
324 
325  if (iparm[IPARM_START_TASK] == PastixTaskAnalyze)
326  {
327  ret = pastix_subtask_blend( pastix_data );
328  if (PASTIX_SUCCESS != ret)
329  {
330  return ret;
331  }
332  iparm[IPARM_START_TASK]++;
333  }
334 
335  /*
336  * Numerical factorization
337  */
338  if (iparm[IPARM_END_TASK] < PastixTaskNumfact) {
339  return PASTIX_SUCCESS;
340  }
341 
342  if (iparm[IPARM_START_TASK] == PastixTaskNumfact)
343  {
344  ret = pastix_task_numfact( pastix_data, spm );
345  if (PASTIX_SUCCESS != ret) {
346  return ret;
347  }
348  iparm[IPARM_START_TASK]++;
349  }
350 
351  /*
352  * Solve
353  */
354  if (iparm[IPARM_END_TASK] < PastixTaskSolve) {
355  return PASTIX_SUCCESS;
356  }
357 
358  if (iparm[IPARM_START_TASK] == PastixTaskSolve) {
359  size = pastix_size_of( spm->flttype ) * spm->n;
360  if ( pastix_data->x0 ) {
361  free(pastix_data->x0);
362  pastix_data->x0 = NULL;
363  }
364  if ( pastix_data->b ) {
365  free(pastix_data->b);
366  pastix_data->b = NULL;
367  }
368 
369  /*
370  * Backup the initial b if we need to perform an iterative
371  * refinement after the solve step
372  */
373  if (iparm[IPARM_END_TASK] > PastixTaskSolve) {
374  pastix_data->b = malloc(size);
375  memcpy(pastix_data->b, b, size);
376  }
377  pastix_task_solve( pastix_data, nrhs, b, spm->n );
378  iparm[IPARM_START_TASK]++;
379 
380  /*
381  * Backup the first solution x0 if the user wants to come back later for
382  * iterative refinement
383  */
384  if (iparm[IPARM_END_TASK] == PastixTaskSolve) {
385  pastix_data->x0 = malloc(size);
386  memcpy(pastix_data->x0, b, size);
387  }
388  }
389 
390  /*
391  * Refinement
392  */
393  if (iparm[IPARM_END_TASK] < PastixTaskRefine) {
394  return PASTIX_SUCCESS;
395  }
396 
397  if (iparm[IPARM_START_TASK] == PastixTaskRefine) {
398  void *refineB = pastix_data->b;
399  void *refineX0 = pastix_data->x0;
400  size = pastix_size_of( spm->flttype ) * spm->n;
401  if ( !refineB ) {
402  if ( !refineX0 ) {
403  /*
404  * Neither b or x0 have been saved.
405  * Then, we need to start with x0 as a null vector. For that, we
406  * backup the original b, and we use the given b as x in the
407  * refinement step to store the solution in place.
408  */
409  /* refineB = malloc(size); */
410  /* memcpy(refineB, b, size); */
411 
412  /* refineX0 = b; */
413  /* memset(refineX0, 0, size); */
414  /* exit(0); */
415  /*
416  * Neither b and x0 have been saved, this should never happen.
417  */
418  fprintf(stderr, "Neither b and x0 have been saved, this should never happen\n");
419  return PASTIX_ERR_UNKNOWN;
420  }
421  else {
422  /*
423  * x0 is saved, but not b. It means that we exit the pastix
424  * function call between the solve and refinemnet
425  * step. Therefor, b holds the original b.
426  */
427  refineB = b;
428  }
429  }
430  else {
431  if ( !refineX0 ) {
432  /*
433  * b is saved, but not x0. It means that we did not exit the
434  * pastix function call between solve and refinement steps.
435  * Therefor, b holds the initial solution x0 from the solve step.
436  */
437  refineX0 = b;
438  }
439  else {
440  /*
441  * Both x0 and b are saved. This should never happen.
442  */
443  fprintf(stderr, "Both b and x0 are defined, this should never happen\n");
444  return PASTIX_ERR_UNKNOWN;
445  }
446  }
447  pastix_task_refine( pastix_data, spm->n, nrhs,
448  refineB, spm->n, refineX0, spm->n );
449  iparm[IPARM_START_TASK]++;
450 
451  /*
452  * Let's return the solution to the user
453  */
454  if ( b != refineX0 ) {
455  memcpy(b, refineB, size);
456  }
457  }
458 
459  if ( pastix_data->x0 != NULL ) {
460  free( pastix_data->x0 );
461  pastix_data->x0 = NULL;
462  }
463  if ( pastix_data->b != NULL ) {
464  free( pastix_data->b );
465  pastix_data->b = NULL;
466  }
467 
468  /*
469  * Cleaning
470  */
471  if (iparm[IPARM_END_TASK] < PastixTaskClean) {
472  return PASTIX_SUCCESS;
473  }
474 
475  if (iparm[IPARM_START_TASK] == PastixTaskClean) {
476  if ( pastix_data->csc != NULL ) {
477  spmatrix_t *spm = (spmatrix_t*)(pastix_data->csc);
478  free( spm );
479  }
480  pastixFinalize( pastix_data_ptr );
481  iparm[IPARM_START_TASK]++;
482  }
483 
484  return PASTIX_SUCCESS;
485 }
486 
487 
488 /**
489  *******************************************************************************
490  *
491  * @brief Expand an spm structure and the already computed data structure
492  * associated if any.
493  *
494  *******************************************************************************
495  *
496  * @param[in] pastix_data
497  * The pastix context in which the spm is used
498  *
499  * @param[in] spm
500  * The multi-dof sparse matrix to expand into a single dof sparse matrix.
501  * On exit, spm contains the expanded matrix. The compressed form of
502  * the matrix is destroyed.
503  *
504  *******************************************************************************/
505 void
506 pastixExpand( const pastix_data_t *pastix_data,
507  spmatrix_t *spm )
508 {
509  spmatrix_t tmp;
510 
511  if ( spm == NULL ) {
512  return;
513  }
514 
515  /* Nothing to expand */
516  if ( spm->dof == 1 ) {
517  return;
518  }
519 
520  if ( ( pastix_data->steps & STEP_ORDERING ) &&
521  ( pastix_data->ordemesh != NULL ) )
522  {
523  pastixOrderExpand( pastix_data->ordemesh, spm );
524 #if !defined(NDEBUG)
525  pastixOrderCheck( pastix_data->ordemesh );
526 #endif
527  }
528  if ( ( pastix_data->steps & STEP_SYMBFACT ) &&
529  ( pastix_data->symbmtx != NULL ) )
530  {
531  pastixSymbolExpand( pastix_data->symbmtx );
532  }
533 
534  spmExpand( spm, &tmp );
535  spmExit( spm );
536  memcpy( spm, &tmp, sizeof(spmatrix_t) );
537 }
pastix_order_s::peritab
pastix_int_t * peritab
Definition: order.h:50
PastixTaskNumfact
@ PastixTaskNumfact
Definition: api.h:189
pastix_order_s::permtab
pastix_int_t * permtab
Definition: order.h:49
pastixOrderExit
void pastixOrderExit(pastix_order_t *const ordeptr)
Free the arrays initialized in the order structure.
Definition: order.c:273
pastix_order_s
Order structure.
Definition: order.h:45
pastix_task_solve
int pastix_task_solve(pastix_data_t *pastix_data, pastix_int_t nrhs, void *b, pastix_int_t ldb)
Solve the given problem.
Definition: pastix_task_solve.c:681
pastix_order_s::vertnbr
pastix_int_t vertnbr
Definition: order.h:47
pastix_subtask_symbfact
int pastix_subtask_symbfact(pastix_data_t *pastix_data)
Computes the symbolic factorization step.
Definition: pastix_subtask_symbfact.c:90
pastixSymbolExpand
void pastixSymbolExpand(symbol_matrix_t *symbptr)
Expand the symbol matrix structure based on the dof information (compressed -> expanded)
Definition: symbol_expand.c:181
pastixOrderCheck
int pastixOrderCheck(const pastix_order_t *const ordeptr)
This routine checks the correctness of the ordering structure.
Definition: order_check.c:39
PastixTaskRefine
@ PastixTaskRefine
Definition: api.h:191
pastix_subtask_blend
int pastix_subtask_blend(pastix_data_t *pastix_data)
Compute the proportional mapping and the final solver structure.
Definition: pastix_subtask_blend.c:121
PastixTaskAnalyze
@ PastixTaskAnalyze
Definition: api.h:188
PastixTaskInit
@ PastixTaskInit
Definition: api.h:185
PASTIX_SUCCESS
@ PASTIX_SUCCESS
Definition: api.h:344
PastixVerboseNo
@ PastixVerboseNo
Definition: api.h:208
PASTIX_ERR_UNKNOWN
@ PASTIX_ERR_UNKNOWN
Definition: api.h:345
PastixTaskSymbfact
@ PastixTaskSymbfact
Definition: api.h:187
pastixExpand
void pastixExpand(const pastix_data_t *pastix_data, spmatrix_t *spm)
Expand an spm structure and the already computed data structure associated if any.
Definition: pastix.c:506
IPARM_MTX_TYPE
@ IPARM_MTX_TYPE
Definition: api.h:143
pastixOrderAlloc
int pastixOrderAlloc(pastix_order_t *const ordeptr, pastix_int_t vertnbr, pastix_int_t cblknbr)
Allocate the order structure.
Definition: order.c:55
pastixInit
void pastixInit(pastix_data_t **pastix_data, PASTIX_Comm pastix_comm, pastix_int_t *iparm, double *dparm)
Initialize the solver instance.
Definition: api.c:857
IPARM_END_TASK
@ IPARM_END_TASK
Definition: api.h:141
IPARM_FLOAT
@ IPARM_FLOAT
Definition: api.h:142
PastixTaskClean
@ PastixTaskClean
Definition: api.h:192
PastixTaskSolve
@ PastixTaskSolve
Definition: api.h:190
IPARM_MODIFY_PARAMETER
@ IPARM_MODIFY_PARAMETER
Definition: api.h:139
pastixOrderExpand
void pastixOrderExpand(pastix_order_t *const ordeptr, spmatrix_t *const spm)
This routine expand the permutation arrays and the rangtab when the spm is using multiple dof per unk...
Definition: order.c:393
PastixTaskOrdering
@ PastixTaskOrdering
Definition: api.h:186
pastix
int pastix(pastix_data_t **pastix_data_ptr, PASTIX_Comm pastix_comm, pastix_int_t n, pastix_int_t *colptr, pastix_int_t *row, void *avals, pastix_int_t *perm, pastix_int_t *invp, void *b, pastix_int_t nrhs, pastix_int_t *iparm, double *dparm)
Main function for compatibility with former releases.
Definition: pastix.c:103
IPARM_VERBOSE
@ IPARM_VERBOSE
Definition: api.h:36
order.h
IPARM_START_TASK
@ IPARM_START_TASK
Definition: api.h:140
pastixInitParam
void pastixInitParam(pastix_int_t *iparm, double *dparm)
Initialize the iparm and dparm arrays to their default values.
Definition: api.c:401
IPARM_DOF_NBR
@ IPARM_DOF_NBR
Definition: api.h:144
pastixFinalize
void pastixFinalize(pastix_data_t **pastix_data)
Finalize the solver instance.
Definition: api.c:880
pastix_subtask_order
int pastix_subtask_order(pastix_data_t *pastix_data, const spmatrix_t *spm, pastix_order_t *myorder)
Computes the ordering of the given graph in parameters.
Definition: pastix_subtask_order.c:97
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
PASTIX_ERR_BADPARAMETER
@ PASTIX_ERR_BADPARAMETER
Definition: api.h:351
pastix_task_numfact
int pastix_task_numfact(pastix_data_t *pastix_data, spmatrix_t *spm)
Perform all the numerical factorization steps: fill the internal block CSC and the solver matrix stru...
Definition: pastix_task_sopalin.c:587