Actual source code: superlu.c
2: /* --------------------------------------------------------------------
4: This file implements a subclass of the SeqAIJ matrix class that uses
5: the SuperLU sparse solver.
6: */
8: /*
9: Defines the data structure for the base matrix type (SeqAIJ)
10: */
11: #include <../src/mat/impls/aij/seq/aij.h>
13: /*
14: SuperLU include files
15: */
16: EXTERN_C_BEGIN
17: #if defined(PETSC_USE_COMPLEX)
18: #if defined(PETSC_USE_REAL_SINGLE)
19: #include <slu_cdefs.h>
20: #else
21: #include <slu_zdefs.h>
22: #endif
23: #else
24: #if defined(PETSC_USE_REAL_SINGLE)
25: #include <slu_sdefs.h>
26: #else
27: #include <slu_ddefs.h>
28: #endif
29: #endif
30: #include <slu_util.h>
31: EXTERN_C_END
33: /*
34: This is the data that defines the SuperLU factored matrix type
35: */
36: typedef struct {
37: SuperMatrix A,L,U,B,X;
38: superlu_options_t options;
39: PetscInt *perm_c; /* column permutation vector */
40: PetscInt *perm_r; /* row permutations from partial pivoting */
41: PetscInt *etree;
42: PetscReal *R, *C;
43: char equed[1];
44: PetscInt lwork;
45: void *work;
46: PetscReal rpg, rcond;
47: mem_usage_t mem_usage;
48: MatStructure flg;
49: SuperLUStat_t stat;
50: Mat A_dup;
51: PetscScalar *rhs_dup;
52: GlobalLU_t Glu;
53: PetscBool needconversion;
55: /* Flag to clean up (non-global) SuperLU objects during Destroy */
56: PetscBool CleanUpSuperLU;
57: } Mat_SuperLU;
59: /*
60: Utility function
61: */
62: static PetscErrorCode MatView_Info_SuperLU(Mat A,PetscViewer viewer)
63: {
64: Mat_SuperLU *lu= (Mat_SuperLU*)A->data;
65: superlu_options_t options;
67: options = lu->options;
69: PetscViewerASCIIPrintf(viewer,"SuperLU run parameters:\n");
70: PetscViewerASCIIPrintf(viewer," Equil: %s\n",(options.Equil != NO) ? "YES" : "NO");
71: PetscViewerASCIIPrintf(viewer," ColPerm: %" PetscInt_FMT "\n",options.ColPerm);
72: PetscViewerASCIIPrintf(viewer," IterRefine: %" PetscInt_FMT "\n",options.IterRefine);
73: PetscViewerASCIIPrintf(viewer," SymmetricMode: %s\n",(options.SymmetricMode != NO) ? "YES" : "NO");
74: PetscViewerASCIIPrintf(viewer," DiagPivotThresh: %g\n",options.DiagPivotThresh);
75: PetscViewerASCIIPrintf(viewer," PivotGrowth: %s\n",(options.PivotGrowth != NO) ? "YES" : "NO");
76: PetscViewerASCIIPrintf(viewer," ConditionNumber: %s\n",(options.ConditionNumber != NO) ? "YES" : "NO");
77: PetscViewerASCIIPrintf(viewer," RowPerm: %" PetscInt_FMT "\n",options.RowPerm);
78: PetscViewerASCIIPrintf(viewer," ReplaceTinyPivot: %s\n",(options.ReplaceTinyPivot != NO) ? "YES" : "NO");
79: PetscViewerASCIIPrintf(viewer," PrintStat: %s\n",(options.PrintStat != NO) ? "YES" : "NO");
80: PetscViewerASCIIPrintf(viewer," lwork: %" PetscInt_FMT "\n",lu->lwork);
81: if (A->factortype == MAT_FACTOR_ILU) {
82: PetscViewerASCIIPrintf(viewer," ILU_DropTol: %g\n",options.ILU_DropTol);
83: PetscViewerASCIIPrintf(viewer," ILU_FillTol: %g\n",options.ILU_FillTol);
84: PetscViewerASCIIPrintf(viewer," ILU_FillFactor: %g\n",options.ILU_FillFactor);
85: PetscViewerASCIIPrintf(viewer," ILU_DropRule: %" PetscInt_FMT "\n",options.ILU_DropRule);
86: PetscViewerASCIIPrintf(viewer," ILU_Norm: %" PetscInt_FMT "\n",options.ILU_Norm);
87: PetscViewerASCIIPrintf(viewer," ILU_MILU: %" PetscInt_FMT "\n",options.ILU_MILU);
88: }
89: return 0;
90: }
92: PetscErrorCode MatSolve_SuperLU_Private(Mat A,Vec b,Vec x)
93: {
94: Mat_SuperLU *lu = (Mat_SuperLU*)A->data;
95: const PetscScalar *barray;
96: PetscScalar *xarray;
97: PetscInt info,i,n;
98: PetscReal ferr,berr;
99: static PetscBool cite = PETSC_FALSE;
101: if (lu->lwork == -1) return 0;
102: PetscCitationsRegister("@article{superlu99,\n author = {James W. Demmel and Stanley C. Eisenstat and\n John R. Gilbert and Xiaoye S. Li and Joseph W. H. Liu},\n title = {A supernodal approach to sparse partial pivoting},\n journal = {SIAM J. Matrix Analysis and Applications},\n year = {1999},\n volume = {20},\n number = {3},\n pages = {720-755}\n}\n",&cite);
104: VecGetLocalSize(x,&n);
105: lu->B.ncol = 1; /* Set the number of right-hand side */
106: if (lu->options.Equil && !lu->rhs_dup) {
107: /* superlu overwrites b when Equil is used, thus create rhs_dup to keep user's b unchanged */
108: PetscMalloc1(n,&lu->rhs_dup);
109: }
110: if (lu->options.Equil) {
111: /* Copy b into rsh_dup */
112: VecGetArrayRead(b,&barray);
113: PetscArraycpy(lu->rhs_dup,barray,n);
114: VecRestoreArrayRead(b,&barray);
115: barray = lu->rhs_dup;
116: } else {
117: VecGetArrayRead(b,&barray);
118: }
119: VecGetArray(x,&xarray);
121: #if defined(PETSC_USE_COMPLEX)
122: #if defined(PETSC_USE_REAL_SINGLE)
123: ((DNformat*)lu->B.Store)->nzval = (singlecomplex*)barray;
124: ((DNformat*)lu->X.Store)->nzval = (singlecomplex*)xarray;
125: #else
126: ((DNformat*)lu->B.Store)->nzval = (doublecomplex*)barray;
127: ((DNformat*)lu->X.Store)->nzval = (doublecomplex*)xarray;
128: #endif
129: #else
130: ((DNformat*)lu->B.Store)->nzval = (void*)barray;
131: ((DNformat*)lu->X.Store)->nzval = xarray;
132: #endif
134: lu->options.Fact = FACTORED; /* Indicate the factored form of A is supplied. */
135: if (A->factortype == MAT_FACTOR_LU) {
136: #if defined(PETSC_USE_COMPLEX)
137: #if defined(PETSC_USE_REAL_SINGLE)
138: PetscStackCall("SuperLU:cgssvx",cgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
139: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
140: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
141: #else
142: PetscStackCall("SuperLU:zgssvx",zgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
143: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
144: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
145: #endif
146: #else
147: #if defined(PETSC_USE_REAL_SINGLE)
148: PetscStackCall("SuperLU:sgssvx",sgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
149: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
150: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
151: #else
152: PetscStackCall("SuperLU:dgssvx",dgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
153: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
154: &lu->Glu,&lu->mem_usage, &lu->stat, &info));
155: #endif
156: #endif
157: } else if (A->factortype == MAT_FACTOR_ILU) {
158: #if defined(PETSC_USE_COMPLEX)
159: #if defined(PETSC_USE_REAL_SINGLE)
160: PetscStackCall("SuperLU:cgsisx",cgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
161: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
162: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
163: #else
164: PetscStackCall("SuperLU:zgsisx",zgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
165: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
166: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
167: #endif
168: #else
169: #if defined(PETSC_USE_REAL_SINGLE)
170: PetscStackCall("SuperLU:sgsisx",sgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
171: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
172: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
173: #else
174: PetscStackCall("SuperLU:dgsisx",dgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
175: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
176: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
177: #endif
178: #endif
179: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
180: if (!lu->options.Equil) {
181: VecRestoreArrayRead(b,&barray);
182: }
183: VecRestoreArray(x,&xarray);
185: if (!info || info == lu->A.ncol+1) {
186: if (lu->options.IterRefine) {
187: PetscPrintf(PETSC_COMM_SELF,"Iterative Refinement:\n");
188: PetscPrintf(PETSC_COMM_SELF," %8s%8s%16s%16s\n", "rhs", "Steps", "FERR", "BERR");
189: for (i = 0; i < 1; ++i) PetscPrintf(PETSC_COMM_SELF," %8d%8d%16e%16e\n", i+1, lu->stat.RefineSteps, ferr, berr);
190: }
191: } else if (info > 0) {
192: if (lu->lwork == -1) {
193: PetscPrintf(PETSC_COMM_SELF," ** Estimated memory: %" PetscInt_FMT " bytes\n", info - lu->A.ncol);
194: } else {
195: PetscPrintf(PETSC_COMM_SELF," Warning: gssvx() returns info %" PetscInt_FMT "\n",info);
196: }
199: if (lu->options.PrintStat) {
200: PetscPrintf(PETSC_COMM_SELF,"MatSolve__SuperLU():\n");
201: PetscStackCall("SuperLU:StatPrint",StatPrint(&lu->stat));
202: }
203: return 0;
204: }
206: PetscErrorCode MatSolve_SuperLU(Mat A,Vec b,Vec x)
207: {
208: Mat_SuperLU *lu = (Mat_SuperLU*)A->data;
210: if (A->factorerrortype) {
211: PetscInfo(A,"MatSolve is called with singular matrix factor, skip\n");
212: VecSetInf(x);
213: return 0;
214: }
216: lu->options.Trans = TRANS;
217: MatSolve_SuperLU_Private(A,b,x);
218: return 0;
219: }
221: PetscErrorCode MatSolveTranspose_SuperLU(Mat A,Vec b,Vec x)
222: {
223: Mat_SuperLU *lu = (Mat_SuperLU*)A->data;
225: if (A->factorerrortype) {
226: PetscInfo(A,"MatSolve is called with singular matrix factor, skip\n");
227: VecSetInf(x);
228: return 0;
229: }
231: lu->options.Trans = NOTRANS;
232: MatSolve_SuperLU_Private(A,b,x);
233: return 0;
234: }
236: static PetscErrorCode MatLUFactorNumeric_SuperLU(Mat F,Mat A,const MatFactorInfo *info)
237: {
238: Mat_SuperLU *lu = (Mat_SuperLU*)F->data;
239: Mat_SeqAIJ *aa;
240: PetscInt sinfo;
241: PetscReal ferr, berr;
242: NCformat *Ustore;
243: SCformat *Lstore;
245: if (lu->flg == SAME_NONZERO_PATTERN) { /* successing numerical factorization */
246: lu->options.Fact = SamePattern;
247: /* Ref: ~SuperLU_3.0/EXAMPLE/dlinsolx2.c */
248: Destroy_SuperMatrix_Store(&lu->A);
249: if (lu->A_dup) {
250: MatCopy_SeqAIJ(A,lu->A_dup,SAME_NONZERO_PATTERN);
251: }
253: if (lu->lwork >= 0) {
254: PetscStackCall("SuperLU:Destroy_SuperNode_Matrix",Destroy_SuperNode_Matrix(&lu->L));
255: PetscStackCall("SuperLU:Destroy_CompCol_Matrix",Destroy_CompCol_Matrix(&lu->U));
256: lu->options.Fact = SamePattern;
257: }
258: }
260: /* Create the SuperMatrix for lu->A=A^T:
261: Since SuperLU likes column-oriented matrices,we pass it the transpose,
262: and then solve A^T X = B in MatSolve(). */
263: if (lu->A_dup) {
264: aa = (Mat_SeqAIJ*)(lu->A_dup)->data;
265: } else {
266: aa = (Mat_SeqAIJ*)(A)->data;
267: }
268: #if defined(PETSC_USE_COMPLEX)
269: #if defined(PETSC_USE_REAL_SINGLE)
270: PetscStackCall("SuperLU:cCreate_CompCol_Matrix",cCreate_CompCol_Matrix(&lu->A,A->cmap->n,A->rmap->n,aa->nz,(singlecomplex*)aa->a,aa->j,aa->i,SLU_NC,SLU_C,SLU_GE));
271: #else
272: PetscStackCall("SuperLU:zCreate_CompCol_Matrix",zCreate_CompCol_Matrix(&lu->A,A->cmap->n,A->rmap->n,aa->nz,(doublecomplex*)aa->a,aa->j,aa->i,SLU_NC,SLU_Z,SLU_GE));
273: #endif
274: #else
275: #if defined(PETSC_USE_REAL_SINGLE)
276: PetscStackCall("SuperLU:sCreate_CompCol_Matrix",sCreate_CompCol_Matrix(&lu->A,A->cmap->n,A->rmap->n,aa->nz,aa->a,aa->j,aa->i,SLU_NC,SLU_S,SLU_GE));
277: #else
278: PetscStackCall("SuperLU:dCreate_CompCol_Matrix",dCreate_CompCol_Matrix(&lu->A,A->cmap->n,A->rmap->n,aa->nz,aa->a,aa->j,aa->i,SLU_NC,SLU_D,SLU_GE));
279: #endif
280: #endif
282: /* Numerical factorization */
283: lu->B.ncol = 0; /* Indicate not to solve the system */
284: if (F->factortype == MAT_FACTOR_LU) {
285: #if defined(PETSC_USE_COMPLEX)
286: #if defined(PETSC_USE_REAL_SINGLE)
287: PetscStackCall("SuperLU:cgssvx",cgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
288: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
289: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
290: #else
291: PetscStackCall("SuperLU:zgssvx",zgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
292: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
293: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
294: #endif
295: #else
296: #if defined(PETSC_USE_REAL_SINGLE)
297: PetscStackCall("SuperLU:sgssvx",sgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
298: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
299: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
300: #else
301: PetscStackCall("SuperLU:dgssvx",dgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
302: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
303: &lu->Glu,&lu->mem_usage, &lu->stat, &sinfo));
304: #endif
305: #endif
306: } else if (F->factortype == MAT_FACTOR_ILU) {
307: /* Compute the incomplete factorization, condition number and pivot growth */
308: #if defined(PETSC_USE_COMPLEX)
309: #if defined(PETSC_USE_REAL_SINGLE)
310: PetscStackCall("SuperLU:cgsisx",cgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r,lu->etree, lu->equed, lu->R, lu->C,
311: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
312: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
313: #else
314: PetscStackCall("SuperLU:zgsisx",zgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r,lu->etree, lu->equed, lu->R, lu->C,
315: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
316: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
317: #endif
318: #else
319: #if defined(PETSC_USE_REAL_SINGLE)
320: PetscStackCall("SuperLU:sgsisx",sgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
321: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
322: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
323: #else
324: PetscStackCall("SuperLU:dgsisx",dgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
325: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
326: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
327: #endif
328: #endif
329: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
330: if (!sinfo || sinfo == lu->A.ncol+1) {
331: if (lu->options.PivotGrowth) {
332: PetscPrintf(PETSC_COMM_SELF," Recip. pivot growth = %e\n", lu->rpg);
333: }
334: if (lu->options.ConditionNumber) {
335: PetscPrintf(PETSC_COMM_SELF," Recip. condition number = %e\n", lu->rcond);
336: }
337: } else if (sinfo > 0) {
338: if (A->erroriffailure) {
339: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot in row %" PetscInt_FMT,sinfo);
340: } else {
341: if (sinfo <= lu->A.ncol) {
342: if (lu->options.ILU_FillTol == 0.0) {
343: F->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
344: }
345: PetscInfo(F,"Number of zero pivots %" PetscInt_FMT ", ILU_FillTol %g\n",sinfo,lu->options.ILU_FillTol);
346: } else if (sinfo == lu->A.ncol + 1) {
347: /*
348: U is nonsingular, but RCOND is less than machine
349: precision, meaning that the matrix is singular to
350: working precision. Nevertheless, the solution and
351: error bounds are computed because there are a number
352: of situations where the computed solution can be more
353: accurate than the value of RCOND would suggest.
354: */
355: PetscInfo(F,"Matrix factor U is nonsingular, but is singular to working precision. The solution is computed. info %" PetscInt_FMT,sinfo);
356: } else { /* sinfo > lu->A.ncol + 1 */
357: F->factorerrortype = MAT_FACTOR_OUTMEMORY;
358: PetscInfo(F,"Number of bytes allocated when memory allocation fails %" PetscInt_FMT "\n",sinfo);
359: }
360: }
361: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB, "info = %" PetscInt_FMT ", the %" PetscInt_FMT "-th argument in gssvx() had an illegal value", sinfo,-sinfo);
363: if (lu->options.PrintStat) {
364: PetscPrintf(PETSC_COMM_SELF,"MatLUFactorNumeric_SuperLU():\n");
365: PetscStackCall("SuperLU:StatPrint",StatPrint(&lu->stat));
366: Lstore = (SCformat*) lu->L.Store;
367: Ustore = (NCformat*) lu->U.Store;
368: PetscPrintf(PETSC_COMM_SELF," No of nonzeros in factor L = %" PetscInt_FMT "\n", Lstore->nnz);
369: PetscPrintf(PETSC_COMM_SELF," No of nonzeros in factor U = %" PetscInt_FMT "\n", Ustore->nnz);
370: PetscPrintf(PETSC_COMM_SELF," No of nonzeros in L+U = %" PetscInt_FMT "\n", Lstore->nnz + Ustore->nnz - lu->A.ncol);
371: PetscPrintf(PETSC_COMM_SELF," L\\U MB %.3f\ttotal MB needed %.3f\n",lu->mem_usage.for_lu/1e6, lu->mem_usage.total_needed/1e6);
372: }
374: lu->flg = SAME_NONZERO_PATTERN;
375: F->ops->solve = MatSolve_SuperLU;
376: F->ops->solvetranspose = MatSolveTranspose_SuperLU;
377: F->ops->matsolve = NULL;
378: return 0;
379: }
381: static PetscErrorCode MatDestroy_SuperLU(Mat A)
382: {
383: Mat_SuperLU *lu=(Mat_SuperLU*)A->data;
385: if (lu->CleanUpSuperLU) { /* Free the SuperLU datastructures */
386: PetscStackCall("SuperLU:Destroy_SuperMatrix_Store",Destroy_SuperMatrix_Store(&lu->A));
387: PetscStackCall("SuperLU:Destroy_SuperMatrix_Store",Destroy_SuperMatrix_Store(&lu->B));
388: PetscStackCall("SuperLU:Destroy_SuperMatrix_Store",Destroy_SuperMatrix_Store(&lu->X));
389: PetscStackCall("SuperLU:StatFree",StatFree(&lu->stat));
390: if (lu->lwork >= 0) {
391: PetscStackCall("SuperLU:Destroy_SuperNode_Matrix",Destroy_SuperNode_Matrix(&lu->L));
392: PetscStackCall("SuperLU:Destroy_CompCol_Matrix",Destroy_CompCol_Matrix(&lu->U));
393: }
394: }
395: PetscFree(lu->etree);
396: PetscFree(lu->perm_r);
397: PetscFree(lu->perm_c);
398: PetscFree(lu->R);
399: PetscFree(lu->C);
400: PetscFree(lu->rhs_dup);
401: MatDestroy(&lu->A_dup);
402: PetscFree(A->data);
404: /* clear composed functions */
405: PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);
406: PetscObjectComposeFunction((PetscObject)A,"MatSuperluSetILUDropTol_C",NULL);
407: return 0;
408: }
410: static PetscErrorCode MatView_SuperLU(Mat A,PetscViewer viewer)
411: {
412: PetscBool iascii;
413: PetscViewerFormat format;
415: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
416: if (iascii) {
417: PetscViewerGetFormat(viewer,&format);
418: if (format == PETSC_VIEWER_ASCII_INFO) {
419: MatView_Info_SuperLU(A,viewer);
420: }
421: }
422: return 0;
423: }
425: PetscErrorCode MatMatSolve_SuperLU(Mat A,Mat B,Mat X)
426: {
427: Mat_SuperLU *lu = (Mat_SuperLU*)A->data;
428: PetscBool flg;
430: PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
432: PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
434: lu->options.Trans = TRANS;
435: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_SuperLU() is not implemented yet");
436: return 0;
437: }
439: static PetscErrorCode MatLUFactorSymbolic_SuperLU(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
440: {
441: Mat_SuperLU *lu = (Mat_SuperLU*)(F->data);
443: lu->flg = DIFFERENT_NONZERO_PATTERN;
444: lu->CleanUpSuperLU = PETSC_TRUE;
445: F->ops->lufactornumeric = MatLUFactorNumeric_SuperLU;
447: /* if we are here, the nonzero pattern has changed unless the user explicitly called MatLUFactorSymbolic */
448: MatDestroy(&lu->A_dup);
449: if (lu->needconversion) {
450: MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&lu->A_dup);
451: }
452: if (lu->options.Equil == YES && !lu->A_dup) { /* superlu overwrites input matrix and rhs when Equil is used, thus create A_dup to keep user's A unchanged */
453: MatDuplicate_SeqAIJ(A,MAT_COPY_VALUES,&lu->A_dup);
454: }
455: return 0;
456: }
458: static PetscErrorCode MatSuperluSetILUDropTol_SuperLU(Mat F,PetscReal dtol)
459: {
460: Mat_SuperLU *lu= (Mat_SuperLU*)F->data;
462: lu->options.ILU_DropTol = dtol;
463: return 0;
464: }
466: /*@
467: MatSuperluSetILUDropTol - Set SuperLU ILU drop tolerance
469: Logically Collective on Mat
471: Input Parameters:
472: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-SuperLU interface
473: - dtol - drop tolerance
475: Options Database:
476: . -mat_superlu_ilu_droptol <dtol> - the drop tolerance
478: Level: beginner
480: References:
481: . * - SuperLU Users' Guide
483: .seealso: MatGetFactor()
484: @*/
485: PetscErrorCode MatSuperluSetILUDropTol(Mat F,PetscReal dtol)
486: {
489: PetscTryMethod(F,"MatSuperluSetILUDropTol_C",(Mat,PetscReal),(F,dtol));
490: return 0;
491: }
493: PetscErrorCode MatFactorGetSolverType_seqaij_superlu(Mat A,MatSolverType *type)
494: {
495: *type = MATSOLVERSUPERLU;
496: return 0;
497: }
499: /*MC
500: MATSOLVERSUPERLU = "superlu" - A solver package providing solvers LU and ILU for sequential matrices
501: via the external package SuperLU.
503: Use ./configure --download-superlu to have PETSc installed with SuperLU
505: Use -pc_type lu -pc_factor_mat_solver_type superlu to use this direct solver
507: Options Database Keys:
508: + -mat_superlu_equil <FALSE> - Equil (None)
509: . -mat_superlu_colperm <COLAMD> - (choose one of) NATURAL MMD_ATA MMD_AT_PLUS_A COLAMD
510: . -mat_superlu_iterrefine <NOREFINE> - (choose one of) NOREFINE SINGLE DOUBLE EXTRA
511: . -mat_superlu_symmetricmode: <FALSE> - SymmetricMode (None)
512: . -mat_superlu_diagpivotthresh <1> - DiagPivotThresh (None)
513: . -mat_superlu_pivotgrowth <FALSE> - PivotGrowth (None)
514: . -mat_superlu_conditionnumber <FALSE> - ConditionNumber (None)
515: . -mat_superlu_rowperm <NOROWPERM> - (choose one of) NOROWPERM LargeDiag
516: . -mat_superlu_replacetinypivot <FALSE> - ReplaceTinyPivot (None)
517: . -mat_superlu_printstat <FALSE> - PrintStat (None)
518: . -mat_superlu_lwork <0> - size of work array in bytes used by factorization (None)
519: . -mat_superlu_ilu_droptol <0> - ILU_DropTol (None)
520: . -mat_superlu_ilu_filltol <0> - ILU_FillTol (None)
521: . -mat_superlu_ilu_fillfactor <0> - ILU_FillFactor (None)
522: . -mat_superlu_ilu_droprull <0> - ILU_DropRule (None)
523: . -mat_superlu_ilu_norm <0> - ILU_Norm (None)
524: - -mat_superlu_ilu_milu <0> - ILU_MILU (None)
526: Notes:
527: Do not confuse this with MATSOLVERSUPERLU_DIST which is for parallel sparse solves
529: Cannot currently use ordering provided by PETSc.
531: Level: beginner
533: .seealso: PCLU, PCILU, MATSOLVERSUPERLU_DIST, MATSOLVERMUMPS, PCFactorSetMatSolverType(), MatSolverType
534: M*/
536: static PetscErrorCode MatGetFactor_seqaij_superlu(Mat A,MatFactorType ftype,Mat *F)
537: {
538: Mat B;
539: Mat_SuperLU *lu;
540: PetscInt indx,m=A->rmap->n,n=A->cmap->n;
541: PetscBool flg,set;
542: PetscReal real_input;
543: const char *colperm[] ={"NATURAL","MMD_ATA","MMD_AT_PLUS_A","COLAMD"}; /* MY_PERMC - not supported by the petsc interface yet */
544: const char *iterrefine[]={"NOREFINE", "SINGLE", "DOUBLE", "EXTRA"};
545: const char *rowperm[] ={"NOROWPERM", "LargeDiag"}; /* MY_PERMC - not supported by the petsc interface yet */
548: MatCreate(PetscObjectComm((PetscObject)A),&B);
549: MatSetSizes(B,A->rmap->n,A->cmap->n,PETSC_DETERMINE,PETSC_DETERMINE);
550: PetscStrallocpy("superlu",&((PetscObject)B)->type_name);
551: MatSetUp(B);
552: B->trivialsymbolic = PETSC_TRUE;
553: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
554: B->ops->lufactorsymbolic = MatLUFactorSymbolic_SuperLU;
555: B->ops->ilufactorsymbolic = MatLUFactorSymbolic_SuperLU;
556: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
558: PetscFree(B->solvertype);
559: PetscStrallocpy(MATSOLVERSUPERLU,&B->solvertype);
561: B->ops->getinfo = MatGetInfo_External;
562: B->ops->destroy = MatDestroy_SuperLU;
563: B->ops->view = MatView_SuperLU;
564: B->factortype = ftype;
565: B->assembled = PETSC_TRUE; /* required by -ksp_view */
566: B->preallocated = PETSC_TRUE;
568: PetscNewLog(B,&lu);
570: if (ftype == MAT_FACTOR_LU) {
571: set_default_options(&lu->options);
572: /* Comments from SuperLU_4.0/SRC/dgssvx.c:
573: "Whether or not the system will be equilibrated depends on the
574: scaling of the matrix A, but if equilibration is used, A is
575: overwritten by diag(R)*A*diag(C) and B by diag(R)*B
576: (if options->Trans=NOTRANS) or diag(C)*B (if options->Trans = TRANS or CONJ)."
577: We set 'options.Equil = NO' as default because additional space is needed for it.
578: */
579: lu->options.Equil = NO;
580: } else if (ftype == MAT_FACTOR_ILU) {
581: /* Set the default input options of ilu: */
582: PetscStackCall("SuperLU:ilu_set_default_options",ilu_set_default_options(&lu->options));
583: }
584: lu->options.PrintStat = NO;
586: /* Initialize the statistics variables. */
587: PetscStackCall("SuperLU:StatInit",StatInit(&lu->stat));
588: lu->lwork = 0; /* allocate space internally by system malloc */
590: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"SuperLU Options","Mat");
591: PetscOptionsBool("-mat_superlu_equil","Equil","None",(PetscBool)lu->options.Equil,(PetscBool*)&lu->options.Equil,NULL);
592: PetscOptionsEList("-mat_superlu_colperm","ColPerm","None",colperm,4,colperm[3],&indx,&flg);
593: if (flg) lu->options.ColPerm = (colperm_t)indx;
594: PetscOptionsEList("-mat_superlu_iterrefine","IterRefine","None",iterrefine,4,iterrefine[0],&indx,&flg);
595: if (flg) lu->options.IterRefine = (IterRefine_t)indx;
596: PetscOptionsBool("-mat_superlu_symmetricmode","SymmetricMode","None",(PetscBool)lu->options.SymmetricMode,&flg,&set);
597: if (set && flg) lu->options.SymmetricMode = YES;
598: PetscOptionsReal("-mat_superlu_diagpivotthresh","DiagPivotThresh","None",lu->options.DiagPivotThresh,&real_input,&flg);
599: if (flg) lu->options.DiagPivotThresh = (double) real_input;
600: PetscOptionsBool("-mat_superlu_pivotgrowth","PivotGrowth","None",(PetscBool)lu->options.PivotGrowth,&flg,&set);
601: if (set && flg) lu->options.PivotGrowth = YES;
602: PetscOptionsBool("-mat_superlu_conditionnumber","ConditionNumber","None",(PetscBool)lu->options.ConditionNumber,&flg,&set);
603: if (set && flg) lu->options.ConditionNumber = YES;
604: PetscOptionsEList("-mat_superlu_rowperm","rowperm","None",rowperm,2,rowperm[lu->options.RowPerm],&indx,&flg);
605: if (flg) lu->options.RowPerm = (rowperm_t)indx;
606: PetscOptionsBool("-mat_superlu_replacetinypivot","ReplaceTinyPivot","None",(PetscBool)lu->options.ReplaceTinyPivot,&flg,&set);
607: if (set && flg) lu->options.ReplaceTinyPivot = YES;
608: PetscOptionsBool("-mat_superlu_printstat","PrintStat","None",(PetscBool)lu->options.PrintStat,&flg,&set);
609: if (set && flg) lu->options.PrintStat = YES;
610: PetscOptionsInt("-mat_superlu_lwork","size of work array in bytes used by factorization","None",lu->lwork,&lu->lwork,NULL);
611: if (lu->lwork > 0) {
612: /* lwork is in bytes, hence PetscMalloc() is used here, not PetscMalloc1()*/
613: PetscMalloc(lu->lwork,&lu->work);
614: } else if (lu->lwork != 0 && lu->lwork != -1) {
615: PetscPrintf(PETSC_COMM_SELF," Warning: lwork %" PetscInt_FMT " is not supported by SUPERLU. The default lwork=0 is used.\n",lu->lwork);
616: lu->lwork = 0;
617: }
618: /* ilu options */
619: PetscOptionsReal("-mat_superlu_ilu_droptol","ILU_DropTol","None",lu->options.ILU_DropTol,&real_input,&flg);
620: if (flg) lu->options.ILU_DropTol = (double) real_input;
621: PetscOptionsReal("-mat_superlu_ilu_filltol","ILU_FillTol","None",lu->options.ILU_FillTol,&real_input,&flg);
622: if (flg) lu->options.ILU_FillTol = (double) real_input;
623: PetscOptionsReal("-mat_superlu_ilu_fillfactor","ILU_FillFactor","None",lu->options.ILU_FillFactor,&real_input,&flg);
624: if (flg) lu->options.ILU_FillFactor = (double) real_input;
625: PetscOptionsInt("-mat_superlu_ilu_droprull","ILU_DropRule","None",lu->options.ILU_DropRule,&lu->options.ILU_DropRule,NULL);
626: PetscOptionsInt("-mat_superlu_ilu_norm","ILU_Norm","None",lu->options.ILU_Norm,&indx,&flg);
627: if (flg) lu->options.ILU_Norm = (norm_t)indx;
628: PetscOptionsInt("-mat_superlu_ilu_milu","ILU_MILU","None",lu->options.ILU_MILU,&indx,&flg);
629: if (flg) lu->options.ILU_MILU = (milu_t)indx;
630: PetscOptionsEnd();
632: /* Allocate spaces (notice sizes are for the transpose) */
633: PetscMalloc1(m,&lu->etree);
634: PetscMalloc1(n,&lu->perm_r);
635: PetscMalloc1(m,&lu->perm_c);
636: PetscMalloc1(n,&lu->R);
637: PetscMalloc1(m,&lu->C);
639: /* create rhs and solution x without allocate space for .Store */
640: #if defined(PETSC_USE_COMPLEX)
641: #if defined(PETSC_USE_REAL_SINGLE)
642: PetscStackCall("SuperLU:cCreate_Dense_Matrix(",cCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_C, SLU_GE));
643: PetscStackCall("SuperLU:cCreate_Dense_Matrix(",cCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_C, SLU_GE));
644: #else
645: PetscStackCall("SuperLU:zCreate_Dense_Matrix",zCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_Z, SLU_GE));
646: PetscStackCall("SuperLU:zCreate_Dense_Matrix",zCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_Z, SLU_GE));
647: #endif
648: #else
649: #if defined(PETSC_USE_REAL_SINGLE)
650: PetscStackCall("SuperLU:sCreate_Dense_Matrix",sCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_S, SLU_GE));
651: PetscStackCall("SuperLU:sCreate_Dense_Matrix",sCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_S, SLU_GE));
652: #else
653: PetscStackCall("SuperLU:dCreate_Dense_Matrix",dCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_D, SLU_GE));
654: PetscStackCall("SuperLU:dCreate_Dense_Matrix",dCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_D, SLU_GE));
655: #endif
656: #endif
658: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqaij_superlu);
659: PetscObjectComposeFunction((PetscObject)B,"MatSuperluSetILUDropTol_C",MatSuperluSetILUDropTol_SuperLU);
660: B->data = lu;
662: *F = B;
663: return 0;
664: }
666: static PetscErrorCode MatGetFactor_seqsell_superlu(Mat A,MatFactorType ftype,Mat *F)
667: {
668: Mat_SuperLU *lu;
670: MatGetFactor_seqaij_superlu(A,ftype,F);
671: lu = (Mat_SuperLU*)((*F)->data);
672: lu->needconversion = PETSC_TRUE;
673: return 0;
674: }
676: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_SuperLU(void)
677: {
678: MatSolverTypeRegister(MATSOLVERSUPERLU,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_seqaij_superlu);
679: MatSolverTypeRegister(MATSOLVERSUPERLU,MATSEQAIJ,MAT_FACTOR_ILU,MatGetFactor_seqaij_superlu);
680: MatSolverTypeRegister(MATSOLVERSUPERLU,MATSEQSELL,MAT_FACTOR_LU,MatGetFactor_seqsell_superlu);
681: MatSolverTypeRegister(MATSOLVERSUPERLU,MATSEQSELL,MAT_FACTOR_ILU,MatGetFactor_seqsell_superlu);
682: return 0;
683: }