18e6d0c30SPeter Brune #include <../src/ksp/pc/impls/gamg/gamg.h> /*I "petscpc.h" I*/ 28e6d0c30SPeter Brune #include <petsc-private/kspimpl.h> 38e6d0c30SPeter Brune 48e6d0c30SPeter Brune typedef struct { 58e6d0c30SPeter Brune PetscReal dummy; /* empty struct; save for later */ 68e6d0c30SPeter Brune } PC_GAMG_Classical; 78e6d0c30SPeter Brune 88e6d0c30SPeter Brune 98e6d0c30SPeter Brune #undef __FUNCT__ 10bfde193fSPeter Brune #define __FUNCT__ "PCGAMGClassicalCreateGhostVector_Private" 11bfde193fSPeter Brune PetscErrorCode PCGAMGClassicalCreateGhostVector_Private(Mat G,Vec *gvec,PetscInt **global) 128e6d0c30SPeter Brune { 138e6d0c30SPeter Brune Mat_MPIAIJ *aij = (Mat_MPIAIJ*)G->data; 148e6d0c30SPeter Brune PetscErrorCode ierr; 158e6d0c30SPeter Brune PetscBool isMPIAIJ; 168e6d0c30SPeter Brune 178e6d0c30SPeter Brune PetscFunctionBegin; 188e6d0c30SPeter Brune ierr = PetscObjectTypeCompare((PetscObject)G, MATMPIAIJ, &isMPIAIJ); CHKERRQ(ierr); 198e6d0c30SPeter Brune if (isMPIAIJ) { 208e6d0c30SPeter Brune if (gvec)ierr = VecDuplicate(aij->lvec,gvec);CHKERRQ(ierr); 218e6d0c30SPeter Brune if (global)*global = aij->garray; 228e6d0c30SPeter Brune } else { 238e6d0c30SPeter Brune /* no off-processor nodes */ 248e6d0c30SPeter Brune if (gvec)*gvec = NULL; 258e6d0c30SPeter Brune if (global)*global = NULL; 268e6d0c30SPeter Brune } 278e6d0c30SPeter Brune PetscFunctionReturn(0); 288e6d0c30SPeter Brune } 298e6d0c30SPeter Brune 308e6d0c30SPeter Brune #undef __FUNCT__ 31bfde193fSPeter Brune #define __FUNCT__ "PCGAMGClassicalGraphSplitting_Private" 328e6d0c30SPeter Brune /* 338e6d0c30SPeter Brune Split the relevant graph into diagonal and off-diagonal parts in local numbering; for now this 348e6d0c30SPeter Brune a roundabout private interface to the mats' internal diag and offdiag mats. 358e6d0c30SPeter Brune */ 36bfde193fSPeter Brune PetscErrorCode PCGAMGClassicalGraphSplitting_Private(Mat G,Mat *Gd, Mat *Go) 378e6d0c30SPeter Brune { 388e6d0c30SPeter Brune Mat_MPIAIJ *aij = (Mat_MPIAIJ*)G->data; 398e6d0c30SPeter Brune PetscErrorCode ierr; 408e6d0c30SPeter Brune PetscBool isMPIAIJ; 418e6d0c30SPeter Brune PetscFunctionBegin; 428e6d0c30SPeter Brune ierr = PetscObjectTypeCompare((PetscObject)G, MATMPIAIJ, &isMPIAIJ ); CHKERRQ(ierr); 438e6d0c30SPeter Brune if (isMPIAIJ) { 448e6d0c30SPeter Brune *Gd = aij->A; 458e6d0c30SPeter Brune *Go = aij->B; 468e6d0c30SPeter Brune } else { 478e6d0c30SPeter Brune *Gd = G; 488e6d0c30SPeter Brune *Go = NULL; 498e6d0c30SPeter Brune } 508e6d0c30SPeter Brune PetscFunctionReturn(0); 518e6d0c30SPeter Brune } 528e6d0c30SPeter Brune 538e6d0c30SPeter Brune #undef __FUNCT__ 548e6d0c30SPeter Brune #define __FUNCT__ "PCGAMGGraph_Classical" 5565b3d5b6SPeter Brune PetscErrorCode PCGAMGGraph_Classical(PC pc,const Mat A,Mat *G) 568e6d0c30SPeter Brune { 57550383edSPeter Brune PetscInt s,f,n,idx,lidx,gidx; 58e5a0faa4SPeter Brune PetscInt r,c,ncols; 598e6d0c30SPeter Brune const PetscInt *rcol; 608e6d0c30SPeter Brune const PetscScalar *rval; 61e5a0faa4SPeter Brune PetscInt *gcol; 628e6d0c30SPeter Brune PetscScalar *gval; 63e5a0faa4SPeter Brune PetscReal rmax; 64550383edSPeter Brune PetscInt cmax = 0; 658e6d0c30SPeter Brune PC_MG *mg; 668e6d0c30SPeter Brune PC_GAMG *gamg; 678e6d0c30SPeter Brune PetscErrorCode ierr; 688e6d0c30SPeter Brune PetscInt *gsparse,*lsparse; 69e5a0faa4SPeter Brune PetscScalar *Amax; 708e6d0c30SPeter Brune MatType mtype; 718e6d0c30SPeter Brune 728e6d0c30SPeter Brune PetscFunctionBegin; 738e6d0c30SPeter Brune mg = (PC_MG *)pc->data; 748e6d0c30SPeter Brune gamg = (PC_GAMG *)mg->innerctx; 758e6d0c30SPeter Brune 768e6d0c30SPeter Brune ierr = MatGetOwnershipRange(A,&s,&f);CHKERRQ(ierr); 77550383edSPeter Brune n=f-s; 78550383edSPeter Brune ierr = PetscMalloc(sizeof(PetscInt)*n,&lsparse);CHKERRQ(ierr); 79550383edSPeter Brune ierr = PetscMalloc(sizeof(PetscInt)*n,&gsparse);CHKERRQ(ierr); 80550383edSPeter Brune ierr = PetscMalloc(sizeof(PetscScalar)*n,&Amax);CHKERRQ(ierr); 818e6d0c30SPeter Brune 82550383edSPeter Brune for (r = 0;r < n;r++) { 838e6d0c30SPeter Brune lsparse[r] = 0; 84550383edSPeter Brune gsparse[r] = 0; 858e6d0c30SPeter Brune } 868e6d0c30SPeter Brune 87550383edSPeter Brune for (r = s;r < f;r++) { 88e5a0faa4SPeter Brune /* determine the maximum off-diagonal in each row */ 89e5a0faa4SPeter Brune rmax = 0.; 90550383edSPeter Brune ierr = MatGetRow(A,r,&ncols,&rcol,&rval);CHKERRQ(ierr); 91e5a0faa4SPeter Brune for (c = 0; c < ncols; c++) { 92e5a0faa4SPeter Brune if (PetscAbsScalar(rval[c]) > rmax && rcol[c] != r) { 93e5a0faa4SPeter Brune rmax = PetscAbsScalar(rval[c]); 94e5a0faa4SPeter Brune } 95e5a0faa4SPeter Brune } 96550383edSPeter Brune Amax[r-s] = rmax; 97550383edSPeter Brune if (ncols > cmax) cmax = ncols; 98550383edSPeter Brune lidx = 0; 99550383edSPeter Brune gidx = 0; 100e5a0faa4SPeter Brune /* create the local and global sparsity patterns */ 1018e6d0c30SPeter Brune for (c = 0; c < ncols; c++) { 102550383edSPeter Brune if (PetscAbsScalar(rval[c]) > gamg->threshold*PetscRealPart(Amax[r-s])) { 103550383edSPeter Brune if (rcol[c] < f && rcol[c] >= s) { 104550383edSPeter Brune lidx++; 105550383edSPeter Brune } else { 106550383edSPeter Brune gidx++; 1078e6d0c30SPeter Brune } 1088e6d0c30SPeter Brune } 1098e6d0c30SPeter Brune } 110550383edSPeter Brune ierr = MatRestoreRow(A,r,&ncols,&rcol,&rval);CHKERRQ(ierr); 111550383edSPeter Brune lsparse[r-s] = lidx; 112550383edSPeter Brune gsparse[r-s] = gidx; 1138e6d0c30SPeter Brune } 114e5a0faa4SPeter Brune ierr = PetscMalloc(sizeof(PetscScalar)*cmax,&gval);CHKERRQ(ierr); 115e5a0faa4SPeter Brune ierr = PetscMalloc(sizeof(PetscInt)*cmax,&gcol);CHKERRQ(ierr); 116e5a0faa4SPeter Brune 1178e6d0c30SPeter Brune ierr = MatCreate(PetscObjectComm((PetscObject)A),G); CHKERRQ(ierr); 1188e6d0c30SPeter Brune ierr = MatGetType(A,&mtype);CHKERRQ(ierr); 1198e6d0c30SPeter Brune ierr = MatSetType(*G,mtype);CHKERRQ(ierr); 120550383edSPeter Brune ierr = MatSetSizes(*G,n,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1218e6d0c30SPeter Brune ierr = MatMPIAIJSetPreallocation(*G,0,lsparse,0,gsparse);CHKERRQ(ierr); 1228e6d0c30SPeter Brune ierr = MatSeqAIJSetPreallocation(*G,0,lsparse);CHKERRQ(ierr); 1238e6d0c30SPeter Brune for (r = s;r < f;r++) { 1248e6d0c30SPeter Brune ierr = MatGetRow(A,r,&ncols,&rcol,&rval);CHKERRQ(ierr); 1258e6d0c30SPeter Brune idx = 0; 1268e6d0c30SPeter Brune for (c = 0; c < ncols; c++) { 1278e6d0c30SPeter Brune /* classical strength of connection */ 1286774142eSPeter Brune if (PetscAbsScalar(rval[c]) > gamg->threshold*PetscRealPart(Amax[r-s])) { 1298e6d0c30SPeter Brune gcol[idx] = rcol[c]; 1308e6d0c30SPeter Brune gval[idx] = rval[c]; 1318e6d0c30SPeter Brune idx++; 1328e6d0c30SPeter Brune } 1338e6d0c30SPeter Brune } 1348e6d0c30SPeter Brune ierr = MatSetValues(*G,1,&r,idx,gcol,gval,INSERT_VALUES);CHKERRQ(ierr); 1358e6d0c30SPeter Brune ierr = MatRestoreRow(A,r,&ncols,&rcol,&rval);CHKERRQ(ierr); 1368e6d0c30SPeter Brune } 1378e6d0c30SPeter Brune ierr = MatAssemblyBegin(*G, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1388e6d0c30SPeter Brune ierr = MatAssemblyEnd(*G, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1398e6d0c30SPeter Brune 1408e6d0c30SPeter Brune ierr = PetscFree(gval);CHKERRQ(ierr); 1418e6d0c30SPeter Brune ierr = PetscFree(gcol);CHKERRQ(ierr); 1428e6d0c30SPeter Brune ierr = PetscFree(lsparse);CHKERRQ(ierr); 1438e6d0c30SPeter Brune ierr = PetscFree(gsparse);CHKERRQ(ierr); 144e5a0faa4SPeter Brune ierr = PetscFree(Amax);CHKERRQ(ierr); 1458e6d0c30SPeter Brune PetscFunctionReturn(0); 1468e6d0c30SPeter Brune } 1478e6d0c30SPeter Brune 1488e6d0c30SPeter Brune 1498e6d0c30SPeter Brune #undef __FUNCT__ 1508e6d0c30SPeter Brune #define __FUNCT__ "PCGAMGCoarsen_Classical" 1518e6d0c30SPeter Brune PetscErrorCode PCGAMGCoarsen_Classical(PC pc,Mat *G,PetscCoarsenData **agg_lists) 1528e6d0c30SPeter Brune { 1538e6d0c30SPeter Brune PetscErrorCode ierr; 1548e6d0c30SPeter Brune MatCoarsen crs; 1558e6d0c30SPeter Brune MPI_Comm fcomm = ((PetscObject)pc)->comm; 1568e6d0c30SPeter Brune 1578e6d0c30SPeter Brune PetscFunctionBegin; 1588e6d0c30SPeter Brune 1598e6d0c30SPeter Brune 1608e6d0c30SPeter Brune /* construct the graph if necessary */ 1618e6d0c30SPeter Brune if (!G) { 1628e6d0c30SPeter Brune SETERRQ(fcomm,PETSC_ERR_ARG_WRONGSTATE,"Must set Graph in PC in PCGAMG before coarsening"); 1638e6d0c30SPeter Brune } 1648e6d0c30SPeter Brune 1658e6d0c30SPeter Brune ierr = MatCoarsenCreate(fcomm,&crs);CHKERRQ(ierr); 1668e6d0c30SPeter Brune ierr = MatCoarsenSetFromOptions(crs);CHKERRQ(ierr); 1678e6d0c30SPeter Brune ierr = MatCoarsenSetAdjacency(crs,*G);CHKERRQ(ierr); 1688e6d0c30SPeter Brune ierr = MatCoarsenSetStrictAggs(crs,PETSC_TRUE);CHKERRQ(ierr); 1698e6d0c30SPeter Brune ierr = MatCoarsenApply(crs);CHKERRQ(ierr); 1708e6d0c30SPeter Brune ierr = MatCoarsenGetData(crs,agg_lists);CHKERRQ(ierr); 1718e6d0c30SPeter Brune ierr = MatCoarsenDestroy(&crs);CHKERRQ(ierr); 1728e6d0c30SPeter Brune 1738e6d0c30SPeter Brune PetscFunctionReturn(0); 1748e6d0c30SPeter Brune } 1758e6d0c30SPeter Brune 1768e6d0c30SPeter Brune #undef __FUNCT__ 177bfde193fSPeter Brune #define __FUNCT__ "PCGAMGClassicalGhost_Private" 1788e6d0c30SPeter Brune /* 1798e6d0c30SPeter Brune Find all ghost nodes that are coarse and output the fine/coarse splitting for those as well 1808e6d0c30SPeter Brune 1818e6d0c30SPeter Brune Input: 1828e6d0c30SPeter Brune G - graph; 1838e6d0c30SPeter Brune gvec - Global Vector 1848e6d0c30SPeter Brune avec - Local part of the scattered vec 1858e6d0c30SPeter Brune bvec - Global part of the scattered vec 1868e6d0c30SPeter Brune 1878e6d0c30SPeter Brune Output: 1888e6d0c30SPeter Brune findx - indirection t 1898e6d0c30SPeter Brune 1908e6d0c30SPeter Brune */ 191bfde193fSPeter Brune PetscErrorCode PCGAMGClassicalGhost_Private(Mat G,Vec v,Vec gv) 1928e6d0c30SPeter Brune { 1938e6d0c30SPeter Brune PetscErrorCode ierr; 1948e6d0c30SPeter Brune Mat_MPIAIJ *aij = (Mat_MPIAIJ*)G->data; 1958e6d0c30SPeter Brune PetscBool isMPIAIJ; 1968e6d0c30SPeter Brune 1978e6d0c30SPeter Brune PetscFunctionBegin; 1988e6d0c30SPeter Brune ierr = PetscObjectTypeCompare((PetscObject)G, MATMPIAIJ, &isMPIAIJ ); CHKERRQ(ierr); 1998e6d0c30SPeter Brune if (isMPIAIJ) { 2008e6d0c30SPeter Brune ierr = VecScatterBegin(aij->Mvctx,v,gv,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2018e6d0c30SPeter Brune ierr = VecScatterEnd(aij->Mvctx,v,gv,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2028e6d0c30SPeter Brune } 2038e6d0c30SPeter Brune PetscFunctionReturn(0); 2048e6d0c30SPeter Brune } 2058e6d0c30SPeter Brune 2068e6d0c30SPeter Brune #undef __FUNCT__ 2078e6d0c30SPeter Brune #define __FUNCT__ "PCGAMGProlongator_Classical" 20865b3d5b6SPeter Brune PetscErrorCode PCGAMGProlongator_Classical(PC pc, const Mat A, const Mat G, PetscCoarsenData *agg_lists,Mat *P) 2098e6d0c30SPeter Brune { 2108e6d0c30SPeter Brune PetscErrorCode ierr; 2118e6d0c30SPeter Brune MPI_Comm comm; 2128e6d0c30SPeter Brune Mat lG,gG,lA,gA; /* on and off diagonal matrices */ 2138e6d0c30SPeter Brune PetscInt fn; /* fine local blocked sizes */ 2148e6d0c30SPeter Brune PetscInt cn; /* coarse local blocked sizes */ 2158e6d0c30SPeter Brune PetscInt gn; /* size of the off-diagonal fine vector */ 2168e6d0c30SPeter Brune PetscInt fs,fe; /* fine (row) ownership range*/ 2178e6d0c30SPeter Brune PetscInt cs,ce; /* coarse (column) ownership range */ 2188e6d0c30SPeter Brune PetscInt i,j,k; /* indices! */ 2198e6d0c30SPeter Brune PetscBool iscoarse; /* flag for determining if a node is coarse */ 2208e6d0c30SPeter Brune PetscInt *lcid,*gcid; /* on and off-processor coarse unknown IDs */ 2218e6d0c30SPeter Brune PetscInt *lsparse,*gsparse; /* on and off-processor sparsity patterns for prolongator */ 2228e6d0c30SPeter Brune PetscScalar pij; 2238e6d0c30SPeter Brune const PetscScalar *rval; 2248e6d0c30SPeter Brune const PetscInt *rcol; 2258e6d0c30SPeter Brune PetscScalar g_pos,g_neg,a_pos,a_neg,diag,invdiag,alpha,beta; 2268e6d0c30SPeter Brune Vec F; /* vec of coarse size */ 2278e6d0c30SPeter Brune Vec C; /* vec of fine size */ 2288e6d0c30SPeter Brune Vec gF; /* vec of off-diagonal fine size */ 2298e6d0c30SPeter Brune MatType mtype; 2308e6d0c30SPeter Brune PetscInt c_indx; 2318e6d0c30SPeter Brune const PetscScalar *vcols; 2328e6d0c30SPeter Brune const PetscInt *icols; 2338e6d0c30SPeter Brune PetscScalar c_scalar; 2348e6d0c30SPeter Brune PetscInt ncols,col; 2358e6d0c30SPeter Brune PetscInt row_f,row_c; 2368e6d0c30SPeter Brune PetscInt cmax=0,ncolstotal,idx; 2378e6d0c30SPeter Brune PetscScalar *pvals; 2388e6d0c30SPeter Brune PetscInt *pcols; 2398e6d0c30SPeter Brune 2408e6d0c30SPeter Brune PetscFunctionBegin; 2418e6d0c30SPeter Brune comm = ((PetscObject)pc)->comm; 2428e6d0c30SPeter Brune ierr = MatGetOwnershipRange(A,&fs,&fe); CHKERRQ(ierr); 2438e6d0c30SPeter Brune fn = (fe - fs); 2448e6d0c30SPeter Brune 2458e6d0c30SPeter Brune ierr = MatGetVecs(A,&F,NULL);CHKERRQ(ierr); 2468e6d0c30SPeter Brune 2478e6d0c30SPeter Brune /* get the number of local unknowns and the indices of the local unknowns */ 2488e6d0c30SPeter Brune 2498e6d0c30SPeter Brune ierr = PetscMalloc(sizeof(PetscInt)*fn,&lsparse);CHKERRQ(ierr); 2508e6d0c30SPeter Brune ierr = PetscMalloc(sizeof(PetscInt)*fn,&gsparse);CHKERRQ(ierr); 2518e6d0c30SPeter Brune ierr = PetscMalloc(sizeof(PetscInt)*fn,&lcid);CHKERRQ(ierr); 2528e6d0c30SPeter Brune 2538e6d0c30SPeter Brune /* count the number of coarse unknowns */ 2548e6d0c30SPeter Brune cn = 0; 2558e6d0c30SPeter Brune for (i=0;i<fn;i++) { 2568e6d0c30SPeter Brune /* filter out singletons */ 2578e6d0c30SPeter Brune ierr = PetscCDEmptyAt(agg_lists,i,&iscoarse); CHKERRQ(ierr); 2588e6d0c30SPeter Brune lcid[i] = -1; 2598e6d0c30SPeter Brune if (!iscoarse) { 2608e6d0c30SPeter Brune cn++; 2618e6d0c30SPeter Brune } 2628e6d0c30SPeter Brune } 2638e6d0c30SPeter Brune 2648e6d0c30SPeter Brune /* create the coarse vector */ 2658e6d0c30SPeter Brune ierr = VecCreateMPI(comm,cn,PETSC_DECIDE,&C);CHKERRQ(ierr); 2668e6d0c30SPeter Brune ierr = VecGetOwnershipRange(C,&cs,&ce);CHKERRQ(ierr); 2678e6d0c30SPeter Brune 2688e6d0c30SPeter Brune /* construct a global vector indicating the global indices of the coarse unknowns */ 2698e6d0c30SPeter Brune cn = 0; 2708e6d0c30SPeter Brune for (i=0;i<fn;i++) { 2718e6d0c30SPeter Brune ierr = PetscCDEmptyAt(agg_lists,i,&iscoarse); CHKERRQ(ierr); 2728e6d0c30SPeter Brune if (!iscoarse) { 2738e6d0c30SPeter Brune lcid[i] = cs+cn; 2748e6d0c30SPeter Brune cn++; 2758e6d0c30SPeter Brune } else { 2768e6d0c30SPeter Brune lcid[i] = -1; 2778e6d0c30SPeter Brune } 278*167fb786SPeter Brune *((PetscInt *)&c_scalar) = lcid[i]; 2798e6d0c30SPeter Brune c_indx = fs+i; 2808e6d0c30SPeter Brune ierr = VecSetValues(F,1,&c_indx,&c_scalar,INSERT_VALUES);CHKERRQ(ierr); 2818e6d0c30SPeter Brune } 2828e6d0c30SPeter Brune 2838e6d0c30SPeter Brune ierr = VecAssemblyBegin(F);CHKERRQ(ierr); 2848e6d0c30SPeter Brune ierr = VecAssemblyEnd(F);CHKERRQ(ierr); 2858e6d0c30SPeter Brune 286*167fb786SPeter Brune /* split the operator into two */ 287bfde193fSPeter Brune ierr = PCGAMGClassicalGraphSplitting_Private(G,&lG,&gG);CHKERRQ(ierr); 288bfde193fSPeter Brune ierr = PCGAMGClassicalGraphSplitting_Private(A,&lA,&gA);CHKERRQ(ierr); 2898e6d0c30SPeter Brune 2908e6d0c30SPeter Brune /* scatter to the ghost vector */ 291bfde193fSPeter Brune ierr = PCGAMGClassicalCreateGhostVector_Private(G,&gF,NULL);CHKERRQ(ierr); 292bfde193fSPeter Brune ierr = PCGAMGClassicalGhost_Private(G,F,gF);CHKERRQ(ierr); 2938e6d0c30SPeter Brune 2948e6d0c30SPeter Brune if (gG) { 2958e6d0c30SPeter Brune ierr = VecGetSize(gF,&gn);CHKERRQ(ierr); 2968e6d0c30SPeter Brune ierr = PetscMalloc(sizeof(PetscInt)*gn,&gcid);CHKERRQ(ierr); 2978e6d0c30SPeter Brune for (i=0;i<gn;i++) { 2988e6d0c30SPeter Brune ierr = VecGetValues(gF,1,&i,&c_scalar);CHKERRQ(ierr); 299*167fb786SPeter Brune gcid[i] = *((PetscInt *)&c_scalar); 3008e6d0c30SPeter Brune } 3018e6d0c30SPeter Brune } 3028e6d0c30SPeter Brune 3038e6d0c30SPeter Brune ierr = VecDestroy(&F);CHKERRQ(ierr); 3048e6d0c30SPeter Brune ierr = VecDestroy(&gF);CHKERRQ(ierr); 3058e6d0c30SPeter Brune ierr = VecDestroy(&C);CHKERRQ(ierr); 3068e6d0c30SPeter Brune 3078e6d0c30SPeter Brune /* count the on and off processor sparsity patterns for the prolongator */ 3088e6d0c30SPeter Brune 3098e6d0c30SPeter Brune for (i=0;i<fn;i++) { 3108e6d0c30SPeter Brune /* on */ 3118e6d0c30SPeter Brune ncolstotal = ncols; 3128e6d0c30SPeter Brune lsparse[i] = 0; 313e5a0faa4SPeter Brune gsparse[i] = 0; 3148e6d0c30SPeter Brune if (lcid[i] >= 0) { 3158e6d0c30SPeter Brune lsparse[i] = 1; 3168e6d0c30SPeter Brune gsparse[i] = 0; 3178e6d0c30SPeter Brune } else { 318*167fb786SPeter Brune ierr = MatGetRow(lG,i,&ncols,&icols,&vcols);CHKERRQ(ierr); 3198e6d0c30SPeter Brune for (j = 0;j < ncols;j++) { 3208e6d0c30SPeter Brune col = icols[j]; 3218e6d0c30SPeter Brune if (lcid[col] >= 0 && vcols[j] != 0.) { 3228e6d0c30SPeter Brune lsparse[i] += 1; 3238e6d0c30SPeter Brune } 3248e6d0c30SPeter Brune } 3258e6d0c30SPeter Brune ierr = MatRestoreRow(lG,i,&ncols,&icols,&vcols);CHKERRQ(ierr); 3268e6d0c30SPeter Brune ncolstotal += ncols; 3278e6d0c30SPeter Brune /* off */ 3288e6d0c30SPeter Brune if (gG) { 3298e6d0c30SPeter Brune ierr = MatGetRow(gG,i,&ncols,&icols,&vcols);CHKERRQ(ierr); 3308e6d0c30SPeter Brune for (j = 0; j < ncols; j++) { 3318e6d0c30SPeter Brune col = icols[j]; 3328e6d0c30SPeter Brune if (gcid[col] >= 0 && vcols[j] != 0.) { 3338e6d0c30SPeter Brune gsparse[i] += 1; 3348e6d0c30SPeter Brune } 3358e6d0c30SPeter Brune } 336*167fb786SPeter Brune ierr = MatRestoreRow(gG,i,&ncols,&icols,&vcols);CHKERRQ(ierr); 3378e6d0c30SPeter Brune } 3388e6d0c30SPeter Brune if (ncolstotal > cmax) cmax = ncolstotal; 3398e6d0c30SPeter Brune } 3408e6d0c30SPeter Brune } 3418e6d0c30SPeter Brune 3428e6d0c30SPeter Brune ierr = PetscMalloc(sizeof(PetscInt)*cmax,&pcols);CHKERRQ(ierr); 3438e6d0c30SPeter Brune ierr = PetscMalloc(sizeof(PetscScalar)*cmax,&pvals);CHKERRQ(ierr); 3448e6d0c30SPeter Brune 3458e6d0c30SPeter Brune /* preallocate and create the prolongator */ 3468e6d0c30SPeter Brune ierr = MatCreate(comm,P); CHKERRQ(ierr); 3478e6d0c30SPeter Brune ierr = MatGetType(G,&mtype);CHKERRQ(ierr); 3488e6d0c30SPeter Brune ierr = MatSetType(*P,mtype);CHKERRQ(ierr); 3498e6d0c30SPeter Brune 3508e6d0c30SPeter Brune ierr = MatSetSizes(*P,fn,cn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 3518e6d0c30SPeter Brune ierr = MatMPIAIJSetPreallocation(*P,0,lsparse,0,gsparse);CHKERRQ(ierr); 3528e6d0c30SPeter Brune ierr = MatSeqAIJSetPreallocation(*P,0,lsparse);CHKERRQ(ierr); 3538e6d0c30SPeter Brune 3548e6d0c30SPeter Brune /* loop over local fine nodes -- get the diagonal, the sum of positive and negative strong and weak weights, and set up the row */ 3558e6d0c30SPeter Brune for (i = 0;i < fn;i++) { 3568e6d0c30SPeter Brune /* determine on or off */ 3578e6d0c30SPeter Brune row_f = i + fs; 3588e6d0c30SPeter Brune row_c = lcid[i]; 3598e6d0c30SPeter Brune if (row_c >= 0) { 3608e6d0c30SPeter Brune pij = 1.; 3618e6d0c30SPeter Brune ierr = MatSetValues(*P,1,&row_f,1,&row_c,&pij,INSERT_VALUES);CHKERRQ(ierr); 3628e6d0c30SPeter Brune } else { 36365b3d5b6SPeter Brune PetscInt nstrong=0,ntotal=0; 3648e6d0c30SPeter Brune g_pos = 0.; 3658e6d0c30SPeter Brune g_neg = 0.; 3668e6d0c30SPeter Brune a_pos = 0.; 3678e6d0c30SPeter Brune a_neg = 0.; 3688e6d0c30SPeter Brune diag = 0.; 3698e6d0c30SPeter Brune 3708e6d0c30SPeter Brune /* local strong connections */ 3718e6d0c30SPeter Brune ierr = MatGetRow(lG,i,&ncols,&rcol,&rval);CHKERRQ(ierr); 3728e6d0c30SPeter Brune for (k = 0; k < ncols; k++) { 373e5a0faa4SPeter Brune if (lcid[rcol[k]] >= 0) { 3748e6d0c30SPeter Brune if (PetscRealPart(rval[k]) > 0) { 3758e6d0c30SPeter Brune g_pos += rval[k]; 3768e6d0c30SPeter Brune } else { 3778e6d0c30SPeter Brune g_neg += rval[k]; 3788e6d0c30SPeter Brune } 379e5a0faa4SPeter Brune nstrong++; 380e5a0faa4SPeter Brune } 3818e6d0c30SPeter Brune } 3828e6d0c30SPeter Brune ierr = MatRestoreRow(lG,i,&ncols,&rcol,&rval);CHKERRQ(ierr); 3838e6d0c30SPeter Brune 3848e6d0c30SPeter Brune /* ghosted strong connections */ 3858e6d0c30SPeter Brune if (gG) { 3868e6d0c30SPeter Brune ierr = MatGetRow(gG,i,&ncols,&rcol,&rval);CHKERRQ(ierr); 3878e6d0c30SPeter Brune for (k = 0; k < ncols; k++) { 388e5a0faa4SPeter Brune if (gcid[rcol[k]] >= 0) { 3898e6d0c30SPeter Brune if (PetscRealPart(rval[k]) > 0.) { 390bfde193fSPeter Brune g_pos += rval[k]; 3918e6d0c30SPeter Brune } else { 392bfde193fSPeter Brune g_neg += rval[k]; 3938e6d0c30SPeter Brune } 394e5a0faa4SPeter Brune nstrong++; 395e5a0faa4SPeter Brune } 3968e6d0c30SPeter Brune } 3978e6d0c30SPeter Brune ierr = MatRestoreRow(gG,i,&ncols,&rcol,&rval);CHKERRQ(ierr); 3988e6d0c30SPeter Brune } 3998e6d0c30SPeter Brune 4008e6d0c30SPeter Brune /* local all connections */ 4018e6d0c30SPeter Brune ierr = MatGetRow(lA,i,&ncols,&rcol,&rval);CHKERRQ(ierr); 4028e6d0c30SPeter Brune for (k = 0; k < ncols; k++) { 403e5a0faa4SPeter Brune if (rcol[k] != i) { 4048e6d0c30SPeter Brune if (PetscRealPart(rval[k]) > 0) { 405bfde193fSPeter Brune a_pos += rval[k]; 4068e6d0c30SPeter Brune } else { 407bfde193fSPeter Brune a_neg += rval[k]; 4088e6d0c30SPeter Brune } 409e5a0faa4SPeter Brune ntotal++; 410e5a0faa4SPeter Brune } else diag = rval[k]; 4118e6d0c30SPeter Brune } 4128e6d0c30SPeter Brune ierr = MatRestoreRow(lA,i,&ncols,&rcol,&rval);CHKERRQ(ierr); 4138e6d0c30SPeter Brune 4148e6d0c30SPeter Brune /* ghosted all connections */ 4158e6d0c30SPeter Brune if (gA) { 4168e6d0c30SPeter Brune ierr = MatGetRow(gA,i,&ncols,&rcol,&rval);CHKERRQ(ierr); 4178e6d0c30SPeter Brune for (k = 0; k < ncols; k++) { 4188e6d0c30SPeter Brune if (PetscRealPart(rval[k]) > 0.) { 4198e6d0c30SPeter Brune a_pos += PetscRealPart(rval[k]); 4208e6d0c30SPeter Brune } else { 4218e6d0c30SPeter Brune a_neg += PetscRealPart(rval[k]); 4228e6d0c30SPeter Brune } 423e5a0faa4SPeter Brune ntotal++; 4248e6d0c30SPeter Brune } 4258e6d0c30SPeter Brune ierr = MatRestoreRow(gA,i,&ncols,&rcol,&rval);CHKERRQ(ierr); 4268e6d0c30SPeter Brune } 4278e6d0c30SPeter Brune 4288e6d0c30SPeter Brune if (g_neg == 0.) { 4298e6d0c30SPeter Brune alpha = 0.; 4308e6d0c30SPeter Brune } else { 4318e6d0c30SPeter Brune alpha = -a_neg/g_neg; 4328e6d0c30SPeter Brune } 4338e6d0c30SPeter Brune 4348e6d0c30SPeter Brune if (g_pos == 0.) { 4358e6d0c30SPeter Brune diag += a_pos; 4368e6d0c30SPeter Brune beta = 0.; 4378e6d0c30SPeter Brune } else { 4388e6d0c30SPeter Brune beta = -a_pos/g_pos; 4398e6d0c30SPeter Brune } 440e5a0faa4SPeter Brune if (diag == 0.) { 441e5a0faa4SPeter Brune invdiag = 0.; 442e5a0faa4SPeter Brune } else invdiag = 1. / diag; 4438e6d0c30SPeter Brune /* on */ 4448e6d0c30SPeter Brune ierr = MatGetRow(lG,i,&ncols,&icols,&vcols);CHKERRQ(ierr); 4458e6d0c30SPeter Brune idx = 0; 4468e6d0c30SPeter Brune for (j = 0;j < ncols;j++) { 4478e6d0c30SPeter Brune col = icols[j]; 4488e6d0c30SPeter Brune if (lcid[col] >= 0 && vcols[j] != 0.) { 4498e6d0c30SPeter Brune row_f = i + fs; 4508e6d0c30SPeter Brune row_c = lcid[col]; 4518e6d0c30SPeter Brune /* set the values for on-processor ones */ 4528e6d0c30SPeter Brune if (PetscRealPart(vcols[j]) < 0.) { 4538e6d0c30SPeter Brune pij = vcols[j]*alpha*invdiag; 4548e6d0c30SPeter Brune } else { 4558e6d0c30SPeter Brune pij = vcols[j]*beta*invdiag; 4568e6d0c30SPeter Brune } 4578e6d0c30SPeter Brune if (PetscAbsScalar(pij) != 0.) { 4588e6d0c30SPeter Brune pvals[idx] = pij; 4598e6d0c30SPeter Brune pcols[idx] = row_c; 4608e6d0c30SPeter Brune idx++; 4618e6d0c30SPeter Brune } 4628e6d0c30SPeter Brune } 4638e6d0c30SPeter Brune } 4648e6d0c30SPeter Brune ierr = MatRestoreRow(lG,i,&ncols,&icols,&vcols);CHKERRQ(ierr); 4658e6d0c30SPeter Brune /* off */ 4668e6d0c30SPeter Brune if (gG) { 4678e6d0c30SPeter Brune ierr = MatGetRow(gG,i,&ncols,&icols,&vcols);CHKERRQ(ierr); 4688e6d0c30SPeter Brune for (j = 0; j < ncols; j++) { 4698e6d0c30SPeter Brune col = icols[j]; 4708e6d0c30SPeter Brune if (gcid[col] >= 0 && vcols[j] != 0.) { 4718e6d0c30SPeter Brune row_f = i + fs; 4728e6d0c30SPeter Brune row_c = gcid[col]; 4738e6d0c30SPeter Brune /* set the values for on-processor ones */ 4748e6d0c30SPeter Brune if (PetscRealPart(vcols[j]) < 0.) { 4758e6d0c30SPeter Brune pij = vcols[j]*alpha*invdiag; 4768e6d0c30SPeter Brune } else { 4778e6d0c30SPeter Brune pij = vcols[j]*beta*invdiag; 4788e6d0c30SPeter Brune } 4798e6d0c30SPeter Brune if (PetscAbsScalar(pij) != 0.) { 4808e6d0c30SPeter Brune pvals[idx] = pij; 4818e6d0c30SPeter Brune pcols[idx] = row_c; 4828e6d0c30SPeter Brune idx++; 4838e6d0c30SPeter Brune } 4848e6d0c30SPeter Brune } 4858e6d0c30SPeter Brune } 4863c9ab2c3SPeter Brune ierr = MatRestoreRow(gG,i,&ncols,&icols,&vcols);CHKERRQ(ierr); 4873c9ab2c3SPeter Brune } 4888e6d0c30SPeter Brune ierr = MatSetValues(*P,1,&row_f,idx,pcols,pvals,INSERT_VALUES);CHKERRQ(ierr); 4898e6d0c30SPeter Brune } 4908e6d0c30SPeter Brune } 4913c9ab2c3SPeter Brune 4928e6d0c30SPeter Brune ierr = MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4938e6d0c30SPeter Brune ierr = MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4948e6d0c30SPeter Brune 4958e6d0c30SPeter Brune ierr = PetscFree(lsparse);CHKERRQ(ierr); 4968e6d0c30SPeter Brune ierr = PetscFree(gsparse);CHKERRQ(ierr); 4978e6d0c30SPeter Brune ierr = PetscFree(pcols);CHKERRQ(ierr); 4988e6d0c30SPeter Brune ierr = PetscFree(pvals);CHKERRQ(ierr); 4998e6d0c30SPeter Brune ierr = PetscFree(lcid);CHKERRQ(ierr); 5008e6d0c30SPeter Brune if (gG) {ierr = PetscFree(gcid);CHKERRQ(ierr);} 5018e6d0c30SPeter Brune 5028e6d0c30SPeter Brune PetscFunctionReturn(0); 5038e6d0c30SPeter Brune } 5048e6d0c30SPeter Brune 5058e6d0c30SPeter Brune #undef __FUNCT__ 5068e6d0c30SPeter Brune #define __FUNCT__ "PCGAMGDestroy_Classical" 5078e6d0c30SPeter Brune PetscErrorCode PCGAMGDestroy_Classical(PC pc) 5088e6d0c30SPeter Brune { 5098e6d0c30SPeter Brune PetscErrorCode ierr; 5108e6d0c30SPeter Brune PC_MG *mg = (PC_MG*)pc->data; 5118e6d0c30SPeter Brune PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 5128e6d0c30SPeter Brune 5138e6d0c30SPeter Brune PetscFunctionBegin; 5148e6d0c30SPeter Brune ierr = PetscFree(pc_gamg->subctx);CHKERRQ(ierr); 5158e6d0c30SPeter Brune PetscFunctionReturn(0); 5168e6d0c30SPeter Brune } 5178e6d0c30SPeter Brune 5188e6d0c30SPeter Brune #undef __FUNCT__ 5198e6d0c30SPeter Brune #define __FUNCT__ "PCGAMGSetFromOptions_Classical" 5208e6d0c30SPeter Brune PetscErrorCode PCGAMGSetFromOptions_Classical(PC pc) 5218e6d0c30SPeter Brune { 5228e6d0c30SPeter Brune PetscFunctionBegin; 5238e6d0c30SPeter Brune PetscFunctionReturn(0); 5248e6d0c30SPeter Brune } 5258e6d0c30SPeter Brune 5268e6d0c30SPeter Brune #undef __FUNCT__ 5278e6d0c30SPeter Brune #define __FUNCT__ "PCGAMGSetData_Classical" 5288e6d0c30SPeter Brune PetscErrorCode PCGAMGSetData_Classical(PC pc, Mat A) 5298e6d0c30SPeter Brune { 5308e6d0c30SPeter Brune PC_MG *mg = (PC_MG*)pc->data; 5318e6d0c30SPeter Brune PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 5328e6d0c30SPeter Brune 5338e6d0c30SPeter Brune PetscFunctionBegin; 5348e6d0c30SPeter Brune /* no data for classical AMG */ 5358e6d0c30SPeter Brune pc_gamg->data = NULL; 5368e6d0c30SPeter Brune pc_gamg->data_cell_cols = 1; 5378e6d0c30SPeter Brune pc_gamg->data_cell_rows = 1; 5388e6d0c30SPeter Brune pc_gamg->data_sz = 0; 5398e6d0c30SPeter Brune PetscFunctionReturn(0); 5408e6d0c30SPeter Brune } 5418e6d0c30SPeter Brune 5428e6d0c30SPeter Brune /* -------------------------------------------------------------------------- */ 5438e6d0c30SPeter Brune /* 5448e6d0c30SPeter Brune PCCreateGAMG_Classical 5458e6d0c30SPeter Brune 5468e6d0c30SPeter Brune */ 5478e6d0c30SPeter Brune #undef __FUNCT__ 5488e6d0c30SPeter Brune #define __FUNCT__ "PCCreateGAMG_Classical" 5498e6d0c30SPeter Brune PetscErrorCode PCCreateGAMG_Classical(PC pc) 5508e6d0c30SPeter Brune { 5518e6d0c30SPeter Brune PetscErrorCode ierr; 5528e6d0c30SPeter Brune PC_MG *mg = (PC_MG*)pc->data; 5538e6d0c30SPeter Brune PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 5548e6d0c30SPeter Brune PC_GAMG_Classical *pc_gamg_classical; 5558e6d0c30SPeter Brune 5568e6d0c30SPeter Brune PetscFunctionBegin; 5578e6d0c30SPeter Brune if (pc_gamg->subctx) { 5588e6d0c30SPeter Brune /* call base class */ 5598e6d0c30SPeter Brune ierr = PCDestroy_GAMG(pc);CHKERRQ(ierr); 5608e6d0c30SPeter Brune } 5618e6d0c30SPeter Brune 5628e6d0c30SPeter Brune /* create sub context for SA */ 5638e6d0c30SPeter Brune ierr = PetscNewLog(pc, PC_GAMG_Classical, &pc_gamg_classical);CHKERRQ(ierr); 5648e6d0c30SPeter Brune pc_gamg->subctx = pc_gamg_classical; 5658e6d0c30SPeter Brune pc->ops->setfromoptions = PCGAMGSetFromOptions_Classical; 5668e6d0c30SPeter Brune /* reset does not do anything; setup not virtual */ 5678e6d0c30SPeter Brune 5688e6d0c30SPeter Brune /* set internal function pointers */ 5698e6d0c30SPeter Brune pc_gamg->ops->destroy = PCGAMGDestroy_Classical; 5708e6d0c30SPeter Brune pc_gamg->ops->graph = PCGAMGGraph_Classical; 5718e6d0c30SPeter Brune pc_gamg->ops->coarsen = PCGAMGCoarsen_Classical; 5728e6d0c30SPeter Brune pc_gamg->ops->prolongator = PCGAMGProlongator_Classical; 5738e6d0c30SPeter Brune pc_gamg->ops->optprol = NULL; 5748e6d0c30SPeter Brune 5758e6d0c30SPeter Brune pc_gamg->ops->createdefaultdata = PCGAMGSetData_Classical; 5768e6d0c30SPeter Brune PetscFunctionReturn(0); 5778e6d0c30SPeter Brune } 578