Actual source code: spooles.c
1: #define PETCSMAT_DLL
3: /*
4: Provides an interface to the Spooles serial sparse solver
5: */
6: #include src/mat/impls/aij/seq/aij.h
7: #include src/mat/impls/sbaij/seq/sbaij.h
8: #include src/mat/impls/aij/seq/spooles/spooles.h
10: /* make sun CC happy */
11: static void (*f)(void);
16: PetscErrorCode MatConvert_Spooles_Base(Mat A,MatType type,MatReuse reuse,Mat *newmat)
17: {
19: Mat B=*newmat;
20: Mat_Spooles *lu=(Mat_Spooles*)A->spptr;
23: if (reuse == MAT_INITIAL_MATRIX) {
24: MatDuplicate(A,MAT_COPY_VALUES,&B);
25: }
26: /* Reset the stashed function pointers set by inherited routines */
27: B->ops->duplicate = lu->MatDuplicate;
28: B->ops->choleskyfactorsymbolic = lu->MatCholeskyFactorSymbolic;
29: B->ops->lufactorsymbolic = lu->MatLUFactorSymbolic;
30: B->ops->view = lu->MatView;
31: B->ops->assemblyend = lu->MatAssemblyEnd;
32: B->ops->destroy = lu->MatDestroy;
34: PetscObjectQueryFunction((PetscObject)B,"MatMPISBAIJSetPreallocation_C",&f);
35: if (f) {
36: PetscObjectComposeFunction((PetscObject)B,"MatMPISBAIJSetPreallocation_C","",(PetscVoidFunction)lu->MatPreallocate);
37: }
39: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaijspooles_seqaij_C","",PETSC_NULL);
40: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijspooles_C","",PETSC_NULL);
41: PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpiaijspooles_mpiaij_C","",PETSC_NULL);
42: PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpiaij_mpiaijspooles_C","",PETSC_NULL);
43: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaijspooles_seqsbaij_C","",PETSC_NULL);
44: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqsbaijspooles_C","",PETSC_NULL);
45: PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpisbaijspooles_mpisbaij_C","",PETSC_NULL);
46: PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpisbaij_mpisbaijspooles_C","",PETSC_NULL);
48: PetscObjectChangeTypeName((PetscObject)B,type);
49: *newmat = B;
50: return(0);
51: }
56: PetscErrorCode MatDestroy_SeqAIJSpooles(Mat A)
57: {
58: Mat_Spooles *lu = (Mat_Spooles*)A->spptr;
60:
62: if (lu->CleanUpSpooles) {
63: FrontMtx_free(lu->frontmtx);
64: IV_free(lu->newToOldIV);
65: IV_free(lu->oldToNewIV);
66: InpMtx_free(lu->mtxA);
67: ETree_free(lu->frontETree);
68: IVL_free(lu->symbfacIVL);
69: SubMtxManager_free(lu->mtxmanager);
70: Graph_free(lu->graph);
71: }
72: MatConvert_Spooles_Base(A,lu->basetype,MAT_REUSE_MATRIX,&A);
73: (*A->ops->destroy)(A);
74: return(0);
75: }
79: PetscErrorCode MatSolve_SeqAIJSpooles(Mat A,Vec b,Vec x)
80: {
81: Mat_Spooles *lu = (Mat_Spooles*)A->spptr;
82: PetscScalar *array;
83: DenseMtx *mtxY, *mtxX ;
84: PetscErrorCode ierr;
85: PetscInt irow,neqns=A->cmap.n,nrow=A->rmap.n,*iv;
86: #if defined(PETSC_USE_COMPLEX)
87: double x_real,x_imag;
88: #else
89: double *entX;
90: #endif
93: mtxY = DenseMtx_new();
94: DenseMtx_init(mtxY, lu->options.typeflag, 0, 0, nrow, 1, 1, nrow); /* column major */
95: VecGetArray(b,&array);
97: if (lu->options.useQR) { /* copy b to mtxY */
98: for ( irow = 0 ; irow < nrow; irow++ )
99: #if !defined(PETSC_USE_COMPLEX)
100: DenseMtx_setRealEntry(mtxY, irow, 0, *array++);
101: #else
102: DenseMtx_setComplexEntry(mtxY, irow, 0, PetscRealPart(array[irow]), PetscImaginaryPart(array[irow]));
103: #endif
104: } else { /* copy permuted b to mtxY */
105: iv = IV_entries(lu->oldToNewIV);
106: for ( irow = 0 ; irow < nrow; irow++ )
107: #if !defined(PETSC_USE_COMPLEX)
108: DenseMtx_setRealEntry(mtxY, *iv++, 0, *array++);
109: #else
110: DenseMtx_setComplexEntry(mtxY,*iv++,0,PetscRealPart(array[irow]),PetscImaginaryPart(array[irow]));
111: #endif
112: }
113: VecRestoreArray(b,&array);
115: mtxX = DenseMtx_new();
116: DenseMtx_init(mtxX, lu->options.typeflag, 0, 0, neqns, 1, 1, neqns);
117: if (lu->options.useQR) {
118: FrontMtx_QR_solve(lu->frontmtx, lu->mtxA, mtxX, mtxY, lu->mtxmanager,
119: lu->cpus, lu->options.msglvl, lu->options.msgFile);
120: } else {
121: FrontMtx_solve(lu->frontmtx, mtxX, mtxY, lu->mtxmanager,
122: lu->cpus, lu->options.msglvl, lu->options.msgFile);
123: }
124: if ( lu->options.msglvl > 2 ) {
125: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n right hand side matrix after permutation");
126: DenseMtx_writeForHumanEye(mtxY, lu->options.msgFile);
127: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n solution matrix in new ordering");
128: DenseMtx_writeForHumanEye(mtxX, lu->options.msgFile);
129: fflush(lu->options.msgFile);
130: }
132: /* permute solution into original ordering, then copy to x */
133: DenseMtx_permuteRows(mtxX, lu->newToOldIV);
134: VecGetArray(x,&array);
136: #if !defined(PETSC_USE_COMPLEX)
137: entX = DenseMtx_entries(mtxX);
138: DVcopy(neqns, array, entX);
139: #else
140: for (irow=0; irow<nrow; irow++){
141: DenseMtx_complexEntry(mtxX,irow,0,&x_real,&x_imag);
142: array[irow] = x_real+x_imag*PETSC_i;
143: }
144: #endif
146: VecRestoreArray(x,&array);
147:
148: /* free memory */
149: DenseMtx_free(mtxX);
150: DenseMtx_free(mtxY);
151: return(0);
152: }
156: PetscErrorCode MatFactorNumeric_SeqAIJSpooles(Mat A,MatFactorInfo *info,Mat *F)
157: {
158: Mat_Spooles *lu = (Mat_Spooles*)(*F)->spptr;
159: ChvManager *chvmanager ;
160: Chv *rootchv ;
161: IVL *adjIVL;
162: PetscErrorCode ierr;
163: PetscInt nz,nrow=A->rmap.n,irow,nedges,neqns=A->cmap.n,*ai,*aj,i,*diag=0,fierr;
164: PetscScalar *av;
165: double cputotal,facops;
166: #if defined(PETSC_USE_COMPLEX)
167: PetscInt nz_row,*aj_tmp;
168: PetscScalar *av_tmp;
169: #else
170: PetscInt *ivec1,*ivec2,j;
171: double *dvec;
172: #endif
173: PetscTruth isAIJ,isSeqAIJ;
174:
176: if (lu->flg == DIFFERENT_NONZERO_PATTERN) { /* first numeric factorization */
177: (*F)->ops->solve = MatSolve_SeqAIJSpooles;
178: (*F)->ops->destroy = MatDestroy_SeqAIJSpooles;
179: (*F)->assembled = PETSC_TRUE;
180:
181: /* set Spooles options */
182: SetSpoolesOptions(A, &lu->options);
184: lu->mtxA = InpMtx_new();
185: }
187: /* copy A to Spooles' InpMtx object */
188: PetscTypeCompare((PetscObject)A,MATSEQAIJSPOOLES,&isSeqAIJ);
189: PetscTypeCompare((PetscObject)A,MATAIJSPOOLES,&isAIJ);
190: if (isSeqAIJ || isAIJ){
191: Mat_SeqAIJ *mat = (Mat_SeqAIJ*)A->data;
192: ai=mat->i; aj=mat->j; av=mat->a;
193: if (lu->options.symflag == SPOOLES_NONSYMMETRIC) {
194: nz=mat->nz;
195: } else { /* SPOOLES_SYMMETRIC || SPOOLES_HERMITIAN */
196: nz=(mat->nz + A->rmap.n)/2;
197: diag=mat->diag;
198: }
199: } else { /* A is SBAIJ */
200: Mat_SeqSBAIJ *mat = (Mat_SeqSBAIJ*)A->data;
201: ai=mat->i; aj=mat->j; av=mat->a;
202: nz=mat->nz;
203: }
204: InpMtx_init(lu->mtxA, INPMTX_BY_ROWS, lu->options.typeflag, nz, 0);
205:
206: #if defined(PETSC_USE_COMPLEX)
207: for (irow=0; irow<nrow; irow++) {
208: if ( lu->options.symflag == SPOOLES_NONSYMMETRIC || !isAIJ){
209: nz_row = ai[irow+1] - ai[irow];
210: aj_tmp = aj + ai[irow];
211: av_tmp = av + ai[irow];
212: } else {
213: nz_row = ai[irow+1] - diag[irow];
214: aj_tmp = aj + diag[irow];
215: av_tmp = av + diag[irow];
216: }
217: for (i=0; i<nz_row; i++){
218: InpMtx_inputComplexEntry(lu->mtxA, irow, *aj_tmp++,PetscRealPart(*av_tmp),PetscImaginaryPart(*av_tmp));
219: av_tmp++;
220: }
221: }
222: #else
223: ivec1 = InpMtx_ivec1(lu->mtxA);
224: ivec2 = InpMtx_ivec2(lu->mtxA);
225: dvec = InpMtx_dvec(lu->mtxA);
226: if ( lu->options.symflag == SPOOLES_NONSYMMETRIC || !isAIJ){
227: for (irow = 0; irow < nrow; irow++){
228: for (i = ai[irow]; i<ai[irow+1]; i++) ivec1[i] = irow;
229: }
230: IVcopy(nz, ivec2, aj);
231: DVcopy(nz, dvec, av);
232: } else {
233: nz = 0;
234: for (irow = 0; irow < nrow; irow++){
235: for (j = diag[irow]; j<ai[irow+1]; j++) {
236: ivec1[nz] = irow;
237: ivec2[nz] = aj[j];
238: dvec[nz] = av[j];
239: nz++;
240: }
241: }
242: }
243: InpMtx_inputRealTriples(lu->mtxA, nz, ivec1, ivec2, dvec);
244: #endif
246: InpMtx_changeStorageMode(lu->mtxA, INPMTX_BY_VECTORS);
247: if ( lu->options.msglvl > 0 ) {
248: printf("\n\n input matrix");
249: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n input matrix");
250: InpMtx_writeForHumanEye(lu->mtxA, lu->options.msgFile);
251: fflush(lu->options.msgFile);
252: }
254: if ( lu->flg == DIFFERENT_NONZERO_PATTERN){ /* first numeric factorization */
255: /*---------------------------------------------------
256: find a low-fill ordering
257: (1) create the Graph object
258: (2) order the graph
259: -------------------------------------------------------*/
260: if (lu->options.useQR){
261: adjIVL = InpMtx_adjForATA(lu->mtxA);
262: } else {
263: adjIVL = InpMtx_fullAdjacency(lu->mtxA);
264: }
265: nedges = IVL_tsize(adjIVL);
267: lu->graph = Graph_new();
268: Graph_init2(lu->graph, 0, neqns, 0, nedges, neqns, nedges, adjIVL, NULL, NULL);
269: if ( lu->options.msglvl > 2 ) {
270: if (lu->options.useQR){
271: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n graph of A^T A");
272: } else {
273: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n graph of the input matrix");
274: }
275: Graph_writeForHumanEye(lu->graph, lu->options.msgFile);
276: fflush(lu->options.msgFile);
277: }
279: switch (lu->options.ordering) {
280: case 0:
281: lu->frontETree = orderViaBestOfNDandMS(lu->graph,
282: lu->options.maxdomainsize, lu->options.maxzeros, lu->options.maxsize,
283: lu->options.seed, lu->options.msglvl, lu->options.msgFile); break;
284: case 1:
285: lu->frontETree = orderViaMMD(lu->graph,lu->options.seed,lu->options.msglvl,lu->options.msgFile); break;
286: case 2:
287: lu->frontETree = orderViaMS(lu->graph, lu->options.maxdomainsize,
288: lu->options.seed,lu->options.msglvl,lu->options.msgFile); break;
289: case 3:
290: lu->frontETree = orderViaND(lu->graph, lu->options.maxdomainsize,
291: lu->options.seed,lu->options.msglvl,lu->options.msgFile); break;
292: default:
293: SETERRQ(PETSC_ERR_ARG_WRONG,"Unknown Spooles's ordering");
294: }
296: if ( lu->options.msglvl > 0 ) {
297: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n front tree from ordering");
298: ETree_writeForHumanEye(lu->frontETree, lu->options.msgFile);
299: fflush(lu->options.msgFile);
300: }
301:
302: /* get the permutation, permute the front tree */
303: lu->oldToNewIV = ETree_oldToNewVtxPerm(lu->frontETree);
304: lu->oldToNew = IV_entries(lu->oldToNewIV);
305: lu->newToOldIV = ETree_newToOldVtxPerm(lu->frontETree);
306: if (!lu->options.useQR) ETree_permuteVertices(lu->frontETree, lu->oldToNewIV);
308: /* permute the matrix */
309: if (lu->options.useQR){
310: InpMtx_permute(lu->mtxA, NULL, lu->oldToNew);
311: } else {
312: InpMtx_permute(lu->mtxA, lu->oldToNew, lu->oldToNew);
313: if ( lu->options.symflag == SPOOLES_SYMMETRIC) {
314: InpMtx_mapToUpperTriangle(lu->mtxA);
315: }
316: #if defined(PETSC_USE_COMPLEX)
317: if ( lu->options.symflag == SPOOLES_HERMITIAN ) {
318: InpMtx_mapToUpperTriangleH(lu->mtxA);
319: }
320: #endif
321: InpMtx_changeCoordType(lu->mtxA, INPMTX_BY_CHEVRONS);
322: }
323: InpMtx_changeStorageMode(lu->mtxA, INPMTX_BY_VECTORS);
325: /* get symbolic factorization */
326: if (lu->options.useQR){
327: lu->symbfacIVL = SymbFac_initFromGraph(lu->frontETree, lu->graph);
328: IVL_overwrite(lu->symbfacIVL, lu->oldToNewIV);
329: IVL_sortUp(lu->symbfacIVL);
330: ETree_permuteVertices(lu->frontETree, lu->oldToNewIV);
331: } else {
332: lu->symbfacIVL = SymbFac_initFromInpMtx(lu->frontETree, lu->mtxA);
333: }
334: if ( lu->options.msglvl > 2 ) {
335: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n old-to-new permutation vector");
336: IV_writeForHumanEye(lu->oldToNewIV, lu->options.msgFile);
337: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n new-to-old permutation vector");
338: IV_writeForHumanEye(lu->newToOldIV, lu->options.msgFile);
339: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n front tree after permutation");
340: ETree_writeForHumanEye(lu->frontETree, lu->options.msgFile);
341: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n input matrix after permutation");
342: InpMtx_writeForHumanEye(lu->mtxA, lu->options.msgFile);
343: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n symbolic factorization");
344: IVL_writeForHumanEye(lu->symbfacIVL, lu->options.msgFile);
345: fflush(lu->options.msgFile);
346: }
348: lu->frontmtx = FrontMtx_new();
349: lu->mtxmanager = SubMtxManager_new();
350: SubMtxManager_init(lu->mtxmanager, NO_LOCK, 0);
352: } else { /* new num factorization using previously computed symbolic factor */
354: if (lu->options.pivotingflag) { /* different FrontMtx is required */
355: FrontMtx_free(lu->frontmtx);
356: lu->frontmtx = FrontMtx_new();
357: } else {
358: FrontMtx_clearData (lu->frontmtx);
359: }
361: SubMtxManager_free(lu->mtxmanager);
362: lu->mtxmanager = SubMtxManager_new();
363: SubMtxManager_init(lu->mtxmanager, NO_LOCK, 0);
365: /* permute mtxA */
366: if (lu->options.useQR){
367: InpMtx_permute(lu->mtxA, NULL, lu->oldToNew);
368: } else {
369: InpMtx_permute(lu->mtxA, lu->oldToNew, lu->oldToNew);
370: if ( lu->options.symflag == SPOOLES_SYMMETRIC ) {
371: InpMtx_mapToUpperTriangle(lu->mtxA);
372: }
373: InpMtx_changeCoordType(lu->mtxA, INPMTX_BY_CHEVRONS);
374: }
375: InpMtx_changeStorageMode(lu->mtxA, INPMTX_BY_VECTORS);
376: if ( lu->options.msglvl > 2 ) {
377: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n input matrix after permutation");
378: InpMtx_writeForHumanEye(lu->mtxA, lu->options.msgFile);
379: }
380: } /* end of if( lu->flg == DIFFERENT_NONZERO_PATTERN) */
381:
382: if (lu->options.useQR){
383: FrontMtx_init(lu->frontmtx, lu->frontETree, lu->symbfacIVL, lu->options.typeflag,
384: SPOOLES_SYMMETRIC, FRONTMTX_DENSE_FRONTS,
385: SPOOLES_NO_PIVOTING, NO_LOCK, 0, NULL,
386: lu->mtxmanager, lu->options.msglvl, lu->options.msgFile);
387: } else {
388: FrontMtx_init(lu->frontmtx, lu->frontETree, lu->symbfacIVL, lu->options.typeflag, lu->options.symflag,
389: FRONTMTX_DENSE_FRONTS, lu->options.pivotingflag, NO_LOCK, 0, NULL,
390: lu->mtxmanager, lu->options.msglvl, lu->options.msgFile);
391: }
393: if ( lu->options.symflag == SPOOLES_SYMMETRIC ) { /* || SPOOLES_HERMITIAN ? */
394: if ( lu->options.patchAndGoFlag == 1 ) {
395: lu->frontmtx->patchinfo = PatchAndGoInfo_new();
396: PatchAndGoInfo_init(lu->frontmtx->patchinfo, 1, lu->options.toosmall, lu->options.fudge,
397: lu->options.storeids, lu->options.storevalues);
398: } else if ( lu->options.patchAndGoFlag == 2 ) {
399: lu->frontmtx->patchinfo = PatchAndGoInfo_new();
400: PatchAndGoInfo_init(lu->frontmtx->patchinfo, 2, lu->options.toosmall, lu->options.fudge,
401: lu->options.storeids, lu->options.storevalues);
402: }
403: }
405: /* numerical factorization */
406: chvmanager = ChvManager_new();
407: ChvManager_init(chvmanager, NO_LOCK, 1);
408: DVfill(10, lu->cpus, 0.0);
409: if (lu->options.useQR){
410: facops = 0.0 ;
411: FrontMtx_QR_factor(lu->frontmtx, lu->mtxA, chvmanager,
412: lu->cpus, &facops, lu->options.msglvl, lu->options.msgFile);
413: if ( lu->options.msglvl > 1 ) {
414: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n factor matrix");
415: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n facops = %9.2f", facops);
416: }
417: } else {
418: IVfill(20, lu->stats, 0);
419: rootchv = FrontMtx_factorInpMtx(lu->frontmtx, lu->mtxA, lu->options.tau, 0.0,
420: chvmanager, &fierr, lu->cpus,lu->stats,lu->options.msglvl,lu->options.msgFile);
421: if (rootchv) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"\n matrix found to be singular");
422: if (fierr >= 0) SETERRQ1(PETSC_ERR_LIB,"\n error encountered at front %D", fierr);
423:
424: if(lu->options.FrontMtxInfo){
425: PetscPrintf(PETSC_COMM_SELF,"\n %8d pivots, %8d pivot tests, %8d delayed rows and columns\n",lu->stats[0], lu->stats[1], lu->stats[2]);
426: cputotal = lu->cpus[8] ;
427: if ( cputotal > 0.0 ) {
428: PetscPrintf(PETSC_COMM_SELF,
429: "\n cpus cpus/totaltime"
430: "\n initialize fronts %8.3f %6.2f"
431: "\n load original entries %8.3f %6.2f"
432: "\n update fronts %8.3f %6.2f"
433: "\n assemble postponed data %8.3f %6.2f"
434: "\n factor fronts %8.3f %6.2f"
435: "\n extract postponed data %8.3f %6.2f"
436: "\n store factor entries %8.3f %6.2f"
437: "\n miscellaneous %8.3f %6.2f"
438: "\n total time %8.3f \n",
439: lu->cpus[0], 100.*lu->cpus[0]/cputotal,
440: lu->cpus[1], 100.*lu->cpus[1]/cputotal,
441: lu->cpus[2], 100.*lu->cpus[2]/cputotal,
442: lu->cpus[3], 100.*lu->cpus[3]/cputotal,
443: lu->cpus[4], 100.*lu->cpus[4]/cputotal,
444: lu->cpus[5], 100.*lu->cpus[5]/cputotal,
445: lu->cpus[6], 100.*lu->cpus[6]/cputotal,
446: lu->cpus[7], 100.*lu->cpus[7]/cputotal, cputotal);
447: }
448: }
449: }
450: ChvManager_free(chvmanager);
452: if ( lu->options.msglvl > 0 ) {
453: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n factor matrix");
454: FrontMtx_writeForHumanEye(lu->frontmtx, lu->options.msgFile);
455: fflush(lu->options.msgFile);
456: }
458: if ( lu->options.symflag == SPOOLES_SYMMETRIC ) { /* || SPOOLES_HERMITIAN ? */
459: if ( lu->options.patchAndGoFlag == 1 ) {
460: if ( lu->frontmtx->patchinfo->fudgeIV != NULL ) {
461: if (lu->options.msglvl > 0 ){
462: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n small pivots found at these locations");
463: IV_writeForHumanEye(lu->frontmtx->patchinfo->fudgeIV, lu->options.msgFile);
464: }
465: }
466: PatchAndGoInfo_free(lu->frontmtx->patchinfo);
467: } else if ( lu->options.patchAndGoFlag == 2 ) {
468: if (lu->options.msglvl > 0 ){
469: if ( lu->frontmtx->patchinfo->fudgeIV != NULL ) {
470: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n small pivots found at these locations");
471: IV_writeForHumanEye(lu->frontmtx->patchinfo->fudgeIV, lu->options.msgFile);
472: }
473: if ( lu->frontmtx->patchinfo->fudgeDV != NULL ) {
474: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n perturbations");
475: DV_writeForHumanEye(lu->frontmtx->patchinfo->fudgeDV, lu->options.msgFile);
476: }
477: }
478: PatchAndGoInfo_free(lu->frontmtx->patchinfo);
479: }
480: }
482: /* post-process the factorization */
483: FrontMtx_postProcess(lu->frontmtx, lu->options.msglvl, lu->options.msgFile);
484: if ( lu->options.msglvl > 2 ) {
485: PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n factor matrix after post-processing");
486: FrontMtx_writeForHumanEye(lu->frontmtx, lu->options.msgFile);
487: fflush(lu->options.msgFile);
488: }
490: lu->flg = SAME_NONZERO_PATTERN;
491: lu->CleanUpSpooles = PETSC_TRUE;
492: return(0);
493: }
498: PetscErrorCode MatConvert_SeqAIJ_SeqAIJSpooles(Mat A,MatType type,MatReuse reuse,Mat *newmat)
499: {
501: Mat B=*newmat;
502: Mat_Spooles *lu;
505: PetscNew(Mat_Spooles,&lu);
506: if (reuse == MAT_INITIAL_MATRIX) {
507: /* This routine is inherited, so we know the type is correct. */
508: MatDuplicate(A,MAT_COPY_VALUES,&B);
509: lu->MatDuplicate = B->ops->duplicate;
510: lu->MatCholeskyFactorSymbolic = B->ops->choleskyfactorsymbolic;
511: lu->MatLUFactorSymbolic = B->ops->lufactorsymbolic;
512: lu->MatView = B->ops->view;
513: lu->MatAssemblyEnd = B->ops->assemblyend;
514: lu->MatDestroy = B->ops->destroy;
515: } else {
516: lu->MatDuplicate = A->ops->duplicate;
517: lu->MatCholeskyFactorSymbolic = A->ops->choleskyfactorsymbolic;
518: lu->MatLUFactorSymbolic = A->ops->lufactorsymbolic;
519: lu->MatView = A->ops->view;
520: lu->MatAssemblyEnd = A->ops->assemblyend;
521: lu->MatDestroy = A->ops->destroy;
522: }
523: B->spptr = (void*)lu;
524: lu->basetype = MATSEQAIJ;
525: lu->useQR = PETSC_FALSE;
526: lu->CleanUpSpooles = PETSC_FALSE;
528: B->ops->duplicate = MatDuplicate_Spooles;
529: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqAIJSpooles;
530: B->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqAIJSpooles;
531: B->ops->view = MatView_SeqAIJSpooles;
532: B->ops->assemblyend = MatAssemblyEnd_SeqAIJSpooles;
533: B->ops->destroy = MatDestroy_SeqAIJSpooles;
535: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaijspooles_seqaij_C",
536: "MatConvert_Spooles_Base",MatConvert_Spooles_Base);
537: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqaijspooles_C",
538: "MatConvert_SeqAIJ_SeqAIJSpooles",MatConvert_SeqAIJ_SeqAIJSpooles);
539: /* PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJSPOOLES); */
540: PetscObjectChangeTypeName((PetscObject)B,type);
541: *newmat = B;
542: return(0);
543: }
548: PetscErrorCode MatDuplicate_Spooles(Mat A, MatDuplicateOption op, Mat *M) {
550: Mat_Spooles *lu=(Mat_Spooles *)A->spptr;
553: (*lu->MatDuplicate)(A,op,M);
554: PetscMemcpy((*M)->spptr,lu,sizeof(Mat_Spooles));
555: return(0);
556: }
558: /*MC
559: MATSEQAIJSPOOLES - MATSEQAIJSPOOLES = "seqaijspooles" - A matrix type providing direct solvers (LU or Cholesky) for sequential matrices
560: via the external package SPOOLES.
562: If SPOOLES is installed (see the manual for
563: instructions on how to declare the existence of external packages),
564: a matrix type can be constructed which invokes SPOOLES solvers.
565: After calling MatCreate(...,A), simply call MatSetType(A,MATSEQAIJSPOOLES).
567: This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is
568: supported for this matrix type. One can also call MatConvert for an inplace conversion to or from
569: the MATSEQAIJ type without data copy.
571: Options Database Keys:
572: + -mat_type seqaijspooles - sets the matrix type to "seqaijspooles" during a call to MatSetFromOptions()
573: . -mat_spooles_tau <tau> - upper bound on the magnitude of the largest element in L or U
574: . -mat_spooles_seed <seed> - random number seed used for ordering
575: . -mat_spooles_msglvl <msglvl> - message output level
576: . -mat_spooles_ordering <BestOfNDandMS,MMD,MS,ND> - ordering used
577: . -mat_spooles_maxdomainsize <n> - maximum subgraph size used by Spooles orderings
578: . -mat_spooles_maxzeros <n> - maximum number of zeros inside a supernode
579: . -mat_spooles_maxsize <n> - maximum size of a supernode
580: . -mat_spooles_FrontMtxInfo <true,fase> - print Spooles information about the computed factorization
581: . -mat_spooles_symmetryflag <0,1,2> - 0: SPOOLES_SYMMETRIC, 1: SPOOLES_HERMITIAN, 2: SPOOLES_NONSYMMETRIC
582: . -mat_spooles_patchAndGoFlag <0,1,2> - 0: no patch, 1: use PatchAndGo strategy 1, 2: use PatchAndGo strategy 2
583: . -mat_spooles_toosmall <dt> - drop tolerance for PatchAndGo strategy 1
584: . -mat_spooles_storeids <bool integer> - if nonzero, stores row and col numbers where patches were applied in an IV object
585: . -mat_spooles_fudge <delta> - fudge factor for rescaling diagonals with PatchAndGo strategy 2
586: - -mat_spooles_storevalues <bool integer> - if nonzero and PatchAndGo strategy 2 is used, store change in diagonal value in a DV object
588: Level: beginner
590: .seealso: PCLU
591: M*/
596: PetscErrorCode MatCreate_SeqAIJSpooles(Mat A)
597: {
601: MatSetType(A,MATSEQAIJ);
602: MatConvert_SeqAIJ_SeqAIJSpooles(A,MATSEQAIJSPOOLES,MAT_REUSE_MATRIX,&A);
603: return(0);
604: }
607: /*MC
608: MATAIJSPOOLES - MATAIJSPOOLES = "aijspooles" - A matrix type providing direct solvers (LU or Cholesky) for sequential and parellel matrices
609: via the external package SPOOLES.
611: If SPOOLES is installed (see the manual for
612: instructions on how to declare the existence of external packages),
613: a matrix type can be constructed which invokes SPOOLES solvers.
614: After calling MatCreate(...,A), simply call MatSetType(A,MATAIJSPOOLES).
615: This matrix type is supported for double precision real and complex.
617: This matrix inherits from MATAIJ. As a result, MatSeqAIJSetPreallocation and MatMPIAIJSetPreallocation are
618: supported for this matrix type. One can also call MatConvert for an inplace conversion to or from
619: the MATAIJ type without data copy.
621: Options Database Keys:
622: + -mat_type aijspooles - sets the matrix type to "aijspooles" during a call to MatSetFromOptions()
623: . -mat_spooles_tau <tau> - upper bound on the magnitude of the largest element in L or U
624: . -mat_spooles_seed <seed> - random number seed used for ordering
625: . -mat_spooles_msglvl <msglvl> - message output level
626: . -mat_spooles_ordering <BestOfNDandMS,MMD,MS,ND> - ordering used
627: . -mat_spooles_maxdomainsize <n> - maximum subgraph size used by Spooles orderings
628: . -mat_spooles_maxzeros <n> - maximum number of zeros inside a supernode
629: . -mat_spooles_maxsize <n> - maximum size of a supernode
630: . -mat_spooles_FrontMtxInfo <true,fase> - print Spooles information about the computed factorization
631: . -mat_spooles_symmetryflag <0,1,2> - 0: SPOOLES_SYMMETRIC, 1: SPOOLES_HERMITIAN, 2: SPOOLES_NONSYMMETRIC
632: . -mat_spooles_patchAndGoFlag <0,1,2> - 0: no patch, 1: use PatchAndGo strategy 1, 2: use PatchAndGo strategy 2
633: . -mat_spooles_toosmall <dt> - drop tolerance for PatchAndGo strategy 1
634: . -mat_spooles_storeids <bool integer> - if nonzero, stores row and col numbers where patches were applied in an IV object
635: . -mat_spooles_fudge <delta> - fudge factor for rescaling diagonals with PatchAndGo strategy 2
636: - -mat_spooles_storevalues <bool integer> - if nonzero and PatchAndGo strategy 2 is used, store change in diagonal value in a DV object
638: Level: beginner
640: .seealso: PCLU
641: M*/
645: PetscErrorCode MatCreate_AIJSpooles(Mat A)
646: {
648: PetscMPIInt size;
651: MPI_Comm_size(A->comm,&size);
652: if (size == 1) {
653: MatSetType(A,MATSEQAIJ);
654: MatConvert_SeqAIJ_SeqAIJSpooles(A,MATSEQAIJSPOOLES,MAT_REUSE_MATRIX,&A);
655: } else {
656: MatSetType(A,MATMPIAIJ);
657: MatConvert_MPIAIJ_MPIAIJSpooles(A,MATMPIAIJSPOOLES,MAT_REUSE_MATRIX,&A);
658: }
659: return(0);
660: }