Actual source code: strumpack.c
petsc-3.8.4 2018-03-24
1: #include <../src/mat/impls/aij/seq/aij.h>
2: #include <../src/mat/impls/aij/mpi/mpiaij.h>
3: #include <StrumpackSparseSolver.h>
5: static PetscErrorCode MatGetDiagonal_STRUMPACK(Mat A,Vec v)
6: {
8: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Mat type: STRUMPACK factor");
9: return(0);
10: }
12: static PetscErrorCode MatDestroy_STRUMPACK(Mat A)
13: {
14: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)A->spptr;
15: PetscErrorCode ierr;
16: PetscBool flg;
19: /* Deallocate STRUMPACK storage */
20: PetscStackCall("STRUMPACK_destroy",STRUMPACK_destroy(S));
21: PetscFree(A->spptr);
22: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
23: if (flg) {
24: MatDestroy_SeqAIJ(A);
25: } else {
26: MatDestroy_MPIAIJ(A);
27: }
29: /* clear composed functions */
30: PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);
31: PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSRelCompTol_C",NULL);
32: PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSMinSize_C",NULL);
33: PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetColPerm_C",NULL);
35: return(0);
36: }
38: static PetscErrorCode MatSTRUMPACKSetHSSRelCompTol_STRUMPACK(Mat F,PetscReal rctol)
39: {
40: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
43: PetscStackCall("STRUMPACK_set_rctol", STRUMPACK_set_rctol(*S,rctol));
44: return(0);
45: }
47: /*@
48: MatSTRUMPACKSetHSSRelCompTol - Set STRUMPACK relative tolerance for HSS compression
49: Logically Collective on Mat
51: Input Parameters:
52: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
53: - rctol - relative compression tolerance
55: Options Database:
56: . -mat_strumpack_rctol <rctol>
58: Level: beginner
60: References:
61: . STRUMPACK manual
63: .seealso: MatGetFactor()
64: @*/
65: PetscErrorCode MatSTRUMPACKSetHSSRelCompTol(Mat F,PetscReal rctol)
66: {
72: PetscTryMethod(F,"MatSTRUMPACKSetHSSRelCompTol_C",(Mat,PetscReal),(F,rctol));
73: return(0);
74: }
76: static PetscErrorCode MatSTRUMPACKSetHSSMinSize_STRUMPACK(Mat F,PetscInt hssminsize)
77: {
78: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
81: PetscStackCall("STRUMPACK_set_minimum_HSS_size", STRUMPACK_set_minimum_HSS_size(*S,hssminsize));
82: return(0);
83: }
85: /*@
86: MatSTRUMPACKSetHSSMinSize - Set STRUMPACK minimum dense matrix size for low-rank approximation
87: Logically Collective on Mat
89: Input Parameters:
90: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
91: - hssminsize - minimum dense matrix size for low-rank approximation
93: Options Database:
94: . -mat_strumpack_hssminsize <hssminsize>
96: Level: beginner
98: References:
99: . STRUMPACK manual
101: .seealso: MatGetFactor()
102: @*/
103: PetscErrorCode MatSTRUMPACKSetHSSMinSize(Mat F,PetscInt hssminsize)
104: {
110: PetscTryMethod(F,"MatSTRUMPACKSetHSSMinSize_C",(Mat,PetscInt),(F,hssminsize));
111: return(0);
112: }
114: static PetscErrorCode MatSTRUMPACKSetColPerm_STRUMPACK(Mat F,PetscBool cperm)
115: {
116: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
119: PetscStackCall("STRUMPACK_set_mc64job",STRUMPACK_set_mc64job(*S,cperm ? 5 : 0));
120: return(0);
121: }
123: /*@
124: MatSTRUMPACKSetColPerm - Set whether STRUMPACK should try to permute the columns of the matrix in order to get a nonzero diagonal
125: Logically Collective on Mat
127: Input Parameters:
128: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
129: - cperm - whether or not to permute (internally) the columns of the matrix
131: Options Database:
132: . -mat_strumpack_colperm <cperm>
134: Level: beginner
136: References:
137: . STRUMPACK manual
139: .seealso: MatGetFactor()
140: @*/
141: PetscErrorCode MatSTRUMPACKSetColPerm(Mat F,PetscBool cperm)
142: {
148: PetscTryMethod(F,"MatSTRUMPACKSetColPerm_C",(Mat,PetscBool),(F,cperm));
149: return(0);
150: }
152: static PetscErrorCode MatSolve_STRUMPACK(Mat A,Vec b_mpi,Vec x)
153: {
154: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)A->spptr;
155: STRUMPACK_RETURN_CODE sp_err;
156: PetscErrorCode ierr;
157: const PetscScalar *bptr;
158: PetscScalar *xptr;
161: VecGetArray(x,&xptr);
162: VecGetArrayRead(b_mpi,&bptr);
164: PetscStackCall("STRUMPACK_solve",sp_err = STRUMPACK_solve(*S,(PetscScalar*)bptr,xptr,0));
165: switch (sp_err) {
166: case STRUMPACK_SUCCESS: break;
167: case STRUMPACK_MATRIX_NOT_SET: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix was not set"); break; }
168: case STRUMPACK_REORDERING_ERROR: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix reordering failed"); break; }
169: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: solve failed");
170: }
171: VecRestoreArray(x,&xptr);
172: VecRestoreArrayRead(b_mpi,&bptr);
173: return(0);
174: }
176: static PetscErrorCode MatMatSolve_STRUMPACK(Mat A,Mat B_mpi,Mat X)
177: {
178: PetscErrorCode ierr;
179: PetscBool flg;
182: PetscObjectTypeCompareAny((PetscObject)B_mpi,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
183: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
184: PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
185: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
186: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_STRUMPACK() is not implemented yet");
187: return(0);
188: }
190: static PetscErrorCode MatFactorInfo_STRUMPACK(Mat A,PetscViewer viewer)
191: {
192: PetscErrorCode ierr;
195: /* check if matrix is strumpack type */
196: if (A->ops->solve != MatSolve_STRUMPACK) return(0);
197: PetscViewerASCIIPrintf(viewer,"STRUMPACK sparse solver!\n");
198: return(0);
199: }
201: static PetscErrorCode MatView_STRUMPACK(Mat A,PetscViewer viewer)
202: {
203: PetscErrorCode ierr;
204: PetscBool iascii;
205: PetscViewerFormat format;
208: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
209: if (iascii) {
210: PetscViewerGetFormat(viewer,&format);
211: if (format == PETSC_VIEWER_ASCII_INFO) {
212: MatFactorInfo_STRUMPACK(A,viewer);
213: }
214: }
215: return(0);
216: }
218: static PetscErrorCode MatLUFactorNumeric_STRUMPACK(Mat F,Mat A,const MatFactorInfo *info)
219: {
220: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
221: STRUMPACK_RETURN_CODE sp_err;
222: Mat_SeqAIJ *A_d,*A_o;
223: Mat_MPIAIJ *mat;
224: PetscErrorCode ierr;
225: PetscInt M=A->rmap->N,m=A->rmap->n;
226: PetscBool flg;
229: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);
230: if (flg) { /* A is MATMPIAIJ */
231: mat = (Mat_MPIAIJ*)A->data;
232: A_d = (Mat_SeqAIJ*)(mat->A)->data;
233: A_o = (Mat_SeqAIJ*)(mat->B)->data;
234: PetscStackCall("STRUMPACK_set_MPIAIJ_matrix",STRUMPACK_set_MPIAIJ_matrix(*S,&m,A_d->i,A_d->j,A_d->a,A_o->i,A_o->j,A_o->a,mat->garray));
235: } else { /* A is MATSEQAIJ */
236: A_d = (Mat_SeqAIJ*)A->data;
237: PetscStackCall("STRUMPACK_set_csr_matrix",STRUMPACK_set_csr_matrix(*S,&M,A_d->i,A_d->j,A_d->a,0));
238: }
240: /* Reorder and Factor the matrix. */
241: /* TODO figure out how to avoid reorder if the matrix values changed, but the pattern remains the same. */
242: PetscStackCall("STRUMPACK_reorder",sp_err = STRUMPACK_reorder(*S));
243: PetscStackCall("STRUMPACK_factor",sp_err = STRUMPACK_factor(*S));
244: switch (sp_err) {
245: case STRUMPACK_SUCCESS: break;
246: case STRUMPACK_MATRIX_NOT_SET: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix was not set"); break; }
247: case STRUMPACK_REORDERING_ERROR: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix reordering failed"); break; }
248: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: factorization failed");
249: }
250: return(0);
251: }
253: static PetscErrorCode MatLUFactorSymbolic_STRUMPACK(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
254: {
256: F->ops->lufactornumeric = MatLUFactorNumeric_STRUMPACK;
257: F->ops->solve = MatSolve_STRUMPACK;
258: F->ops->matsolve = MatMatSolve_STRUMPACK;
259: return(0);
260: }
262: static PetscErrorCode MatFactorGetSolverPackage_aij_strumpack(Mat A,const MatSolverPackage *type)
263: {
265: *type = MATSOLVERSTRUMPACK;
266: return(0);
267: }
269: /*MC
270: MATSOLVERSSTRUMPACK = "strumpack" - A solver package providing a direct sparse solver (PCLU)
271: and a preconditioner (PCILU) using low-rank compression via the external package STRUMPACK.
273: Consult the STRUMPACK-sparse manual for more info.
275: Use
276: ./configure --download-strumpack
277: to have PETSc installed with STRUMPACK
279: Use
280: -pc_type lu -pc_factor_mat_solver_package strumpack
281: to use this as an exact (direct) solver, use
282: -pc_type ilu -pc_factor_mat_solver_package strumpack
283: to enable low-rank compression (i.e, use as a preconditioner).
285: Works with AIJ matrices
287: Options Database Keys:
288: + -mat_strumpack_rctol <1e-4> - Relative compression tolerance (None)
289: . -mat_strumpack_hssminsize <512> - Minimum size of dense block for HSS compression (None)
290: . -mat_strumpack_colperm <TRUE> - Permute matrix to make diagonal nonzeros (None)
292: Level: beginner
294: .seealso: PCLU, PCILU, MATSOLVERSUPERLU_DIST, MATSOLVERMUMPS, PCFactorSetMatSolverPackage(), MatSolverPackage
295: M*/
296: static PetscErrorCode MatGetFactor_aij_strumpack(Mat A,MatFactorType ftype,Mat *F)
297: {
298: Mat B;
299: STRUMPACK_SparseSolver *S;
300: PetscErrorCode ierr;
301: PetscInt M=A->rmap->N,N=A->cmap->N;
302: STRUMPACK_INTERFACE iface;
303: PetscBool verb,flg,set;
304: PetscReal rctol;
305: PetscInt hssminsize;
306: int argc;
307: char **args,*copts,*pname;
308: size_t len;
309: const STRUMPACK_PRECISION table[2][2][2] = {{{STRUMPACK_FLOATCOMPLEX_64,STRUMPACK_DOUBLECOMPLEX_64},
310: {STRUMPACK_FLOAT_64,STRUMPACK_DOUBLE_64}},
311: {{STRUMPACK_FLOATCOMPLEX,STRUMPACK_DOUBLECOMPLEX},
312: {STRUMPACK_FLOAT,STRUMPACK_DOUBLE}}};
313: const STRUMPACK_PRECISION prec = table[(sizeof(PetscInt)==8)?0:1][(PETSC_SCALAR==PETSC_COMPLEX)?0:1][(PETSC_REAL==PETSC_FLOAT)?0:1];
316: /* Create the factorization matrix */
317: MatCreate(PetscObjectComm((PetscObject)A),&B);
318: MatSetSizes(B,A->rmap->n,A->cmap->n,M,N);
319: MatSetType(B,((PetscObject)A)->type_name);
320: MatSeqAIJSetPreallocation(B,0,NULL);
321: MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
322: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
323: B->ops->lufactorsymbolic = MatLUFactorSymbolic_STRUMPACK;
324: B->ops->ilufactorsymbolic = MatLUFactorSymbolic_STRUMPACK;
325: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
326: B->ops->view = MatView_STRUMPACK;
327: B->ops->destroy = MatDestroy_STRUMPACK;
328: B->ops->getdiagonal = MatGetDiagonal_STRUMPACK;
329: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_aij_strumpack);
330: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSRelCompTol_C",MatSTRUMPACKSetHSSRelCompTol_STRUMPACK);
331: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSMinSize_C",MatSTRUMPACKSetHSSMinSize_STRUMPACK);
332: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetColPerm_C",MatSTRUMPACKSetColPerm_STRUMPACK);
333: B->factortype = ftype;
334: B->assembled = PETSC_TRUE; /* required by -ksp_view */
335: PetscNewLog(B,&S);
336: B->spptr = S;
338: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
339: iface = flg ? STRUMPACK_MT : STRUMPACK_MPI_DIST;
341: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"STRUMPACK Options","Mat");
343: verb = PetscLogPrintInfo ? PETSC_TRUE : PETSC_FALSE;
344: PetscOptionsBool("-mat_strumpack_verbose","Print STRUMPACK information","None",verb,&verb,NULL);
346: PetscOptionsGetAll(NULL,&copts);
347: PetscStrlen(copts,&len);
348: len += PETSC_MAX_PATH_LEN+1;
349: PetscMalloc1(len,&pname);
350: /* first string is assumed to be the program name, so add program name to options string */
351: PetscGetProgramName(pname,len);
352: PetscStrcat(pname," ");
353: PetscStrcat(pname,copts);
354: PetscStrToArray(pname,' ',&argc,&args);
355: PetscFree(copts);
356: PetscFree(pname);
358: PetscStackCall("STRUMPACK_init",STRUMPACK_init(S,PetscObjectComm((PetscObject)A),prec,iface,argc,args,verb));
360: PetscStackCall("STRUMPACK_get_rctol",rctol = (PetscReal)STRUMPACK_get_rctol(*S));
361: PetscOptionsReal("-mat_strumpack_rctol","Relative compression tolerance","None",rctol,&rctol,&set);
362: if (set) PetscStackCall("STRUMPACK_set_rctol",STRUMPACK_set_rctol(*S,(double)rctol));
364: PetscStackCall("STRUMPACK_get_mc64job",flg = (STRUMPACK_get_mc64job(*S) == 0) ? PETSC_FALSE : PETSC_TRUE);
365: PetscOptionsBool("-mat_strumpack_colperm","Find a col perm to get nonzero diagonal","None",flg,&flg,&set);
366: if (set) PetscStackCall("STRUMPACK_set_mc64job",STRUMPACK_set_mc64job(*S,flg ? 5 : 0));
368: PetscStackCall("STRUMPACK_get_minimum_HSS_size",hssminsize = (PetscInt)STRUMPACK_get_minimum_HSS_size(*S));
369: PetscOptionsInt("-mat_strumpack_hssminsize","Minimum size of dense block for HSS compression","None",hssminsize,&hssminsize,&set);
370: if (set) PetscStackCall("STRUMPACK_set_minimum_HSS_size",STRUMPACK_set_minimum_HSS_size(*S,(int)hssminsize));
372: PetscOptionsEnd();
374: if (ftype == MAT_FACTOR_ILU) {
375: /* When enabling HSS compression, the STRUMPACK solver becomes an incomplete */
376: /* (or approximate) LU factorization. */
377: PetscStackCall("STRUMPACK_use_HSS",STRUMPACK_use_HSS(*S,1));
378: /* Disable the outer iterative solver from STRUMPACK. */
379: /* When STRUMPACK is used as a direct solver, it will by default do iterative refinement. */
380: /* When STRUMPACK is used with as an approximate factorization preconditioner (by enabling */
381: /* low-rank compression), it will use it's own GMRES. Here we can disable the */
382: /* outer iterative solver, as PETSc uses STRUMPACK from within a KSP. */
383: PetscStackCall("STRUMPACK_set_Krylov_solver", STRUMPACK_set_Krylov_solver(*S, STRUMPACK_DIRECT));
384: }
386: /* Allow the user to set or overwrite the above options from the command line */
387: PetscStackCall("STRUMPACK_set_from_options",STRUMPACK_set_from_options(*S));
388: PetscStrToArrayDestroy(argc,args);
390: *F = B;
391: return(0);
392: }
394: PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_STRUMPACK(void)
395: {
399: MatSolverPackageRegister(MATSOLVERSTRUMPACK,MATMPIAIJ,MAT_FACTOR_LU,MatGetFactor_aij_strumpack);
400: MatSolverPackageRegister(MATSOLVERSTRUMPACK,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_aij_strumpack);
401: MatSolverPackageRegister(MATSOLVERSTRUMPACK,MATMPIAIJ,MAT_FACTOR_ILU,MatGetFactor_aij_strumpack);
402: MatSolverPackageRegister(MATSOLVERSTRUMPACK,MATSEQAIJ,MAT_FACTOR_ILU,MatGetFactor_aij_strumpack);
403: return(0);
404: }