PaStiX Handbook 6.4.0
Loading...
Searching...
No Matches
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-2024 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
8 * Univ. Bordeaux. All rights reserved.
9 *
10 * @version 6.4.0
11 * @author Mathieu Faverge
12 * @author Pierre Ramet
13 * @author Xavier Lacoste
14 * @author Gregoire Pichon
15 * @author Matias Hastaran
16 * @date 2024-07-05
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 */
102int
103pastix( 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 if ( pastix_data == NULL ) {
170 printf("Pastix old interface: pastixTaskInit step not called before calling following steps\n");
172 }
173
174 /*
175 * Initialize the internal spm structure.
176 * We perform if only if starting step is lower than numerical
177 * factorization, because further steps are using the internal bcsc for
178 * computations with A.
179 */
180 if (iparm[IPARM_START_TASK] <= PastixTaskNumfact) {
181 if ( (pastix_data->csc != NULL) &&
182 ((pastix_data->csc->n != n) ||
183 (pastix_data->csc->nnz != (colptr[n] - colptr[0])) ||
184 (pastix_data->csc->colptr != colptr) ||
185 (pastix_data->csc->rowptr != row)) )
186 {
187 /*
188 * This is a new csc, we need to delete the old one stored in pastix_data, and create a new one
189 * We do not use spmExit, because the user allocated the fields.
190 */
191 memFree_null( pastix_data->csc );
192 }
193
194 if ( pastix_data->csc == NULL )
195 {
196 /*
197 * Check and set the matrix type
198 */
199 if (iparm[IPARM_FLOAT] == -1) {
200 printf("Pastix old interface: you have to set iparm[IPARM_FLOAT]\n");
202 }
203 if (iparm[IPARM_MTX_TYPE] == -1) {
204 printf("Pastix old interface: you have to set iparm[IPARM_MTX_TYPE]\n");
206 }
207 if (iparm[IPARM_DOF_NBR] < 1) {
208 fprintf(stderr,
209 "pastix: Variadic dofs are not supported in old pastix interface.\n"
210 " Please switch to the new interface to use this feature, \n"
211 " or set to a positive value\n");
213 }
214
215 spm = malloc(sizeof( spmatrix_t ));
216 spmInit( spm );
217
218 spm->mtxtype = iparm[IPARM_MTX_TYPE];
219 spm->flttype = iparm[IPARM_FLOAT];
220 spm->fmttype = SpmCSC;
221
222 spm->n = n;
223 spm->nnz = colptr[n] - colptr[0];
224 spm->dof = iparm[IPARM_DOF_NBR];
225
226 spm->colptr = colptr;
227 spm->rowptr = row;
228 spm->values = avals;
229
230 spm->baseval = spmFindBase(spm);
231 spmUpdateComputedFields( spm );
232
233 pastix_data->csc = spm;
234 }
235 else {
236 /* Cast to overwrite the const */
237 spm = (spmatrix_t*)(pastix_data->csc);
238 }
239
240 /*
241 * Update value field if given only at numerical steps
242 */
243 if ( spm->values == NULL ) {
244 spm->values = avals;
245 }
246 }
247 else {
248 /* Cast to overwrite the const */
249 spm = (spmatrix_t*)(pastix_data->csc);
250 }
251
252 /*
253 * Ordering
254 */
256 {
257 pastix_order_t *o = NULL;
258
259 if ( (perm != NULL) || (invp != NULL) ) {
260 o = malloc( sizeof(pastix_order_t) );
261 ret = pastixOrderAlloc( o, 0, 0 );
262 if (PASTIX_SUCCESS != ret) {
263 free(o);
264 return ret;
265 }
266
267 if ( perm != NULL ) {
268 MALLOC_INTERN(o->permtab, n, pastix_int_t);
269 memcpy( o->permtab, perm, n * sizeof(pastix_int_t) );
270 o->vertnbr = n;
271 }
272 if ( invp != NULL ) {
273 MALLOC_INTERN(o->peritab, n, pastix_int_t);
274 memcpy( o->peritab, invp, n * sizeof(pastix_int_t) );
275 o->vertnbr = n;
276 }
277 }
278
279 ret = pastix_subtask_order( pastix_data, spm, o );
280 if (PASTIX_SUCCESS != ret) {
281 if ( o != NULL ) {
283 free(o);
284 }
285 return ret;
286 }
287
288 if ( o != NULL ) {
289 if ( perm != NULL ) {
290 assert( o->permtab != NULL );
291 assert( o->vertnbr == n );
292 memcpy( perm, o->permtab, n * sizeof(pastix_int_t));
293 }
294 if ( invp != NULL ) {
295 assert( o->peritab != NULL );
296 assert( o->vertnbr == n );
297 memcpy( invp, o->peritab, n * sizeof(pastix_int_t));
298 }
300 free(o);
301 }
302 iparm[IPARM_START_TASK]++;
303 }
304
305 /*
306 * Symbolic factorization
307 */
308 if (iparm[IPARM_END_TASK] < PastixTaskSymbfact) {
309 return PASTIX_SUCCESS;
310 }
311
313 {
314 ret = pastix_subtask_symbfact( pastix_data );
315 if (PASTIX_SUCCESS != ret)
316 {
317 return ret;
318 }
319 iparm[IPARM_START_TASK]++;
320 }
321
322 /*
323 * Analyze step
324 */
325 if (iparm[IPARM_END_TASK] < PastixTaskAnalyze) {
326 return PASTIX_SUCCESS;
327 }
328
330 {
331 ret = pastix_subtask_blend( pastix_data );
332 if (PASTIX_SUCCESS != ret)
333 {
334 return ret;
335 }
336 iparm[IPARM_START_TASK]++;
337 }
338
339 /*
340 * Numerical factorization
341 */
342 if (iparm[IPARM_END_TASK] < PastixTaskNumfact) {
343 return PASTIX_SUCCESS;
344 }
345
347 {
348 ret = pastix_task_numfact( pastix_data, spm );
349 if (PASTIX_SUCCESS != ret) {
350 return ret;
351 }
352 iparm[IPARM_START_TASK]++;
353 }
354
355 /*
356 * Solve
357 */
358 if (iparm[IPARM_END_TASK] < PastixTaskSolve) {
359 return PASTIX_SUCCESS;
360 }
361
362 if (iparm[IPARM_START_TASK] == PastixTaskSolve) {
363 size = pastix_size_of( spm->flttype ) * spm->n;
364 if ( pastix_data->x0 ) {
365 free(pastix_data->x0);
366 pastix_data->x0 = NULL;
367 }
368 if ( pastix_data->b ) {
369 free(pastix_data->b);
370 pastix_data->b = NULL;
371 }
372
373 /*
374 * Backup the initial b if we need to perform an iterative
375 * refinement after the solve step
376 */
377 if (iparm[IPARM_END_TASK] > PastixTaskSolve) {
378 pastix_data->b = malloc(size);
379 memcpy(pastix_data->b, b, size);
380 }
381 pastix_task_solve( pastix_data, spm->nexp, nrhs, b, spm->nexp );
382 iparm[IPARM_START_TASK]++;
383
384 /*
385 * Backup the first solution x0 if the user wants to come back later for
386 * iterative refinement
387 */
388 if (iparm[IPARM_END_TASK] == PastixTaskSolve) {
389 pastix_data->x0 = malloc(size);
390 memcpy(pastix_data->x0, b, size);
391 }
392 }
393
394 /*
395 * Refinement
396 */
397 if (iparm[IPARM_END_TASK] < PastixTaskRefine) {
398 return PASTIX_SUCCESS;
399 }
400
401 if (iparm[IPARM_START_TASK] == PastixTaskRefine) {
402 void *refineB = pastix_data->b;
403 void *refineX0 = pastix_data->x0;
404 size = pastix_size_of( spm->flttype ) * spm->n;
405 if ( !refineB ) {
406 if ( !refineX0 ) {
407 /*
408 * Neither b or x0 have been saved.
409 * Then, we need to start with x0 as a null vector. For that, we
410 * backup the original b, and we use the given b as x in the
411 * refinement step to store the solution in place.
412 */
413 /* refineB = malloc(size); */
414 /* memcpy(refineB, b, size); */
415
416 /* refineX0 = b; */
417 /* memset(refineX0, 0, size); */
418 /* exit(0); */
419 /*
420 * Neither b and x0 have been saved, this should never happen.
421 */
422 fprintf(stderr, "Neither b and x0 have been saved, this should never happen\n");
423 return PASTIX_ERR_UNKNOWN;
424 }
425 else {
426 /*
427 * x0 is saved, but not b. It means that we exit the pastix
428 * function call between the solve and refinemnet
429 * step. Therefor, b holds the original b.
430 */
431 refineB = b;
432 }
433 }
434 else {
435 if ( !refineX0 ) {
436 /*
437 * b is saved, but not x0. It means that we did not exit the
438 * pastix function call between solve and refinement steps.
439 * Therefor, b holds the initial solution x0 from the solve step.
440 */
441 refineX0 = b;
442 }
443 else {
444 /*
445 * Both x0 and b are saved. This should never happen.
446 */
447 fprintf(stderr, "Both b and x0 are defined, this should never happen\n");
448 return PASTIX_ERR_UNKNOWN;
449 }
450 }
451 pastix_task_refine( pastix_data, spm->n, nrhs,
452 refineB, spm->n, refineX0, spm->n );
453 iparm[IPARM_START_TASK]++;
454
455 /*
456 * Let's return the solution to the user
457 */
458 if ( b != refineX0 ) {
459 memcpy(b, refineB, size);
460 }
461 }
462
463 if ( pastix_data->x0 != NULL ) {
464 free( pastix_data->x0 );
465 pastix_data->x0 = NULL;
466 }
467 if ( pastix_data->b != NULL ) {
468 free( pastix_data->b );
469 pastix_data->b = NULL;
470 }
471
472 /*
473 * Cleaning
474 */
475 if (iparm[IPARM_END_TASK] < PastixTaskClean) {
476 return PASTIX_SUCCESS;
477 }
478
479 if (iparm[IPARM_START_TASK] == PastixTaskClean) {
480 if ( pastix_data->csc != NULL ) {
481 spmatrix_t *spm = (spmatrix_t*)(pastix_data->csc);
482 free( spm );
483 }
484 pastixFinalize( pastix_data_ptr );
485 iparm[IPARM_START_TASK]++;
486 }
487
488 return PASTIX_SUCCESS;
489}
490
491
492/**
493 *******************************************************************************
494 *
495 * @brief Expand an spm structure and the already computed data structure
496 * associated if any.
497 *
498 *******************************************************************************
499 *
500 * @param[in] pastix_data
501 * The pastix context in which the spm is used
502 *
503 * @param[in] spm
504 * The multi-dof sparse matrix to expand into a single dof sparse matrix.
505 * On exit, spm contains the expanded matrix. The compressed form of
506 * the matrix is destroyed.
507 *
508 *******************************************************************************/
509void
510pastixExpand( const pastix_data_t *pastix_data,
511 spmatrix_t *spm )
512{
513 spmatrix_t tmp;
514
515 if ( spm == NULL ) {
516 return;
517 }
518
519 /* Nothing to expand */
520 if ( spm->dof == 1 ) {
521 return;
522 }
523
524 if ( ( pastix_data->steps & STEP_ORDERING ) &&
525 ( pastix_data->ordemesh != NULL ) )
526 {
527 pastixOrderExpand( pastix_data->ordemesh, spm );
528#if !defined(NDEBUG)
529 pastixOrderCheck( pastix_data->ordemesh );
530#endif
531 }
532 if ( ( pastix_data->steps & STEP_SYMBFACT ) &&
533 ( pastix_data->symbmtx != NULL ) )
534 {
535 pastixSymbolExpand( pastix_data->symbmtx );
536 }
537
538 spmExpand( spm, &tmp );
539 spmExit( spm );
540 memcpy( spm, &tmp, sizeof(spmatrix_t) );
541}
BEGIN_C_DECLS typedef int pastix_int_t
Definition datatypes.h:51
int pastix_subtask_symbfact(pastix_data_t *pastix_data)
Computes the symbolic factorization step.
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.
int pastix_subtask_blend(pastix_data_t *pastix_data)
Compute the proportional mapping and the final solver structure.
void pastixFinalize(pastix_data_t **pastix_data)
Finalize the solver instance.
Definition api.c:928
void pastixInitParam(pastix_int_t *iparm, double *dparm)
Initialize the iparm and dparm arrays to their default values.
Definition api.c:420
void pastixInit(pastix_data_t **pastix_data, PASTIX_Comm pastix_comm, pastix_int_t *iparm, double *dparm)
Initialize the solver instance.
Definition api.c:905
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
@ PastixTaskRefine
Definition api.h:202
@ PastixTaskInit
Definition api.h:196
@ PastixTaskSymbfact
Definition api.h:198
@ PastixTaskOrdering
Definition api.h:197
@ PastixTaskSolve
Definition api.h:201
@ PastixTaskAnalyze
Definition api.h:199
@ PastixTaskClean
Definition api.h:203
@ PastixTaskNumfact
Definition api.h:200
@ IPARM_MODIFY_PARAMETER
Definition api.h:146
@ IPARM_MTX_TYPE
Definition api.h:150
@ IPARM_START_TASK
Definition api.h:147
@ IPARM_DOF_NBR
Definition api.h:151
@ IPARM_FLOAT
Definition api.h:149
@ IPARM_VERBOSE
Definition api.h:36
@ IPARM_END_TASK
Definition api.h:148
@ PastixVerboseNo
Definition api.h:221
@ PASTIX_ERR_UNKNOWN
Definition api.h:368
@ PASTIX_SUCCESS
Definition api.h:367
@ PASTIX_ERR_BADPARAMETER
Definition api.h:374
pastix_int_t * permtab
Definition order.h:51
pastix_int_t * peritab
Definition order.h:52
pastix_int_t vertnbr
Definition order.h:49
void pastixOrderExpand(pastix_order_t *ordeptr, const spmatrix_t *spm)
This routine expand the permutation arrays and the rangtab when the spm is using multiple dof per unk...
Definition order.c:396
int pastixOrderAlloc(pastix_order_t *ordeptr, pastix_int_t vertnbr, pastix_int_t cblknbr)
Allocate the order structure.
Definition order.c:55
int pastixOrderCheck(const pastix_order_t *ordeptr)
This routine checks the correctness of the ordering structure.
Definition order_check.c:39
void pastixOrderExit(pastix_order_t *ordeptr)
Free the arrays initialized in the order structure.
Definition order.c:273
Order structure.
Definition order.h:47
void pastixSymbolExpand(symbol_matrix_t *symbptr)
Expand the symbol matrix structure based on the dof information (compressed -> expanded)
pastix_order_t * ordemesh
Definition pastixdata.h:98
const spmatrix_t * csc
Definition pastixdata.h:90
symbol_matrix_t * symbmtx
Definition pastixdata.h:100
pastix_int_t steps
Definition pastixdata.h:73
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.
int pastix_task_solve(pastix_data_t *pastix_data, pastix_int_t m, pastix_int_t nrhs, void *B, pastix_int_t ldb)
Solve the given problem.
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...
Main PaStiX data structure.
Definition pastixdata.h:68
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:510