Actual source code: util.c
petsc-3.8.4 2018-03-24
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>
7: /*
8: Produces a set of block column indices of the matrix row, one for each block represented in the original row
10: n - the number of block indices in cc[]
11: cc - the block indices (must be large enough to contain the indices)
12: */
13: PETSC_STATIC_INLINE PetscErrorCode MatCollapseRow(Mat Amat,PetscInt row,PetscInt bs,PetscInt *n,PetscInt *cc)
14: {
15: PetscInt cnt = -1,nidx,j;
16: const PetscInt *idx;
20: MatGetRow(Amat,row,&nidx,&idx,NULL);
21: if (nidx) {
22: cnt = 0;
23: cc[cnt] = idx[0]/bs;
24: for (j=1; j<nidx; j++) {
25: if (cc[cnt] < idx[j]/bs) cc[++cnt] = idx[j]/bs;
26: }
27: }
28: MatRestoreRow(Amat,row,&nidx,&idx,NULL);
29: *n = cnt+1;
30: return(0);
31: }
33: /*
34: Produces a set of block column indices of the matrix block row, one for each block represented in the original set of rows
36: ncollapsed - the number of block indices
37: collapsed - the block indices (must be large enough to contain the indices)
38: */
39: PETSC_STATIC_INLINE PetscErrorCode MatCollapseRows(Mat Amat,PetscInt start,PetscInt bs,PetscInt *w0,PetscInt *w1,PetscInt *w2,PetscInt *ncollapsed,PetscInt **collapsed)
40: {
41: PetscInt i,nprev,*cprev = w0,ncur = 0,*ccur = w1,*merged = w2,*cprevtmp;
45: MatCollapseRow(Amat,start,bs,&nprev,cprev);
46: for (i=start+1; i<start+bs; i++) {
47: MatCollapseRow(Amat,i,bs,&ncur,ccur);
48: PetscMergeIntArray(nprev,cprev,ncur,ccur,&nprev,&merged);
49: cprevtmp = cprev; cprev = merged; merged = cprevtmp;
50: }
51: *ncollapsed = nprev;
52: if (collapsed) *collapsed = cprev;
53: return(0);
54: }
57: /* -------------------------------------------------------------------------- */
58: /*
59: PCGAMGCreateGraph - create simple scaled scalar graph from matrix
61: Input Parameter:
62: . Amat - matrix
63: Output Parameter:
64: . a_Gmaat - eoutput scalar graph (symmetric?)
65: */
66: PetscErrorCode PCGAMGCreateGraph(Mat Amat, Mat *a_Gmat)
67: {
69: PetscInt Istart,Iend,Ii,jj,kk,ncols,nloc,NN,MM,bs;
70: MPI_Comm comm;
71: Mat Gmat;
72: MatType mtype;
75: PetscObjectGetComm((PetscObject)Amat,&comm);
76: MatGetOwnershipRange(Amat, &Istart, &Iend);
77: MatGetSize(Amat, &MM, &NN);
78: MatGetBlockSize(Amat, &bs);
79: nloc = (Iend-Istart)/bs;
81: #if defined PETSC_GAMG_USE_LOG
82: PetscLogEventBegin(petsc_gamg_setup_events[GRAPH],0,0,0,0);
83: #endif
85: if (bs > 1) {
86: const PetscScalar *vals;
87: const PetscInt *idx;
88: PetscInt *d_nnz, *o_nnz,*w0,*w1,*w2;
89: PetscBool ismpiaij,isseqaij;
91: /*
92: Determine the preallocation needed for the scalar matrix derived from the vector matrix.
93: */
95: PetscObjectBaseTypeCompare((PetscObject)Amat,MATSEQAIJ,&isseqaij);
96: PetscObjectBaseTypeCompare((PetscObject)Amat,MATMPIAIJ,&ismpiaij);
98: PetscMalloc2(nloc, &d_nnz,isseqaij ? 0 : nloc, &o_nnz);
100: if (isseqaij) {
101: PetscInt max_d_nnz;
103: /*
104: Determine exact preallocation count for (sequential) scalar matrix
105: */
106: MatSeqAIJGetMaxRowNonzeros(Amat,&max_d_nnz);
107: max_d_nnz = PetscMin(nloc,bs*max_d_nnz);
108: PetscMalloc3(max_d_nnz, &w0,max_d_nnz, &w1,max_d_nnz, &w2);
109: for (Ii = 0, jj = 0; Ii < Iend; Ii += bs, jj++) {
110: MatCollapseRows(Amat,Ii,bs,w0,w1,w2,&d_nnz[jj],NULL);
111: }
112: PetscFree3(w0,w1,w2);
114: } else if (ismpiaij) {
115: Mat Daij,Oaij;
116: const PetscInt *garray;
117: PetscInt max_d_nnz;
119: MatMPIAIJGetSeqAIJ(Amat,&Daij,&Oaij,&garray);
121: /*
122: Determine exact preallocation count for diagonal block portion of scalar matrix
123: */
124: MatSeqAIJGetMaxRowNonzeros(Daij,&max_d_nnz);
125: max_d_nnz = PetscMin(nloc,bs*max_d_nnz);
126: PetscMalloc3(max_d_nnz, &w0,max_d_nnz, &w1,max_d_nnz, &w2);
127: for (Ii = 0, jj = 0; Ii < Iend - Istart; Ii += bs, jj++) {
128: MatCollapseRows(Daij,Ii,bs,w0,w1,w2,&d_nnz[jj],NULL);
129: }
130: PetscFree3(w0,w1,w2);
132: /*
133: Over estimate (usually grossly over), preallocation count for off-diagonal portion of scalar matrix
134: */
135: for (Ii = 0, jj = 0; Ii < Iend - Istart; Ii += bs, jj++) {
136: o_nnz[jj] = 0;
137: for (kk=0; kk<bs; kk++) { /* rows that get collapsed to a single row */
138: MatGetRow(Oaij,Ii+kk,&ncols,0,0);
139: o_nnz[jj] += ncols;
140: MatRestoreRow(Oaij,Ii+kk,&ncols,0,0);
141: }
142: if (o_nnz[jj] > (NN/bs-nloc)) o_nnz[jj] = NN/bs-nloc;
143: }
145: } else {
146: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Require AIJ matrix type");
147: }
149: /* get scalar copy (norms) of matrix */
150: MatGetType(Amat,&mtype);
151: MatCreate(comm, &Gmat);
152: MatSetSizes(Gmat,nloc,nloc,PETSC_DETERMINE,PETSC_DETERMINE);
153: MatSetBlockSizes(Gmat, 1, 1);
154: MatSetType(Gmat, mtype);
155: MatSeqAIJSetPreallocation(Gmat,0,d_nnz);
156: MatMPIAIJSetPreallocation(Gmat,0,d_nnz,0,o_nnz);
157: PetscFree2(d_nnz,o_nnz);
159: for (Ii = Istart; Ii < Iend; Ii++) {
160: PetscInt dest_row = Ii/bs;
161: MatGetRow(Amat,Ii,&ncols,&idx,&vals);
162: for (jj=0; jj<ncols; jj++) {
163: PetscInt dest_col = idx[jj]/bs;
164: PetscScalar sv = PetscAbs(PetscRealPart(vals[jj]));
165: MatSetValues(Gmat,1,&dest_row,1,&dest_col,&sv,ADD_VALUES);
166: }
167: MatRestoreRow(Amat,Ii,&ncols,&idx,&vals);
168: }
169: MatAssemblyBegin(Gmat,MAT_FINAL_ASSEMBLY);
170: MatAssemblyEnd(Gmat,MAT_FINAL_ASSEMBLY);
171: } else {
172: /* just copy scalar matrix - abs() not taken here but scaled later */
173: MatDuplicate(Amat, MAT_COPY_VALUES, &Gmat);
174: }
176: #if defined PETSC_GAMG_USE_LOG
177: PetscLogEventEnd(petsc_gamg_setup_events[GRAPH],0,0,0,0);
178: #endif
180: *a_Gmat = Gmat;
181: return(0);
182: }
184: /* -------------------------------------------------------------------------- */
185: /*@C
186: PCGAMGFilterGraph - filter (remove zero and possibly small values from the) graph and make it symmetric if requested
188: Collective on Mat
190: Input Parameter:
191: + a_Gmat - the graph
192: . vfilter - threshold paramter [0,1)
193: - symm - make the result symmetric
195: Level: developer
197: Notes: This is called before graph coarsers are called.
199: .seealso: PCGAMGSetThreshold()
200: @*/
201: PetscErrorCode PCGAMGFilterGraph(Mat *a_Gmat,PetscReal vfilter,PetscBool symm)
202: {
203: PetscErrorCode ierr;
204: PetscInt Istart,Iend,Ii,jj,ncols,nnz0,nnz1, NN, MM, nloc;
205: PetscMPIInt rank;
206: Mat Gmat = *a_Gmat, tGmat, matTrans;
207: MPI_Comm comm;
208: const PetscScalar *vals;
209: const PetscInt *idx;
210: PetscInt *d_nnz, *o_nnz;
211: Vec diag;
212: MatType mtype;
215: #if defined PETSC_GAMG_USE_LOG
216: PetscLogEventBegin(petsc_gamg_setup_events[GRAPH],0,0,0,0);
217: #endif
218: /* scale Gmat for all values between -1 and 1 */
219: MatCreateVecs(Gmat, &diag, 0);
220: MatGetDiagonal(Gmat, diag);
221: VecReciprocal(diag);
222: VecSqrtAbs(diag);
223: MatDiagonalScale(Gmat, diag, diag);
224: VecDestroy(&diag);
226: if (vfilter < 0.0 && !symm) {
227: /* Just use the provided matrix as the graph but make all values positive */
228: MatInfo info;
229: PetscScalar *avals;
230: PetscBool isaij,ismpiaij;
231: PetscObjectBaseTypeCompare((PetscObject)Gmat,MATSEQAIJ,&isaij);
232: PetscObjectBaseTypeCompare((PetscObject)Gmat,MATMPIAIJ,&ismpiaij);
233: if (!isaij && !ismpiaij) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Require (MPI)AIJ matrix type");
234: if (isaij) {
235: MatGetInfo(Gmat,MAT_LOCAL,&info);
236: MatSeqAIJGetArray(Gmat,&avals);
237: for (jj = 0; jj<info.nz_used; jj++) avals[jj] = PetscAbsScalar(avals[jj]);
238: MatSeqAIJRestoreArray(Gmat,&avals);
239: } else {
240: Mat_MPIAIJ *aij = (Mat_MPIAIJ*)Gmat->data;
241: MatGetInfo(aij->A,MAT_LOCAL,&info);
242: MatSeqAIJGetArray(aij->A,&avals);
243: for (jj = 0; jj<info.nz_used; jj++) avals[jj] = PetscAbsScalar(avals[jj]);
244: MatSeqAIJRestoreArray(aij->A,&avals);
245: MatGetInfo(aij->B,MAT_LOCAL,&info);
246: MatSeqAIJGetArray(aij->B,&avals);
247: for (jj = 0; jj<info.nz_used; jj++) avals[jj] = PetscAbsScalar(avals[jj]);
248: MatSeqAIJRestoreArray(aij->B,&avals);
249: }
250: #if defined PETSC_GAMG_USE_LOG
251: PetscLogEventEnd(petsc_gamg_setup_events[GRAPH],0,0,0,0);
252: #endif
253: return(0);
254: }
256: PetscObjectGetComm((PetscObject)Gmat,&comm);
257: MPI_Comm_rank(comm,&rank);
258: MatGetOwnershipRange(Gmat, &Istart, &Iend);
259: nloc = Iend - Istart;
260: MatGetSize(Gmat, &MM, &NN);
262: if (symm) {
263: MatTranspose(Gmat, MAT_INITIAL_MATRIX, &matTrans);
264: }
266: /* Determine upper bound on nonzeros needed in new filtered matrix */
267: PetscMalloc2(nloc, &d_nnz,nloc, &o_nnz);
268: for (Ii = Istart, jj = 0; Ii < Iend; Ii++, jj++) {
269: MatGetRow(Gmat,Ii,&ncols,NULL,NULL);
270: d_nnz[jj] = ncols;
271: o_nnz[jj] = ncols;
272: MatRestoreRow(Gmat,Ii,&ncols,NULL,NULL);
273: if (symm) {
274: MatGetRow(matTrans,Ii,&ncols,NULL,NULL);
275: d_nnz[jj] += ncols;
276: o_nnz[jj] += ncols;
277: MatRestoreRow(matTrans,Ii,&ncols,NULL,NULL);
278: }
279: if (d_nnz[jj] > nloc) d_nnz[jj] = nloc;
280: if (o_nnz[jj] > (MM-nloc)) o_nnz[jj] = MM - nloc;
281: }
282: MatGetType(Gmat,&mtype);
283: MatCreate(comm, &tGmat);
284: MatSetSizes(tGmat,nloc,nloc,MM,MM);
285: MatSetBlockSizes(tGmat, 1, 1);
286: MatSetType(tGmat, mtype);
287: MatSeqAIJSetPreallocation(tGmat,0,d_nnz);
288: MatMPIAIJSetPreallocation(tGmat,0,d_nnz,0,o_nnz);
289: PetscFree2(d_nnz,o_nnz);
290: if (symm) {
291: MatDestroy(&matTrans);
292: } else {
293: /* all entries are generated locally so MatAssembly will be slightly faster for large process counts */
294: MatSetOption(tGmat,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
295: }
297: for (Ii = Istart, nnz0 = nnz1 = 0; Ii < Iend; Ii++) {
298: MatGetRow(Gmat,Ii,&ncols,&idx,&vals);
299: for (jj=0; jj<ncols; jj++,nnz0++) {
300: PetscScalar sv = PetscAbs(PetscRealPart(vals[jj]));
301: if (PetscRealPart(sv) > vfilter) {
302: nnz1++;
303: if (symm) {
304: sv *= 0.5;
305: MatSetValues(tGmat,1,&Ii,1,&idx[jj],&sv,ADD_VALUES);
306: MatSetValues(tGmat,1,&idx[jj],1,&Ii,&sv,ADD_VALUES);
307: } else {
308: MatSetValues(tGmat,1,&Ii,1,&idx[jj],&sv,ADD_VALUES);
309: }
310: }
311: }
312: MatRestoreRow(Gmat,Ii,&ncols,&idx,&vals);
313: }
314: MatAssemblyBegin(tGmat,MAT_FINAL_ASSEMBLY);
315: MatAssemblyEnd(tGmat,MAT_FINAL_ASSEMBLY);
317: #if defined PETSC_GAMG_USE_LOG
318: PetscLogEventEnd(petsc_gamg_setup_events[GRAPH],0,0,0,0);
319: #endif
321: #if defined(PETSC_USE_INFO)
322: {
323: double t1 = (!nnz0) ? 1. : 100.*(double)nnz1/(double)nnz0, t2 = (!nloc) ? 1. : (double)nnz0/(double)nloc;
324: PetscInfo4(*a_Gmat,"\t %g%% nnz after filtering, with threshold %g, %g nnz ave. (N=%D)\n",t1,vfilter,t2,MM);
325: }
326: #endif
327: MatDestroy(&Gmat);
328: *a_Gmat = tGmat;
329: return(0);
330: }
332: /* -------------------------------------------------------------------------- */
333: /*
334: PCGAMGGetDataWithGhosts - hacks into Mat MPIAIJ so this must have npe > 1
336: Input Parameter:
337: . Gmat - MPIAIJ matrix for scattters
338: . data_sz - number of data terms per node (# cols in output)
339: . data_in[nloc*data_sz] - column oriented data
340: Output Parameter:
341: . a_stride - numbrt of rows of output
342: . a_data_out[stride*data_sz] - output data with ghosts
343: */
344: PetscErrorCode PCGAMGGetDataWithGhosts(Mat Gmat,PetscInt data_sz,PetscReal data_in[],PetscInt *a_stride,PetscReal **a_data_out)
345: {
347: MPI_Comm comm;
348: Vec tmp_crds;
349: Mat_MPIAIJ *mpimat = (Mat_MPIAIJ*)Gmat->data;
350: PetscInt nnodes,num_ghosts,dir,kk,jj,my0,Iend,nloc;
351: PetscScalar *data_arr;
352: PetscReal *datas;
353: PetscBool isMPIAIJ;
356: PetscObjectGetComm((PetscObject)Gmat,&comm);
357: PetscObjectBaseTypeCompare((PetscObject)Gmat, MATMPIAIJ, &isMPIAIJ);
358: MatGetOwnershipRange(Gmat, &my0, &Iend);
359: nloc = Iend - my0;
360: VecGetLocalSize(mpimat->lvec, &num_ghosts);
361: nnodes = num_ghosts + nloc;
362: *a_stride = nnodes;
363: MatCreateVecs(Gmat, &tmp_crds, 0);
365: PetscMalloc1(data_sz*nnodes, &datas);
366: for (dir=0; dir<data_sz; dir++) {
367: /* set local, and global */
368: for (kk=0; kk<nloc; kk++) {
369: PetscInt gid = my0 + kk;
370: PetscScalar crd = (PetscScalar)data_in[dir*nloc + kk]; /* col oriented */
371: datas[dir*nnodes + kk] = PetscRealPart(crd);
373: VecSetValues(tmp_crds, 1, &gid, &crd, INSERT_VALUES);
374: }
375: VecAssemblyBegin(tmp_crds);
376: VecAssemblyEnd(tmp_crds);
377: /* get ghost datas */
378: VecScatterBegin(mpimat->Mvctx,tmp_crds,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
379: VecScatterEnd(mpimat->Mvctx,tmp_crds,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
380: VecGetArray(mpimat->lvec, &data_arr);
381: for (kk=nloc,jj=0;jj<num_ghosts;kk++,jj++) datas[dir*nnodes + kk] = PetscRealPart(data_arr[jj]);
382: VecRestoreArray(mpimat->lvec, &data_arr);
383: }
384: VecDestroy(&tmp_crds);
385: *a_data_out = datas;
386: return(0);
387: }
390: /*
391: *
392: * PCGAMGHashTableCreate
393: */
395: PetscErrorCode PCGAMGHashTableCreate(PetscInt a_size, PCGAMGHashTable *a_tab)
396: {
398: PetscInt kk;
401: a_tab->size = a_size;
402: PetscMalloc1(a_size, &a_tab->table);
403: PetscMalloc1(a_size, &a_tab->data);
404: for (kk=0; kk<a_size; kk++) a_tab->table[kk] = -1;
405: return(0);
406: }
408: PetscErrorCode PCGAMGHashTableDestroy(PCGAMGHashTable *a_tab)
409: {
413: PetscFree(a_tab->table);
414: PetscFree(a_tab->data);
415: return(0);
416: }
418: PetscErrorCode PCGAMGHashTableAdd(PCGAMGHashTable *a_tab, PetscInt a_key, PetscInt a_data)
419: {
420: PetscInt kk,idx;
423: if (a_key<0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Negative key %d.",a_key);
424: for (kk = 0, idx = GAMG_HASH(a_key);
425: kk < a_tab->size;
426: kk++, idx = (idx==(a_tab->size-1)) ? 0 : idx + 1) {
428: if (a_tab->table[idx] == a_key) {
429: /* exists */
430: a_tab->data[idx] = a_data;
431: break;
432: } else if (a_tab->table[idx] == -1) {
433: /* add */
434: a_tab->table[idx] = a_key;
435: a_tab->data[idx] = a_data;
436: break;
437: }
438: }
439: if (kk==a_tab->size) {
440: /* this is not to efficient, waiting until completely full */
441: PetscInt oldsize = a_tab->size, new_size = 2*a_tab->size + 5, *oldtable = a_tab->table, *olddata = a_tab->data;
444: a_tab->size = new_size;
446: PetscMalloc1(a_tab->size, &a_tab->table);
447: PetscMalloc1(a_tab->size, &a_tab->data);
449: for (kk=0;kk<a_tab->size;kk++) a_tab->table[kk] = -1;
450: for (kk=0;kk<oldsize;kk++) {
451: if (oldtable[kk] != -1) {
452: PCGAMGHashTableAdd(a_tab, oldtable[kk], olddata[kk]);
453: }
454: }
455: PetscFree(oldtable);
456: PetscFree(olddata);
457: PCGAMGHashTableAdd(a_tab, a_key, a_data);
458: }
459: return(0);
460: }