xref: /petsc/src/mat/impls/aij/mpi/strumpack/strumpack.c (revision bf88b097425cd9825fae466b4b414b1007061480) 
1 #include <../src/mat/impls/aij/seq/aij.h>
2 #include <../src/mat/impls/aij/mpi/mpiaij.h>
3 #include <StrumpackSparseSolver.h>
4 
5 /*
6     CENTRALIZED - STRUMPACK expects the entire sparse matrix and right-hand side on every process.
7     DISTRIBUTED - STRUMPACK expects the sparse matrix and right-hand side to be distributed across the entire MPI communicator.
8 */
9 typedef enum {CENTRALIZED, DISTRIBUTED} STRUMPACK_MatInputMode;
10 const char *STRUMPACK_MatInputModes[] = {"CENTRALIZED","DISTRIBUTED","STRUMPACK_MatInputMode","PETSC_",0};
11 
12 typedef struct {
13   STRUMPACK_SparseSolver S;
14   STRUMPACK_MatInputMode MatInputMode;
15 } STRUMPACK_data;
16 
17 extern PetscErrorCode MatFactorInfo_STRUMPACK_MPI(Mat,PetscViewer);
18 extern PetscErrorCode MatLUFactorNumeric_STRUMPACK_MPI(Mat,Mat,const MatFactorInfo*);
19 extern PetscErrorCode MatDestroy_STRUMPACK_MPI(Mat);
20 extern PetscErrorCode MatView_STRUMPACK_MPI(Mat,PetscViewer);
21 extern PetscErrorCode MatSolve_STRUMPACK_MPI(Mat,Vec,Vec);
22 extern PetscErrorCode MatLUFactorSymbolic_STRUMPACK_MPI(Mat,Mat,IS,IS,const MatFactorInfo*);
23 extern PetscErrorCode MatDestroy_MPIAIJ(Mat);
24 
25 #undef __FUNCT__
26 #define __FUNCT__ "MatGetDiagonal_STRUMPACK_MPI"
27 PetscErrorCode MatGetDiagonal_STRUMPACK_MPI(Mat A,Vec v)
28 {
29   PetscFunctionBegin;
30   SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Mat type: STRUMPACK_MPI factor");
31   PetscFunctionReturn(0);
32 }
33 
34 #undef __FUNCT__
35 #define __FUNCT__ "MatDestroy_STRUMPACK_MPI"
36 PetscErrorCode MatDestroy_STRUMPACK_MPI(Mat A)
37 {
38   STRUMPACK_data *sp = (STRUMPACK_data*)A->spptr;
39   PetscErrorCode ierr;
40   PetscBool      flg;
41 
42   PetscFunctionBegin;
43   /* Deallocate STRUMPACK_MPI storage */
44   PetscStackCall("STRUMPACK_destroy",STRUMPACK_destroy(&(sp->S)));
45   ierr = PetscFree(A->spptr);CHKERRQ(ierr);
46   ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr);
47   if (flg) {
48     ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr);
49   } else {
50     ierr = MatDestroy_MPIAIJ(A);CHKERRQ(ierr);
51   }
52   PetscFunctionReturn(0);
53 }
54 
55 #undef __FUNCT__
56 #define __FUNCT__ "MatSolve_STRUMPACK_MPI"
57 PetscErrorCode MatSolve_STRUMPACK_MPI(Mat A,Vec b_mpi,Vec x)
58 {
59   STRUMPACK_data        *sp = (STRUMPACK_data*)A->spptr;
60   STRUMPACK_RETURN_CODE sp_err;
61   PetscErrorCode        ierr;
62   PetscMPIInt           size;
63   PetscInt              N=A->cmap->N;
64   const PetscScalar     *bptr;
65   PetscScalar           *xptr;
66   Vec                   x_seq,b_seq;
67   IS                    iden;
68   VecScatter            scat;
69 
70   PetscFunctionBegin;
71   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr);
72   if (size > 1 && sp->MatInputMode == CENTRALIZED) {
73     ierr = VecCreateSeq(PETSC_COMM_SELF,N,&x_seq);CHKERRQ(ierr);
74     ierr = VecGetArray(x_seq,&xptr);CHKERRQ(ierr);
75     /* global mat input, convert b to b_seq */
76     ierr = VecCreateSeq(PETSC_COMM_SELF,N,&b_seq);CHKERRQ(ierr);
77     ierr = ISCreateStride(PETSC_COMM_SELF,N,0,1,&iden);CHKERRQ(ierr);
78     ierr = VecScatterCreate(b_mpi,iden,b_seq,iden,&scat);CHKERRQ(ierr);
79     ierr = ISDestroy(&iden);CHKERRQ(ierr);
80     ierr = VecScatterBegin(scat,b_mpi,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
81     ierr = VecScatterEnd(scat,b_mpi,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
82     ierr = VecGetArrayRead(b_seq,&bptr);CHKERRQ(ierr);
83   } else { /* size==1 || distributed mat input */
84     ierr = VecGetArray(x,&xptr);CHKERRQ(ierr);
85     ierr = VecGetArrayRead(b_mpi,&bptr);CHKERRQ(ierr);
86   }
87 
88   PetscStackCall("STRUMPACK_solve",sp_err = STRUMPACK_solve(sp->S,(PetscScalar*)bptr,xptr,0));
89 
90   if (sp_err != STRUMPACK_SUCCESS) {
91     if (sp_err == STRUMPACK_MATRIX_NOT_SET)        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix was not set");
92     else if (sp_err == STRUMPACK_REORDERING_ERROR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix reordering failed");
93     else                                           SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: solve failed");
94   }
95 
96   if (size > 1 && sp->MatInputMode == CENTRALIZED) {
97     ierr = VecRestoreArrayRead(b_seq,&bptr);CHKERRQ(ierr);
98     ierr = VecDestroy(&b_seq);CHKERRQ(ierr);
99     /* convert seq x to mpi x */
100     ierr = VecRestoreArray(x_seq,&xptr);CHKERRQ(ierr);
101     ierr = VecScatterBegin(scat,x_seq,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
102     ierr = VecScatterEnd(scat,x_seq,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
103     ierr = VecScatterDestroy(&scat);CHKERRQ(ierr);
104     ierr = VecDestroy(&x_seq);CHKERRQ(ierr);
105   } else {
106     ierr = VecRestoreArray(x,&xptr);CHKERRQ(ierr);
107     ierr = VecRestoreArrayRead(b_mpi,&bptr);CHKERRQ(ierr);
108   }
109 
110   PetscFunctionReturn(0);
111 }
112 
113 #undef __FUNCT__
114 #define __FUNCT__ "MatMatSolve_STRUMPACK_MPI"
115 PetscErrorCode MatMatSolve_STRUMPACK_MPI(Mat A,Mat B_mpi,Mat X)
116 {
117   PetscErrorCode   ierr;
118   PetscBool        flg;
119 
120   PetscFunctionBegin;
121   ierr = PetscObjectTypeCompareAny((PetscObject)B_mpi,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr);
122   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
123   ierr = PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr);
124   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
125   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_STRUMPACK_MPI() is not implemented yet");
126   PetscFunctionReturn(0);
127 }
128 
129 
130 #undef __FUNCT__
131 #define __FUNCT__ "MatLUFactorNumeric_STRUMPACK_MPI"
132 PetscErrorCode MatLUFactorNumeric_STRUMPACK_MPI(Mat F,Mat A,const MatFactorInfo *info)
133 {
134   STRUMPACK_data        *sp = (STRUMPACK_data*)F->spptr;
135   STRUMPACK_RETURN_CODE sp_err;
136   Mat                   *tseq,A_seq = NULL;
137   Mat_SeqAIJ            *A_d,*A_o;
138   PetscErrorCode        ierr;
139   PetscInt              M=A->rmap->N,m=A->rmap->n,N=A->cmap->N;
140   PetscMPIInt           size;
141   IS                    isrow;
142 
143   PetscFunctionBegin;
144   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr);
145 
146   if (sp->MatInputMode == CENTRALIZED) { /* global mat input */
147     if (size > 1) { /* convert mpi A to seq mat A */
148       ierr = ISCreateStride(PETSC_COMM_SELF,M,0,1,&isrow);CHKERRQ(ierr);
149       ierr = MatGetSubMatrices(A,1,&isrow,&isrow,MAT_INITIAL_MATRIX,&tseq);CHKERRQ(ierr);
150       ierr = ISDestroy(&isrow);CHKERRQ(ierr);
151       A_seq = *tseq;
152       ierr  = PetscFree(tseq);CHKERRQ(ierr);
153       A_d   = (Mat_SeqAIJ*)A_seq->data;
154     } else {
155       PetscBool flg;
156       ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);CHKERRQ(ierr);
157       if (flg) {
158         Mat_MPIAIJ *At = (Mat_MPIAIJ*)A->data;
159         A = At->A;
160       }
161       A_d = (Mat_SeqAIJ*)A->data;
162     }
163     PetscStackCall("STRUMPACK_set_csr_matrix",STRUMPACK_set_csr_matrix(sp->S,&N,A_d->i,A_d->j,A_d->a,0));
164   } else { /* distributed mat input */
165     Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data;
166     A_d = (Mat_SeqAIJ*)(mat->A)->data;
167     A_o = (Mat_SeqAIJ*)(mat->B)->data;
168     PetscStackCall("STRUMPACK_set_MPIAIJ_matrix",STRUMPACK_set_MPIAIJ_matrix(sp->S,&m,A_d->i,A_d->j,A_d->a,A_o->i,A_o->j,A_o->a,mat->garray));
169   }
170 
171   /* Reorder and Factor the matrix. */
172   /* TODO figure out how to avoid reorder if the matrix values changed, but the pattern remains the same. */
173   PetscStackCall("STRUMPACK_reorder",sp_err = STRUMPACK_reorder(sp->S));
174   PetscStackCall("STRUMPACK_factor",sp_err = STRUMPACK_factor(sp->S));
175   if (sp_err != STRUMPACK_SUCCESS) {
176     if (sp_err == STRUMPACK_MATRIX_NOT_SET)        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix was not set");
177     else if (sp_err == STRUMPACK_REORDERING_ERROR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix reordering failed");
178     else                                           SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: factorization failed");
179   }
180 
181   if (sp->MatInputMode == CENTRALIZED && size > 1) {
182     ierr = MatDestroy(&A_seq);CHKERRQ(ierr);
183   }
184 
185   PetscFunctionReturn(0);
186 }
187 
188 #undef __FUNCT__
189 #define __FUNCT__ "MatLUFactorSymbolic_STRUMPACK_MPI"
190 PetscErrorCode MatLUFactorSymbolic_STRUMPACK_MPI(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
191 {
192   PetscFunctionBegin;
193   F->ops->lufactornumeric = MatLUFactorNumeric_STRUMPACK_MPI;
194   F->ops->solve           = MatSolve_STRUMPACK_MPI;
195   F->ops->matsolve        = MatMatSolve_STRUMPACK_MPI;
196   PetscFunctionReturn(0);
197 }
198 
199 #undef __FUNCT__
200 #define __FUNCT__ "MatFactorGetSolverPackage_aij_strumpack_mpi"
201 PetscErrorCode MatFactorGetSolverPackage_aij_strumpack_mpi(Mat A,const MatSolverPackage *type)
202 {
203   PetscFunctionBegin;
204   *type = MATSOLVERSTRUMPACK_MPI;
205   PetscFunctionReturn(0);
206 }
207 
208 #undef __FUNCT__
209 #define __FUNCT__ "MatGetFactor_aij_strumpack_mpi"
210 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_strumpack_mpi(Mat A,MatFactorType ftype,Mat *F)
211 {
212   Mat                 B;
213   STRUMPACK_data      *sp;
214   PetscErrorCode      ierr;
215   PetscInt            M=A->rmap->N,N=A->cmap->N;
216   PetscMPIInt         size;
217   STRUMPACK_INTERFACE iface;
218   PetscBool           verb;
219   PetscInt            argc;
220   char                **args;
221 
222   PetscFunctionBegin;
223   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr);
224 
225   /* Create the factorization matrix */
226   ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr);
227   ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,M,N);CHKERRQ(ierr);
228   ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
229   ierr = MatSeqAIJSetPreallocation(B,0,NULL);
230   ierr = MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);CHKERRQ(ierr);
231 
232   B->ops->lufactorsymbolic = MatLUFactorSymbolic_STRUMPACK_MPI;
233   B->ops->view             = MatView_STRUMPACK_MPI;
234   B->ops->destroy          = MatDestroy_STRUMPACK_MPI;
235   B->ops->getdiagonal      = MatGetDiagonal_STRUMPACK_MPI;
236 
237   ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_aij_strumpack_mpi);CHKERRQ(ierr);
238 
239   B->factortype = MAT_FACTOR_LU;
240 
241   ierr     = PetscNewLog(B,&sp);CHKERRQ(ierr);
242   B->spptr = sp;
243 
244   ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"STRUMPACK_MPI Options","Mat");CHKERRQ(ierr);
245 
246   sp->MatInputMode = DISTRIBUTED;
247   iface = STRUMPACK_MPI_DIST;
248 
249   ierr = PetscOptionsEnum("-mat_strumpack_mpi_matinput","Matrix input mode (global or distributed)","None",STRUMPACK_MatInputModes,
250                           (PetscEnum)sp->MatInputMode,(PetscEnum*)&sp->MatInputMode,NULL);CHKERRQ(ierr);
251   if (sp->MatInputMode == DISTRIBUTED && size == 1) sp->MatInputMode = CENTRALIZED;
252   if (sp->MatInputMode == DISTRIBUTED) {
253     iface = STRUMPACK_MPI_DIST;
254   } else if (sp->MatInputMode == CENTRALIZED) {
255     iface = STRUMPACK_MPI;
256   }
257   if (PetscLogPrintInfo) verb = PETSC_TRUE;
258   else verb = PETSC_FALSE;
259   ierr = PetscOptionsBool("-mat_strumpack_verbose","Print STRUMPACK information","None",verb,&verb,NULL);CHKERRQ(ierr);
260 
261   PetscOptionsEnd();
262 
263   ierr = PetscGetArgs(&argc,&args);CHKERRQ(ierr);
264 
265 #if defined(PETSC_USE_64BIT_INDICES)
266 #if defined(PETSC_USE_COMPLEX)
267 #if defined(PETSC_USE_REAL_SINGLE)
268   PetscStackCall("STRUMPACK_init",STRUMPACK_init(&(sp->S),PetscObjectComm((PetscObject)A),STRUMPACK_FLOATCOMPLEX_64,iface,argc,args,verb));
269 #else
270   PetscStackCall("STRUMPACK_init",STRUMPACK_init(&(sp->S),PetscObjectComm((PetscObject)A),STRUMPACK_DOUBLECOMPLEX_64,iface,argc,args,verb));
271 #endif
272 #else
273 #if defined(PETSC_USE_REAL_SINGLE)
274   PetscStackCall("STRUMPACK_init",STRUMPACK_init(&(sp->S),PetscObjectComm((PetscObject)A),STRUMPACK_FLOAT_64,iface,argc,args,verb));
275 #else
276   PetscStackCall("STRUMPACK_init",STRUMPACK_init(&(sp->S),PetscObjectComm((PetscObject)A),STRUMPACK_DOUBLE_64,iface,argc,args,verb));
277 #endif
278 #endif
279 #else
280 #if defined(PETSC_USE_COMPLEX)
281 #if defined(PETSC_USE_REAL_SINGLE)
282   PetscStackCall("STRUMPACK_init",STRUMPACK_init(&(sp->S),PetscObjectComm((PetscObject)A),STRUMPACK_FLOATCOMPLEX,iface,argc,args,verb));
283 #else
284   PetscStackCall("STRUMPACK_init",STRUMPACK_init(&(sp->S),PetscObjectComm((PetscObject)A),STRUMPACK_DOUBLECOMPLEX,iface,argc,args,verb));
285 #endif
286 #else
287 #if defined(PETSC_USE_REAL_SINGLE)
288   PetscStackCall("STRUMPACK_init",STRUMPACK_init(&(sp->S),PetscObjectComm((PetscObject)A),STRUMPACK_FLOAT,iface,argc,args,verb));
289 #else
290   PetscStackCall("STRUMPACK_init",STRUMPACK_init(&(sp->S),PetscObjectComm((PetscObject)A),STRUMPACK_DOUBLE,iface,argc,args,verb));
291 #endif
292 #endif
293 #endif
294   PetscStackCall("STRUMPACK_set_from_options",STRUMPACK_set_from_options(sp->S));
295 
296   *F = B;
297   PetscFunctionReturn(0);
298 }
299 
300 #undef __FUNCT__
301 #define __FUNCT__ "MatSolverPackageRegister_STRUMPACK_MPI"
302 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_STRUMPACK_MPI(void)
303 {
304   PetscErrorCode ierr;
305 
306   PetscFunctionBegin;
307   ierr = MatSolverPackageRegister(MATSOLVERSTRUMPACK_MPI,MATMPIAIJ,MAT_FACTOR_LU,MatGetFactor_aij_strumpack_mpi);CHKERRQ(ierr);
308   ierr = MatSolverPackageRegister(MATSOLVERSTRUMPACK_MPI,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_aij_strumpack_mpi);CHKERRQ(ierr);
309   PetscFunctionReturn(0);
310 }
311 
312 #undef __FUNCT__
313 #define __FUNCT__ "MatFactorInfo_STRUMPACK_MPI"
314 PetscErrorCode MatFactorInfo_STRUMPACK_MPI(Mat A,PetscViewer viewer)
315 {
316   PetscErrorCode  ierr;
317 
318   PetscFunctionBegin;
319   /* check if matrix is strumpack type */
320   if (A->ops->solve != MatSolve_STRUMPACK_MPI) PetscFunctionReturn(0);
321   ierr = PetscViewerASCIIPrintf(viewer,"STRUMPACK_MPI sparse solver!\n");CHKERRQ(ierr);
322   PetscFunctionReturn(0);
323 }
324 
325 #undef __FUNCT__
326 #define __FUNCT__ "MatView_STRUMPACK_MPI"
327 PetscErrorCode MatView_STRUMPACK_MPI(Mat A,PetscViewer viewer)
328 {
329   PetscErrorCode    ierr;
330   PetscBool         iascii;
331   PetscViewerFormat format;
332 
333   PetscFunctionBegin;
334   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
335   if (iascii) {
336     ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
337     if (format == PETSC_VIEWER_ASCII_INFO) {
338       ierr = MatFactorInfo_STRUMPACK_MPI(A,viewer);CHKERRQ(ierr);
339     }
340   }
341   PetscFunctionReturn(0);
342 }
343 
344 
345 /*MC
346   MATSOLVERSSTRUMPACK_MPI - Parallel direct solver package for LU factorization
347 
348   Use ./configure --download-strumpack --download-parmetis --download-metis --download-ptscotch (and ScaLAPACK?/BLACS?) to have PETSc installed with SuperLU_DIST
349 
350   Use -pc_type lu -pc_factor_mat_solver_package strumpack to us this direct solver
351 
352    Works with AIJ matrices
353 
354 .seealso: PCLU
355 
356 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage
357 
358 M*/
359 
360