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