xref: /petsc/src/ksp/pc/impls/gamg/util.c (revision 8f3cd775ab5449ffc0e0df3829b3fa5c4fb1eb7d)
16618991cSMark Adams /*
26618991cSMark Adams  GAMG geometric-algebric multigrid PC - Mark Adams 2011
36618991cSMark Adams  */
46618991cSMark Adams #include <petsc/private/matimpl.h>
56618991cSMark Adams #include <../src/ksp/pc/impls/gamg/gamg.h>           /*I "petscpc.h" I*/
66618991cSMark Adams 
76618991cSMark Adams /*
86618991cSMark Adams    Produces a set of block column indices of the matrix row, one for each block represented in the original row
96618991cSMark Adams 
106618991cSMark Adams    n - the number of block indices in cc[]
116618991cSMark Adams    cc - the block indices (must be large enough to contain the indices)
126618991cSMark Adams */
136618991cSMark Adams PETSC_STATIC_INLINE PetscErrorCode MatCollapseRow(Mat Amat,PetscInt row,PetscInt bs,PetscInt *n,PetscInt *cc)
146618991cSMark Adams {
156618991cSMark Adams   PetscInt       cnt = -1,nidx,j;
166618991cSMark Adams   const PetscInt *idx;
176618991cSMark Adams   PetscErrorCode ierr;
186618991cSMark Adams 
196618991cSMark Adams   PetscFunctionBegin;
206618991cSMark Adams   ierr = MatGetRow(Amat,row,&nidx,&idx,NULL);CHKERRQ(ierr);
216618991cSMark Adams   if (nidx) {
226618991cSMark Adams     cnt = 0;
236618991cSMark Adams     cc[cnt] = idx[0]/bs;
246618991cSMark Adams     for (j=1; j<nidx; j++) {
256618991cSMark Adams       if (cc[cnt] < idx[j]/bs) cc[++cnt] = idx[j]/bs;
266618991cSMark Adams     }
276618991cSMark Adams   }
286618991cSMark Adams   ierr = MatRestoreRow(Amat,row,&nidx,&idx,NULL);CHKERRQ(ierr);
296618991cSMark Adams   *n = cnt+1;
306618991cSMark Adams   PetscFunctionReturn(0);
316618991cSMark Adams }
326618991cSMark Adams 
336618991cSMark Adams /*
346618991cSMark Adams     Produces a set of block column indices of the matrix block row, one for each block represented in the original set of rows
356618991cSMark Adams 
366618991cSMark Adams     ncollapsed - the number of block indices
376618991cSMark Adams     collapsed - the block indices (must be large enough to contain the indices)
386618991cSMark Adams */
396618991cSMark Adams PETSC_STATIC_INLINE PetscErrorCode MatCollapseRows(Mat Amat,PetscInt start,PetscInt bs,PetscInt *w0,PetscInt *w1,PetscInt *w2,PetscInt *ncollapsed,PetscInt **collapsed)
406618991cSMark Adams {
416618991cSMark Adams   PetscInt       i,nprev,*cprev = w0,ncur = 0,*ccur = w1,*merged = w2,*cprevtmp;
426618991cSMark Adams   PetscErrorCode ierr;
436618991cSMark Adams 
446618991cSMark Adams   PetscFunctionBegin;
456618991cSMark Adams   ierr = MatCollapseRow(Amat,start,bs,&nprev,cprev);CHKERRQ(ierr);
466618991cSMark Adams   for (i=start+1; i<start+bs; i++) {
476618991cSMark Adams     ierr  = MatCollapseRow(Amat,i,bs,&ncur,ccur);CHKERRQ(ierr);
486618991cSMark Adams     ierr  = PetscMergeIntArray(nprev,cprev,ncur,ccur,&nprev,&merged);CHKERRQ(ierr);
496618991cSMark Adams     cprevtmp = cprev; cprev = merged; merged = cprevtmp;
506618991cSMark Adams   }
516618991cSMark Adams   *ncollapsed = nprev;
526618991cSMark Adams   if (collapsed) *collapsed  = cprev;
536618991cSMark Adams   PetscFunctionReturn(0);
546618991cSMark Adams }
556618991cSMark Adams 
566618991cSMark Adams 
576618991cSMark Adams /* -------------------------------------------------------------------------- */
586618991cSMark Adams /*
596618991cSMark Adams    PCGAMGCreateGraph - create simple scaled scalar graph from matrix
606618991cSMark Adams 
616618991cSMark Adams  Input Parameter:
626618991cSMark Adams  . Amat - matrix
636618991cSMark Adams  Output Parameter:
646618991cSMark Adams  . a_Gmaat - eoutput scalar graph (symmetric?)
656618991cSMark Adams  */
666618991cSMark Adams PetscErrorCode PCGAMGCreateGraph(Mat Amat, Mat *a_Gmat)
676618991cSMark Adams {
686618991cSMark Adams   PetscErrorCode ierr;
69f42dcbb3SMark Adams   PetscInt       Istart,Iend,Ii,jj,kk,ncols,nloc,NN,MM,bs;
706618991cSMark Adams   MPI_Comm       comm;
716618991cSMark Adams   Mat            Gmat;
726618991cSMark Adams   MatType        mtype;
736618991cSMark Adams 
746618991cSMark Adams   PetscFunctionBegin;
756618991cSMark Adams   ierr = PetscObjectGetComm((PetscObject)Amat,&comm);CHKERRQ(ierr);
766618991cSMark Adams   ierr = MatGetOwnershipRange(Amat, &Istart, &Iend);CHKERRQ(ierr);
776618991cSMark Adams   ierr = MatGetSize(Amat, &MM, &NN);CHKERRQ(ierr);
786618991cSMark Adams   ierr = MatGetBlockSize(Amat, &bs);CHKERRQ(ierr);
796618991cSMark Adams   nloc = (Iend-Istart)/bs;
806618991cSMark Adams 
816618991cSMark Adams #if defined PETSC_GAMG_USE_LOG
826618991cSMark Adams   ierr = PetscLogEventBegin(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr);
836618991cSMark Adams #endif
846618991cSMark Adams 
856618991cSMark Adams   if (bs > 1) {
866618991cSMark Adams     const PetscScalar *vals;
876618991cSMark Adams     const PetscInt    *idx;
88f42dcbb3SMark Adams     PetscInt          *d_nnz, *o_nnz,*w0,*w1,*w2;
896618991cSMark Adams     PetscBool         ismpiaij,isseqaij;
906618991cSMark Adams 
916618991cSMark Adams     /*
926618991cSMark Adams        Determine the preallocation needed for the scalar matrix derived from the vector matrix.
936618991cSMark Adams     */
946618991cSMark Adams 
954099cc6bSBarry Smith     ierr = PetscObjectBaseTypeCompare((PetscObject)Amat,MATSEQAIJ,&isseqaij);CHKERRQ(ierr);
964099cc6bSBarry Smith     ierr = PetscObjectBaseTypeCompare((PetscObject)Amat,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr);
976618991cSMark Adams 
986618991cSMark Adams     ierr = PetscMalloc2(nloc, &d_nnz,isseqaij ? 0 : nloc, &o_nnz);CHKERRQ(ierr);
996618991cSMark Adams 
1006618991cSMark Adams     if (isseqaij) {
1016618991cSMark Adams       PetscInt       max_d_nnz;
1026618991cSMark Adams 
1036618991cSMark Adams       /*
1046618991cSMark Adams           Determine exact preallocation count for (sequential) scalar matrix
1056618991cSMark Adams       */
1066618991cSMark Adams       ierr = MatSeqAIJGetMaxRowNonzeros(Amat,&max_d_nnz);CHKERRQ(ierr);
1076618991cSMark Adams       max_d_nnz = PetscMin(nloc,bs*max_d_nnz);CHKERRQ(ierr);
1086618991cSMark Adams       ierr = PetscMalloc3(max_d_nnz, &w0,max_d_nnz, &w1,max_d_nnz, &w2);CHKERRQ(ierr);
1096618991cSMark Adams       for (Ii = 0, jj = 0; Ii < Iend; Ii += bs, jj++) {
1106618991cSMark Adams         ierr = MatCollapseRows(Amat,Ii,bs,w0,w1,w2,&d_nnz[jj],NULL);CHKERRQ(ierr);
1116618991cSMark Adams       }
1126618991cSMark Adams       ierr = PetscFree3(w0,w1,w2);CHKERRQ(ierr);
1136618991cSMark Adams 
1146618991cSMark Adams     } else if (ismpiaij) {
1156618991cSMark Adams       Mat            Daij,Oaij;
1166618991cSMark Adams       const PetscInt *garray;
1176618991cSMark Adams       PetscInt       max_d_nnz;
1186618991cSMark Adams 
1196618991cSMark Adams       ierr = MatMPIAIJGetSeqAIJ(Amat,&Daij,&Oaij,&garray);CHKERRQ(ierr);
1206618991cSMark Adams 
1216618991cSMark Adams       /*
1226618991cSMark Adams           Determine exact preallocation count for diagonal block portion of scalar matrix
1236618991cSMark Adams       */
1246618991cSMark Adams       ierr = MatSeqAIJGetMaxRowNonzeros(Daij,&max_d_nnz);CHKERRQ(ierr);
1256618991cSMark Adams       max_d_nnz = PetscMin(nloc,bs*max_d_nnz);CHKERRQ(ierr);
1266618991cSMark Adams       ierr = PetscMalloc3(max_d_nnz, &w0,max_d_nnz, &w1,max_d_nnz, &w2);CHKERRQ(ierr);
1276618991cSMark Adams       for (Ii = 0, jj = 0; Ii < Iend - Istart; Ii += bs, jj++) {
1286618991cSMark Adams         ierr = MatCollapseRows(Daij,Ii,bs,w0,w1,w2,&d_nnz[jj],NULL);CHKERRQ(ierr);
1296618991cSMark Adams       }
1306618991cSMark Adams       ierr = PetscFree3(w0,w1,w2);CHKERRQ(ierr);
1316618991cSMark Adams 
1326618991cSMark Adams       /*
1336618991cSMark Adams          Over estimate (usually grossly over), preallocation count for off-diagonal portion of scalar matrix
1346618991cSMark Adams       */
1356618991cSMark Adams       for (Ii = 0, jj = 0; Ii < Iend - Istart; Ii += bs, jj++) {
1366618991cSMark Adams         o_nnz[jj] = 0;
1376618991cSMark Adams         for (kk=0; kk<bs; kk++) { /* rows that get collapsed to a single row */
1386618991cSMark Adams           ierr = MatGetRow(Oaij,Ii+kk,&ncols,0,0);CHKERRQ(ierr);
1396618991cSMark Adams           o_nnz[jj] += ncols;
1406618991cSMark Adams           ierr = MatRestoreRow(Oaij,Ii+kk,&ncols,0,0);CHKERRQ(ierr);
1416618991cSMark Adams         }
1426618991cSMark Adams         if (o_nnz[jj] > (NN/bs-nloc)) o_nnz[jj] = NN/bs-nloc;
1436618991cSMark Adams       }
1446618991cSMark Adams 
1456618991cSMark Adams     } else {
146359038b3SMark Adams       SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Require AIJ matrix type");
1476618991cSMark Adams     }
1486618991cSMark Adams 
149359038b3SMark Adams     /* get scalar copy (norms) of matrix */
1506618991cSMark Adams     ierr = MatGetType(Amat,&mtype);CHKERRQ(ierr);
1516618991cSMark Adams     ierr = MatCreate(comm, &Gmat);CHKERRQ(ierr);
1526618991cSMark Adams     ierr = MatSetSizes(Gmat,nloc,nloc,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
1536618991cSMark Adams     ierr = MatSetBlockSizes(Gmat, 1, 1);CHKERRQ(ierr);
1546618991cSMark Adams     ierr = MatSetType(Gmat, mtype);CHKERRQ(ierr);
1556618991cSMark Adams     ierr = MatSeqAIJSetPreallocation(Gmat,0,d_nnz);CHKERRQ(ierr);
1566618991cSMark Adams     ierr = MatMPIAIJSetPreallocation(Gmat,0,d_nnz,0,o_nnz);CHKERRQ(ierr);
1576618991cSMark Adams     ierr = PetscFree2(d_nnz,o_nnz);CHKERRQ(ierr);
1586618991cSMark Adams 
1596618991cSMark Adams     for (Ii = Istart; Ii < Iend; Ii++) {
1606618991cSMark Adams       PetscInt dest_row = Ii/bs;
1616618991cSMark Adams       ierr = MatGetRow(Amat,Ii,&ncols,&idx,&vals);CHKERRQ(ierr);
1626618991cSMark Adams       for (jj=0; jj<ncols; jj++) {
1636618991cSMark Adams         PetscInt    dest_col = idx[jj]/bs;
1646618991cSMark Adams         PetscScalar sv       = PetscAbs(PetscRealPart(vals[jj]));
1656618991cSMark Adams         ierr = MatSetValues(Gmat,1,&dest_row,1,&dest_col,&sv,ADD_VALUES);CHKERRQ(ierr);
1666618991cSMark Adams       }
1676618991cSMark Adams       ierr = MatRestoreRow(Amat,Ii,&ncols,&idx,&vals);CHKERRQ(ierr);
1686618991cSMark Adams     }
1696618991cSMark Adams     ierr = MatAssemblyBegin(Gmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1706618991cSMark Adams     ierr = MatAssemblyEnd(Gmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1716618991cSMark Adams   } else {
1726618991cSMark Adams     /* just copy scalar matrix - abs() not taken here but scaled later */
1736618991cSMark Adams     ierr = MatDuplicate(Amat, MAT_COPY_VALUES, &Gmat);CHKERRQ(ierr);
1746618991cSMark Adams   }
1756618991cSMark Adams 
1766618991cSMark Adams #if defined PETSC_GAMG_USE_LOG
1776618991cSMark Adams   ierr = PetscLogEventEnd(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr);
1786618991cSMark Adams #endif
1796618991cSMark Adams 
1806618991cSMark Adams   *a_Gmat = Gmat;
1816618991cSMark Adams   PetscFunctionReturn(0);
1826618991cSMark Adams }
1836618991cSMark Adams 
1846618991cSMark Adams /* -------------------------------------------------------------------------- */
185a37438d7SBarry Smith /*@C
186a37438d7SBarry Smith    PCGAMGFilterGraph - filter (remove zero and possibly small values from the) graph and make it symmetric if requested
187a37438d7SBarry Smith 
188a37438d7SBarry Smith    Collective on Mat
1896618991cSMark Adams 
1906618991cSMark Adams    Input Parameter:
191a37438d7SBarry Smith +   a_Gmat - the graph
1926618991cSMark Adams .   vfilter - threshold paramter [0,1)
193a37438d7SBarry Smith -   symm - make the result symmetric
194a37438d7SBarry Smith 
195a37438d7SBarry Smith    Level: developer
196a37438d7SBarry Smith 
19795452b02SPatrick Sanan    Notes:
19895452b02SPatrick Sanan     This is called before graph coarsers are called.
199a37438d7SBarry Smith 
200a37438d7SBarry Smith .seealso: PCGAMGSetThreshold()
201a37438d7SBarry Smith @*/
2026618991cSMark Adams PetscErrorCode PCGAMGFilterGraph(Mat *a_Gmat,PetscReal vfilter,PetscBool symm)
2036618991cSMark Adams {
2046618991cSMark Adams   PetscErrorCode    ierr;
2056618991cSMark Adams   PetscInt          Istart,Iend,Ii,jj,ncols,nnz0,nnz1, NN, MM, nloc;
2066618991cSMark Adams   PetscMPIInt       rank;
2076618991cSMark Adams   Mat               Gmat  = *a_Gmat, tGmat, matTrans;
2086618991cSMark Adams   MPI_Comm          comm;
2096618991cSMark Adams   const PetscScalar *vals;
2106618991cSMark Adams   const PetscInt    *idx;
2116618991cSMark Adams   PetscInt          *d_nnz, *o_nnz;
2126618991cSMark Adams   Vec               diag;
2136618991cSMark Adams   MatType           mtype;
2146618991cSMark Adams 
2156618991cSMark Adams   PetscFunctionBegin;
2166618991cSMark Adams #if defined PETSC_GAMG_USE_LOG
2176618991cSMark Adams   ierr = PetscLogEventBegin(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr);
2186618991cSMark Adams #endif
2196618991cSMark Adams   /* scale Gmat for all values between -1 and 1 */
2206618991cSMark Adams   ierr = MatCreateVecs(Gmat, &diag, 0);CHKERRQ(ierr);
2216618991cSMark Adams   ierr = MatGetDiagonal(Gmat, diag);CHKERRQ(ierr);
2226618991cSMark Adams   ierr = VecReciprocal(diag);CHKERRQ(ierr);
2236618991cSMark Adams   ierr = VecSqrtAbs(diag);CHKERRQ(ierr);
2246618991cSMark Adams   ierr = MatDiagonalScale(Gmat, diag, diag);CHKERRQ(ierr);
2256618991cSMark Adams   ierr = VecDestroy(&diag);CHKERRQ(ierr);
2266618991cSMark Adams 
2276618991cSMark Adams   if (vfilter < 0.0 && !symm) {
2286618991cSMark Adams     /* Just use the provided matrix as the graph but make all values positive */
2296618991cSMark Adams     MatInfo     info;
2306618991cSMark Adams     PetscScalar *avals;
231359038b3SMark Adams     PetscBool isaij,ismpiaij;
2324099cc6bSBarry Smith     ierr = PetscObjectBaseTypeCompare((PetscObject)Gmat,MATSEQAIJ,&isaij);CHKERRQ(ierr);
2334099cc6bSBarry Smith     ierr = PetscObjectBaseTypeCompare((PetscObject)Gmat,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr);
234359038b3SMark Adams     if (!isaij && !ismpiaij) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Require (MPI)AIJ matrix type");
235359038b3SMark Adams     if (isaij) {
2366618991cSMark Adams       ierr = MatGetInfo(Gmat,MAT_LOCAL,&info);CHKERRQ(ierr);
2376618991cSMark Adams       ierr = MatSeqAIJGetArray(Gmat,&avals);CHKERRQ(ierr);
2386618991cSMark Adams       for (jj = 0; jj<info.nz_used; jj++) avals[jj] = PetscAbsScalar(avals[jj]);
2396618991cSMark Adams       ierr = MatSeqAIJRestoreArray(Gmat,&avals);CHKERRQ(ierr);
2406618991cSMark Adams     } else {
241359038b3SMark Adams       Mat_MPIAIJ  *aij = (Mat_MPIAIJ*)Gmat->data;
2426618991cSMark Adams       ierr = MatGetInfo(aij->A,MAT_LOCAL,&info);CHKERRQ(ierr);
2436618991cSMark Adams       ierr = MatSeqAIJGetArray(aij->A,&avals);CHKERRQ(ierr);
2446618991cSMark Adams       for (jj = 0; jj<info.nz_used; jj++) avals[jj] = PetscAbsScalar(avals[jj]);
2456618991cSMark Adams       ierr = MatSeqAIJRestoreArray(aij->A,&avals);CHKERRQ(ierr);
2466618991cSMark Adams       ierr = MatGetInfo(aij->B,MAT_LOCAL,&info);CHKERRQ(ierr);
2476618991cSMark Adams       ierr = MatSeqAIJGetArray(aij->B,&avals);CHKERRQ(ierr);
2486618991cSMark Adams       for (jj = 0; jj<info.nz_used; jj++) avals[jj] = PetscAbsScalar(avals[jj]);
2496618991cSMark Adams       ierr = MatSeqAIJRestoreArray(aij->B,&avals);CHKERRQ(ierr);
2506618991cSMark Adams     }
2516618991cSMark Adams #if defined PETSC_GAMG_USE_LOG
2526618991cSMark Adams     ierr = PetscLogEventEnd(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr);
2536618991cSMark Adams #endif
2546618991cSMark Adams     PetscFunctionReturn(0);
2556618991cSMark Adams   }
2566618991cSMark Adams 
2576618991cSMark Adams   ierr = PetscObjectGetComm((PetscObject)Gmat,&comm);CHKERRQ(ierr);
2586618991cSMark Adams   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
2596618991cSMark Adams   ierr = MatGetOwnershipRange(Gmat, &Istart, &Iend);CHKERRQ(ierr);
2606618991cSMark Adams   nloc = Iend - Istart;
2616618991cSMark Adams   ierr = MatGetSize(Gmat, &MM, &NN);CHKERRQ(ierr);
2626618991cSMark Adams 
2636618991cSMark Adams   if (symm) {
2646618991cSMark Adams     ierr = MatTranspose(Gmat, MAT_INITIAL_MATRIX, &matTrans);CHKERRQ(ierr);
2656618991cSMark Adams   }
2666618991cSMark Adams 
2676618991cSMark Adams   /* Determine upper bound on nonzeros needed in new filtered matrix */
2686618991cSMark Adams   ierr = PetscMalloc2(nloc, &d_nnz,nloc, &o_nnz);CHKERRQ(ierr);
2696618991cSMark Adams   for (Ii = Istart, jj = 0; Ii < Iend; Ii++, jj++) {
2706618991cSMark Adams     ierr      = MatGetRow(Gmat,Ii,&ncols,NULL,NULL);CHKERRQ(ierr);
2716618991cSMark Adams     d_nnz[jj] = ncols;
2726618991cSMark Adams     o_nnz[jj] = ncols;
2736618991cSMark Adams     ierr      = MatRestoreRow(Gmat,Ii,&ncols,NULL,NULL);CHKERRQ(ierr);
2746618991cSMark Adams     if (symm) {
2756618991cSMark Adams       ierr       = MatGetRow(matTrans,Ii,&ncols,NULL,NULL);CHKERRQ(ierr);
2766618991cSMark Adams       d_nnz[jj] += ncols;
2776618991cSMark Adams       o_nnz[jj] += ncols;
2786618991cSMark Adams       ierr       = MatRestoreRow(matTrans,Ii,&ncols,NULL,NULL);CHKERRQ(ierr);
2796618991cSMark Adams     }
2806618991cSMark Adams     if (d_nnz[jj] > nloc) d_nnz[jj] = nloc;
2816618991cSMark Adams     if (o_nnz[jj] > (MM-nloc)) o_nnz[jj] = MM - nloc;
2826618991cSMark Adams   }
2836618991cSMark Adams   ierr = MatGetType(Gmat,&mtype);CHKERRQ(ierr);
2846618991cSMark Adams   ierr = MatCreate(comm, &tGmat);CHKERRQ(ierr);
2856618991cSMark Adams   ierr = MatSetSizes(tGmat,nloc,nloc,MM,MM);CHKERRQ(ierr);
2866618991cSMark Adams   ierr = MatSetBlockSizes(tGmat, 1, 1);CHKERRQ(ierr);
2876618991cSMark Adams   ierr = MatSetType(tGmat, mtype);CHKERRQ(ierr);
2886618991cSMark Adams   ierr = MatSeqAIJSetPreallocation(tGmat,0,d_nnz);CHKERRQ(ierr);
2896618991cSMark Adams   ierr = MatMPIAIJSetPreallocation(tGmat,0,d_nnz,0,o_nnz);CHKERRQ(ierr);
2906618991cSMark Adams   ierr = PetscFree2(d_nnz,o_nnz);CHKERRQ(ierr);
2916618991cSMark Adams   if (symm) {
2926618991cSMark Adams     ierr = MatDestroy(&matTrans);CHKERRQ(ierr);
2936618991cSMark Adams   } else {
2946618991cSMark Adams     /* all entries are generated locally so MatAssembly will be slightly faster for large process counts */
2956618991cSMark Adams     ierr = MatSetOption(tGmat,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
2966618991cSMark Adams   }
2976618991cSMark Adams 
2986618991cSMark Adams   for (Ii = Istart, nnz0 = nnz1 = 0; Ii < Iend; Ii++) {
2996618991cSMark Adams     ierr = MatGetRow(Gmat,Ii,&ncols,&idx,&vals);CHKERRQ(ierr);
3006618991cSMark Adams     for (jj=0; jj<ncols; jj++,nnz0++) {
3016618991cSMark Adams       PetscScalar sv = PetscAbs(PetscRealPart(vals[jj]));
3026618991cSMark Adams       if (PetscRealPart(sv) > vfilter) {
3036618991cSMark Adams         nnz1++;
3046618991cSMark Adams         if (symm) {
3056618991cSMark Adams           sv  *= 0.5;
3066618991cSMark Adams           ierr = MatSetValues(tGmat,1,&Ii,1,&idx[jj],&sv,ADD_VALUES);CHKERRQ(ierr);
3076618991cSMark Adams           ierr = MatSetValues(tGmat,1,&idx[jj],1,&Ii,&sv,ADD_VALUES);CHKERRQ(ierr);
3086618991cSMark Adams         } else {
3096618991cSMark Adams           ierr = MatSetValues(tGmat,1,&Ii,1,&idx[jj],&sv,ADD_VALUES);CHKERRQ(ierr);
3106618991cSMark Adams         }
3116618991cSMark Adams       }
3126618991cSMark Adams     }
3136618991cSMark Adams     ierr = MatRestoreRow(Gmat,Ii,&ncols,&idx,&vals);CHKERRQ(ierr);
3146618991cSMark Adams   }
3156618991cSMark Adams   ierr = MatAssemblyBegin(tGmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
3166618991cSMark Adams   ierr = MatAssemblyEnd(tGmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
3176618991cSMark Adams 
3186618991cSMark Adams #if defined PETSC_GAMG_USE_LOG
3196618991cSMark Adams   ierr = PetscLogEventEnd(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr);
3206618991cSMark Adams #endif
3216618991cSMark Adams 
3226618991cSMark Adams #if defined(PETSC_USE_INFO)
3236618991cSMark Adams   {
3246618991cSMark Adams     double t1 = (!nnz0) ? 1. : 100.*(double)nnz1/(double)nnz0, t2 = (!nloc) ? 1. : (double)nnz0/(double)nloc;
3256618991cSMark Adams     ierr = PetscInfo4(*a_Gmat,"\t %g%% nnz after filtering, with threshold %g, %g nnz ave. (N=%D)\n",t1,vfilter,t2,MM);CHKERRQ(ierr);
3266618991cSMark Adams   }
3276618991cSMark Adams #endif
3286618991cSMark Adams   ierr    = MatDestroy(&Gmat);CHKERRQ(ierr);
3296618991cSMark Adams   *a_Gmat = tGmat;
3306618991cSMark Adams   PetscFunctionReturn(0);
3316618991cSMark Adams }
3326618991cSMark Adams 
3336618991cSMark Adams /* -------------------------------------------------------------------------- */
3346618991cSMark Adams /*
335359038b3SMark Adams    PCGAMGGetDataWithGhosts - hacks into Mat MPIAIJ so this must have npe > 1
3366618991cSMark Adams 
3376618991cSMark Adams    Input Parameter:
3386618991cSMark Adams    . Gmat - MPIAIJ matrix for scattters
3396618991cSMark Adams    . data_sz - number of data terms per node (# cols in output)
3406618991cSMark Adams    . data_in[nloc*data_sz] - column oriented data
3416618991cSMark Adams    Output Parameter:
3426618991cSMark Adams    . a_stride - numbrt of rows of output
3436618991cSMark Adams    . a_data_out[stride*data_sz] - output data with ghosts
3446618991cSMark Adams */
3456618991cSMark Adams PetscErrorCode PCGAMGGetDataWithGhosts(Mat Gmat,PetscInt data_sz,PetscReal data_in[],PetscInt *a_stride,PetscReal **a_data_out)
3466618991cSMark Adams {
3476618991cSMark Adams   PetscErrorCode ierr;
3486618991cSMark Adams   MPI_Comm       comm;
3496618991cSMark Adams   Vec            tmp_crds;
3506618991cSMark Adams   Mat_MPIAIJ     *mpimat = (Mat_MPIAIJ*)Gmat->data;
3516618991cSMark Adams   PetscInt       nnodes,num_ghosts,dir,kk,jj,my0,Iend,nloc;
3526618991cSMark Adams   PetscScalar    *data_arr;
3536618991cSMark Adams   PetscReal      *datas;
3546618991cSMark Adams   PetscBool      isMPIAIJ;
3556618991cSMark Adams 
3566618991cSMark Adams   PetscFunctionBegin;
3576618991cSMark Adams   ierr = PetscObjectGetComm((PetscObject)Gmat,&comm);CHKERRQ(ierr);
3584099cc6bSBarry Smith   ierr      = PetscObjectBaseTypeCompare((PetscObject)Gmat, MATMPIAIJ, &isMPIAIJ);CHKERRQ(ierr);
3596618991cSMark Adams   ierr      = MatGetOwnershipRange(Gmat, &my0, &Iend);CHKERRQ(ierr);
3606618991cSMark Adams   nloc      = Iend - my0;
3616618991cSMark Adams   ierr      = VecGetLocalSize(mpimat->lvec, &num_ghosts);CHKERRQ(ierr);
3626618991cSMark Adams   nnodes    = num_ghosts + nloc;
3636618991cSMark Adams   *a_stride = nnodes;
3646618991cSMark Adams   ierr      = MatCreateVecs(Gmat, &tmp_crds, 0);CHKERRQ(ierr);
3656618991cSMark Adams 
3666618991cSMark Adams   ierr = PetscMalloc1(data_sz*nnodes, &datas);CHKERRQ(ierr);
3676618991cSMark Adams   for (dir=0; dir<data_sz; dir++) {
3686618991cSMark Adams     /* set local, and global */
3696618991cSMark Adams     for (kk=0; kk<nloc; kk++) {
3706618991cSMark Adams       PetscInt    gid = my0 + kk;
3716618991cSMark Adams       PetscScalar crd = (PetscScalar)data_in[dir*nloc + kk]; /* col oriented */
3726618991cSMark Adams       datas[dir*nnodes + kk] = PetscRealPart(crd);
3736618991cSMark Adams 
3746618991cSMark Adams       ierr = VecSetValues(tmp_crds, 1, &gid, &crd, INSERT_VALUES);CHKERRQ(ierr);
3756618991cSMark Adams     }
3766618991cSMark Adams     ierr = VecAssemblyBegin(tmp_crds);CHKERRQ(ierr);
3776618991cSMark Adams     ierr = VecAssemblyEnd(tmp_crds);CHKERRQ(ierr);
3786618991cSMark Adams     /* get ghost datas */
3796618991cSMark Adams     ierr = VecScatterBegin(mpimat->Mvctx,tmp_crds,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
3806618991cSMark Adams     ierr = VecScatterEnd(mpimat->Mvctx,tmp_crds,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
3816618991cSMark Adams     ierr = VecGetArray(mpimat->lvec, &data_arr);CHKERRQ(ierr);
3826618991cSMark Adams     for (kk=nloc,jj=0;jj<num_ghosts;kk++,jj++) datas[dir*nnodes + kk] = PetscRealPart(data_arr[jj]);
3836618991cSMark Adams     ierr = VecRestoreArray(mpimat->lvec, &data_arr);CHKERRQ(ierr);
3846618991cSMark Adams   }
3856618991cSMark Adams   ierr        = VecDestroy(&tmp_crds);CHKERRQ(ierr);
3866618991cSMark Adams   *a_data_out = datas;
3876618991cSMark Adams   PetscFunctionReturn(0);
3886618991cSMark Adams }
3896618991cSMark Adams 
3901943db53SBarry Smith PetscErrorCode PCGAMGHashTableCreate(PetscInt a_size, PCGAMGHashTable *a_tab)
3916618991cSMark Adams {
3926618991cSMark Adams   PetscErrorCode ierr;
3936618991cSMark Adams   PetscInt       kk;
3946618991cSMark Adams 
3956618991cSMark Adams   PetscFunctionBegin;
3966618991cSMark Adams   a_tab->size = a_size;
397*8f3cd775SBarry Smith   ierr = PetscMalloc2(a_size, &a_tab->table,a_size, &a_tab->data);CHKERRQ(ierr);
3986618991cSMark Adams   for (kk=0; kk<a_size; kk++) a_tab->table[kk] = -1;
3996618991cSMark Adams   PetscFunctionReturn(0);
4006618991cSMark Adams }
4016618991cSMark Adams 
4021943db53SBarry Smith PetscErrorCode PCGAMGHashTableDestroy(PCGAMGHashTable *a_tab)
4036618991cSMark Adams {
4046618991cSMark Adams   PetscErrorCode ierr;
4056618991cSMark Adams 
4066618991cSMark Adams   PetscFunctionBegin;
407*8f3cd775SBarry Smith   ierr = PetscFree2(a_tab->table,a_tab->data);CHKERRQ(ierr);
4086618991cSMark Adams   PetscFunctionReturn(0);
4096618991cSMark Adams }
4106618991cSMark Adams 
4111943db53SBarry Smith PetscErrorCode PCGAMGHashTableAdd(PCGAMGHashTable *a_tab, PetscInt a_key, PetscInt a_data)
4126618991cSMark Adams {
4136618991cSMark Adams   PetscInt kk,idx;
4146618991cSMark Adams 
4156618991cSMark Adams   PetscFunctionBegin;
416*8f3cd775SBarry Smith   if (a_key<0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Negative key %D.",a_key);
417*8f3cd775SBarry Smith   for (kk = 0, idx = GAMG_HASH(a_key); kk < a_tab->size; kk++, idx = (idx==(a_tab->size-1)) ? 0 : idx + 1) {
4186618991cSMark Adams     if (a_tab->table[idx] == a_key) {
4196618991cSMark Adams       /* exists */
4206618991cSMark Adams       a_tab->data[idx] = a_data;
4216618991cSMark Adams       break;
4226618991cSMark Adams     } else if (a_tab->table[idx] == -1) {
4236618991cSMark Adams       /* add */
4246618991cSMark Adams       a_tab->table[idx] = a_key;
4256618991cSMark Adams       a_tab->data[idx]  = a_data;
4266618991cSMark Adams       break;
4276618991cSMark Adams     }
4286618991cSMark Adams   }
4296618991cSMark Adams   if (kk==a_tab->size) {
4306618991cSMark Adams     /* this is not to efficient, waiting until completely full */
4316618991cSMark Adams     PetscInt       oldsize = a_tab->size, new_size = 2*a_tab->size + 5, *oldtable = a_tab->table, *olddata = a_tab->data;
4326618991cSMark Adams     PetscErrorCode ierr;
4336618991cSMark Adams 
4346618991cSMark Adams     a_tab->size = new_size;
435*8f3cd775SBarry Smith     ierr = PetscMalloc2(a_tab->size, &a_tab->table,a_tab->size, &a_tab->data);CHKERRQ(ierr);
4366618991cSMark Adams     for (kk=0;kk<a_tab->size;kk++) a_tab->table[kk] = -1;
4376618991cSMark Adams     for (kk=0;kk<oldsize;kk++) {
4386618991cSMark Adams       if (oldtable[kk] != -1) {
4391943db53SBarry Smith         ierr = PCGAMGHashTableAdd(a_tab, oldtable[kk], olddata[kk]);CHKERRQ(ierr);
4406618991cSMark Adams        }
4416618991cSMark Adams     }
442*8f3cd775SBarry Smith     ierr = PetscFree2(oldtable,olddata);CHKERRQ(ierr);
4431943db53SBarry Smith     ierr = PCGAMGHashTableAdd(a_tab, a_key, a_data);CHKERRQ(ierr);
4446618991cSMark Adams   }
4456618991cSMark Adams   PetscFunctionReturn(0);
4466618991cSMark Adams }
4476618991cSMark Adams 
448