xref: /petsc/src/mat/impls/aij/mpi/pastix/pastix.c (revision ef46b1a67e276116c83b5d4ce8efc2932ea4fc0a) 
1 /*
2  Provides an interface to the PaStiX sparse solver
3  */
4 #include <../src/mat/impls/aij/seq/aij.h>
5 #include <../src/mat/impls/aij/mpi/mpiaij.h>
6 #include <../src/mat/impls/sbaij/seq/sbaij.h>
7 #include <../src/mat/impls/sbaij/mpi/mpisbaij.h>
8 
9 #if defined(PETSC_USE_COMPLEX)
10 #define _H_COMPLEX
11 #endif
12 
13 EXTERN_C_BEGIN
14 #include <pastix.h>
15 EXTERN_C_END
16 
17 #if defined(PETSC_USE_COMPLEX)
18 #if defined(PETSC_USE_REAL_SINGLE)
19 #define PASTIX_CALL c_pastix
20 #else
21 #define PASTIX_CALL z_pastix
22 #endif
23 
24 #else /* PETSC_USE_COMPLEX */
25 
26 #if defined(PETSC_USE_REAL_SINGLE)
27 #define PASTIX_CALL s_pastix
28 #else
29 #define PASTIX_CALL d_pastix
30 #endif
31 
32 #endif /* PETSC_USE_COMPLEX */
33 
34 typedef PetscScalar PastixScalar;
35 
36 typedef struct Mat_Pastix_ {
37   pastix_data_t *pastix_data;    /* Pastix data storage structure                        */
38   MatStructure  matstruc;
39   PetscInt      n;               /* Number of columns in the matrix                      */
40   PetscInt      *colptr;         /* Index of first element of each column in row and val */
41   PetscInt      *row;            /* Row of each element of the matrix                    */
42   PetscScalar   *val;            /* Value of each element of the matrix                  */
43   PetscInt      *perm;           /* Permutation tabular                                  */
44   PetscInt      *invp;           /* Reverse permutation tabular                          */
45   PetscScalar   *rhs;            /* Rhight-hand-side member                              */
46   PetscInt      rhsnbr;          /* Rhight-hand-side number (must be 1)                  */
47   PetscInt      iparm[IPARM_SIZE];       /* Integer parameters                                   */
48   double        dparm[DPARM_SIZE];       /* Floating point parameters                            */
49   MPI_Comm      pastix_comm;     /* PaStiX MPI communicator                              */
50   PetscMPIInt   commRank;        /* MPI rank                                             */
51   PetscMPIInt   commSize;        /* MPI communicator size                                */
52   PetscBool     CleanUpPastix;   /* Boolean indicating if we call PaStiX clean step      */
53   VecScatter    scat_rhs;
54   VecScatter    scat_sol;
55   Vec           b_seq;
56 } Mat_Pastix;
57 
58 extern PetscErrorCode MatDuplicate_Pastix(Mat,MatDuplicateOption,Mat*);
59 
60 /*
61    convert Petsc seqaij matrix to CSC: colptr[n], row[nz], val[nz]
62 
63   input:
64     A       - matrix in seqaij or mpisbaij (bs=1) format
65     valOnly - FALSE: spaces are allocated and values are set for the CSC
66               TRUE:  Only fill values
67   output:
68     n       - Size of the matrix
69     colptr  - Index of first element of each column in row and val
70     row     - Row of each element of the matrix
71     values  - Value of each element of the matrix
72  */
73 PetscErrorCode MatConvertToCSC(Mat A,PetscBool valOnly,PetscInt *n,PetscInt **colptr,PetscInt **row,PetscScalar **values)
74 {
75   Mat_SeqAIJ  *aa      = (Mat_SeqAIJ*)A->data;
76   PetscInt    *rowptr  = aa->i;
77   PetscInt    *col     = aa->j;
78   PetscScalar *rvalues = aa->a;
79   PetscInt     m       = A->rmap->N;
80   PetscInt     nnz;
81   PetscInt     i,j, k;
82   PetscInt     base    = 1;
83   PetscInt     idx;
84   PetscInt     colidx;
85   PetscInt    *colcount;
86   PetscBool    isSBAIJ;
87   PetscBool    isSeqSBAIJ;
88   PetscBool    isMpiSBAIJ;
89   PetscBool    isSym;
90 
91   PetscFunctionBegin;
92   PetscCall(MatIsSymmetric(A,0.0,&isSym));
93   PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSBAIJ,&isSBAIJ));
94   PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ));
95   PetscCall(PetscObjectTypeCompare((PetscObject)A,MATMPISBAIJ,&isMpiSBAIJ));
96 
97   *n = A->cmap->N;
98 
99   /* PaStiX only needs triangular matrix if matrix is symmetric
100    */
101   if (isSym && !(isSBAIJ || isSeqSBAIJ || isMpiSBAIJ)) nnz = (aa->nz - *n)/2 + *n;
102   else nnz = aa->nz;
103 
104   if (!valOnly) {
105     PetscCall(PetscMalloc1((*n)+1,colptr));
106     PetscCall(PetscMalloc1(nnz,row));
107     PetscCall(PetscMalloc1(nnz,values));
108 
109     if (isSBAIJ || isSeqSBAIJ || isMpiSBAIJ) {
110       PetscCall(PetscArraycpy (*colptr, rowptr, (*n)+1));
111       for (i = 0; i < *n+1; i++) (*colptr)[i] += base;
112       PetscCall(PetscArraycpy (*row, col, nnz));
113       for (i = 0; i < nnz; i++) (*row)[i] += base;
114       PetscCall(PetscArraycpy (*values, rvalues, nnz));
115     } else {
116       PetscCall(PetscMalloc1(*n,&colcount));
117 
118       for (i = 0; i < m; i++) colcount[i] = 0;
119       /* Fill-in colptr */
120       for (i = 0; i < m; i++) {
121         for (j = rowptr[i]; j < rowptr[i+1]; j++) {
122           if (!isSym || col[j] <= i)  colcount[col[j]]++;
123         }
124       }
125 
126       (*colptr)[0] = base;
127       for (j = 0; j < *n; j++) {
128         (*colptr)[j+1] = (*colptr)[j] + colcount[j];
129         /* in next loop we fill starting from (*colptr)[colidx] - base */
130         colcount[j] = -base;
131       }
132 
133       /* Fill-in rows and values */
134       for (i = 0; i < m; i++) {
135         for (j = rowptr[i]; j < rowptr[i+1]; j++) {
136           if (!isSym || col[j] <= i) {
137             colidx         = col[j];
138             idx            = (*colptr)[colidx] + colcount[colidx];
139             (*row)[idx]    = i + base;
140             (*values)[idx] = rvalues[j];
141             colcount[colidx]++;
142           }
143         }
144       }
145       PetscCall(PetscFree(colcount));
146     }
147   } else {
148     /* Fill-in only values */
149     for (i = 0; i < m; i++) {
150       for (j = rowptr[i]; j < rowptr[i+1]; j++) {
151         colidx = col[j];
152         if ((isSBAIJ || isSeqSBAIJ || isMpiSBAIJ) ||!isSym || col[j] <= i) {
153           /* look for the value to fill */
154           for (k = (*colptr)[colidx] - base; k < (*colptr)[colidx + 1] - base; k++) {
155             if (((*row)[k]-base) == i) {
156               (*values)[k] = rvalues[j];
157               break;
158             }
159           }
160           /* data structure of sparse matrix has changed */
161           PetscCheckFalse(k == (*colptr)[colidx + 1] - base,PETSC_COMM_SELF,PETSC_ERR_PLIB,"overflow on k %" PetscInt_FMT,k);
162         }
163       }
164     }
165   }
166   PetscFunctionReturn(0);
167 }
168 
169 /*
170   Call clean step of PaStiX if lu->CleanUpPastix == true.
171   Free the CSC matrix.
172  */
173 PetscErrorCode MatDestroy_Pastix(Mat A)
174 {
175   Mat_Pastix *lu = (Mat_Pastix*)A->data;
176 
177   PetscFunctionBegin;
178   if (lu->CleanUpPastix) {
179     /* Terminate instance, deallocate memories */
180     PetscCall(VecScatterDestroy(&lu->scat_rhs));
181     PetscCall(VecDestroy(&lu->b_seq));
182     PetscCall(VecScatterDestroy(&lu->scat_sol));
183 
184     lu->iparm[IPARM_START_TASK]=API_TASK_CLEAN;
185     lu->iparm[IPARM_END_TASK]  =API_TASK_CLEAN;
186 
187     PASTIX_CALL(&(lu->pastix_data),
188                 lu->pastix_comm,
189                 lu->n,
190                 lu->colptr,
191                 lu->row,
192                 (PastixScalar*)lu->val,
193                 lu->perm,
194                 lu->invp,
195                 (PastixScalar*)lu->rhs,
196                 lu->rhsnbr,
197                 lu->iparm,
198                 lu->dparm);
199     PetscCheck(lu->iparm[IPARM_ERROR_NUMBER] == 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by PaStiX in destroy: iparm(IPARM_ERROR_NUMBER)=%" PetscInt_FMT,lu->iparm[IPARM_ERROR_NUMBER]);
200     PetscCall(PetscFree(lu->colptr));
201     PetscCall(PetscFree(lu->row));
202     PetscCall(PetscFree(lu->val));
203     PetscCall(PetscFree(lu->perm));
204     PetscCall(PetscFree(lu->invp));
205     PetscCallMPI(MPI_Comm_free(&(lu->pastix_comm)));
206   }
207   PetscCall(PetscFree(A->data));
208   PetscFunctionReturn(0);
209 }
210 
211 /*
212   Gather right-hand-side.
213   Call for Solve step.
214   Scatter solution.
215  */
216 PetscErrorCode MatSolve_PaStiX(Mat A,Vec b,Vec x)
217 {
218   Mat_Pastix  *lu = (Mat_Pastix*)A->data;
219   PetscScalar *array;
220   Vec          x_seq;
221 
222   PetscFunctionBegin;
223   lu->rhsnbr = 1;
224   x_seq      = lu->b_seq;
225   if (lu->commSize > 1) {
226     /* PaStiX only supports centralized rhs. Scatter b into a sequential rhs vector */
227     PetscCall(VecScatterBegin(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD));
228     PetscCall(VecScatterEnd(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD));
229     PetscCall(VecGetArray(x_seq,&array));
230   } else {  /* size == 1 */
231     PetscCall(VecCopy(b,x));
232     PetscCall(VecGetArray(x,&array));
233   }
234   lu->rhs = array;
235   if (lu->commSize == 1) {
236     PetscCall(VecRestoreArray(x,&array));
237   } else {
238     PetscCall(VecRestoreArray(x_seq,&array));
239   }
240 
241   /* solve phase */
242   /*-------------*/
243   lu->iparm[IPARM_START_TASK] = API_TASK_SOLVE;
244   lu->iparm[IPARM_END_TASK]   = API_TASK_REFINE;
245   lu->iparm[IPARM_RHS_MAKING] = API_RHS_B;
246 
247   PASTIX_CALL(&(lu->pastix_data),
248               lu->pastix_comm,
249               lu->n,
250               lu->colptr,
251               lu->row,
252               (PastixScalar*)lu->val,
253               lu->perm,
254               lu->invp,
255               (PastixScalar*)lu->rhs,
256               lu->rhsnbr,
257               lu->iparm,
258               lu->dparm);
259   PetscCheck(lu->iparm[IPARM_ERROR_NUMBER] == 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by PaStiX in solve phase: lu->iparm[IPARM_ERROR_NUMBER] = %" PetscInt_FMT,lu->iparm[IPARM_ERROR_NUMBER]);
260 
261   if (lu->commSize == 1) {
262     PetscCall(VecRestoreArray(x,&(lu->rhs)));
263   } else {
264     PetscCall(VecRestoreArray(x_seq,&(lu->rhs)));
265   }
266 
267   if (lu->commSize > 1) { /* convert PaStiX centralized solution to petsc mpi x */
268     PetscCall(VecScatterBegin(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD));
269     PetscCall(VecScatterEnd(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD));
270   }
271   PetscFunctionReturn(0);
272 }
273 
274 /*
275   Numeric factorisation using PaStiX solver.
276 
277  */
278 PetscErrorCode MatFactorNumeric_PaStiX(Mat F,Mat A,const MatFactorInfo *info)
279 {
280   Mat_Pastix     *lu =(Mat_Pastix*)(F)->data;
281   Mat            *tseq;
282   PetscInt       icntl;
283   PetscInt       M=A->rmap->N;
284   PetscBool      valOnly,flg, isSym;
285   IS             is_iden;
286   Vec            b;
287   IS             isrow;
288   PetscBool      isSeqAIJ,isSeqSBAIJ,isMPIAIJ;
289 
290   PetscFunctionBegin;
291   PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ));
292   PetscCall(PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ));
293   PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ));
294   if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
295     (F)->ops->solve = MatSolve_PaStiX;
296 
297     /* Initialize a PASTIX instance */
298     PetscCallMPI(MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(lu->pastix_comm)));
299     PetscCallMPI(MPI_Comm_rank(lu->pastix_comm, &lu->commRank));
300     PetscCallMPI(MPI_Comm_size(lu->pastix_comm, &lu->commSize));
301 
302     /* Set pastix options */
303     lu->iparm[IPARM_MODIFY_PARAMETER] = API_NO;
304     lu->iparm[IPARM_START_TASK]       = API_TASK_INIT;
305     lu->iparm[IPARM_END_TASK]         = API_TASK_INIT;
306 
307     lu->rhsnbr = 1;
308 
309     /* Call to set default pastix options */
310     PASTIX_CALL(&(lu->pastix_data),
311                 lu->pastix_comm,
312                 lu->n,
313                 lu->colptr,
314                 lu->row,
315                 (PastixScalar*)lu->val,
316                 lu->perm,
317                 lu->invp,
318                 (PastixScalar*)lu->rhs,
319                 lu->rhsnbr,
320                 lu->iparm,
321                 lu->dparm);
322     PetscCheck(lu->iparm[IPARM_ERROR_NUMBER] == 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by PaStiX in MatFactorNumeric: iparm(IPARM_ERROR_NUMBER)=%" PetscInt_FMT,lu->iparm[IPARM_ERROR_NUMBER]);
323 
324     PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"PaStiX Options","Mat");
325     icntl = -1;
326     lu->iparm[IPARM_VERBOSE] = API_VERBOSE_NOT;
327     PetscCall(PetscOptionsInt("-mat_pastix_verbose","iparm[IPARM_VERBOSE] : level of printing (0 to 2)","None",lu->iparm[IPARM_VERBOSE],&icntl,&flg));
328     if ((flg && icntl >= 0) || PetscLogPrintInfo) {
329       lu->iparm[IPARM_VERBOSE] =  icntl;
330     }
331     icntl=-1;
332     PetscCall(PetscOptionsInt("-mat_pastix_threadnbr","iparm[IPARM_THREAD_NBR] : Number of thread by MPI node","None",lu->iparm[IPARM_THREAD_NBR],&icntl,&flg));
333     if ((flg && icntl > 0)) {
334       lu->iparm[IPARM_THREAD_NBR] = icntl;
335     }
336     PetscOptionsEnd();
337     valOnly = PETSC_FALSE;
338   } else {
339     if (isSeqAIJ || isMPIAIJ) {
340       PetscCall(PetscFree(lu->colptr));
341       PetscCall(PetscFree(lu->row));
342       PetscCall(PetscFree(lu->val));
343       valOnly = PETSC_FALSE;
344     } else valOnly = PETSC_TRUE;
345   }
346 
347   lu->iparm[IPARM_MATRIX_VERIFICATION] = API_YES;
348 
349   /* convert mpi A to seq mat A */
350   PetscCall(ISCreateStride(PETSC_COMM_SELF,M,0,1,&isrow));
351   PetscCall(MatCreateSubMatrices(A,1,&isrow,&isrow,MAT_INITIAL_MATRIX,&tseq));
352   PetscCall(ISDestroy(&isrow));
353 
354   PetscCall(MatConvertToCSC(*tseq,valOnly, &lu->n, &lu->colptr, &lu->row, &lu->val));
355   PetscCall(MatIsSymmetric(*tseq,0.0,&isSym));
356   PetscCall(MatDestroyMatrices(1,&tseq));
357 
358   if (!lu->perm) {
359     PetscCall(PetscMalloc1(lu->n,&(lu->perm)));
360     PetscCall(PetscMalloc1(lu->n,&(lu->invp)));
361   }
362 
363   if (isSym) {
364     /* On symmetric matrix, LLT */
365     lu->iparm[IPARM_SYM]           = API_SYM_YES;
366     lu->iparm[IPARM_FACTORIZATION] = API_FACT_LDLT;
367   } else {
368     /* On unsymmetric matrix, LU */
369     lu->iparm[IPARM_SYM]           = API_SYM_NO;
370     lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
371   }
372 
373   /*----------------*/
374   if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
375     if (!(isSeqAIJ || isSeqSBAIJ) && !lu->b_seq) {
376       /* PaStiX only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
377       PetscCall(VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq));
378       PetscCall(ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden));
379       PetscCall(MatCreateVecs(A,NULL,&b));
380       PetscCall(VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs));
381       PetscCall(VecScatterCreate(lu->b_seq,is_iden,b,is_iden,&lu->scat_sol));
382       PetscCall(ISDestroy(&is_iden));
383       PetscCall(VecDestroy(&b));
384     }
385     lu->iparm[IPARM_START_TASK] = API_TASK_ORDERING;
386     lu->iparm[IPARM_END_TASK]   = API_TASK_NUMFACT;
387 
388     PASTIX_CALL(&(lu->pastix_data),
389                 lu->pastix_comm,
390                 lu->n,
391                 lu->colptr,
392                 lu->row,
393                 (PastixScalar*)lu->val,
394                 lu->perm,
395                 lu->invp,
396                 (PastixScalar*)lu->rhs,
397                 lu->rhsnbr,
398                 lu->iparm,
399                 lu->dparm);
400     PetscCheck(lu->iparm[IPARM_ERROR_NUMBER] == 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by PaStiX in analysis phase: iparm(IPARM_ERROR_NUMBER)=%" PetscInt_FMT,lu->iparm[IPARM_ERROR_NUMBER]);
401   } else {
402     lu->iparm[IPARM_START_TASK] = API_TASK_NUMFACT;
403     lu->iparm[IPARM_END_TASK]   = API_TASK_NUMFACT;
404     PASTIX_CALL(&(lu->pastix_data),
405                 lu->pastix_comm,
406                 lu->n,
407                 lu->colptr,
408                 lu->row,
409                 (PastixScalar*)lu->val,
410                 lu->perm,
411                 lu->invp,
412                 (PastixScalar*)lu->rhs,
413                 lu->rhsnbr,
414                 lu->iparm,
415                 lu->dparm);
416     PetscCheck(lu->iparm[IPARM_ERROR_NUMBER] == 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by PaStiX in analysis phase: iparm(IPARM_ERROR_NUMBER)=%" PetscInt_FMT,lu->iparm[IPARM_ERROR_NUMBER]);
417   }
418 
419   (F)->assembled    = PETSC_TRUE;
420   lu->matstruc      = SAME_NONZERO_PATTERN;
421   lu->CleanUpPastix = PETSC_TRUE;
422   PetscFunctionReturn(0);
423 }
424 
425 /* Note the Petsc r and c permutations are ignored */
426 PetscErrorCode MatLUFactorSymbolic_AIJPASTIX(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
427 {
428   Mat_Pastix *lu = (Mat_Pastix*)F->data;
429 
430   PetscFunctionBegin;
431   lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
432   lu->iparm[IPARM_SYM]           = API_SYM_YES;
433   lu->matstruc                   = DIFFERENT_NONZERO_PATTERN;
434   F->ops->lufactornumeric        = MatFactorNumeric_PaStiX;
435   PetscFunctionReturn(0);
436 }
437 
438 PetscErrorCode MatCholeskyFactorSymbolic_SBAIJPASTIX(Mat F,Mat A,IS r,const MatFactorInfo *info)
439 {
440   Mat_Pastix *lu = (Mat_Pastix*)(F)->data;
441 
442   PetscFunctionBegin;
443   lu->iparm[IPARM_FACTORIZATION]  = API_FACT_LLT;
444   lu->iparm[IPARM_SYM]            = API_SYM_NO;
445   lu->matstruc                    = DIFFERENT_NONZERO_PATTERN;
446   (F)->ops->choleskyfactornumeric = MatFactorNumeric_PaStiX;
447   PetscFunctionReturn(0);
448 }
449 
450 PetscErrorCode MatView_PaStiX(Mat A,PetscViewer viewer)
451 {
452   PetscBool         iascii;
453   PetscViewerFormat format;
454 
455   PetscFunctionBegin;
456   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii));
457   if (iascii) {
458     PetscCall(PetscViewerGetFormat(viewer,&format));
459     if (format == PETSC_VIEWER_ASCII_INFO) {
460       Mat_Pastix *lu=(Mat_Pastix*)A->data;
461 
462       PetscCall(PetscViewerASCIIPrintf(viewer,"PaStiX run parameters:\n"));
463       PetscCall(PetscViewerASCIIPrintf(viewer,"  Matrix type :                      %s \n",((lu->iparm[IPARM_SYM] == API_SYM_YES) ? "Symmetric" : "Unsymmetric")));
464       PetscCall(PetscViewerASCIIPrintf(viewer,"  Level of printing (0,1,2):         %" PetscInt_FMT " \n",lu->iparm[IPARM_VERBOSE]));
465       PetscCall(PetscViewerASCIIPrintf(viewer,"  Number of refinements iterations : %" PetscInt_FMT " \n",lu->iparm[IPARM_NBITER]));
466       PetscCall(PetscPrintf(PETSC_COMM_SELF,"  Error :                        %g \n",lu->dparm[DPARM_RELATIVE_ERROR]));
467     }
468   }
469   PetscFunctionReturn(0);
470 }
471 
472 /*MC
473      MATSOLVERPASTIX  - A solver package providing direct solvers (LU) for distributed
474   and sequential matrices via the external package PaStiX.
475 
476   Use ./configure --download-pastix --download-ptscotch  to have PETSc installed with PasTiX
477 
478   Use -pc_type lu -pc_factor_mat_solver_type pastix to use this direct solver
479 
480   Options Database Keys:
481 + -mat_pastix_verbose   <0,1,2>   - print level
482 - -mat_pastix_threadnbr <integer> - Set the thread number by MPI task.
483 
484   Notes:
485     This only works for matrices with symmetric nonzero structure, if you pass it a matrix with
486    nonsymmetric structure PasTiX and hence PETSc return with an error.
487 
488   Level: beginner
489 
490 .seealso: PCFactorSetMatSolverType(), MatSolverType
491 
492 M*/
493 
494 PetscErrorCode MatGetInfo_PaStiX(Mat A,MatInfoType flag,MatInfo *info)
495 {
496   Mat_Pastix *lu =(Mat_Pastix*)A->data;
497 
498   PetscFunctionBegin;
499   info->block_size        = 1.0;
500   info->nz_allocated      = lu->iparm[IPARM_NNZEROS];
501   info->nz_used           = lu->iparm[IPARM_NNZEROS];
502   info->nz_unneeded       = 0.0;
503   info->assemblies        = 0.0;
504   info->mallocs           = 0.0;
505   info->memory            = 0.0;
506   info->fill_ratio_given  = 0;
507   info->fill_ratio_needed = 0;
508   info->factor_mallocs    = 0;
509   PetscFunctionReturn(0);
510 }
511 
512 static PetscErrorCode MatFactorGetSolverType_pastix(Mat A,MatSolverType *type)
513 {
514   PetscFunctionBegin;
515   *type = MATSOLVERPASTIX;
516   PetscFunctionReturn(0);
517 }
518 
519 /*
520     The seq and mpi versions of this function are the same
521 */
522 static PetscErrorCode MatGetFactor_seqaij_pastix(Mat A,MatFactorType ftype,Mat *F)
523 {
524   Mat         B;
525   Mat_Pastix *pastix;
526 
527   PetscFunctionBegin;
528   PetscCheck(ftype == MAT_FACTOR_LU,PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix");
529   /* Create the factorization matrix */
530   PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B));
531   PetscCall(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N));
532   PetscCall(PetscStrallocpy("pastix",&((PetscObject)B)->type_name));
533   PetscCall(MatSetUp(B));
534 
535   B->trivialsymbolic       = PETSC_TRUE;
536   B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
537   B->ops->view             = MatView_PaStiX;
538   B->ops->getinfo          = MatGetInfo_PaStiX;
539 
540   PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_pastix));
541 
542   B->factortype = MAT_FACTOR_LU;
543 
544   /* set solvertype */
545   PetscCall(PetscFree(B->solvertype));
546   PetscCall(PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype));
547 
548   PetscCall(PetscNewLog(B,&pastix));
549 
550   pastix->CleanUpPastix = PETSC_FALSE;
551   pastix->scat_rhs      = NULL;
552   pastix->scat_sol      = NULL;
553   B->ops->getinfo       = MatGetInfo_External;
554   B->ops->destroy       = MatDestroy_Pastix;
555   B->data               = (void*)pastix;
556 
557   *F = B;
558   PetscFunctionReturn(0);
559 }
560 
561 static PetscErrorCode MatGetFactor_mpiaij_pastix(Mat A,MatFactorType ftype,Mat *F)
562 {
563   Mat         B;
564   Mat_Pastix *pastix;
565 
566   PetscFunctionBegin;
567   PetscCheck(ftype == MAT_FACTOR_LU,PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix");
568   /* Create the factorization matrix */
569   PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B));
570   PetscCall(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N));
571   PetscCall(PetscStrallocpy("pastix",&((PetscObject)B)->type_name));
572   PetscCall(MatSetUp(B));
573 
574   B->trivialsymbolic       = PETSC_TRUE;
575   B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
576   B->ops->view             = MatView_PaStiX;
577   B->ops->getinfo          = MatGetInfo_PaStiX;
578   PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_pastix));
579 
580   B->factortype = MAT_FACTOR_LU;
581 
582   /* set solvertype */
583   PetscCall(PetscFree(B->solvertype));
584   PetscCall(PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype));
585 
586   PetscCall(PetscNewLog(B,&pastix));
587 
588   pastix->CleanUpPastix = PETSC_FALSE;
589   pastix->scat_rhs      = NULL;
590   pastix->scat_sol      = NULL;
591   B->ops->getinfo       = MatGetInfo_External;
592   B->ops->destroy       = MatDestroy_Pastix;
593   B->data               = (void*)pastix;
594 
595   *F = B;
596   PetscFunctionReturn(0);
597 }
598 
599 static PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
600 {
601   Mat         B;
602   Mat_Pastix *pastix;
603 
604   PetscFunctionBegin;
605   PetscCheck(ftype == MAT_FACTOR_CHOLESKY,PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix");
606   /* Create the factorization matrix */
607   PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B));
608   PetscCall(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N));
609   PetscCall(PetscStrallocpy("pastix",&((PetscObject)B)->type_name));
610   PetscCall(MatSetUp(B));
611 
612   B->trivialsymbolic             = PETSC_TRUE;
613   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
614   B->ops->view                   = MatView_PaStiX;
615   B->ops->getinfo                = MatGetInfo_PaStiX;
616   PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_pastix));
617 
618   B->factortype = MAT_FACTOR_CHOLESKY;
619 
620   /* set solvertype */
621   PetscCall(PetscFree(B->solvertype));
622   PetscCall(PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype));
623 
624   PetscCall(PetscNewLog(B,&pastix));
625 
626   pastix->CleanUpPastix = PETSC_FALSE;
627   pastix->scat_rhs      = NULL;
628   pastix->scat_sol      = NULL;
629   B->ops->getinfo       = MatGetInfo_External;
630   B->ops->destroy       = MatDestroy_Pastix;
631   B->data               = (void*)pastix;
632   *F = B;
633   PetscFunctionReturn(0);
634 }
635 
636 static PetscErrorCode MatGetFactor_mpisbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
637 {
638   Mat         B;
639   Mat_Pastix *pastix;
640 
641   PetscFunctionBegin;
642   PetscCheck(ftype == MAT_FACTOR_CHOLESKY,PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix");
643 
644   /* Create the factorization matrix */
645   PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B));
646   PetscCall(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N));
647   PetscCall(PetscStrallocpy("pastix",&((PetscObject)B)->type_name));
648   PetscCall(MatSetUp(B));
649 
650   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
651   B->ops->view                   = MatView_PaStiX;
652   B->ops->getinfo                = MatGetInfo_PaStiX;
653   B->ops->destroy                = MatDestroy_Pastix;
654   PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_pastix));
655 
656   B->factortype = MAT_FACTOR_CHOLESKY;
657 
658   /* set solvertype */
659   PetscCall(PetscFree(B->solvertype));
660   PetscCall(PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype));
661 
662   PetscCall(PetscNewLog(B,&pastix));
663 
664   pastix->CleanUpPastix = PETSC_FALSE;
665   pastix->scat_rhs      = NULL;
666   pastix->scat_sol      = NULL;
667   B->data               = (void*)pastix;
668 
669   *F = B;
670   PetscFunctionReturn(0);
671 }
672 
673 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Pastix(void)
674 {
675   PetscFunctionBegin;
676   PetscCall(MatSolverTypeRegister(MATSOLVERPASTIX,MATMPIAIJ,        MAT_FACTOR_LU,MatGetFactor_mpiaij_pastix));
677   PetscCall(MatSolverTypeRegister(MATSOLVERPASTIX,MATSEQAIJ,        MAT_FACTOR_LU,MatGetFactor_seqaij_pastix));
678   PetscCall(MatSolverTypeRegister(MATSOLVERPASTIX,MATMPISBAIJ,      MAT_FACTOR_CHOLESKY,MatGetFactor_mpisbaij_pastix));
679   PetscCall(MatSolverTypeRegister(MATSOLVERPASTIX,MATSEQSBAIJ,      MAT_FACTOR_CHOLESKY,MatGetFactor_seqsbaij_pastix));
680   PetscFunctionReturn(0);
681 }
682