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