xref: /petsc/src/ksp/pc/impls/gamg/classical.c (revision 167fb786fdb731ab03cfde3016f1ed5bc75cc9ed)
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