Actual source code: util.c

petsc-3.8.4 2018-03-24
Report Typos and Errors
  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: }