Actual source code: gamg.c
petsc-3.4.5 2014-06-29
1: /*
2: GAMG geometric-algebric multigrid PC - Mark Adams 2011
3: */
4: #include "petsc-private/matimpl.h"
5: #include <../src/ksp/pc/impls/gamg/gamg.h> /*I "petscpc.h" I*/
6: #include <petsc-private/kspimpl.h>
7: #include <../src/ksp/pc/impls/bjacobi/bjacobi.h> /* Hack to access same_local_solves */
9: #if defined PETSC_GAMG_USE_LOG
10: PetscLogEvent petsc_gamg_setup_events[NUM_SET];
11: #endif
13: #if defined PETSC_USE_LOG
14: PetscLogEvent PC_GAMGGgraph_AGG;
15: PetscLogEvent PC_GAMGGgraph_GEO;
16: PetscLogEvent PC_GAMGCoarsen_AGG;
17: PetscLogEvent PC_GAMGCoarsen_GEO;
18: PetscLogEvent PC_GAMGProlongator_AGG;
19: PetscLogEvent PC_GAMGProlongator_GEO;
20: PetscLogEvent PC_GAMGOptprol_AGG;
21: PetscLogEvent PC_GAMGKKTProl_AGG;
22: #endif
24: #define GAMG_MAXLEVELS 30
26: /* #define GAMG_STAGES */
27: #if (defined PETSC_GAMG_USE_LOG && defined GAMG_STAGES)
28: static PetscLogStage gamg_stages[GAMG_MAXLEVELS];
29: #endif
31: static PetscFunctionList GAMGList = 0;
32: static PetscBool PCGAMGPackageInitialized;
34: /* ----------------------------------------------------------------------------- */
37: PetscErrorCode PCReset_GAMG(PC pc)
38: {
40: PC_MG *mg = (PC_MG*)pc->data;
41: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
44: if (pc_gamg->data) { /* this should not happen, cleaned up in SetUp */
45: PetscPrintf(PetscObjectComm((PetscObject)pc),"***[%d]%s this should not happen, cleaned up in SetUp\n",0,__FUNCT__);
46: PetscFree(pc_gamg->data);
47: }
48: pc_gamg->data = NULL; pc_gamg->data_sz = 0;
50: if (pc_gamg->orig_data) {
51: PetscFree(pc_gamg->orig_data);
52: }
53: return(0);
54: }
56: /* private 2x2 Mat Nest for Stokes */
57: typedef struct {
58: Mat A11,A21,A12,Amat;
59: IS prim_is,constr_is;
60: } GAMGKKTMat;
64: static PetscErrorCode GAMGKKTMatCreate(Mat A, PetscBool iskkt, GAMGKKTMat *out)
65: {
67: out->Amat = A;
68: if (iskkt) {
70: IS is_constraint, is_prime;
71: PetscInt nmin,nmax;
73: MatGetOwnershipRange(A, &nmin, &nmax);
74: MatFindZeroDiagonals(A, &is_constraint);
75: ISComplement(is_constraint, nmin, nmax, &is_prime);
77: out->prim_is = is_prime;
78: out->constr_is = is_constraint;
80: MatGetSubMatrix(A, is_prime, is_prime, MAT_INITIAL_MATRIX, &out->A11);
81: MatGetSubMatrix(A, is_prime, is_constraint, MAT_INITIAL_MATRIX, &out->A12);
82: MatGetSubMatrix(A, is_constraint, is_prime, MAT_INITIAL_MATRIX, &out->A21);
83: } else {
84: out->A11 = A;
85: out->A21 = NULL;
86: out->A12 = NULL;
87: out->prim_is = NULL;
88: out->constr_is = NULL;
89: }
90: return(0);
91: }
95: static PetscErrorCode GAMGKKTMatDestroy(GAMGKKTMat *mat)
96: {
100: if (mat->A11 && mat->A11 != mat->Amat) {
101: MatDestroy(&mat->A11);
102: }
103: MatDestroy(&mat->A21);
104: MatDestroy(&mat->A12);
106: ISDestroy(&mat->prim_is);
107: ISDestroy(&mat->constr_is);
108: return(0);
109: }
111: /* -------------------------------------------------------------------------- */
112: /*
113: createLevel: create coarse op with RAP. repartition and/or reduce number
114: of active processors.
116: Input Parameter:
117: . pc - parameters + side effect: coarse data in 'pc_gamg->data' and
118: 'pc_gamg->data_sz' are changed via repartitioning/reduction.
119: . Amat_fine - matrix on this fine (k) level
120: . cr_bs - coarse block size
121: . isLast -
122: . stokes -
123: In/Output Parameter:
124: . a_P_inout - prolongation operator to the next level (k-->k-1)
125: . a_nactive_proc - number of active procs
126: Output Parameter:
127: . a_Amat_crs - coarse matrix that is created (k-1)
128: */
132: static PetscErrorCode createLevel(const PC pc,const Mat Amat_fine,const PetscInt cr_bs,const PetscBool isLast,const PetscBool stokes,Mat *a_P_inout,Mat *a_Amat_crs,PetscMPIInt *a_nactive_proc)
133: {
134: PetscErrorCode ierr;
135: PC_MG *mg = (PC_MG*)pc->data;
136: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
137: const PetscBool repart = pc_gamg->repart;
138: const PetscInt min_eq_proc = pc_gamg->min_eq_proc, coarse_max = pc_gamg->coarse_eq_limit;
139: Mat Cmat,Pold=*a_P_inout;
140: MPI_Comm comm;
141: PetscMPIInt rank,size,new_size,nactive=*a_nactive_proc;
142: PetscInt ncrs_eq,ncrs_prim,f_bs;
145: PetscObjectGetComm((PetscObject)Amat_fine,&comm);
146: MPI_Comm_rank(comm, &rank);
147: MPI_Comm_size(comm, &size);
148: MatGetBlockSize(Amat_fine, &f_bs);
149: /* RAP */
150: MatPtAP(Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, &Cmat);
152: /* set 'ncrs_prim' (nodes), 'ncrs_eq' (equations)*/
153: ncrs_prim = pc_gamg->data_sz/pc_gamg->data_cell_cols/pc_gamg->data_cell_rows;
154: if (pc_gamg->data_sz % (pc_gamg->data_cell_cols*pc_gamg->data_cell_rows)) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"pc_gamg->data_sz %D not divisible by (pc_gamg->data_cell_cols %D *pc_gamg->data_cell_rows %D)",pc_gamg->data_sz,pc_gamg->data_cell_cols,pc_gamg->data_cell_rows);
155: MatGetLocalSize(Cmat, &ncrs_eq, NULL);
157: /* get number of PEs to make active 'new_size', reduce, can be any integer 1-P */
158: {
159: PetscInt ncrs_eq_glob;
160: MatGetSize(Cmat, &ncrs_eq_glob, NULL);
161: new_size = (PetscMPIInt)((float)ncrs_eq_glob/(float)min_eq_proc + 0.5); /* hardwire min. number of eq/proc */
162: if (new_size == 0 || ncrs_eq_glob < coarse_max) new_size = 1;
163: else if (new_size >= nactive) new_size = nactive; /* no change, rare */
164: if (isLast) new_size = 1;
165: }
167: if (!repart && new_size==nactive) *a_Amat_crs = Cmat; /* output - no repartitioning or reduction - could bail here */
168: else {
169: const PetscInt *idx,ndata_rows=pc_gamg->data_cell_rows,ndata_cols=pc_gamg->data_cell_cols,node_data_sz=ndata_rows*ndata_cols;
170: PetscInt *counts,*newproc_idx,ii,jj,kk,strideNew,*tidx,ncrs_prim_new,ncrs_eq_new,nloc_old;
171: IS is_eq_newproc,is_eq_newproc_prim,is_eq_num,is_eq_num_prim,isscat,new_eq_indices;
172: VecScatter vecscat;
173: PetscScalar *array;
174: Vec src_crd, dest_crd;
176: nloc_old = ncrs_eq/cr_bs;
177: if (ncrs_eq % cr_bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ncrs_eq %D not divisible by cr_bs %D",ncrs_eq,cr_bs);
178: #if defined PETSC_GAMG_USE_LOG
179: PetscLogEventBegin(petsc_gamg_setup_events[SET12],0,0,0,0);
180: #endif
181: /* make 'is_eq_newproc' */
182: PetscMalloc(size*sizeof(PetscInt), &counts);
183: if (repart && !stokes) {
184: /* Repartition Cmat_{k} and move colums of P^{k}_{k-1} and coordinates of primal part accordingly */
185: Mat adj;
187: if (pc_gamg->verbose>0) {
188: if (pc_gamg->verbose==1) PetscPrintf(comm,"\t[%d]%s repartition: size (active): %d --> %d, neq = %d\n",rank,__FUNCT__,*a_nactive_proc,new_size,ncrs_eq);
189: else {
190: PetscInt n;
191: MPI_Allreduce(&ncrs_eq, &n, 1, MPIU_INT, MPI_SUM, comm);
192: PetscPrintf(comm,"\t[%d]%s repartition: size (active): %d --> %d, neq = %d\n",rank,__FUNCT__,*a_nactive_proc,new_size,n);
193: }
194: }
196: /* get 'adj' */
197: if (cr_bs == 1) {
198: MatConvert(Cmat, MATMPIADJ, MAT_INITIAL_MATRIX, &adj);
199: } else {
200: /* make a scalar matrix to partition (no Stokes here) */
201: Mat tMat;
202: PetscInt Istart_crs,Iend_crs,ncols,jj,Ii;
203: const PetscScalar *vals;
204: const PetscInt *idx;
205: PetscInt *d_nnz, *o_nnz, M, N;
206: static PetscInt llev = 0;
208: PetscMalloc(ncrs_prim*sizeof(PetscInt), &d_nnz);
209: PetscMalloc(ncrs_prim*sizeof(PetscInt), &o_nnz);
210: MatGetOwnershipRange(Cmat, &Istart_crs, &Iend_crs);
211: MatGetSize(Cmat, &M, &N);
212: for (Ii = Istart_crs, jj = 0; Ii < Iend_crs; Ii += cr_bs, jj++) {
213: MatGetRow(Cmat,Ii,&ncols,0,0);
214: d_nnz[jj] = ncols/cr_bs;
215: o_nnz[jj] = ncols/cr_bs;
216: MatRestoreRow(Cmat,Ii,&ncols,0,0);
217: if (d_nnz[jj] > ncrs_prim) d_nnz[jj] = ncrs_prim;
218: if (o_nnz[jj] > (M/cr_bs-ncrs_prim)) o_nnz[jj] = M/cr_bs-ncrs_prim;
219: }
221: MatCreate(comm, &tMat);
222: MatSetSizes(tMat, ncrs_prim, ncrs_prim,PETSC_DETERMINE, PETSC_DETERMINE);
223: MatSetType(tMat,MATAIJ);
224: MatSeqAIJSetPreallocation(tMat,0,d_nnz);
225: MatMPIAIJSetPreallocation(tMat,0,d_nnz,0,o_nnz);
226: PetscFree(d_nnz);
227: PetscFree(o_nnz);
229: for (ii = Istart_crs; ii < Iend_crs; ii++) {
230: PetscInt dest_row = ii/cr_bs;
231: MatGetRow(Cmat,ii,&ncols,&idx,&vals);
232: for (jj = 0; jj < ncols; jj++) {
233: PetscInt dest_col = idx[jj]/cr_bs;
234: PetscScalar v = 1.0;
235: MatSetValues(tMat,1,&dest_row,1,&dest_col,&v,ADD_VALUES);
236: }
237: MatRestoreRow(Cmat,ii,&ncols,&idx,&vals);
238: }
239: MatAssemblyBegin(tMat,MAT_FINAL_ASSEMBLY);
240: MatAssemblyEnd(tMat,MAT_FINAL_ASSEMBLY);
242: if (llev++ == -1) {
243: PetscViewer viewer; char fname[32];
244: PetscSNPrintf(fname,sizeof(fname),"part_mat_%D.mat",llev);
245: PetscViewerBinaryOpen(comm,fname,FILE_MODE_WRITE,&viewer);
246: MatView(tMat, viewer);
247: PetscViewerDestroy(&viewer);
248: }
250: MatConvert(tMat, MATMPIADJ, MAT_INITIAL_MATRIX, &adj);
252: MatDestroy(&tMat);
253: } /* create 'adj' */
255: { /* partition: get newproc_idx */
256: char prefix[256];
257: const char *pcpre;
258: const PetscInt *is_idx;
259: MatPartitioning mpart;
260: IS proc_is;
261: PetscInt targetPE;
263: MatPartitioningCreate(comm, &mpart);
264: MatPartitioningSetAdjacency(mpart, adj);
265: PCGetOptionsPrefix(pc, &pcpre);
266: PetscSNPrintf(prefix,sizeof(prefix),"%spc_gamg_",pcpre ? pcpre : "");
267: PetscObjectSetOptionsPrefix((PetscObject)mpart,prefix);
268: MatPartitioningSetFromOptions(mpart);
269: MatPartitioningSetNParts(mpart, new_size);
270: MatPartitioningApply(mpart, &proc_is);
271: MatPartitioningDestroy(&mpart);
273: /* collect IS info */
274: PetscMalloc(ncrs_eq*sizeof(PetscInt), &newproc_idx);
275: ISGetIndices(proc_is, &is_idx);
276: targetPE = 1; /* bring to "front" of machine */
277: /*targetPE = size/new_size;*/ /* spread partitioning across machine */
278: for (kk = jj = 0 ; kk < nloc_old ; kk++) {
279: for (ii = 0 ; ii < cr_bs ; ii++, jj++) {
280: newproc_idx[jj] = is_idx[kk] * targetPE; /* distribution */
281: }
282: }
283: ISRestoreIndices(proc_is, &is_idx);
284: ISDestroy(&proc_is);
285: }
286: MatDestroy(&adj);
288: ISCreateGeneral(comm, ncrs_eq, newproc_idx, PETSC_COPY_VALUES, &is_eq_newproc);
289: if (newproc_idx != 0) {
290: PetscFree(newproc_idx);
291: }
292: } else { /* simple aggreagtion of parts -- 'is_eq_newproc' */
294: PetscInt rfactor,targetPE;
295: /* find factor */
296: if (new_size == 1) rfactor = size; /* easy */
297: else {
298: PetscReal best_fact = 0.;
299: jj = -1;
300: for (kk = 1 ; kk <= size ; kk++) {
301: if (size%kk==0) { /* a candidate */
302: PetscReal nactpe = (PetscReal)size/(PetscReal)kk, fact = nactpe/(PetscReal)new_size;
303: if (fact > 1.0) fact = 1./fact; /* keep fact < 1 */
304: if (fact > best_fact) {
305: best_fact = fact; jj = kk;
306: }
307: }
308: }
309: if (jj != -1) rfactor = jj;
310: else rfactor = 1; /* does this happen .. a prime */
311: }
312: new_size = size/rfactor;
314: if (new_size==nactive) {
315: *a_Amat_crs = Cmat; /* output - no repartitioning or reduction, bail out because nested here */
316: PetscFree(counts);
317: if (pc_gamg->verbose>0) {
318: PetscPrintf(comm,"\t[%d]%s aggregate processors noop: new_size=%d, neq(loc)=%d\n",rank,__FUNCT__,new_size,ncrs_eq);
319: }
320: return(0);
321: }
323: if (pc_gamg->verbose) PetscPrintf(comm,"\t[%d]%s number of equations (loc) %d with simple aggregation\n",rank,__FUNCT__,ncrs_eq);
324: targetPE = rank/rfactor;
325: ISCreateStride(comm, ncrs_eq, targetPE, 0, &is_eq_newproc);
327: if (stokes) {
328: ISCreateStride(comm, ncrs_prim*cr_bs, targetPE, 0, &is_eq_newproc_prim);
329: }
330: } /* end simple 'is_eq_newproc' */
332: /*
333: Create an index set from the is_eq_newproc index set to indicate the mapping TO
334: */
335: ISPartitioningToNumbering(is_eq_newproc, &is_eq_num);
336: if (stokes) {
337: ISPartitioningToNumbering(is_eq_newproc_prim, &is_eq_num_prim);
338: } else is_eq_num_prim = is_eq_num;
339: /*
340: Determine how many equations/vertices are assigned to each processor
341: */
342: ISPartitioningCount(is_eq_newproc, size, counts);
343: ncrs_eq_new = counts[rank];
344: ISDestroy(&is_eq_newproc);
345: if (stokes) {
346: ISPartitioningCount(is_eq_newproc_prim, size, counts);
347: ISDestroy(&is_eq_newproc_prim);
348: ncrs_prim_new = counts[rank]/cr_bs; /* nodes */
349: } else ncrs_prim_new = ncrs_eq_new/cr_bs; /* eqs */
351: PetscFree(counts);
352: #if defined PETSC_GAMG_USE_LOG
353: PetscLogEventEnd(petsc_gamg_setup_events[SET12],0,0,0,0);
354: #endif
356: /* move data (for primal equations only) */
357: /* Create a vector to contain the newly ordered element information */
358: VecCreate(comm, &dest_crd);
359: VecSetSizes(dest_crd, node_data_sz*ncrs_prim_new, PETSC_DECIDE);
360: VecSetFromOptions(dest_crd); /* this is needed! */
361: /*
362: There are 'ndata_rows*ndata_cols' data items per node, (one can think of the vectors of having
363: a block size of ...). Note, ISs are expanded into equation space by 'cr_bs'.
364: */
365: PetscMalloc((ncrs_prim*node_data_sz)*sizeof(PetscInt), &tidx);
366: ISGetIndices(is_eq_num_prim, &idx);
367: for (ii=0,jj=0; ii<ncrs_prim; ii++) {
368: PetscInt id = idx[ii*cr_bs]/cr_bs; /* get node back */
369: for (kk=0; kk<node_data_sz; kk++, jj++) tidx[jj] = id*node_data_sz + kk;
370: }
371: ISRestoreIndices(is_eq_num_prim, &idx);
372: ISCreateGeneral(comm, node_data_sz*ncrs_prim, tidx, PETSC_COPY_VALUES, &isscat);
373: PetscFree(tidx);
374: /*
375: Create a vector to contain the original vertex information for each element
376: */
377: VecCreateSeq(PETSC_COMM_SELF, node_data_sz*ncrs_prim, &src_crd);
378: for (jj=0; jj<ndata_cols; jj++) {
379: const PetscInt stride0=ncrs_prim*pc_gamg->data_cell_rows;
380: for (ii=0; ii<ncrs_prim; ii++) {
381: for (kk=0; kk<ndata_rows; kk++) {
382: PetscInt ix = ii*ndata_rows + kk + jj*stride0, jx = ii*node_data_sz + kk*ndata_cols + jj;
383: PetscScalar tt = (PetscScalar)pc_gamg->data[ix];
384: VecSetValues(src_crd, 1, &jx, &tt, INSERT_VALUES);
385: }
386: }
387: }
388: VecAssemblyBegin(src_crd);
389: VecAssemblyEnd(src_crd);
390: /*
391: Scatter the element vertex information (still in the original vertex ordering)
392: to the correct processor
393: */
394: VecScatterCreate(src_crd, NULL, dest_crd, isscat, &vecscat);
395: ISDestroy(&isscat);
396: VecScatterBegin(vecscat,src_crd,dest_crd,INSERT_VALUES,SCATTER_FORWARD);
397: VecScatterEnd(vecscat,src_crd,dest_crd,INSERT_VALUES,SCATTER_FORWARD);
398: VecScatterDestroy(&vecscat);
399: VecDestroy(&src_crd);
400: /*
401: Put the element vertex data into a new allocation of the gdata->ele
402: */
403: PetscFree(pc_gamg->data);
404: PetscMalloc(node_data_sz*ncrs_prim_new*sizeof(PetscReal), &pc_gamg->data);
406: pc_gamg->data_sz = node_data_sz*ncrs_prim_new;
407: strideNew = ncrs_prim_new*ndata_rows;
409: VecGetArray(dest_crd, &array);
410: for (jj=0; jj<ndata_cols; jj++) {
411: for (ii=0; ii<ncrs_prim_new; ii++) {
412: for (kk=0; kk<ndata_rows; kk++) {
413: PetscInt ix = ii*ndata_rows + kk + jj*strideNew, jx = ii*node_data_sz + kk*ndata_cols + jj;
414: pc_gamg->data[ix] = PetscRealPart(array[jx]);
415: }
416: }
417: }
418: VecRestoreArray(dest_crd, &array);
419: VecDestroy(&dest_crd);
421: /* move A and P (columns) with new layout */
422: #if defined PETSC_GAMG_USE_LOG
423: PetscLogEventBegin(petsc_gamg_setup_events[SET13],0,0,0,0);
424: #endif
426: /*
427: Invert for MatGetSubMatrix
428: */
429: ISInvertPermutation(is_eq_num, ncrs_eq_new, &new_eq_indices);
430: ISSort(new_eq_indices); /* is this needed? */
431: ISSetBlockSize(new_eq_indices, cr_bs);
432: if (is_eq_num != is_eq_num_prim) {
433: ISDestroy(&is_eq_num_prim); /* could be same as 'is_eq_num' */
434: }
435: ISDestroy(&is_eq_num);
436: #if defined PETSC_GAMG_USE_LOG
437: PetscLogEventEnd(petsc_gamg_setup_events[SET13],0,0,0,0);
438: PetscLogEventBegin(petsc_gamg_setup_events[SET14],0,0,0,0);
439: #endif
440: /* 'a_Amat_crs' output */
441: {
442: Mat mat;
443: MatGetSubMatrix(Cmat, new_eq_indices, new_eq_indices, MAT_INITIAL_MATRIX, &mat);
444: *a_Amat_crs = mat;
446: if (!PETSC_TRUE) {
447: PetscInt cbs, rbs;
448: MatGetBlockSizes(Cmat, &rbs, &cbs);
449: PetscPrintf(MPI_COMM_SELF,"[%d]%s Old Mat rbs=%d cbs=%d\n",rank,__FUNCT__,rbs,cbs);
450: MatGetBlockSizes(mat, &rbs, &cbs);
451: PetscPrintf(MPI_COMM_SELF,"[%d]%s New Mat rbs=%d cbs=%d cr_bs=%d\n",rank,__FUNCT__,rbs,cbs,cr_bs);
452: }
453: }
454: MatDestroy(&Cmat);
456: #if defined PETSC_GAMG_USE_LOG
457: PetscLogEventEnd(petsc_gamg_setup_events[SET14],0,0,0,0);
458: #endif
459: /* prolongator */
460: {
461: IS findices;
462: PetscInt Istart,Iend;
463: Mat Pnew;
464: MatGetOwnershipRange(Pold, &Istart, &Iend);
465: #if defined PETSC_GAMG_USE_LOG
466: PetscLogEventBegin(petsc_gamg_setup_events[SET15],0,0,0,0);
467: #endif
468: ISCreateStride(comm,Iend-Istart,Istart,1,&findices);
469: ISSetBlockSize(findices,f_bs);
470: MatGetSubMatrix(Pold, findices, new_eq_indices, MAT_INITIAL_MATRIX, &Pnew);
471: ISDestroy(&findices);
473: if (!PETSC_TRUE) {
474: PetscInt cbs, rbs;
475: MatGetBlockSizes(Pold, &rbs, &cbs);
476: PetscPrintf(MPI_COMM_SELF,"[%d]%s Pold rbs=%d cbs=%d\n",rank,__FUNCT__,rbs,cbs);
477: MatGetBlockSizes(Pnew, &rbs, &cbs);
478: PetscPrintf(MPI_COMM_SELF,"[%d]%s Pnew rbs=%d cbs=%d\n",rank,__FUNCT__,rbs,cbs);
479: }
480: #if defined PETSC_GAMG_USE_LOG
481: PetscLogEventEnd(petsc_gamg_setup_events[SET15],0,0,0,0);
482: #endif
483: MatDestroy(a_P_inout);
485: /* output - repartitioned */
486: *a_P_inout = Pnew;
487: }
488: ISDestroy(&new_eq_indices);
490: *a_nactive_proc = new_size; /* output */
491: }
493: /* outout matrix data */
494: if (!PETSC_TRUE) {
495: PetscViewer viewer; char fname[32]; static int llev=0; Cmat = *a_Amat_crs;
496: if (llev==0) {
497: sprintf(fname,"Cmat_%d.m",llev++);
498: PetscViewerASCIIOpen(comm,fname,&viewer);
499: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
500: MatView(Amat_fine, viewer);
501: PetscViewerDestroy(&viewer);
502: }
503: sprintf(fname,"Cmat_%d.m",llev++);
504: PetscViewerASCIIOpen(comm,fname,&viewer);
505: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
506: MatView(Cmat, viewer);
507: PetscViewerDestroy(&viewer);
508: }
509: return(0);
510: }
512: /* -------------------------------------------------------------------------- */
513: /*
514: PCSetUp_GAMG - Prepares for the use of the GAMG preconditioner
515: by setting data structures and options.
517: Input Parameter:
518: . pc - the preconditioner context
520: Application Interface Routine: PCSetUp()
522: Notes:
523: The interface routine PCSetUp() is not usually called directly by
524: the user, but instead is called by PCApply() if necessary.
525: */
528: PetscErrorCode PCSetUp_GAMG(PC pc)
529: {
531: PC_MG *mg = (PC_MG*)pc->data;
532: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
533: Mat Pmat = pc->pmat;
534: PetscInt fine_level,level,level1,bs,M,qq,lidx,nASMBlocksArr[GAMG_MAXLEVELS];
535: MPI_Comm comm;
536: PetscMPIInt rank,size,nactivepe;
537: Mat Aarr[GAMG_MAXLEVELS],Parr[GAMG_MAXLEVELS];
538: PetscReal emaxs[GAMG_MAXLEVELS];
539: IS *ASMLocalIDsArr[GAMG_MAXLEVELS];
540: GAMGKKTMat kktMatsArr[GAMG_MAXLEVELS];
541: PetscLogDouble nnz0=0.,nnztot=0.;
542: MatInfo info;
543: PetscBool stokes = PETSC_FALSE, redo_mesh_setup = (PetscBool)(!pc_gamg->reuse_prol);
546: PetscObjectGetComm((PetscObject)pc,&comm);
547: MPI_Comm_rank(comm,&rank);
548: MPI_Comm_size(comm,&size);
550: if (pc_gamg->verbose>2) PetscPrintf(comm,"[%d]%s pc_gamg->setup_count=%d pc->setupcalled=%d\n",rank,__FUNCT__,pc_gamg->setup_count,pc->setupcalled);
552: if (pc_gamg->setup_count++ > 0) {
553: if (redo_mesh_setup) {
554: /* reset everything */
555: PCReset_MG(pc);
556: pc->setupcalled = 0;
557: } else {
558: PC_MG_Levels **mglevels = mg->levels;
559: /* just do Galerkin grids */
560: Mat B,dA,dB;
562: if (!pc->setupcalled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCSetUp() has not been called yet");
563: if (pc_gamg->Nlevels > 1) {
564: /* currently only handle case where mat and pmat are the same on coarser levels */
565: KSPGetOperators(mglevels[pc_gamg->Nlevels-1]->smoothd,&dA,&dB,NULL);
566: /* (re)set to get dirty flag */
567: KSPSetOperators(mglevels[pc_gamg->Nlevels-1]->smoothd,dA,dB,SAME_NONZERO_PATTERN);
569: for (level=pc_gamg->Nlevels-2; level>=0; level--) {
570: /* the first time through the matrix structure has changed from repartitioning */
571: if (pc_gamg->setup_count==2) {
572: MatPtAP(dB,mglevels[level+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);
573: MatDestroy(&mglevels[level]->A);
575: mglevels[level]->A = B;
576: } else {
577: KSPGetOperators(mglevels[level]->smoothd,NULL,&B,NULL);
578: MatPtAP(dB,mglevels[level+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);
579: }
580: KSPSetOperators(mglevels[level]->smoothd,B,B,SAME_NONZERO_PATTERN);
581: dB = B;
582: }
583: }
585: PCSetUp_MG(pc);
587: /* PCSetUp_MG seems to insists on setting this to GMRES */
588: KSPSetType(mglevels[0]->smoothd, KSPPREONLY);
589: return(0);
590: }
591: }
593: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);
595: GAMGKKTMatCreate(Pmat, stokes, &kktMatsArr[0]);
597: if (!pc_gamg->data) {
598: if (pc_gamg->orig_data) {
599: MatGetBlockSize(Pmat, &bs);
600: MatGetLocalSize(Pmat, &qq, NULL);
602: pc_gamg->data_sz = (qq/bs)*pc_gamg->orig_data_cell_rows*pc_gamg->orig_data_cell_cols;
603: pc_gamg->data_cell_rows = pc_gamg->orig_data_cell_rows;
604: pc_gamg->data_cell_cols = pc_gamg->orig_data_cell_cols;
606: PetscMalloc(pc_gamg->data_sz*sizeof(PetscReal), &pc_gamg->data);
607: for (qq=0; qq<pc_gamg->data_sz; qq++) pc_gamg->data[qq] = pc_gamg->orig_data[qq];
608: } else {
609: if (!pc_gamg->ops->createdefaultdata) SETERRQ(comm,PETSC_ERR_PLIB,"'createdefaultdata' not set(?) need to support NULL data");
610: if (stokes) SETERRQ(comm,PETSC_ERR_PLIB,"Need data (eg, PCSetCoordinates) for Stokes problems");
611: pc_gamg->ops->createdefaultdata(pc, kktMatsArr[0].A11);
612: }
613: }
615: /* cache original data for reuse */
616: if (!pc_gamg->orig_data && redo_mesh_setup) {
617: PetscMalloc(pc_gamg->data_sz*sizeof(PetscReal), &pc_gamg->orig_data);
618: for (qq=0; qq<pc_gamg->data_sz; qq++) pc_gamg->orig_data[qq] = pc_gamg->data[qq];
619: pc_gamg->orig_data_cell_rows = pc_gamg->data_cell_rows;
620: pc_gamg->orig_data_cell_cols = pc_gamg->data_cell_cols;
621: }
623: /* get basic dims */
624: if (stokes) bs = pc_gamg->data_cell_rows; /* this is agg-mg specific */
625: else {
626: MatGetBlockSize(Pmat, &bs);
627: }
629: MatGetSize(Pmat, &M, &qq);
630: if (pc_gamg->verbose) {
631: PetscInt NN = M;
632: if (pc_gamg->verbose==1) {
633: MatGetInfo(Pmat,MAT_LOCAL,&info);
634: MatGetLocalSize(Pmat, &NN, &qq);
635: } else {
636: MatGetInfo(Pmat,MAT_GLOBAL_SUM,&info);
637: }
638: nnz0 = info.nz_used;
639: nnztot = info.nz_used;
640: PetscPrintf(comm,"\t[%d]%s level %d N=%d, n data rows=%d, n data cols=%d, nnz/row (ave)=%d, np=%d\n",
641: rank,__FUNCT__,0,M,pc_gamg->data_cell_rows,pc_gamg->data_cell_cols,
642: (int)(nnz0/(PetscReal)NN),size);
643: }
645: /* Get A_i and R_i */
646: for (level=0, Aarr[0]=Pmat, nactivepe = size; /* hard wired stopping logic */
647: level < (pc_gamg->Nlevels-1) && (level==0 || M>pc_gamg->coarse_eq_limit); /* && (size==1 || nactivepe>1); */
648: level++) {
649: level1 = level + 1;
650: #if defined PETSC_GAMG_USE_LOG
651: PetscLogEventBegin(petsc_gamg_setup_events[SET1],0,0,0,0);
652: #if (defined GAMG_STAGES)
653: PetscLogStagePush(gamg_stages[level]);
654: #endif
655: #endif
656: /* deal with Stokes, get sub matrices */
657: if (level > 0) {
658: GAMGKKTMatCreate(Aarr[level], stokes, &kktMatsArr[level]);
659: }
660: { /* construct prolongator */
661: Mat Gmat;
662: PetscCoarsenData *agg_lists;
663: Mat Prol11,Prol22;
665: pc_gamg->ops->graph(pc,kktMatsArr[level].A11, &Gmat);
666: pc_gamg->ops->coarsen(pc, &Gmat, &agg_lists);
667: pc_gamg->ops->prolongator(pc, kktMatsArr[level].A11, Gmat, agg_lists, &Prol11);
669: /* could have failed to create new level */
670: if (Prol11) {
671: /* get new block size of coarse matrices */
672: MatGetBlockSizes(Prol11, NULL, &bs);
674: if (stokes) {
675: if (!pc_gamg->ops->formkktprol) SETERRQ(comm,PETSC_ERR_USER,"Stokes not supportd by AMG method.");
676: /* R A12 == (T = A21 P)'; G = T' T; coarsen G; form plain agg with G */
677: pc_gamg->ops->formkktprol(pc, Prol11, kktMatsArr[level].A21, &Prol22);
678: }
680: if (pc_gamg->ops->optprol) {
681: /* smooth */
682: pc_gamg->ops->optprol(pc, kktMatsArr[level].A11, &Prol11);
683: }
685: if (stokes) {
686: IS is_row[2];
687: Mat a[4];
689: is_row[0] = kktMatsArr[level].prim_is; is_row[1] = kktMatsArr[level].constr_is;
690: a[0] = Prol11; a[1] = NULL; a[2] = NULL; a[3] = Prol22;
691: MatCreateNest(comm,2,is_row, 2, is_row, a, &Parr[level1]);
692: } else Parr[level1] = Prol11;
693: } else Parr[level1] = NULL;
695: if (pc_gamg->use_aggs_in_gasm) {
696: PetscCDGetASMBlocks(agg_lists, bs, &nASMBlocksArr[level], &ASMLocalIDsArr[level]);
697: }
699: MatDestroy(&Gmat);
700: PetscCDDestroy(agg_lists);
701: } /* construct prolongator scope */
702: #if defined PETSC_GAMG_USE_LOG
703: PetscLogEventEnd(petsc_gamg_setup_events[SET1],0,0,0,0);
704: #endif
705: /* cache eigen estimate */
706: if (pc_gamg->emax_id != -1) {
707: PetscBool flag;
708: PetscObjectComposedDataGetReal((PetscObject)kktMatsArr[level].A11, pc_gamg->emax_id, emaxs[level], flag);
709: if (!flag) emaxs[level] = -1.;
710: } else emaxs[level] = -1.;
711: if (level==0) Aarr[0] = Pmat; /* use Pmat for finest level setup */
712: if (!Parr[level1]) {
713: if (pc_gamg->verbose) {
714: PetscPrintf(comm,"\t[%d]%s stop gridding, level %d\n",rank,__FUNCT__,level);
715: }
716: break;
717: }
718: #if defined PETSC_GAMG_USE_LOG
719: PetscLogEventBegin(petsc_gamg_setup_events[SET2],0,0,0,0);
720: #endif
722: createLevel(pc, Aarr[level], bs, (PetscBool)(level==pc_gamg->Nlevels-2),
723: stokes, &Parr[level1], &Aarr[level1], &nactivepe);
725: #if defined PETSC_GAMG_USE_LOG
726: PetscLogEventEnd(petsc_gamg_setup_events[SET2],0,0,0,0);
727: #endif
728: MatGetSize(Aarr[level1], &M, &qq);
730: if (pc_gamg->verbose > 0) {
731: PetscInt NN = M;
732: if (pc_gamg->verbose==1) {
733: MatGetInfo(Aarr[level1],MAT_LOCAL,&info);
734: MatGetLocalSize(Aarr[level1], &NN, &qq);
735: } else {
736: MatGetInfo(Aarr[level1], MAT_GLOBAL_SUM, &info);
737: }
739: nnztot += info.nz_used;
740: PetscPrintf(comm,"\t\t[%d]%s %d) N=%d, n data cols=%d, nnz/row (ave)=%d, %d active pes\n",
741: rank,__FUNCT__,(int)level1,M,pc_gamg->data_cell_cols,
742: (int)(info.nz_used/(PetscReal)NN), nactivepe);
743: }
745: /* stop if one node -- could pull back for singular problems */
746: if (M/pc_gamg->data_cell_cols < 2) {
747: level++;
748: break;
749: }
750: #if (defined PETSC_GAMG_USE_LOG && defined GAMG_STAGES)
751: PetscLogStagePop();
752: #endif
753: } /* levels */
755: if (pc_gamg->data) {
756: PetscFree(pc_gamg->data);
757: pc_gamg->data = NULL;
758: }
760: if (pc_gamg->verbose) PetscPrintf(comm,"\t[%d]%s %d levels, grid complexity = %g\n",0,__FUNCT__,level+1,nnztot/nnz0);
761: pc_gamg->Nlevels = level + 1;
762: fine_level = level;
763: PCMGSetLevels(pc,pc_gamg->Nlevels,NULL);
765: /* simple setup */
766: if (!PETSC_TRUE) {
767: PC_MG_Levels **mglevels = mg->levels;
768: for (lidx=0,level=pc_gamg->Nlevels-1;
769: lidx<fine_level;
770: lidx++, level--) {
771: PCMGSetInterpolation(pc, lidx+1, Parr[level]);
772: KSPSetOperators(mglevels[lidx]->smoothd, Aarr[level], Aarr[level], SAME_NONZERO_PATTERN);
773: MatDestroy(&Parr[level]);
774: MatDestroy(&Aarr[level]);
775: }
776: KSPSetOperators(mglevels[fine_level]->smoothd, Aarr[0], Aarr[0], SAME_NONZERO_PATTERN);
778: PCSetUp_MG(pc);
779: } else if (pc_gamg->Nlevels > 1) { /* don't setup MG if one level */
780: /* set default smoothers & set operators */
781: for (lidx = 1, level = pc_gamg->Nlevels-2;
782: lidx <= fine_level;
783: lidx++, level--) {
784: KSP smoother;
785: PC subpc;
787: PCMGGetSmoother(pc, lidx, &smoother);
788: KSPGetPC(smoother, &subpc);
790: KSPSetNormType(smoother, KSP_NORM_NONE);
791: /* set ops */
792: KSPSetOperators(smoother, Aarr[level], Aarr[level], SAME_NONZERO_PATTERN);
793: PCMGSetInterpolation(pc, lidx, Parr[level+1]);
795: /* create field split PC, get subsmoother */
796: if (stokes) {
797: KSP *ksps;
798: PetscInt nn;
799: PCFieldSplitSetIS(subpc,"0",kktMatsArr[level].prim_is);
800: PCFieldSplitSetIS(subpc,"1",kktMatsArr[level].constr_is);
801: PCFieldSplitGetSubKSP(subpc,&nn,&ksps);
802: smoother = ksps[0];
803: KSPGetPC(smoother, &subpc);
804: PetscFree(ksps);
805: }
806: GAMGKKTMatDestroy(&kktMatsArr[level]);
808: /* set defaults */
809: KSPSetType(smoother, KSPCHEBYSHEV);
811: /* override defaults and command line args (!) */
812: if (pc_gamg->use_aggs_in_gasm) {
813: PetscInt sz;
814: IS *is;
816: sz = nASMBlocksArr[level];
817: is = ASMLocalIDsArr[level];
818: PCSetType(subpc, PCGASM);
819: if (sz==0) {
820: IS is;
821: PetscInt my0,kk;
822: MatGetOwnershipRange(Aarr[level], &my0, &kk);
823: ISCreateGeneral(PETSC_COMM_SELF, 1, &my0, PETSC_COPY_VALUES, &is);
824: PCGASMSetSubdomains(subpc, 1, &is, NULL);
825: ISDestroy(&is);
826: } else {
827: PetscInt kk;
828: PCGASMSetSubdomains(subpc, sz, is, NULL);
829: for (kk=0; kk<sz; kk++) {
830: ISDestroy(&is[kk]);
831: }
832: PetscFree(is);
833: }
834: PCGASMSetOverlap(subpc, 0);
836: ASMLocalIDsArr[level] = NULL;
837: nASMBlocksArr[level] = 0;
838: PCGASMSetType(subpc, PC_GASM_BASIC);
839: } else {
840: PCSetType(subpc, PCJACOBI);
841: }
842: }
843: {
844: /* coarse grid */
845: KSP smoother,*k2; PC subpc,pc2; PetscInt ii,first;
846: Mat Lmat = Aarr[(level=pc_gamg->Nlevels-1)]; lidx = 0;
847: PCMGGetSmoother(pc, lidx, &smoother);
848: KSPSetOperators(smoother, Lmat, Lmat, SAME_NONZERO_PATTERN);
849: KSPSetNormType(smoother, KSP_NORM_NONE);
850: KSPGetPC(smoother, &subpc);
851: PCSetType(subpc, PCBJACOBI);
852: PCSetUp(subpc);
853: PCBJacobiGetSubKSP(subpc,&ii,&first,&k2);
854: if (ii != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ii %D is not one",ii);
855: KSPGetPC(k2[0],&pc2);
856: PCSetType(pc2, PCLU);
857: PCFactorSetShiftType(pc2,MAT_SHIFT_INBLOCKS);
858: KSPSetTolerances(k2[0],PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,1);
859: /* This flag gets reset by PCBJacobiGetSubKSP(), but our BJacobi really does the same algorithm everywhere (and in
860: * fact, all but one process will have zero dofs), so we reset the flag to avoid having PCView_BJacobi attempt to
861: * view every subdomain as though they were different. */
862: ((PC_BJacobi*)subpc->data)->same_local_solves = PETSC_TRUE;
863: }
865: /* should be called in PCSetFromOptions_GAMG(), but cannot be called prior to PCMGSetLevels() */
866: PetscObjectOptionsBegin((PetscObject)pc);
867: PCSetFromOptions_MG(pc);
868: PetscOptionsEnd();
869: if (mg->galerkin != 2) SETERRQ(comm,PETSC_ERR_USER,"GAMG does Galerkin manually so the -pc_mg_galerkin option must not be used.");
871: /* create cheby smoothers */
872: for (lidx = 1, level = pc_gamg->Nlevels-2;
873: lidx <= fine_level;
874: lidx++, level--) {
875: KSP smoother;
876: PetscBool flag;
877: PC subpc;
879: PCMGGetSmoother(pc, lidx, &smoother);
880: KSPGetPC(smoother, &subpc);
882: /* create field split PC, get subsmoother */
883: if (stokes) {
884: KSP *ksps;
885: PetscInt nn;
886: PCFieldSplitGetSubKSP(subpc,&nn,&ksps);
887: smoother = ksps[0];
888: KSPGetPC(smoother, &subpc);
889: PetscFree(ksps);
890: }
892: /* do my own cheby */
893: PetscObjectTypeCompare((PetscObject)smoother, KSPCHEBYSHEV, &flag);
894: if (flag) {
895: PetscReal emax, emin;
896: PetscObjectTypeCompare((PetscObject)subpc, PCJACOBI, &flag);
897: if (flag && emaxs[level] > 0.0) emax=emaxs[level]; /* eigen estimate only for diagnal PC */
898: else { /* eigen estimate 'emax' */
899: KSP eksp;
900: Mat Lmat = Aarr[level];
901: Vec bb, xx;
903: MatGetVecs(Lmat, &bb, 0);
904: MatGetVecs(Lmat, &xx, 0);
905: {
906: PetscRandom rctx;
907: PetscRandomCreate(comm,&rctx);
908: PetscRandomSetFromOptions(rctx);
909: VecSetRandom(bb,rctx);
910: PetscRandomDestroy(&rctx);
911: }
913: /* zeroing out BC rows -- needed for crazy matrices */
914: {
915: PetscInt Istart,Iend,ncols,jj,Ii;
916: PetscScalar zero = 0.0;
917: MatGetOwnershipRange(Lmat, &Istart, &Iend);
918: for (Ii = Istart, jj = 0; Ii < Iend; Ii++, jj++) {
919: MatGetRow(Lmat,Ii,&ncols,0,0);
920: if (ncols <= 1) {
921: VecSetValues(bb, 1, &Ii, &zero, INSERT_VALUES);
922: }
923: MatRestoreRow(Lmat,Ii,&ncols,0,0);
924: }
925: VecAssemblyBegin(bb);
926: VecAssemblyEnd(bb);
927: }
929: KSPCreate(comm, &eksp);
930: KSPSetTolerances(eksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 10);
931: KSPSetNormType(eksp, KSP_NORM_NONE);
932: KSPSetOptionsPrefix(eksp,((PetscObject)pc)->prefix);
933: KSPAppendOptionsPrefix(eksp, "gamg_est_");
934: KSPSetFromOptions(eksp);
936: KSPSetInitialGuessNonzero(eksp, PETSC_FALSE);
937: KSPSetOperators(eksp, Lmat, Lmat, SAME_NONZERO_PATTERN);
938: KSPSetComputeSingularValues(eksp,PETSC_TRUE);
940: /* set PC type to be same as smoother */
941: KSPSetPC(eksp, subpc);
943: /* solve - keep stuff out of logging */
944: PetscLogEventDeactivate(KSP_Solve);
945: PetscLogEventDeactivate(PC_Apply);
946: KSPSolve(eksp, bb, xx);
947: PetscLogEventActivate(KSP_Solve);
948: PetscLogEventActivate(PC_Apply);
950: KSPComputeExtremeSingularValues(eksp, &emax, &emin);
952: VecDestroy(&xx);
953: VecDestroy(&bb);
954: KSPDestroy(&eksp);
956: if (pc_gamg->verbose > 0) {
957: PetscInt N1, tt;
958: MatGetSize(Aarr[level], &N1, &tt);
959: PetscPrintf(comm,"\t\t\t%s PC setup max eigen=%e min=%e on level %d (N=%d)\n",__FUNCT__,emax,emin,lidx,N1);
960: }
961: }
962: {
963: PetscInt N1, N0;
964: MatGetSize(Aarr[level], &N1, NULL);
965: MatGetSize(Aarr[level+1], &N0, NULL);
966: /* heuristic - is this crap? */
967: /* emin = 1.*emax/((PetscReal)N1/(PetscReal)N0); */
968: emin = emax * pc_gamg->eigtarget[0];
969: emax *= pc_gamg->eigtarget[1];
970: }
971: KSPChebyshevSetEigenvalues(smoother, emax, emin);
972: } /* setup checby flag */
973: } /* non-coarse levels */
975: /* clean up */
976: for (level=1; level<pc_gamg->Nlevels; level++) {
977: MatDestroy(&Parr[level]);
978: MatDestroy(&Aarr[level]);
979: }
981: PCSetUp_MG(pc);
983: if (PETSC_TRUE) {
984: KSP smoother; /* PCSetUp_MG seems to insists on setting this to GMRES on coarse grid */
985: PCMGGetSmoother(pc, 0, &smoother);
986: KSPSetType(smoother, KSPPREONLY);
987: }
988: } else {
989: KSP smoother;
990: if (pc_gamg->verbose) PetscPrintf(comm,"\t[%d]%s one level solver used (system is seen as DD). Using default solver.\n",rank,__FUNCT__);
991: PCMGGetSmoother(pc, 0, &smoother);
992: KSPSetOperators(smoother, Aarr[0], Aarr[0], SAME_NONZERO_PATTERN);
993: KSPSetType(smoother, KSPPREONLY);
994: PCSetUp_MG(pc);
995: }
996: return(0);
997: }
999: /* ------------------------------------------------------------------------- */
1000: /*
1001: PCDestroy_GAMG - Destroys the private context for the GAMG preconditioner
1002: that was created with PCCreate_GAMG().
1004: Input Parameter:
1005: . pc - the preconditioner context
1007: Application Interface Routine: PCDestroy()
1008: */
1011: PetscErrorCode PCDestroy_GAMG(PC pc)
1012: {
1014: PC_MG *mg = (PC_MG*)pc->data;
1015: PC_GAMG *pc_gamg= (PC_GAMG*)mg->innerctx;
1018: PCReset_GAMG(pc);
1019: if (pc_gamg->ops->destroy) {
1020: (*pc_gamg->ops->destroy)(pc);
1021: }
1022: PetscFree(pc_gamg->ops);
1023: PetscFree(pc_gamg->gamg_type_name);
1024: PetscFree(pc_gamg);
1025: PCDestroy_MG(pc);
1026: return(0);
1027: }
1032: /*@
1033: PCGAMGSetProcEqLim - Set number of equations to aim for on coarse grids via
1034: processor reduction.
1036: Not Collective on PC
1038: Input Parameters:
1039: . pc - the preconditioner context
1042: Options Database Key:
1043: . -pc_gamg_process_eq_limit
1045: Level: intermediate
1047: Concepts: Unstructured multrigrid preconditioner
1049: .seealso: ()
1050: @*/
1051: PetscErrorCode PCGAMGSetProcEqLim(PC pc, PetscInt n)
1052: {
1057: PetscTryMethod(pc,"PCGAMGSetProcEqLim_C",(PC,PetscInt),(pc,n));
1058: return(0);
1059: }
1063: static PetscErrorCode PCGAMGSetProcEqLim_GAMG(PC pc, PetscInt n)
1064: {
1065: PC_MG *mg = (PC_MG*)pc->data;
1066: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1069: if (n>0) pc_gamg->min_eq_proc = n;
1070: return(0);
1071: }
1075: /*@
1076: PCGAMGSetCoarseEqLim - Set max number of equations on coarse grids.
1078: Collective on PC
1080: Input Parameters:
1081: . pc - the preconditioner context
1084: Options Database Key:
1085: . -pc_gamg_coarse_eq_limit
1087: Level: intermediate
1089: Concepts: Unstructured multrigrid preconditioner
1091: .seealso: ()
1092: @*/
1093: PetscErrorCode PCGAMGSetCoarseEqLim(PC pc, PetscInt n)
1094: {
1099: PetscTryMethod(pc,"PCGAMGSetCoarseEqLim_C",(PC,PetscInt),(pc,n));
1100: return(0);
1101: }
1105: static PetscErrorCode PCGAMGSetCoarseEqLim_GAMG(PC pc, PetscInt n)
1106: {
1107: PC_MG *mg = (PC_MG*)pc->data;
1108: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1111: if (n>0) pc_gamg->coarse_eq_limit = n;
1112: return(0);
1113: }
1117: /*@
1118: PCGAMGSetRepartitioning - Repartition the coarse grids
1120: Collective on PC
1122: Input Parameters:
1123: . pc - the preconditioner context
1126: Options Database Key:
1127: . -pc_gamg_repartition
1129: Level: intermediate
1131: Concepts: Unstructured multrigrid preconditioner
1133: .seealso: ()
1134: @*/
1135: PetscErrorCode PCGAMGSetRepartitioning(PC pc, PetscBool n)
1136: {
1141: PetscTryMethod(pc,"PCGAMGSetRepartitioning_C",(PC,PetscBool),(pc,n));
1142: return(0);
1143: }
1147: static PetscErrorCode PCGAMGSetRepartitioning_GAMG(PC pc, PetscBool n)
1148: {
1149: PC_MG *mg = (PC_MG*)pc->data;
1150: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1153: pc_gamg->repart = n;
1154: return(0);
1155: }
1159: /*@
1160: PCGAMGSetReuseProl - Reuse prlongation
1162: Collective on PC
1164: Input Parameters:
1165: . pc - the preconditioner context
1168: Options Database Key:
1169: . -pc_gamg_reuse_interpolation
1171: Level: intermediate
1173: Concepts: Unstructured multrigrid preconditioner
1175: .seealso: ()
1176: @*/
1177: PetscErrorCode PCGAMGSetReuseProl(PC pc, PetscBool n)
1178: {
1183: PetscTryMethod(pc,"PCGAMGSetReuseProl_C",(PC,PetscBool),(pc,n));
1184: return(0);
1185: }
1189: static PetscErrorCode PCGAMGSetReuseProl_GAMG(PC pc, PetscBool n)
1190: {
1191: PC_MG *mg = (PC_MG*)pc->data;
1192: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1195: pc_gamg->reuse_prol = n;
1196: return(0);
1197: }
1201: /*@
1202: PCGAMGSetUseASMAggs -
1204: Collective on PC
1206: Input Parameters:
1207: . pc - the preconditioner context
1210: Options Database Key:
1211: . -pc_gamg_use_agg_gasm
1213: Level: intermediate
1215: Concepts: Unstructured multrigrid preconditioner
1217: .seealso: ()
1218: @*/
1219: PetscErrorCode PCGAMGSetUseASMAggs(PC pc, PetscBool n)
1220: {
1225: PetscTryMethod(pc,"PCGAMGSetUseASMAggs_C",(PC,PetscBool),(pc,n));
1226: return(0);
1227: }
1231: static PetscErrorCode PCGAMGSetUseASMAggs_GAMG(PC pc, PetscBool n)
1232: {
1233: PC_MG *mg = (PC_MG*)pc->data;
1234: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1237: pc_gamg->use_aggs_in_gasm = n;
1238: return(0);
1239: }
1243: /*@
1244: PCGAMGSetNlevels -
1246: Not collective on PC
1248: Input Parameters:
1249: . pc - the preconditioner context
1252: Options Database Key:
1253: . -pc_mg_levels
1255: Level: intermediate
1257: Concepts: Unstructured multrigrid preconditioner
1259: .seealso: ()
1260: @*/
1261: PetscErrorCode PCGAMGSetNlevels(PC pc, PetscInt n)
1262: {
1267: PetscTryMethod(pc,"PCGAMGSetNlevels_C",(PC,PetscInt),(pc,n));
1268: return(0);
1269: }
1273: static PetscErrorCode PCGAMGSetNlevels_GAMG(PC pc, PetscInt n)
1274: {
1275: PC_MG *mg = (PC_MG*)pc->data;
1276: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1279: pc_gamg->Nlevels = n;
1280: return(0);
1281: }
1285: /*@
1286: PCGAMGSetThreshold - Relative threshold to use for dropping edges in aggregation graph
1288: Not collective on PC
1290: Input Parameters:
1291: . pc - the preconditioner context
1294: Options Database Key:
1295: . -pc_gamg_threshold
1297: Level: intermediate
1299: Concepts: Unstructured multrigrid preconditioner
1301: .seealso: ()
1302: @*/
1303: PetscErrorCode PCGAMGSetThreshold(PC pc, PetscReal n)
1304: {
1309: PetscTryMethod(pc,"PCGAMGSetThreshold_C",(PC,PetscReal),(pc,n));
1310: return(0);
1311: }
1315: static PetscErrorCode PCGAMGSetThreshold_GAMG(PC pc, PetscReal n)
1316: {
1317: PC_MG *mg = (PC_MG*)pc->data;
1318: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1321: pc_gamg->threshold = n;
1322: return(0);
1323: }
1327: /*@
1328: PCGAMGSetType - Set solution method - calls sub create method
1330: Collective on PC
1332: Input Parameters:
1333: . pc - the preconditioner context
1336: Options Database Key:
1337: . -pc_gamg_type
1339: Level: intermediate
1341: Concepts: Unstructured multrigrid preconditioner
1343: .seealso: ()
1344: @*/
1345: PetscErrorCode PCGAMGSetType(PC pc, PCGAMGType type)
1346: {
1351: PetscTryMethod(pc,"PCGAMGSetType_C",(PC,PCGAMGType),(pc,type));
1352: return(0);
1353: }
1357: static PetscErrorCode PCGAMGSetType_GAMG(PC pc, PCGAMGType type)
1358: {
1359: PetscErrorCode ierr,(*r)(PC);
1360: PC_MG *mg = (PC_MG*)pc->data;
1361: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1364: PetscFunctionListFind(GAMGList,type,&r);
1365: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown GAMG type %s given",type);
1366: if (pc_gamg->ops->destroy) {
1367: (*pc_gamg->ops->destroy)(pc);
1368: PetscMemzero(pc_gamg->ops,sizeof(struct _PCGAMGOps));
1369: }
1370: PetscFree(pc_gamg->gamg_type_name);
1371: PetscStrallocpy(type,&pc_gamg->gamg_type_name);
1372: (*r)(pc);
1373: return(0);
1374: }
1378: PetscErrorCode PCSetFromOptions_GAMG(PC pc)
1379: {
1381: PC_MG *mg = (PC_MG*)pc->data;
1382: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1383: PetscBool flag;
1384: PetscInt two = 2;
1385: MPI_Comm comm;
1388: PetscObjectGetComm((PetscObject)pc,&comm);
1389: PetscOptionsHead("GAMG options");
1390: {
1391: /* -pc_gamg_type */
1392: {
1393: char tname[256] = PCGAMGAGG;
1394: const char *deftype = pc_gamg->gamg_type_name ? pc_gamg->gamg_type_name : tname;
1395: PetscOptionsList("-pc_gamg_type","Type of AMG method","PCGAMGSetType",GAMGList, tname, tname, sizeof(tname), &flag);
1396: /* call PCCreateGAMG_XYZ */
1397: if (flag || !pc_gamg->gamg_type_name) {
1398: PCGAMGSetType(pc, flag ? tname : deftype);
1399: }
1400: }
1401: /* -pc_gamg_verbose */
1402: PetscOptionsInt("-pc_gamg_verbose","Verbose (debugging) output for PCGAMG",
1403: "none", pc_gamg->verbose,
1404: &pc_gamg->verbose, NULL);
1405: /* -pc_gamg_repartition */
1406: PetscOptionsBool("-pc_gamg_repartition",
1407: "Repartion coarse grids (false)",
1408: "PCGAMGRepartitioning",
1409: pc_gamg->repart,
1410: &pc_gamg->repart,
1411: &flag);
1412: /* -pc_gamg_reuse_interpolation */
1413: PetscOptionsBool("-pc_gamg_reuse_interpolation",
1414: "Reuse prolongation operator (true)",
1415: "PCGAMGReuseProl",
1416: pc_gamg->reuse_prol,
1417: &pc_gamg->reuse_prol,
1418: &flag);
1419: /* -pc_gamg_use_agg_gasm */
1420: PetscOptionsBool("-pc_gamg_use_agg_gasm",
1421: "Use aggregation agragates for GASM smoother (false)",
1422: "PCGAMGUseASMAggs",
1423: pc_gamg->use_aggs_in_gasm,
1424: &pc_gamg->use_aggs_in_gasm,
1425: &flag);
1426: /* -pc_gamg_process_eq_limit */
1427: PetscOptionsInt("-pc_gamg_process_eq_limit",
1428: "Limit (goal) on number of equations per process on coarse grids",
1429: "PCGAMGSetProcEqLim",
1430: pc_gamg->min_eq_proc,
1431: &pc_gamg->min_eq_proc,
1432: &flag);
1433: /* -pc_gamg_coarse_eq_limit */
1434: PetscOptionsInt("-pc_gamg_coarse_eq_limit",
1435: "Limit on number of equations for the coarse grid",
1436: "PCGAMGSetCoarseEqLim",
1437: pc_gamg->coarse_eq_limit,
1438: &pc_gamg->coarse_eq_limit,
1439: &flag);
1440: /* -pc_gamg_threshold */
1441: PetscOptionsReal("-pc_gamg_threshold",
1442: "Relative threshold to use for dropping edges in aggregation graph",
1443: "PCGAMGSetThreshold",
1444: pc_gamg->threshold,
1445: &pc_gamg->threshold,
1446: &flag);
1447: if (flag && pc_gamg->verbose) {
1448: PetscPrintf(comm,"\t[%d]%s threshold set %e\n",0,__FUNCT__,pc_gamg->threshold);
1449: }
1450: /* -pc_gamg_eigtarget */
1451: PetscOptionsRealArray("-pc_gamg_eigtarget","Target eigenvalue range as fraction of estimated maximum eigenvalue","PCGAMGSetEigTarget",pc_gamg->eigtarget,&two,NULL);
1452: PetscOptionsInt("-pc_mg_levels",
1453: "Set number of MG levels",
1454: "PCGAMGSetNlevels",
1455: pc_gamg->Nlevels,
1456: &pc_gamg->Nlevels,
1457: &flag);
1459: /* set options for subtype */
1460: if (pc_gamg->ops->setfromoptions) {(*pc_gamg->ops->setfromoptions)(pc);}
1461: }
1462: PetscOptionsTail();
1463: return(0);
1464: }
1466: /* -------------------------------------------------------------------------- */
1467: /*MC
1468: PCGAMG - Geometric algebraic multigrid (AMG) preconditioning framework.
1469: - This is the entry point to GAMG, registered in pcregis.c
1471: Options Database Keys:
1472: Multigrid options(inherited)
1473: + -pc_mg_cycles <1>: 1 for V cycle, 2 for W-cycle (PCMGSetCycleType)
1474: . -pc_mg_smoothup <1>: Number of post-smoothing steps (PCMGSetNumberSmoothUp)
1475: . -pc_mg_smoothdown <1>: Number of pre-smoothing steps (PCMGSetNumberSmoothDown)
1476: - -pc_mg_type <multiplicative>: (one of) additive multiplicative full kascade
1478: Level: intermediate
1480: Concepts: multigrid
1482: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType,
1483: PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(), PCMGSetNumberSmoothDown(),
1484: PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(),
1485: PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(),
1486: PCMGSetCyclesOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR()
1487: M*/
1491: PETSC_EXTERN PetscErrorCode PCCreate_GAMG(PC pc)
1492: {
1494: PC_GAMG *pc_gamg;
1495: PC_MG *mg;
1496: #if defined PETSC_GAMG_USE_LOG
1497: static long count = 0;
1498: #endif
1501: /* PCGAMG is an inherited class of PCMG. Initialize pc as PCMG */
1502: PCSetType(pc, PCMG); /* calls PCCreate_MG() and MGCreate_Private() */
1503: PetscObjectChangeTypeName((PetscObject)pc, PCGAMG);
1505: /* create a supporting struct and attach it to pc */
1506: PetscNewLog(pc, PC_GAMG, &pc_gamg);
1507: mg = (PC_MG*)pc->data;
1508: mg->galerkin = 2; /* Use Galerkin, but it is computed externally */
1509: mg->innerctx = pc_gamg;
1511: PetscNewLog(pc,struct _PCGAMGOps,&pc_gamg->ops);
1513: pc_gamg->setup_count = 0;
1514: /* these should be in subctx but repartitioning needs simple arrays */
1515: pc_gamg->data_sz = 0;
1516: pc_gamg->data = 0;
1518: /* register AMG type */
1519: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
1520: PCGAMGInitializePackage();
1521: #endif
1523: /* overwrite the pointers of PCMG by the functions of base class PCGAMG */
1524: pc->ops->setfromoptions = PCSetFromOptions_GAMG;
1525: pc->ops->setup = PCSetUp_GAMG;
1526: pc->ops->reset = PCReset_GAMG;
1527: pc->ops->destroy = PCDestroy_GAMG;
1529: PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetProcEqLim_C",PCGAMGSetProcEqLim_GAMG);
1530: PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetCoarseEqLim_C",PCGAMGSetCoarseEqLim_GAMG);
1531: PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetRepartitioning_C",PCGAMGSetRepartitioning_GAMG);
1532: PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetReuseProl_C",PCGAMGSetReuseProl_GAMG);
1533: PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetUseASMAggs_C",PCGAMGSetUseASMAggs_GAMG);
1534: PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetThreshold_C",PCGAMGSetThreshold_GAMG);
1535: PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetType_C",PCGAMGSetType_GAMG);
1536: PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetNlevels_C",PCGAMGSetNlevels_GAMG);
1537: pc_gamg->repart = PETSC_FALSE;
1538: pc_gamg->reuse_prol = PETSC_FALSE;
1539: pc_gamg->use_aggs_in_gasm = PETSC_FALSE;
1540: pc_gamg->min_eq_proc = 50;
1541: pc_gamg->coarse_eq_limit = 800;
1542: pc_gamg->threshold = 0.;
1543: pc_gamg->Nlevels = GAMG_MAXLEVELS;
1544: pc_gamg->verbose = 0;
1545: pc_gamg->emax_id = -1;
1546: pc_gamg->eigtarget[0] = 0.05;
1547: pc_gamg->eigtarget[1] = 1.05;
1549: /* private events */
1550: #if defined PETSC_GAMG_USE_LOG
1551: if (count++ == 0) {
1552: PetscLogEventRegister("GAMG: createProl", PC_CLASSID, &petsc_gamg_setup_events[SET1]);
1553: PetscLogEventRegister(" Graph", PC_CLASSID, &petsc_gamg_setup_events[GRAPH]);
1554: /* PetscLogEventRegister(" G.Mat", PC_CLASSID, &petsc_gamg_setup_events[GRAPH_MAT]); */
1555: /* PetscLogEventRegister(" G.Filter", PC_CLASSID, &petsc_gamg_setup_events[GRAPH_FILTER]); */
1556: /* PetscLogEventRegister(" G.Square", PC_CLASSID, &petsc_gamg_setup_events[GRAPH_SQR]); */
1557: PetscLogEventRegister(" MIS/Agg", PC_CLASSID, &petsc_gamg_setup_events[SET4]);
1558: PetscLogEventRegister(" geo: growSupp", PC_CLASSID, &petsc_gamg_setup_events[SET5]);
1559: PetscLogEventRegister(" geo: triangle", PC_CLASSID, &petsc_gamg_setup_events[SET6]);
1560: PetscLogEventRegister(" search&set", PC_CLASSID, &petsc_gamg_setup_events[FIND_V]);
1561: PetscLogEventRegister(" SA: col data", PC_CLASSID, &petsc_gamg_setup_events[SET7]);
1562: PetscLogEventRegister(" SA: frmProl0", PC_CLASSID, &petsc_gamg_setup_events[SET8]);
1563: PetscLogEventRegister(" SA: smooth", PC_CLASSID, &petsc_gamg_setup_events[SET9]);
1564: PetscLogEventRegister("GAMG: partLevel", PC_CLASSID, &petsc_gamg_setup_events[SET2]);
1565: PetscLogEventRegister(" repartition", PC_CLASSID, &petsc_gamg_setup_events[SET12]);
1566: PetscLogEventRegister(" Invert-Sort", PC_CLASSID, &petsc_gamg_setup_events[SET13]);
1567: PetscLogEventRegister(" Move A", PC_CLASSID, &petsc_gamg_setup_events[SET14]);
1568: PetscLogEventRegister(" Move P", PC_CLASSID, &petsc_gamg_setup_events[SET15]);
1570: /* PetscLogEventRegister(" PL move data", PC_CLASSID, &petsc_gamg_setup_events[SET13]); */
1571: /* PetscLogEventRegister("GAMG: fix", PC_CLASSID, &petsc_gamg_setup_events[SET10]); */
1572: /* PetscLogEventRegister("GAMG: set levels", PC_CLASSID, &petsc_gamg_setup_events[SET11]); */
1573: /* create timer stages */
1574: #if defined GAMG_STAGES
1575: {
1576: char str[32];
1577: PetscInt lidx;
1578: sprintf(str,"MG Level %d (finest)",0);
1579: PetscLogStageRegister(str, &gamg_stages[0]);
1580: for (lidx=1; lidx<9; lidx++) {
1581: sprintf(str,"MG Level %d",lidx);
1582: PetscLogStageRegister(str, &gamg_stages[lidx]);
1583: }
1584: }
1585: #endif
1586: }
1587: #endif
1588: /* general events */
1589: #if defined PETSC_USE_LOG
1590: PetscLogEventRegister("PCGAMGgraph_AGG", 0, &PC_GAMGGgraph_AGG);
1591: PetscLogEventRegister("PCGAMGgraph_GEO", PC_CLASSID, &PC_GAMGGgraph_GEO);
1592: PetscLogEventRegister("PCGAMGcoarse_AGG", PC_CLASSID, &PC_GAMGCoarsen_AGG);
1593: PetscLogEventRegister("PCGAMGcoarse_GEO", PC_CLASSID, &PC_GAMGCoarsen_GEO);
1594: PetscLogEventRegister("PCGAMGProl_AGG", PC_CLASSID, &PC_GAMGProlongator_AGG);
1595: PetscLogEventRegister("PCGAMGProl_GEO", PC_CLASSID, &PC_GAMGProlongator_GEO);
1596: PetscLogEventRegister("PCGAMGPOpt_AGG", PC_CLASSID, &PC_GAMGOptprol_AGG);
1597: PetscLogEventRegister("GAMGKKTProl_AGG", PC_CLASSID, &PC_GAMGKKTProl_AGG);
1598: #endif
1600: return(0);
1601: }
1605: /*@C
1606: PCGAMGInitializePackage - This function initializes everything in the PCGAMG package. It is called
1607: from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to PCCreate_GAMG()
1608: when using static libraries.
1610: Level: developer
1612: .keywords: PC, PCGAMG, initialize, package
1613: .seealso: PetscInitialize()
1614: @*/
1615: PetscErrorCode PCGAMGInitializePackage(void)
1616: {
1620: if (PCGAMGPackageInitialized) return(0);
1621: PCGAMGPackageInitialized = PETSC_TRUE;
1622: PetscFunctionListAdd(&GAMGList,PCGAMGGEO,PCCreateGAMG_GEO);
1623: PetscFunctionListAdd(&GAMGList,PCGAMGAGG,PCCreateGAMG_AGG);
1624: PetscRegisterFinalize(PCGAMGFinalizePackage);
1625: return(0);
1626: }
1630: /*@C
1631: PCGAMGFinalizePackage - This function destroys everything in the PCGAMG package. It is
1632: called from PetscFinalize().
1634: Level: developer
1636: .keywords: Petsc, destroy, package
1637: .seealso: PetscFinalize()
1638: @*/
1639: PetscErrorCode PCGAMGFinalizePackage(void)
1640: {
1644: PCGAMGPackageInitialized = PETSC_FALSE;
1645: PetscFunctionListDestroy(&GAMGList);
1646: return(0);
1647: }