Actual source code: tools.c
petsc-3.6.4 2016-04-12
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>
10: /*
11: Produces a set of block column indices of the matrix row, one for each block represented in the original row
13: n - the number of block indices in cc[]
14: cc - the block indices (must be large enough to contain the indices)
15: */
16: PETSC_STATIC_INLINE PetscErrorCode MatCollapseRow(Mat Amat,PetscInt row,PetscInt bs,PetscInt *n,PetscInt *cc)
17: {
18: PetscInt cnt = -1,nidx,j;
19: const PetscInt *idx;
23: MatGetRow(Amat,row,&nidx,&idx,NULL);
24: if (nidx) {
25: cnt = 0;
26: cc[cnt] = idx[0]/bs;
27: for (j=1; j<nidx; j++) {
28: if (cc[cnt] < idx[j]/bs) cc[++cnt] = idx[j]/bs;
29: }
30: }
31: MatRestoreRow(Amat,row,&nidx,&idx,NULL);
32: *n = cnt+1;
33: return(0);
34: }
38: /*
39: Produces a set of block column indices of the matrix block row, one for each block represented in the original set of rows
41: ncollapsed - the number of block indices
42: collapsed - the block indices (must be large enough to contain the indices)
43: */
44: PETSC_STATIC_INLINE PetscErrorCode MatCollapseRows(Mat Amat,PetscInt start,PetscInt bs,PetscInt *w0,PetscInt *w1,PetscInt *w2,PetscInt *ncollapsed,PetscInt **collapsed)
45: {
46: PetscInt i,nprev,*cprev = w0,ncur = 0,*ccur = w1,*merged = w2,*cprevtmp;
50: MatCollapseRow(Amat,start,bs,&nprev,cprev);
51: for (i=start+1; i<start+bs; i++) {
52: MatCollapseRow(Amat,i,bs,&ncur,ccur);
53: PetscMergeIntArray(nprev,cprev,ncur,ccur,&nprev,&merged);
54: cprevtmp = cprev; cprev = merged; merged = cprevtmp;
55: }
56: *ncollapsed = nprev;
57: if (collapsed) *collapsed = cprev;
58: return(0);
59: }
62: /* -------------------------------------------------------------------------- */
63: /*
64: PCGAMGCreateGraph - create simple scaled scalar graph from matrix
66: Input Parameter:
67: . Amat - matrix
68: Output Parameter:
69: . a_Gmaat - eoutput scalar graph (symmetric?)
70: */
73: PetscErrorCode PCGAMGCreateGraph(Mat Amat, Mat *a_Gmat)
74: {
76: PetscInt Istart,Iend,Ii,i,jj,kk,ncols,nloc,NN,MM,bs;
77: MPI_Comm comm;
78: Mat Gmat;
79: MatType mtype;
82: PetscObjectGetComm((PetscObject)Amat,&comm);
83: MatGetOwnershipRange(Amat, &Istart, &Iend);
84: MatGetSize(Amat, &MM, &NN);
85: MatGetBlockSize(Amat, &bs);
86: nloc = (Iend-Istart)/bs;
88: #if defined PETSC_GAMG_USE_LOG
89: PetscLogEventBegin(petsc_gamg_setup_events[GRAPH],0,0,0,0);
90: #endif
92: if (bs > 1) {
93: const PetscScalar *vals;
94: const PetscInt *idx;
95: PetscInt *d_nnz, *o_nnz,*blockmask = NULL,maskcnt,*w0,*w1,*w2;
96: PetscBool ismpiaij,isseqaij;
98: /*
99: Determine the preallocation needed for the scalar matrix derived from the vector matrix.
100: */
102: PetscObjectTypeCompare((PetscObject)Amat,MATSEQAIJ,&isseqaij);
103: PetscObjectTypeCompare((PetscObject)Amat,MATMPIAIJ,&ismpiaij);
105: PetscMalloc2(nloc, &d_nnz,isseqaij ? 0 : nloc, &o_nnz);
107: if (isseqaij) {
108: PetscInt max_d_nnz;
110: /*
111: Determine exact preallocation count for (sequential) scalar matrix
112: */
113: MatSeqAIJGetMaxRowNonzeros(Amat,&max_d_nnz);
114: max_d_nnz = PetscMin(nloc,bs*max_d_nnz);
115: PetscMalloc3(max_d_nnz, &w0,max_d_nnz, &w1,max_d_nnz, &w2);
116: for (Ii = 0, jj = 0; Ii < Iend; Ii += bs, jj++) {
117: MatCollapseRows(Amat,Ii,bs,w0,w1,w2,&d_nnz[jj],NULL);
118: }
119: PetscFree3(w0,w1,w2);
121: } else if (ismpiaij) {
122: Mat Daij,Oaij;
123: const PetscInt *garray;
124: PetscInt max_d_nnz;
126: MatMPIAIJGetSeqAIJ(Amat,&Daij,&Oaij,&garray);
128: /*
129: Determine exact preallocation count for diagonal block portion of scalar matrix
130: */
131: MatSeqAIJGetMaxRowNonzeros(Daij,&max_d_nnz);
132: max_d_nnz = PetscMin(nloc,bs*max_d_nnz);
133: PetscMalloc3(max_d_nnz, &w0,max_d_nnz, &w1,max_d_nnz, &w2);
134: for (Ii = 0, jj = 0; Ii < Iend - Istart; Ii += bs, jj++) {
135: MatCollapseRows(Daij,Ii,bs,w0,w1,w2,&d_nnz[jj],NULL);
136: }
137: PetscFree3(w0,w1,w2);
139: /*
140: Over estimate (usually grossly over), preallocation count for off-diagonal portion of scalar matrix
141: */
142: for (Ii = 0, jj = 0; Ii < Iend - Istart; Ii += bs, jj++) {
143: o_nnz[jj] = 0;
144: for (kk=0; kk<bs; kk++) { /* rows that get collapsed to a single row */
145: MatGetRow(Oaij,Ii+kk,&ncols,0,0);
146: o_nnz[jj] += ncols;
147: MatRestoreRow(Oaij,Ii+kk,&ncols,0,0);
148: }
149: if (o_nnz[jj] > (NN/bs-nloc)) o_nnz[jj] = NN/bs-nloc;
150: }
152: } else {
153: /*
155: This is O(nloc*nloc/bs) work!
157: This is accurate for the "diagonal" block of the matrix but will be grossly high for the
158: off diagonal block most of the time but could be too low for the off-diagonal.
160: This should be fixed to be accurate for the off-diagonal portion. Cannot just use a mask
161: for the off-diagonal portion since for huge matrices that would require too much memory per
162: MPI process.
163: */
164: PetscMalloc1(nloc, &blockmask);
165: for (Ii = Istart, jj = 0; Ii < Iend; Ii += bs, jj++) {
166: o_nnz[jj] = 0;
167: PetscMemzero(blockmask,nloc*sizeof(PetscInt));
168: for (kk=0; kk<bs; kk++) { /* rows that get collapsed to a single row */
169: MatGetRow(Amat,Ii+kk,&ncols,&idx,0);
170: for (i=0; i<ncols; i++) {
171: if (idx[i] >= Istart && idx[i] < Iend) {
172: blockmask[(idx[i] - Istart)/bs] = 1;
173: }
174: }
175: if (ncols > o_nnz[jj]) {
176: o_nnz[jj] = ncols;
177: if (o_nnz[jj] > (NN/bs-nloc)) o_nnz[jj] = NN/bs-nloc;
178: }
179: MatRestoreRow(Amat,Ii+kk,&ncols,&idx,0);
180: }
181: maskcnt = 0;
182: for (i=0; i<nloc; i++) {
183: if (blockmask[i]) maskcnt++;
184: }
185: d_nnz[jj] = maskcnt;
186: }
187: PetscFree(blockmask);
188: }
190: /* get scalar copy (norms) of matrix -- AIJ specific!!! */
191: MatGetType(Amat,&mtype);
192: MatCreate(comm, &Gmat);
193: MatSetSizes(Gmat,nloc,nloc,PETSC_DETERMINE,PETSC_DETERMINE);
194: MatSetBlockSizes(Gmat, 1, 1);
195: MatSetType(Gmat, mtype);
196: MatSeqAIJSetPreallocation(Gmat,0,d_nnz);
197: MatMPIAIJSetPreallocation(Gmat,0,d_nnz,0,o_nnz);
198: PetscFree2(d_nnz,o_nnz);
200: for (Ii = Istart; Ii < Iend; Ii++) {
201: PetscInt dest_row = Ii/bs;
202: MatGetRow(Amat,Ii,&ncols,&idx,&vals);
203: for (jj=0; jj<ncols; jj++) {
204: PetscInt dest_col = idx[jj]/bs;
205: PetscScalar sv = PetscAbs(PetscRealPart(vals[jj]));
206: MatSetValues(Gmat,1,&dest_row,1,&dest_col,&sv,ADD_VALUES);
207: }
208: MatRestoreRow(Amat,Ii,&ncols,&idx,&vals);
209: }
210: MatAssemblyBegin(Gmat,MAT_FINAL_ASSEMBLY);
211: MatAssemblyEnd(Gmat,MAT_FINAL_ASSEMBLY);
212: } else {
213: /* just copy scalar matrix - abs() not taken here but scaled later */
214: MatDuplicate(Amat, MAT_COPY_VALUES, &Gmat);
215: }
217: #if defined PETSC_GAMG_USE_LOG
218: PetscLogEventEnd(petsc_gamg_setup_events[GRAPH],0,0,0,0);
219: #endif
221: *a_Gmat = Gmat;
222: return(0);
223: }
225: /* -------------------------------------------------------------------------- */
226: /*
227: PCGAMGFilterGraph - filter (remove zero and possibly small values from the) graph and symetrize if needed
229: Input Parameter:
230: . vfilter - threshold paramter [0,1)
231: . symm - symetrize?
232: In/Output Parameter:
233: . a_Gmat - original graph
234: */
237: PetscErrorCode PCGAMGFilterGraph(Mat *a_Gmat,PetscReal vfilter,PetscBool symm)
238: {
239: PetscErrorCode ierr;
240: PetscInt Istart,Iend,Ii,jj,ncols,nnz0,nnz1, NN, MM, nloc;
241: PetscMPIInt rank;
242: Mat Gmat = *a_Gmat, tGmat, matTrans;
243: MPI_Comm comm;
244: const PetscScalar *vals;
245: const PetscInt *idx;
246: PetscInt *d_nnz, *o_nnz;
247: Vec diag;
248: MatType mtype;
251: #if defined PETSC_GAMG_USE_LOG
252: PetscLogEventBegin(petsc_gamg_setup_events[GRAPH],0,0,0,0);
253: #endif
254: /* scale Gmat for all values between -1 and 1 */
255: MatCreateVecs(Gmat, &diag, 0);
256: MatGetDiagonal(Gmat, diag);
257: VecReciprocal(diag);
258: VecSqrtAbs(diag);
259: MatDiagonalScale(Gmat, diag, diag);
260: VecDestroy(&diag);
262: if (vfilter < 0.0 && !symm) {
263: /* Just use the provided matrix as the graph but make all values positive */
264: Mat_MPIAIJ *aij = (Mat_MPIAIJ*)Gmat->data;
265: MatInfo info;
266: PetscScalar *avals;
267: PetscMPIInt size;
269: MPI_Comm_size(PetscObjectComm((PetscObject)Gmat),&size);
270: if (size == 1) {
271: MatGetInfo(Gmat,MAT_LOCAL,&info);
272: MatSeqAIJGetArray(Gmat,&avals);
273: for (jj = 0; jj<info.nz_used; jj++) avals[jj] = PetscAbsScalar(avals[jj]);
274: MatSeqAIJRestoreArray(Gmat,&avals);
275: } else {
276: MatGetInfo(aij->A,MAT_LOCAL,&info);
277: MatSeqAIJGetArray(aij->A,&avals);
278: for (jj = 0; jj<info.nz_used; jj++) avals[jj] = PetscAbsScalar(avals[jj]);
279: MatSeqAIJRestoreArray(aij->A,&avals);
280: MatGetInfo(aij->B,MAT_LOCAL,&info);
281: MatSeqAIJGetArray(aij->B,&avals);
282: for (jj = 0; jj<info.nz_used; jj++) avals[jj] = PetscAbsScalar(avals[jj]);
283: MatSeqAIJRestoreArray(aij->B,&avals);
284: }
285: #if defined PETSC_GAMG_USE_LOG
286: PetscLogEventEnd(petsc_gamg_setup_events[GRAPH],0,0,0,0);
287: #endif
288: return(0);
289: }
291: PetscObjectGetComm((PetscObject)Gmat,&comm);
292: MPI_Comm_rank(comm,&rank);
293: MatGetOwnershipRange(Gmat, &Istart, &Iend);
294: nloc = Iend - Istart;
295: MatGetSize(Gmat, &MM, &NN);
297: if (symm) {
298: MatTranspose(Gmat, MAT_INITIAL_MATRIX, &matTrans);
299: }
301: /* Determine upper bound on nonzeros needed in new filtered matrix */
302: PetscMalloc2(nloc, &d_nnz,nloc, &o_nnz);
303: for (Ii = Istart, jj = 0; Ii < Iend; Ii++, jj++) {
304: MatGetRow(Gmat,Ii,&ncols,NULL,NULL);
305: d_nnz[jj] = ncols;
306: o_nnz[jj] = ncols;
307: MatRestoreRow(Gmat,Ii,&ncols,NULL,NULL);
308: if (symm) {
309: MatGetRow(matTrans,Ii,&ncols,NULL,NULL);
310: d_nnz[jj] += ncols;
311: o_nnz[jj] += ncols;
312: MatRestoreRow(matTrans,Ii,&ncols,NULL,NULL);
313: }
314: if (d_nnz[jj] > nloc) d_nnz[jj] = nloc;
315: if (o_nnz[jj] > (MM-nloc)) o_nnz[jj] = MM - nloc;
316: }
317: MatGetType(Gmat,&mtype);
318: MatCreate(comm, &tGmat);
319: MatSetSizes(tGmat,nloc,nloc,MM,MM);
320: MatSetBlockSizes(tGmat, 1, 1);
321: MatSetType(tGmat, mtype);
322: MatSeqAIJSetPreallocation(tGmat,0,d_nnz);
323: MatMPIAIJSetPreallocation(tGmat,0,d_nnz,0,o_nnz);
324: PetscFree2(d_nnz,o_nnz);
325: if (symm) {
326: MatDestroy(&matTrans);
327: } else {
328: /* all entries are generated locally so MatAssembly will be slightly faster for large process counts */
329: MatSetOption(tGmat,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
330: }
332: for (Ii = Istart, nnz0 = nnz1 = 0; Ii < Iend; Ii++) {
333: MatGetRow(Gmat,Ii,&ncols,&idx,&vals);
334: for (jj=0; jj<ncols; jj++,nnz0++) {
335: PetscScalar sv = PetscAbs(PetscRealPart(vals[jj]));
336: if (PetscRealPart(sv) > vfilter) {
337: nnz1++;
338: if (symm) {
339: sv *= 0.5;
340: MatSetValues(tGmat,1,&Ii,1,&idx[jj],&sv,ADD_VALUES);
341: MatSetValues(tGmat,1,&idx[jj],1,&Ii,&sv,ADD_VALUES);
342: } else {
343: MatSetValues(tGmat,1,&Ii,1,&idx[jj],&sv,ADD_VALUES);
344: }
345: }
346: }
347: MatRestoreRow(Gmat,Ii,&ncols,&idx,&vals);
348: }
349: MatAssemblyBegin(tGmat,MAT_FINAL_ASSEMBLY);
350: MatAssemblyEnd(tGmat,MAT_FINAL_ASSEMBLY);
352: #if defined PETSC_GAMG_USE_LOG
353: PetscLogEventEnd(petsc_gamg_setup_events[GRAPH],0,0,0,0);
354: #endif
356: #if defined(PETSC_USE_INFO)
357: {
358: double t1 = (!nnz0) ? 1. : 100.*(double)nnz1/(double)nnz0, t2 = (!nloc) ? 1. : (double)nnz0/(double)nloc;
359: PetscInfo4(*a_Gmat,"\t %g%% nnz after filtering, with threshold %g, %g nnz ave. (N=%D)\n",t1,vfilter,t2,MM);
360: }
361: #endif
362: MatDestroy(&Gmat);
363: *a_Gmat = tGmat;
364: return(0);
365: }
367: /* -------------------------------------------------------------------------- */
368: /*
369: PCGAMGGetDataWithGhosts - hacks into Mat MPIAIJ so this must have > 1 pe
371: Input Parameter:
372: . Gmat - MPIAIJ matrix for scattters
373: . data_sz - number of data terms per node (# cols in output)
374: . data_in[nloc*data_sz] - column oriented data
375: Output Parameter:
376: . a_stride - numbrt of rows of output
377: . a_data_out[stride*data_sz] - output data with ghosts
378: */
381: PetscErrorCode PCGAMGGetDataWithGhosts(Mat Gmat,PetscInt data_sz,PetscReal data_in[],PetscInt *a_stride,PetscReal **a_data_out)
382: {
384: MPI_Comm comm;
385: Vec tmp_crds;
386: Mat_MPIAIJ *mpimat = (Mat_MPIAIJ*)Gmat->data;
387: PetscInt nnodes,num_ghosts,dir,kk,jj,my0,Iend,nloc;
388: PetscScalar *data_arr;
389: PetscReal *datas;
390: PetscBool isMPIAIJ;
393: PetscObjectGetComm((PetscObject)Gmat,&comm);
394: PetscObjectTypeCompare((PetscObject)Gmat, MATMPIAIJ, &isMPIAIJ);
395: MatGetOwnershipRange(Gmat, &my0, &Iend);
396: nloc = Iend - my0;
397: VecGetLocalSize(mpimat->lvec, &num_ghosts);
398: nnodes = num_ghosts + nloc;
399: *a_stride = nnodes;
400: MatCreateVecs(Gmat, &tmp_crds, 0);
402: PetscMalloc1(data_sz*nnodes, &datas);
403: for (dir=0; dir<data_sz; dir++) {
404: /* set local, and global */
405: for (kk=0; kk<nloc; kk++) {
406: PetscInt gid = my0 + kk;
407: PetscScalar crd = (PetscScalar)data_in[dir*nloc + kk]; /* col oriented */
408: datas[dir*nnodes + kk] = PetscRealPart(crd);
410: VecSetValues(tmp_crds, 1, &gid, &crd, INSERT_VALUES);
411: }
412: VecAssemblyBegin(tmp_crds);
413: VecAssemblyEnd(tmp_crds);
414: /* get ghost datas */
415: VecScatterBegin(mpimat->Mvctx,tmp_crds,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
416: VecScatterEnd(mpimat->Mvctx,tmp_crds,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
417: VecGetArray(mpimat->lvec, &data_arr);
418: for (kk=nloc,jj=0;jj<num_ghosts;kk++,jj++) datas[dir*nnodes + kk] = PetscRealPart(data_arr[jj]);
419: VecRestoreArray(mpimat->lvec, &data_arr);
420: }
421: VecDestroy(&tmp_crds);
422: *a_data_out = datas;
423: return(0);
424: }
427: /*
428: *
429: * GAMGTableCreate
430: */
434: PetscErrorCode GAMGTableCreate(PetscInt a_size, GAMGHashTable *a_tab)
435: {
437: PetscInt kk;
440: a_tab->size = a_size;
441: PetscMalloc1(a_size, &a_tab->table);
442: PetscMalloc1(a_size, &a_tab->data);
443: for (kk=0; kk<a_size; kk++) a_tab->table[kk] = -1;
444: return(0);
445: }
449: PetscErrorCode GAMGTableDestroy(GAMGHashTable *a_tab)
450: {
454: PetscFree(a_tab->table);
455: PetscFree(a_tab->data);
456: return(0);
457: }
461: PetscErrorCode GAMGTableAdd(GAMGHashTable *a_tab, PetscInt a_key, PetscInt a_data)
462: {
463: PetscInt kk,idx;
466: if (a_key<0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Negative key %d.",a_key);
467: for (kk = 0, idx = GAMG_HASH(a_key);
468: kk < a_tab->size;
469: kk++, idx = (idx==(a_tab->size-1)) ? 0 : idx + 1) {
471: if (a_tab->table[idx] == a_key) {
472: /* exists */
473: a_tab->data[idx] = a_data;
474: break;
475: } else if (a_tab->table[idx] == -1) {
476: /* add */
477: a_tab->table[idx] = a_key;
478: a_tab->data[idx] = a_data;
479: break;
480: }
481: }
482: if (kk==a_tab->size) {
483: /* this is not to efficient, waiting until completely full */
484: PetscInt oldsize = a_tab->size, new_size = 2*a_tab->size + 5, *oldtable = a_tab->table, *olddata = a_tab->data;
487: a_tab->size = new_size;
489: PetscMalloc1(a_tab->size, &a_tab->table);
490: PetscMalloc1(a_tab->size, &a_tab->data);
492: for (kk=0;kk<a_tab->size;kk++) a_tab->table[kk] = -1;
493: for (kk=0;kk<oldsize;kk++) {
494: if (oldtable[kk] != -1) {
495: GAMGTableAdd(a_tab, oldtable[kk], olddata[kk]);
496: }
497: }
498: PetscFree(oldtable);
499: PetscFree(olddata);
500: GAMGTableAdd(a_tab, a_key, a_data);
501: }
502: return(0);
503: }