1 /* 2 GAMG geometric-algebric multigrid PC - Mark Adams 2011 3 */ 4 5 #include <../src/ksp/pc/impls/gamg/gamg.h> /*I "petscpc.h" I*/ 6 #include <petscblaslapack.h> 7 #include <petscdm.h> 8 #include <petsc/private/kspimpl.h> 9 10 typedef struct { 11 PetscInt nsmooths; 12 PetscBool sym_graph; 13 PetscInt square_graph; 14 } PC_GAMG_AGG; 15 16 /*@ 17 PCGAMGSetNSmooths - Set number of smoothing steps (1 is typical) 18 19 Logically Collective on PC 20 21 Input Parameters: 22 . pc - the preconditioner context 23 24 Options Database Key: 25 . -pc_gamg_agg_nsmooths <nsmooth, default=1> - number of smoothing steps to use with smooth aggregation 26 27 Level: intermediate 28 29 .seealso: () 30 @*/ 31 PetscErrorCode PCGAMGSetNSmooths(PC pc, PetscInt n) 32 { 33 PetscErrorCode ierr; 34 35 PetscFunctionBegin; 36 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 37 PetscValidLogicalCollectiveInt(pc,n,2); 38 ierr = PetscTryMethod(pc,"PCGAMGSetNSmooths_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr); 39 PetscFunctionReturn(0); 40 } 41 42 static PetscErrorCode PCGAMGSetNSmooths_AGG(PC pc, PetscInt n) 43 { 44 PC_MG *mg = (PC_MG*)pc->data; 45 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 46 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 47 48 PetscFunctionBegin; 49 pc_gamg_agg->nsmooths = n; 50 PetscFunctionReturn(0); 51 } 52 53 /*@ 54 PCGAMGSetSymGraph - Symmetrize the graph before computing the aggregation. Some algorithms require the graph be symmetric 55 56 Logically Collective on PC 57 58 Input Parameters: 59 + pc - the preconditioner context 60 - n - PETSC_TRUE or PETSC_FALSE 61 62 Options Database Key: 63 . -pc_gamg_sym_graph <true,default=false> - symmetrize the graph before computing the aggregation 64 65 Level: intermediate 66 67 .seealso: PCGAMGSetSquareGraph() 68 @*/ 69 PetscErrorCode PCGAMGSetSymGraph(PC pc, PetscBool n) 70 { 71 PetscErrorCode ierr; 72 73 PetscFunctionBegin; 74 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 75 PetscValidLogicalCollectiveBool(pc,n,2); 76 ierr = PetscTryMethod(pc,"PCGAMGSetSymGraph_C",(PC,PetscBool),(pc,n));CHKERRQ(ierr); 77 PetscFunctionReturn(0); 78 } 79 80 static PetscErrorCode PCGAMGSetSymGraph_AGG(PC pc, PetscBool n) 81 { 82 PC_MG *mg = (PC_MG*)pc->data; 83 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 84 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 85 86 PetscFunctionBegin; 87 pc_gamg_agg->sym_graph = n; 88 PetscFunctionReturn(0); 89 } 90 91 /*@ 92 PCGAMGSetSquareGraph - Square the graph, ie. compute A'*A before aggregating it 93 94 Logically Collective on PC 95 96 Input Parameters: 97 + pc - the preconditioner context 98 - n - 0, 1 or more 99 100 Options Database Key: 101 . -pc_gamg_square_graph <n,default = 1> - number of levels to square the graph on before aggregating it 102 103 Notes: 104 Squaring the graph increases the rate of coarsening (aggressive coarsening) and thereby reduces the complexity of the coarse grids, and generally results in slower solver converge rates. Reducing coarse grid complexity reduced the complexity of Galerkin coarse grid construction considerably. 105 106 Level: intermediate 107 108 .seealso: PCGAMGSetSymGraph(), PCGAMGSetThreshold() 109 @*/ 110 PetscErrorCode PCGAMGSetSquareGraph(PC pc, PetscInt n) 111 { 112 PetscErrorCode ierr; 113 114 PetscFunctionBegin; 115 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 116 PetscValidLogicalCollectiveInt(pc,n,2); 117 ierr = PetscTryMethod(pc,"PCGAMGSetSquareGraph_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr); 118 PetscFunctionReturn(0); 119 } 120 121 static PetscErrorCode PCGAMGSetSquareGraph_AGG(PC pc, PetscInt n) 122 { 123 PC_MG *mg = (PC_MG*)pc->data; 124 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 125 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 126 127 PetscFunctionBegin; 128 pc_gamg_agg->square_graph = n; 129 PetscFunctionReturn(0); 130 } 131 132 static PetscErrorCode PCSetFromOptions_GAMG_AGG(PetscOptionItems *PetscOptionsObject,PC pc) 133 { 134 PetscErrorCode ierr; 135 PC_MG *mg = (PC_MG*)pc->data; 136 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 137 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 138 139 PetscFunctionBegin; 140 ierr = PetscOptionsHead(PetscOptionsObject,"GAMG-AGG options");CHKERRQ(ierr); 141 { 142 ierr = PetscOptionsInt("-pc_gamg_agg_nsmooths","smoothing steps for smoothed aggregation, usually 1","PCGAMGSetNSmooths",pc_gamg_agg->nsmooths,&pc_gamg_agg->nsmooths,NULL);CHKERRQ(ierr); 143 ierr = PetscOptionsBool("-pc_gamg_sym_graph","Set for asymmetric matrices","PCGAMGSetSymGraph",pc_gamg_agg->sym_graph,&pc_gamg_agg->sym_graph,NULL);CHKERRQ(ierr); 144 ierr = PetscOptionsInt("-pc_gamg_square_graph","Number of levels to square graph for faster coarsening and lower coarse grid complexity","PCGAMGSetSquareGraph",pc_gamg_agg->square_graph,&pc_gamg_agg->square_graph,NULL);CHKERRQ(ierr); 145 } 146 ierr = PetscOptionsTail();CHKERRQ(ierr); 147 PetscFunctionReturn(0); 148 } 149 150 /* -------------------------------------------------------------------------- */ 151 static PetscErrorCode PCDestroy_GAMG_AGG(PC pc) 152 { 153 PetscErrorCode ierr; 154 PC_MG *mg = (PC_MG*)pc->data; 155 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 156 157 PetscFunctionBegin; 158 ierr = PetscFree(pc_gamg->subctx);CHKERRQ(ierr); 159 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",NULL);CHKERRQ(ierr); 160 PetscFunctionReturn(0); 161 } 162 163 /* -------------------------------------------------------------------------- */ 164 /* 165 PCSetCoordinates_AGG 166 - collective 167 168 Input Parameter: 169 . pc - the preconditioner context 170 . ndm - dimesion of data (used for dof/vertex for Stokes) 171 . a_nloc - number of vertices local 172 . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...} 173 */ 174 175 static PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords) 176 { 177 PC_MG *mg = (PC_MG*)pc->data; 178 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 179 PetscErrorCode ierr; 180 PetscInt arrsz,kk,ii,jj,nloc,ndatarows,ndf; 181 Mat mat = pc->pmat; 182 183 PetscFunctionBegin; 184 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 185 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 186 nloc = a_nloc; 187 188 /* SA: null space vectors */ 189 ierr = MatGetBlockSize(mat, &ndf);CHKERRQ(ierr); /* this does not work for Stokes */ 190 if (coords && ndf==1) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */ 191 else if (coords) { 192 PetscCheckFalse(ndm > ndf,PETSC_COMM_SELF,PETSC_ERR_PLIB,"degrees of motion %D > block size %D",ndm,ndf); 193 pc_gamg->data_cell_cols = (ndm==2 ? 3 : 6); /* displacement elasticity */ 194 if (ndm != ndf) { 195 PetscCheckFalse(pc_gamg->data_cell_cols != ndf,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Don't know how to create null space for ndm=%D, ndf=%D. Use MatSetNearNullSpace().",ndm,ndf); 196 } 197 } else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */ 198 pc_gamg->data_cell_rows = ndatarows = ndf; 199 PetscCheckFalse(pc_gamg->data_cell_cols <= 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"pc_gamg->data_cell_cols %D <= 0",pc_gamg->data_cell_cols); 200 arrsz = nloc*pc_gamg->data_cell_rows*pc_gamg->data_cell_cols; 201 202 if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) { 203 ierr = PetscFree(pc_gamg->data);CHKERRQ(ierr); 204 ierr = PetscMalloc1(arrsz+1, &pc_gamg->data);CHKERRQ(ierr); 205 } 206 /* copy data in - column oriented */ 207 for (kk=0; kk<nloc; kk++) { 208 const PetscInt M = nloc*pc_gamg->data_cell_rows; /* stride into data */ 209 PetscReal *data = &pc_gamg->data[kk*ndatarows]; /* start of cell */ 210 if (pc_gamg->data_cell_cols==1) *data = 1.0; 211 else { 212 /* translational modes */ 213 for (ii=0;ii<ndatarows;ii++) { 214 for (jj=0;jj<ndatarows;jj++) { 215 if (ii==jj)data[ii*M + jj] = 1.0; 216 else data[ii*M + jj] = 0.0; 217 } 218 } 219 220 /* rotational modes */ 221 if (coords) { 222 if (ndm == 2) { 223 data += 2*M; 224 data[0] = -coords[2*kk+1]; 225 data[1] = coords[2*kk]; 226 } else { 227 data += 3*M; 228 data[0] = 0.0; data[M+0] = coords[3*kk+2]; data[2*M+0] = -coords[3*kk+1]; 229 data[1] = -coords[3*kk+2]; data[M+1] = 0.0; data[2*M+1] = coords[3*kk]; 230 data[2] = coords[3*kk+1]; data[M+2] = -coords[3*kk]; data[2*M+2] = 0.0; 231 } 232 } 233 } 234 } 235 pc_gamg->data_sz = arrsz; 236 PetscFunctionReturn(0); 237 } 238 239 typedef PetscInt NState; 240 static const NState NOT_DONE=-2; 241 static const NState DELETED =-1; 242 static const NState REMOVED =-3; 243 #define IS_SELECTED(s) (s!=DELETED && s!=NOT_DONE && s!=REMOVED) 244 245 /* -------------------------------------------------------------------------- */ 246 /* 247 smoothAggs - greedy grab of with G1 (unsquared graph) -- AIJ specific 248 - AGG-MG specific: clears singletons out of 'selected_2' 249 250 Input Parameter: 251 . Gmat_2 - global matrix of graph (data not defined) base (squared) graph 252 . Gmat_1 - base graph to grab with base graph 253 Input/Output Parameter: 254 . aggs_2 - linked list of aggs with gids) 255 */ 256 static PetscErrorCode smoothAggs(PC pc,Mat Gmat_2, Mat Gmat_1,PetscCoarsenData *aggs_2) 257 { 258 PetscErrorCode ierr; 259 PetscBool isMPI; 260 Mat_SeqAIJ *matA_1, *matB_1=NULL; 261 MPI_Comm comm; 262 PetscInt lid,*ii,*idx,ix,Iend,my0,kk,n,j; 263 Mat_MPIAIJ *mpimat_2 = NULL, *mpimat_1=NULL; 264 const PetscInt nloc = Gmat_2->rmap->n; 265 PetscScalar *cpcol_1_state,*cpcol_2_state,*cpcol_2_par_orig,*lid_parent_gid; 266 PetscInt *lid_cprowID_1; 267 NState *lid_state; 268 Vec ghost_par_orig2; 269 270 PetscFunctionBegin; 271 ierr = PetscObjectGetComm((PetscObject)Gmat_2,&comm);CHKERRQ(ierr); 272 ierr = MatGetOwnershipRange(Gmat_1,&my0,&Iend);CHKERRQ(ierr); 273 274 /* get submatrices */ 275 ierr = PetscStrbeginswith(((PetscObject)Gmat_1)->type_name,MATMPIAIJ,&isMPI);CHKERRQ(ierr); 276 if (isMPI) { 277 /* grab matrix objects */ 278 mpimat_2 = (Mat_MPIAIJ*)Gmat_2->data; 279 mpimat_1 = (Mat_MPIAIJ*)Gmat_1->data; 280 matA_1 = (Mat_SeqAIJ*)mpimat_1->A->data; 281 matB_1 = (Mat_SeqAIJ*)mpimat_1->B->data; 282 283 /* force compressed row storage for B matrix in AuxMat */ 284 ierr = MatCheckCompressedRow(mpimat_1->B,matB_1->nonzerorowcnt,&matB_1->compressedrow,matB_1->i,Gmat_1->rmap->n,-1.0);CHKERRQ(ierr); 285 286 ierr = PetscMalloc1(nloc, &lid_cprowID_1);CHKERRQ(ierr); 287 for (lid = 0; lid < nloc; lid++) lid_cprowID_1[lid] = -1; 288 for (ix=0; ix<matB_1->compressedrow.nrows; ix++) { 289 PetscInt lid = matB_1->compressedrow.rindex[ix]; 290 lid_cprowID_1[lid] = ix; 291 } 292 } else { 293 PetscBool isAIJ; 294 ierr = PetscStrbeginswith(((PetscObject)Gmat_1)->type_name,MATSEQAIJ,&isAIJ);CHKERRQ(ierr); 295 PetscCheckFalse(!isAIJ,PETSC_COMM_SELF,PETSC_ERR_USER,"Require AIJ matrix."); 296 matA_1 = (Mat_SeqAIJ*)Gmat_1->data; 297 lid_cprowID_1 = NULL; 298 } 299 if (nloc>0) { 300 PetscCheckFalse(matB_1 && !matB_1->compressedrow.use,PETSC_COMM_SELF,PETSC_ERR_PLIB,"matB_1 && !matB_1->compressedrow.use: PETSc bug???"); 301 } 302 /* get state of locals and selected gid for deleted */ 303 ierr = PetscMalloc2(nloc, &lid_state,nloc, &lid_parent_gid);CHKERRQ(ierr); 304 for (lid = 0; lid < nloc; lid++) { 305 lid_parent_gid[lid] = -1.0; 306 lid_state[lid] = DELETED; 307 } 308 309 /* set lid_state */ 310 for (lid = 0; lid < nloc; lid++) { 311 PetscCDIntNd *pos; 312 ierr = PetscCDGetHeadPos(aggs_2,lid,&pos);CHKERRQ(ierr); 313 if (pos) { 314 PetscInt gid1; 315 316 ierr = PetscCDIntNdGetID(pos, &gid1);CHKERRQ(ierr); 317 PetscCheckFalse(gid1 != lid+my0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"gid1 %D != lid %D + my0 %D",gid1,lid,my0); 318 lid_state[lid] = gid1; 319 } 320 } 321 322 /* map local to selected local, DELETED means a ghost owns it */ 323 for (lid=kk=0; lid<nloc; lid++) { 324 NState state = lid_state[lid]; 325 if (IS_SELECTED(state)) { 326 PetscCDIntNd *pos; 327 ierr = PetscCDGetHeadPos(aggs_2,lid,&pos);CHKERRQ(ierr); 328 while (pos) { 329 PetscInt gid1; 330 ierr = PetscCDIntNdGetID(pos, &gid1);CHKERRQ(ierr); 331 ierr = PetscCDGetNextPos(aggs_2,lid,&pos);CHKERRQ(ierr); 332 if (gid1 >= my0 && gid1 < Iend) lid_parent_gid[gid1-my0] = (PetscScalar)(lid + my0); 333 } 334 } 335 } 336 /* get 'cpcol_1/2_state' & cpcol_2_par_orig - uses mpimat_1/2->lvec for temp space */ 337 if (isMPI) { 338 Vec tempVec; 339 /* get 'cpcol_1_state' */ 340 ierr = MatCreateVecs(Gmat_1, &tempVec, NULL);CHKERRQ(ierr); 341 for (kk=0,j=my0; kk<nloc; kk++,j++) { 342 PetscScalar v = (PetscScalar)lid_state[kk]; 343 ierr = VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES);CHKERRQ(ierr); 344 } 345 ierr = VecAssemblyBegin(tempVec);CHKERRQ(ierr); 346 ierr = VecAssemblyEnd(tempVec);CHKERRQ(ierr); 347 ierr = VecScatterBegin(mpimat_1->Mvctx,tempVec, mpimat_1->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 348 ierr = VecScatterEnd(mpimat_1->Mvctx,tempVec, mpimat_1->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 349 ierr = VecGetArray(mpimat_1->lvec, &cpcol_1_state);CHKERRQ(ierr); 350 /* get 'cpcol_2_state' */ 351 ierr = VecScatterBegin(mpimat_2->Mvctx,tempVec, mpimat_2->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 352 ierr = VecScatterEnd(mpimat_2->Mvctx,tempVec, mpimat_2->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 353 ierr = VecGetArray(mpimat_2->lvec, &cpcol_2_state);CHKERRQ(ierr); 354 /* get 'cpcol_2_par_orig' */ 355 for (kk=0,j=my0; kk<nloc; kk++,j++) { 356 PetscScalar v = (PetscScalar)lid_parent_gid[kk]; 357 ierr = VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES);CHKERRQ(ierr); 358 } 359 ierr = VecAssemblyBegin(tempVec);CHKERRQ(ierr); 360 ierr = VecAssemblyEnd(tempVec);CHKERRQ(ierr); 361 ierr = VecDuplicate(mpimat_2->lvec, &ghost_par_orig2);CHKERRQ(ierr); 362 ierr = VecScatterBegin(mpimat_2->Mvctx,tempVec, ghost_par_orig2,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 363 ierr = VecScatterEnd(mpimat_2->Mvctx,tempVec, ghost_par_orig2,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 364 ierr = VecGetArray(ghost_par_orig2, &cpcol_2_par_orig);CHKERRQ(ierr); 365 366 ierr = VecDestroy(&tempVec);CHKERRQ(ierr); 367 } /* ismpi */ 368 for (lid=0; lid<nloc; lid++) { 369 NState state = lid_state[lid]; 370 if (IS_SELECTED(state)) { 371 /* steal locals */ 372 ii = matA_1->i; n = ii[lid+1] - ii[lid]; 373 idx = matA_1->j + ii[lid]; 374 for (j=0; j<n; j++) { 375 PetscInt lidj = idx[j], sgid; 376 NState statej = lid_state[lidj]; 377 if (statej==DELETED && (sgid=(PetscInt)PetscRealPart(lid_parent_gid[lidj])) != lid+my0) { /* steal local */ 378 lid_parent_gid[lidj] = (PetscScalar)(lid+my0); /* send this if sgid is not local */ 379 if (sgid >= my0 && sgid < Iend) { /* I'm stealing this local from a local sgid */ 380 PetscInt hav=0,slid=sgid-my0,gidj=lidj+my0; 381 PetscCDIntNd *pos,*last=NULL; 382 /* looking for local from local so id_llist_2 works */ 383 ierr = PetscCDGetHeadPos(aggs_2,slid,&pos);CHKERRQ(ierr); 384 while (pos) { 385 PetscInt gid; 386 ierr = PetscCDIntNdGetID(pos, &gid);CHKERRQ(ierr); 387 if (gid == gidj) { 388 PetscCheckFalse(!last,PETSC_COMM_SELF,PETSC_ERR_PLIB,"last cannot be null"); 389 ierr = PetscCDRemoveNextNode(aggs_2, slid, last);CHKERRQ(ierr); 390 ierr = PetscCDAppendNode(aggs_2, lid, pos);CHKERRQ(ierr); 391 hav = 1; 392 break; 393 } else last = pos; 394 ierr = PetscCDGetNextPos(aggs_2,slid,&pos);CHKERRQ(ierr); 395 } 396 if (hav != 1) { 397 PetscCheckFalse(!hav,PETSC_COMM_SELF,PETSC_ERR_PLIB,"failed to find adj in 'selected' lists - structurally unsymmetric matrix"); 398 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"found node %D times???",hav); 399 } 400 } else { /* I'm stealing this local, owned by a ghost */ 401 PetscCheckFalse(sgid != -1,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mat has an un-symmetric graph. Use '-%spc_gamg_sym_graph true' to symmetrize the graph or '-%spc_gamg_threshold -1' if the matrix is structurally symmetric.",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : ""); 402 ierr = PetscCDAppendID(aggs_2, lid, lidj+my0);CHKERRQ(ierr); 403 } 404 } 405 } /* local neighbors */ 406 } else if (state == DELETED && lid_cprowID_1) { 407 PetscInt sgidold = (PetscInt)PetscRealPart(lid_parent_gid[lid]); 408 /* see if I have a selected ghost neighbor that will steal me */ 409 if ((ix=lid_cprowID_1[lid]) != -1) { 410 ii = matB_1->compressedrow.i; n = ii[ix+1] - ii[ix]; 411 idx = matB_1->j + ii[ix]; 412 for (j=0; j<n; j++) { 413 PetscInt cpid = idx[j]; 414 NState statej = (NState)PetscRealPart(cpcol_1_state[cpid]); 415 if (IS_SELECTED(statej) && sgidold != (PetscInt)statej) { /* ghost will steal this, remove from my list */ 416 lid_parent_gid[lid] = (PetscScalar)statej; /* send who selected */ 417 if (sgidold>=my0 && sgidold<Iend) { /* this was mine */ 418 PetscInt hav=0,oldslidj=sgidold-my0; 419 PetscCDIntNd *pos,*last=NULL; 420 /* remove from 'oldslidj' list */ 421 ierr = PetscCDGetHeadPos(aggs_2,oldslidj,&pos);CHKERRQ(ierr); 422 while (pos) { 423 PetscInt gid; 424 ierr = PetscCDIntNdGetID(pos, &gid);CHKERRQ(ierr); 425 if (lid+my0 == gid) { 426 /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove lid from oldslidj list *\/ */ 427 PetscCheckFalse(!last,PETSC_COMM_SELF,PETSC_ERR_PLIB,"last cannot be null"); 428 ierr = PetscCDRemoveNextNode(aggs_2, oldslidj, last);CHKERRQ(ierr); 429 /* ghost (PetscScalar)statej will add this later */ 430 hav = 1; 431 break; 432 } else last = pos; 433 ierr = PetscCDGetNextPos(aggs_2,oldslidj,&pos);CHKERRQ(ierr); 434 } 435 if (hav != 1) { 436 PetscCheckFalse(!hav,PETSC_COMM_SELF,PETSC_ERR_PLIB,"failed to find adj in 'selected' lists - structurally unsymmetric matrix"); 437 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"found node %D times???",hav); 438 } 439 } else { 440 /* TODO: ghosts remove this later */ 441 } 442 } 443 } 444 } 445 } /* selected/deleted */ 446 } /* node loop */ 447 448 if (isMPI) { 449 PetscScalar *cpcol_2_parent,*cpcol_2_gid; 450 Vec tempVec,ghostgids2,ghostparents2; 451 PetscInt cpid,nghost_2; 452 PCGAMGHashTable gid_cpid; 453 454 ierr = VecGetSize(mpimat_2->lvec, &nghost_2);CHKERRQ(ierr); 455 ierr = MatCreateVecs(Gmat_2, &tempVec, NULL);CHKERRQ(ierr); 456 457 /* get 'cpcol_2_parent' */ 458 for (kk=0,j=my0; kk<nloc; kk++,j++) { 459 ierr = VecSetValues(tempVec, 1, &j, &lid_parent_gid[kk], INSERT_VALUES);CHKERRQ(ierr); 460 } 461 ierr = VecAssemblyBegin(tempVec);CHKERRQ(ierr); 462 ierr = VecAssemblyEnd(tempVec);CHKERRQ(ierr); 463 ierr = VecDuplicate(mpimat_2->lvec, &ghostparents2);CHKERRQ(ierr); 464 ierr = VecScatterBegin(mpimat_2->Mvctx,tempVec, ghostparents2,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 465 ierr = VecScatterEnd(mpimat_2->Mvctx,tempVec, ghostparents2,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 466 ierr = VecGetArray(ghostparents2, &cpcol_2_parent);CHKERRQ(ierr); 467 468 /* get 'cpcol_2_gid' */ 469 for (kk=0,j=my0; kk<nloc; kk++,j++) { 470 PetscScalar v = (PetscScalar)j; 471 ierr = VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES);CHKERRQ(ierr); 472 } 473 ierr = VecAssemblyBegin(tempVec);CHKERRQ(ierr); 474 ierr = VecAssemblyEnd(tempVec);CHKERRQ(ierr); 475 ierr = VecDuplicate(mpimat_2->lvec, &ghostgids2);CHKERRQ(ierr); 476 ierr = VecScatterBegin(mpimat_2->Mvctx,tempVec, ghostgids2,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 477 ierr = VecScatterEnd(mpimat_2->Mvctx,tempVec, ghostgids2,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 478 ierr = VecGetArray(ghostgids2, &cpcol_2_gid);CHKERRQ(ierr); 479 ierr = VecDestroy(&tempVec);CHKERRQ(ierr); 480 481 /* look for deleted ghosts and add to table */ 482 ierr = PCGAMGHashTableCreate(2*nghost_2+1, &gid_cpid);CHKERRQ(ierr); 483 for (cpid = 0; cpid < nghost_2; cpid++) { 484 NState state = (NState)PetscRealPart(cpcol_2_state[cpid]); 485 if (state==DELETED) { 486 PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]); 487 PetscInt sgid_old = (PetscInt)PetscRealPart(cpcol_2_par_orig[cpid]); 488 if (sgid_old == -1 && sgid_new != -1) { 489 PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]); 490 ierr = PCGAMGHashTableAdd(&gid_cpid, gid, cpid);CHKERRQ(ierr); 491 } 492 } 493 } 494 495 /* look for deleted ghosts and see if they moved - remove it */ 496 for (lid=0; lid<nloc; lid++) { 497 NState state = lid_state[lid]; 498 if (IS_SELECTED(state)) { 499 PetscCDIntNd *pos,*last=NULL; 500 /* look for deleted ghosts and see if they moved */ 501 ierr = PetscCDGetHeadPos(aggs_2,lid,&pos);CHKERRQ(ierr); 502 while (pos) { 503 PetscInt gid; 504 ierr = PetscCDIntNdGetID(pos, &gid);CHKERRQ(ierr); 505 506 if (gid < my0 || gid >= Iend) { 507 ierr = PCGAMGHashTableFind(&gid_cpid, gid, &cpid);CHKERRQ(ierr); 508 if (cpid != -1) { 509 /* a moved ghost - */ 510 /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove 'flid' from list *\/ */ 511 ierr = PetscCDRemoveNextNode(aggs_2, lid, last);CHKERRQ(ierr); 512 } else last = pos; 513 } else last = pos; 514 515 ierr = PetscCDGetNextPos(aggs_2,lid,&pos);CHKERRQ(ierr); 516 } /* loop over list of deleted */ 517 } /* selected */ 518 } 519 ierr = PCGAMGHashTableDestroy(&gid_cpid);CHKERRQ(ierr); 520 521 /* look at ghosts, see if they changed - and it */ 522 for (cpid = 0; cpid < nghost_2; cpid++) { 523 PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]); 524 if (sgid_new >= my0 && sgid_new < Iend) { /* this is mine */ 525 PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]); 526 PetscInt slid_new=sgid_new-my0,hav=0; 527 PetscCDIntNd *pos; 528 529 /* search for this gid to see if I have it */ 530 ierr = PetscCDGetHeadPos(aggs_2,slid_new,&pos);CHKERRQ(ierr); 531 while (pos) { 532 PetscInt gidj; 533 ierr = PetscCDIntNdGetID(pos, &gidj);CHKERRQ(ierr); 534 ierr = PetscCDGetNextPos(aggs_2,slid_new,&pos);CHKERRQ(ierr); 535 536 if (gidj == gid) { hav = 1; break; } 537 } 538 if (hav != 1) { 539 /* insert 'flidj' into head of llist */ 540 ierr = PetscCDAppendID(aggs_2, slid_new, gid);CHKERRQ(ierr); 541 } 542 } 543 } 544 545 ierr = VecRestoreArray(mpimat_1->lvec, &cpcol_1_state);CHKERRQ(ierr); 546 ierr = VecRestoreArray(mpimat_2->lvec, &cpcol_2_state);CHKERRQ(ierr); 547 ierr = VecRestoreArray(ghostparents2, &cpcol_2_parent);CHKERRQ(ierr); 548 ierr = VecRestoreArray(ghostgids2, &cpcol_2_gid);CHKERRQ(ierr); 549 ierr = PetscFree(lid_cprowID_1);CHKERRQ(ierr); 550 ierr = VecDestroy(&ghostgids2);CHKERRQ(ierr); 551 ierr = VecDestroy(&ghostparents2);CHKERRQ(ierr); 552 ierr = VecDestroy(&ghost_par_orig2);CHKERRQ(ierr); 553 } 554 555 ierr = PetscFree2(lid_state,lid_parent_gid);CHKERRQ(ierr); 556 PetscFunctionReturn(0); 557 } 558 559 /* -------------------------------------------------------------------------- */ 560 /* 561 PCSetData_AGG - called if data is not set with PCSetCoordinates. 562 Looks in Mat for near null space. 563 Does not work for Stokes 564 565 Input Parameter: 566 . pc - 567 . a_A - matrix to get (near) null space out of. 568 */ 569 static PetscErrorCode PCSetData_AGG(PC pc, Mat a_A) 570 { 571 PetscErrorCode ierr; 572 PC_MG *mg = (PC_MG*)pc->data; 573 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 574 MatNullSpace mnull; 575 576 PetscFunctionBegin; 577 ierr = MatGetNearNullSpace(a_A, &mnull);CHKERRQ(ierr); 578 if (!mnull) { 579 DM dm; 580 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 581 if (!dm) { 582 ierr = MatGetDM(a_A, &dm);CHKERRQ(ierr); 583 } 584 if (dm) { 585 PetscObject deformation; 586 PetscInt Nf; 587 588 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 589 if (Nf) { 590 ierr = DMGetField(dm, 0, NULL, &deformation);CHKERRQ(ierr); 591 ierr = PetscObjectQuery((PetscObject)deformation,"nearnullspace",(PetscObject*)&mnull);CHKERRQ(ierr); 592 if (!mnull) { 593 ierr = PetscObjectQuery((PetscObject)deformation,"nullspace",(PetscObject*)&mnull);CHKERRQ(ierr); 594 } 595 } 596 } 597 } 598 599 if (!mnull) { 600 PetscInt bs,NN,MM; 601 ierr = MatGetBlockSize(a_A, &bs);CHKERRQ(ierr); 602 ierr = MatGetLocalSize(a_A, &MM, &NN);CHKERRQ(ierr); 603 PetscCheckFalse(MM % bs,PETSC_COMM_SELF,PETSC_ERR_PLIB,"MM %D must be divisible by bs %D",MM,bs); 604 ierr = PCSetCoordinates_AGG(pc, bs, MM/bs, NULL);CHKERRQ(ierr); 605 } else { 606 PetscReal *nullvec; 607 PetscBool has_const; 608 PetscInt i,j,mlocal,nvec,bs; 609 const Vec *vecs; 610 const PetscScalar *v; 611 612 ierr = MatGetLocalSize(a_A,&mlocal,NULL);CHKERRQ(ierr); 613 ierr = MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);CHKERRQ(ierr); 614 pc_gamg->data_sz = (nvec+!!has_const)*mlocal; 615 ierr = PetscMalloc1((nvec+!!has_const)*mlocal,&nullvec);CHKERRQ(ierr); 616 if (has_const) for (i=0; i<mlocal; i++) nullvec[i] = 1.0; 617 for (i=0; i<nvec; i++) { 618 ierr = VecGetArrayRead(vecs[i],&v);CHKERRQ(ierr); 619 for (j=0; j<mlocal; j++) nullvec[(i+!!has_const)*mlocal + j] = PetscRealPart(v[j]); 620 ierr = VecRestoreArrayRead(vecs[i],&v);CHKERRQ(ierr); 621 } 622 pc_gamg->data = nullvec; 623 pc_gamg->data_cell_cols = (nvec+!!has_const); 624 ierr = MatGetBlockSize(a_A, &bs);CHKERRQ(ierr); 625 pc_gamg->data_cell_rows = bs; 626 } 627 PetscFunctionReturn(0); 628 } 629 630 /* -------------------------------------------------------------------------- */ 631 /* 632 formProl0 633 634 Input Parameter: 635 . agg_llists - list of arrays with aggregates -- list from selected vertices of aggregate unselected vertices 636 . bs - row block size 637 . nSAvec - column bs of new P 638 . my0crs - global index of start of locals 639 . data_stride - bs*(nloc nodes + ghost nodes) [data_stride][nSAvec] 640 . data_in[data_stride*nSAvec] - local data on fine grid 641 . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist' 642 Output Parameter: 643 . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data 644 . a_Prol - prolongation operator 645 */ 646 static PetscErrorCode formProl0(PetscCoarsenData *agg_llists,PetscInt bs,PetscInt nSAvec,PetscInt my0crs,PetscInt data_stride,PetscReal data_in[],const PetscInt flid_fgid[],PetscReal **a_data_out,Mat a_Prol) 647 { 648 PetscErrorCode ierr; 649 PetscInt Istart,my0,Iend,nloc,clid,flid = 0,aggID,kk,jj,ii,mm,nSelected,minsz,nghosts,out_data_stride; 650 MPI_Comm comm; 651 PetscReal *out_data; 652 PetscCDIntNd *pos; 653 PCGAMGHashTable fgid_flid; 654 655 PetscFunctionBegin; 656 ierr = PetscObjectGetComm((PetscObject)a_Prol,&comm);CHKERRQ(ierr); 657 ierr = MatGetOwnershipRange(a_Prol, &Istart, &Iend);CHKERRQ(ierr); 658 nloc = (Iend-Istart)/bs; my0 = Istart/bs; 659 PetscCheckFalse((Iend-Istart) % bs,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Iend %D - Istart %D must be divisible by bs %D",Iend,Istart,bs); 660 Iend /= bs; 661 nghosts = data_stride/bs - nloc; 662 663 ierr = PCGAMGHashTableCreate(2*nghosts+1, &fgid_flid);CHKERRQ(ierr); 664 for (kk=0; kk<nghosts; kk++) { 665 ierr = PCGAMGHashTableAdd(&fgid_flid, flid_fgid[nloc+kk], nloc+kk);CHKERRQ(ierr); 666 } 667 668 /* count selected -- same as number of cols of P */ 669 for (nSelected=mm=0; mm<nloc; mm++) { 670 PetscBool ise; 671 ierr = PetscCDEmptyAt(agg_llists, mm, &ise);CHKERRQ(ierr); 672 if (!ise) nSelected++; 673 } 674 ierr = MatGetOwnershipRangeColumn(a_Prol, &ii, &jj);CHKERRQ(ierr); 675 PetscCheckFalse((ii/nSAvec) != my0crs,PETSC_COMM_SELF,PETSC_ERR_PLIB,"ii %D /nSAvec %D != my0crs %D",ii,nSAvec,my0crs); 676 PetscCheckFalse(nSelected != (jj-ii)/nSAvec,PETSC_COMM_SELF,PETSC_ERR_PLIB,"nSelected %D != (jj %D - ii %D)/nSAvec %D",nSelected,jj,ii,nSAvec); 677 678 /* aloc space for coarse point data (output) */ 679 out_data_stride = nSelected*nSAvec; 680 681 ierr = PetscMalloc1(out_data_stride*nSAvec, &out_data);CHKERRQ(ierr); 682 for (ii=0;ii<out_data_stride*nSAvec;ii++) out_data[ii]=PETSC_MAX_REAL; 683 *a_data_out = out_data; /* output - stride nSelected*nSAvec */ 684 685 /* find points and set prolongation */ 686 minsz = 100; 687 for (mm = clid = 0; mm < nloc; mm++) { 688 ierr = PetscCDSizeAt(agg_llists, mm, &jj);CHKERRQ(ierr); 689 if (jj > 0) { 690 const PetscInt lid = mm, cgid = my0crs + clid; 691 PetscInt cids[100]; /* max bs */ 692 PetscBLASInt asz =jj,M=asz*bs,N=nSAvec,INFO; 693 PetscBLASInt Mdata=M+((N-M>0) ? N-M : 0),LDA=Mdata,LWORK=N*bs; 694 PetscScalar *qqc,*qqr,*TAU,*WORK; 695 PetscInt *fids; 696 PetscReal *data; 697 698 /* count agg */ 699 if (asz<minsz) minsz = asz; 700 701 /* get block */ 702 ierr = PetscMalloc5(Mdata*N, &qqc,M*N, &qqr,N, &TAU,LWORK, &WORK,M, &fids);CHKERRQ(ierr); 703 704 aggID = 0; 705 ierr = PetscCDGetHeadPos(agg_llists,lid,&pos);CHKERRQ(ierr); 706 while (pos) { 707 PetscInt gid1; 708 ierr = PetscCDIntNdGetID(pos, &gid1);CHKERRQ(ierr); 709 ierr = PetscCDGetNextPos(agg_llists,lid,&pos);CHKERRQ(ierr); 710 711 if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0; 712 else { 713 ierr = PCGAMGHashTableFind(&fgid_flid, gid1, &flid);CHKERRQ(ierr); 714 PetscCheckFalse(flid < 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Cannot find gid1 in table"); 715 } 716 /* copy in B_i matrix - column oriented */ 717 data = &data_in[flid*bs]; 718 for (ii = 0; ii < bs; ii++) { 719 for (jj = 0; jj < N; jj++) { 720 PetscReal d = data[jj*data_stride + ii]; 721 qqc[jj*Mdata + aggID*bs + ii] = d; 722 } 723 } 724 /* set fine IDs */ 725 for (kk=0; kk<bs; kk++) fids[aggID*bs + kk] = flid_fgid[flid]*bs + kk; 726 aggID++; 727 } 728 729 /* pad with zeros */ 730 for (ii = asz*bs; ii < Mdata; ii++) { 731 for (jj = 0; jj < N; jj++, kk++) { 732 qqc[jj*Mdata + ii] = .0; 733 } 734 } 735 736 /* QR */ 737 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 738 PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO)); 739 ierr = PetscFPTrapPop();CHKERRQ(ierr); 740 PetscCheckFalse(INFO != 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"xGEQRF error"); 741 /* get R - column oriented - output B_{i+1} */ 742 { 743 PetscReal *data = &out_data[clid*nSAvec]; 744 for (jj = 0; jj < nSAvec; jj++) { 745 for (ii = 0; ii < nSAvec; ii++) { 746 PetscCheckFalse(data[jj*out_data_stride + ii] != PETSC_MAX_REAL,PETSC_COMM_SELF,PETSC_ERR_PLIB,"data[jj*out_data_stride + ii] != %e",(double)PETSC_MAX_REAL); 747 if (ii <= jj) data[jj*out_data_stride + ii] = PetscRealPart(qqc[jj*Mdata + ii]); 748 else data[jj*out_data_stride + ii] = 0.; 749 } 750 } 751 } 752 753 /* get Q - row oriented */ 754 PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO)); 755 PetscCheckFalse(INFO != 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"xORGQR error arg %d",-INFO); 756 757 for (ii = 0; ii < M; ii++) { 758 for (jj = 0; jj < N; jj++) { 759 qqr[N*ii + jj] = qqc[jj*Mdata + ii]; 760 } 761 } 762 763 /* add diagonal block of P0 */ 764 for (kk=0; kk<N; kk++) { 765 cids[kk] = N*cgid + kk; /* global col IDs in P0 */ 766 } 767 ierr = MatSetValues(a_Prol,M,fids,N,cids,qqr,INSERT_VALUES);CHKERRQ(ierr); 768 ierr = PetscFree5(qqc,qqr,TAU,WORK,fids);CHKERRQ(ierr); 769 clid++; 770 } /* coarse agg */ 771 } /* for all fine nodes */ 772 ierr = MatAssemblyBegin(a_Prol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 773 ierr = MatAssemblyEnd(a_Prol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 774 ierr = PCGAMGHashTableDestroy(&fgid_flid);CHKERRQ(ierr); 775 PetscFunctionReturn(0); 776 } 777 778 static PetscErrorCode PCView_GAMG_AGG(PC pc,PetscViewer viewer) 779 { 780 PetscErrorCode ierr; 781 PC_MG *mg = (PC_MG*)pc->data; 782 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 783 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 784 785 PetscFunctionBegin; 786 ierr = PetscViewerASCIIPrintf(viewer," AGG specific options\n");CHKERRQ(ierr); 787 ierr = PetscViewerASCIIPrintf(viewer," Symmetric graph %s\n",pc_gamg_agg->sym_graph ? "true" : "false");CHKERRQ(ierr); 788 ierr = PetscViewerASCIIPrintf(viewer," Number of levels to square graph %D\n",pc_gamg_agg->square_graph);CHKERRQ(ierr); 789 ierr = PetscViewerASCIIPrintf(viewer," Number smoothing steps %D\n",pc_gamg_agg->nsmooths);CHKERRQ(ierr); 790 PetscFunctionReturn(0); 791 } 792 793 /* -------------------------------------------------------------------------- */ 794 /* 795 PCGAMGGraph_AGG 796 797 Input Parameter: 798 . pc - this 799 . Amat - matrix on this fine level 800 Output Parameter: 801 . a_Gmat - 802 */ 803 static PetscErrorCode PCGAMGGraph_AGG(PC pc,Mat Amat,Mat *a_Gmat) 804 { 805 PetscErrorCode ierr; 806 PC_MG *mg = (PC_MG*)pc->data; 807 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 808 const PetscReal vfilter = pc_gamg->threshold[pc_gamg->current_level]; 809 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 810 Mat Gmat; 811 MPI_Comm comm; 812 PetscBool /* set,flg , */ symm; 813 814 PetscFunctionBegin; 815 ierr = PetscObjectGetComm((PetscObject)Amat,&comm);CHKERRQ(ierr); 816 ierr = PetscLogEventBegin(PC_GAMGGraph_AGG,0,0,0,0);CHKERRQ(ierr); 817 818 /* ierr = MatIsSymmetricKnown(Amat, &set, &flg);CHKERRQ(ierr); || !(set && flg) -- this causes lot of symm calls */ 819 symm = (PetscBool)(pc_gamg_agg->sym_graph); /* && !pc_gamg_agg->square_graph; */ 820 821 ierr = PCGAMGCreateGraph(Amat, &Gmat);CHKERRQ(ierr); 822 ierr = PCGAMGFilterGraph(&Gmat, vfilter, symm);CHKERRQ(ierr); 823 *a_Gmat = Gmat; 824 ierr = PetscLogEventEnd(PC_GAMGGraph_AGG,0,0,0,0);CHKERRQ(ierr); 825 PetscFunctionReturn(0); 826 } 827 828 /* -------------------------------------------------------------------------- */ 829 /* 830 PCGAMGCoarsen_AGG 831 832 Input Parameter: 833 . a_pc - this 834 Input/Output Parameter: 835 . a_Gmat1 - graph on this fine level - coarsening can change this (squares it) 836 Output Parameter: 837 . agg_lists - list of aggregates 838 */ 839 static PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc,Mat *a_Gmat1,PetscCoarsenData **agg_lists) 840 { 841 PetscErrorCode ierr; 842 PC_MG *mg = (PC_MG*)a_pc->data; 843 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 844 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 845 Mat mat,Gmat2, Gmat1 = *a_Gmat1; /* squared graph */ 846 IS perm; 847 PetscInt Istart,Iend,Ii,nloc,bs,n,m; 848 PetscInt *permute; 849 PetscBool *bIndexSet; 850 MatCoarsen crs; 851 MPI_Comm comm; 852 PetscReal hashfact; 853 PetscInt iSwapIndex; 854 PetscRandom random; 855 856 PetscFunctionBegin; 857 ierr = PetscLogEventBegin(PC_GAMGCoarsen_AGG,0,0,0,0);CHKERRQ(ierr); 858 ierr = PetscObjectGetComm((PetscObject)Gmat1,&comm);CHKERRQ(ierr); 859 ierr = MatGetLocalSize(Gmat1, &n, &m);CHKERRQ(ierr); 860 ierr = MatGetBlockSize(Gmat1, &bs);CHKERRQ(ierr); 861 PetscCheckFalse(bs != 1,PETSC_COMM_SELF,PETSC_ERR_PLIB,"bs %D must be 1",bs); 862 nloc = n/bs; 863 864 if (pc_gamg->current_level < pc_gamg_agg->square_graph) { 865 ierr = PCGAMGSquareGraph_GAMG(a_pc,Gmat1,&Gmat2);CHKERRQ(ierr); 866 } else Gmat2 = Gmat1; 867 868 /* get MIS aggs - randomize */ 869 ierr = PetscMalloc1(nloc, &permute);CHKERRQ(ierr); 870 ierr = PetscCalloc1(nloc, &bIndexSet);CHKERRQ(ierr); 871 for (Ii = 0; Ii < nloc; Ii++) permute[Ii] = Ii; 872 ierr = PetscRandomCreate(PETSC_COMM_SELF,&random);CHKERRQ(ierr); 873 ierr = MatGetOwnershipRange(Gmat1, &Istart, &Iend);CHKERRQ(ierr); 874 for (Ii = 0; Ii < nloc; Ii++) { 875 ierr = PetscRandomGetValueReal(random,&hashfact);CHKERRQ(ierr); 876 iSwapIndex = (PetscInt) (hashfact*nloc)%nloc; 877 if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) { 878 PetscInt iTemp = permute[iSwapIndex]; 879 permute[iSwapIndex] = permute[Ii]; 880 permute[Ii] = iTemp; 881 bIndexSet[iSwapIndex] = PETSC_TRUE; 882 } 883 } 884 ierr = PetscFree(bIndexSet);CHKERRQ(ierr); 885 ierr = PetscRandomDestroy(&random);CHKERRQ(ierr); 886 ierr = ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm);CHKERRQ(ierr); 887 ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET4],0,0,0,0);CHKERRQ(ierr); 888 ierr = MatCoarsenCreate(comm, &crs);CHKERRQ(ierr); 889 ierr = MatCoarsenSetFromOptions(crs);CHKERRQ(ierr); 890 ierr = MatCoarsenSetGreedyOrdering(crs, perm);CHKERRQ(ierr); 891 ierr = MatCoarsenSetAdjacency(crs, Gmat2);CHKERRQ(ierr); 892 ierr = MatCoarsenSetStrictAggs(crs, PETSC_TRUE);CHKERRQ(ierr); 893 ierr = MatCoarsenApply(crs);CHKERRQ(ierr); 894 ierr = MatCoarsenGetData(crs, agg_lists);CHKERRQ(ierr); /* output */ 895 ierr = MatCoarsenDestroy(&crs);CHKERRQ(ierr); 896 897 ierr = ISDestroy(&perm);CHKERRQ(ierr); 898 ierr = PetscFree(permute);CHKERRQ(ierr); 899 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET4],0,0,0,0);CHKERRQ(ierr); 900 901 /* smooth aggs */ 902 if (Gmat2 != Gmat1) { 903 const PetscCoarsenData *llist = *agg_lists; 904 ierr = smoothAggs(a_pc,Gmat2, Gmat1, *agg_lists);CHKERRQ(ierr); 905 ierr = MatDestroy(&Gmat1);CHKERRQ(ierr); 906 *a_Gmat1 = Gmat2; /* output */ 907 ierr = PetscCDGetMat(llist, &mat);CHKERRQ(ierr); 908 PetscCheckFalse(mat,comm,PETSC_ERR_ARG_WRONG, "Auxilary matrix with squared graph????"); 909 } else { 910 const PetscCoarsenData *llist = *agg_lists; 911 /* see if we have a matrix that takes precedence (returned from MatCoarsenApply) */ 912 ierr = PetscCDGetMat(llist, &mat);CHKERRQ(ierr); 913 if (mat) { 914 ierr = MatDestroy(&Gmat1);CHKERRQ(ierr); 915 *a_Gmat1 = mat; /* output */ 916 } 917 } 918 ierr = PetscLogEventEnd(PC_GAMGCoarsen_AGG,0,0,0,0);CHKERRQ(ierr); 919 PetscFunctionReturn(0); 920 } 921 922 /* -------------------------------------------------------------------------- */ 923 /* 924 PCGAMGProlongator_AGG 925 926 Input Parameter: 927 . pc - this 928 . Amat - matrix on this fine level 929 . Graph - used to get ghost data for nodes in 930 . agg_lists - list of aggregates 931 Output Parameter: 932 . a_P_out - prolongation operator to the next level 933 */ 934 static PetscErrorCode PCGAMGProlongator_AGG(PC pc,Mat Amat,Mat Gmat,PetscCoarsenData *agg_lists,Mat *a_P_out) 935 { 936 PC_MG *mg = (PC_MG*)pc->data; 937 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 938 const PetscInt col_bs = pc_gamg->data_cell_cols; 939 PetscErrorCode ierr; 940 PetscInt Istart,Iend,nloc,ii,jj,kk,my0,nLocalSelected,bs; 941 Mat Prol; 942 PetscMPIInt size; 943 MPI_Comm comm; 944 PetscReal *data_w_ghost; 945 PetscInt myCrs0, nbnodes=0, *flid_fgid; 946 MatType mtype; 947 948 PetscFunctionBegin; 949 ierr = PetscObjectGetComm((PetscObject)Amat,&comm);CHKERRQ(ierr); 950 PetscCheckFalse(col_bs < 1,comm,PETSC_ERR_PLIB,"Column bs cannot be less than 1"); 951 ierr = PetscLogEventBegin(PC_GAMGProlongator_AGG,0,0,0,0);CHKERRQ(ierr); 952 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 953 ierr = MatGetOwnershipRange(Amat, &Istart, &Iend);CHKERRQ(ierr); 954 ierr = MatGetBlockSize(Amat, &bs);CHKERRQ(ierr); 955 nloc = (Iend-Istart)/bs; my0 = Istart/bs; 956 PetscCheckFalse((Iend-Istart) % bs,PETSC_COMM_SELF,PETSC_ERR_PLIB,"(Iend %D - Istart %D) not divisible by bs %D",Iend,Istart,bs); 957 958 /* get 'nLocalSelected' */ 959 for (ii=0, nLocalSelected = 0; ii < nloc; ii++) { 960 PetscBool ise; 961 /* filter out singletons 0 or 1? */ 962 ierr = PetscCDEmptyAt(agg_lists, ii, &ise);CHKERRQ(ierr); 963 if (!ise) nLocalSelected++; 964 } 965 966 /* create prolongator, create P matrix */ 967 ierr = MatGetType(Amat,&mtype);CHKERRQ(ierr); 968 ierr = MatCreate(comm, &Prol);CHKERRQ(ierr); 969 ierr = MatSetSizes(Prol,nloc*bs,nLocalSelected*col_bs,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 970 ierr = MatSetBlockSizes(Prol, bs, col_bs);CHKERRQ(ierr); 971 ierr = MatSetType(Prol, mtype);CHKERRQ(ierr); 972 ierr = MatSeqAIJSetPreallocation(Prol,col_bs, NULL);CHKERRQ(ierr); 973 ierr = MatMPIAIJSetPreallocation(Prol,col_bs, NULL, col_bs, NULL);CHKERRQ(ierr); 974 975 /* can get all points "removed" */ 976 ierr = MatGetSize(Prol, &kk, &ii);CHKERRQ(ierr); 977 if (!ii) { 978 ierr = PetscInfo(pc,"%s: No selected points on coarse grid\n");CHKERRQ(ierr); 979 ierr = MatDestroy(&Prol);CHKERRQ(ierr); 980 *a_P_out = NULL; /* out */ 981 ierr = PetscLogEventEnd(PC_GAMGProlongator_AGG,0,0,0,0);CHKERRQ(ierr); 982 PetscFunctionReturn(0); 983 } 984 ierr = PetscInfo(pc,"%s: New grid %D nodes\n",((PetscObject)pc)->prefix,ii/col_bs);CHKERRQ(ierr); 985 ierr = MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk);CHKERRQ(ierr); 986 987 PetscCheckFalse((kk-myCrs0) % col_bs,PETSC_COMM_SELF,PETSC_ERR_PLIB,"(kk %D -myCrs0 %D) not divisible by col_bs %D",kk,myCrs0,col_bs); 988 myCrs0 = myCrs0/col_bs; 989 PetscCheckFalse((kk/col_bs-myCrs0) != nLocalSelected,PETSC_COMM_SELF,PETSC_ERR_PLIB,"(kk %D/col_bs %D - myCrs0 %D) != nLocalSelected %D)",kk,col_bs,myCrs0,nLocalSelected); 990 991 /* create global vector of data in 'data_w_ghost' */ 992 ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET7],0,0,0,0);CHKERRQ(ierr); 993 if (size > 1) { /* */ 994 PetscReal *tmp_gdata,*tmp_ldata,*tp2; 995 ierr = PetscMalloc1(nloc, &tmp_ldata);CHKERRQ(ierr); 996 for (jj = 0; jj < col_bs; jj++) { 997 for (kk = 0; kk < bs; kk++) { 998 PetscInt ii,stride; 999 const PetscReal *tp = pc_gamg->data + jj*bs*nloc + kk; 1000 for (ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp; 1001 1002 ierr = PCGAMGGetDataWithGhosts(Gmat, 1, tmp_ldata, &stride, &tmp_gdata);CHKERRQ(ierr); 1003 1004 if (!jj && !kk) { /* now I know how many todal nodes - allocate */ 1005 ierr = PetscMalloc1(stride*bs*col_bs, &data_w_ghost);CHKERRQ(ierr); 1006 nbnodes = bs*stride; 1007 } 1008 tp2 = data_w_ghost + jj*bs*stride + kk; 1009 for (ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii]; 1010 ierr = PetscFree(tmp_gdata);CHKERRQ(ierr); 1011 } 1012 } 1013 ierr = PetscFree(tmp_ldata);CHKERRQ(ierr); 1014 } else { 1015 nbnodes = bs*nloc; 1016 data_w_ghost = (PetscReal*)pc_gamg->data; 1017 } 1018 1019 /* get P0 */ 1020 if (size > 1) { 1021 PetscReal *fid_glid_loc,*fiddata; 1022 PetscInt stride; 1023 1024 ierr = PetscMalloc1(nloc, &fid_glid_loc);CHKERRQ(ierr); 1025 for (kk=0; kk<nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0+kk); 1026 ierr = PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata);CHKERRQ(ierr); 1027 ierr = PetscMalloc1(stride, &flid_fgid);CHKERRQ(ierr); 1028 for (kk=0; kk<stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk]; 1029 ierr = PetscFree(fiddata);CHKERRQ(ierr); 1030 1031 PetscCheckFalse(stride != nbnodes/bs,PETSC_COMM_SELF,PETSC_ERR_PLIB,"stride %D != nbnodes %D/bs %D",stride,nbnodes,bs); 1032 ierr = PetscFree(fid_glid_loc);CHKERRQ(ierr); 1033 } else { 1034 ierr = PetscMalloc1(nloc, &flid_fgid);CHKERRQ(ierr); 1035 for (kk=0; kk<nloc; kk++) flid_fgid[kk] = my0 + kk; 1036 } 1037 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET7],0,0,0,0);CHKERRQ(ierr); 1038 ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET8],0,0,0,0);CHKERRQ(ierr); 1039 { 1040 PetscReal *data_out = NULL; 1041 ierr = formProl0(agg_lists, bs, col_bs, myCrs0, nbnodes,data_w_ghost, flid_fgid, &data_out, Prol);CHKERRQ(ierr); 1042 ierr = PetscFree(pc_gamg->data);CHKERRQ(ierr); 1043 1044 pc_gamg->data = data_out; 1045 pc_gamg->data_cell_rows = col_bs; 1046 pc_gamg->data_sz = col_bs*col_bs*nLocalSelected; 1047 } 1048 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET8],0,0,0,0);CHKERRQ(ierr); 1049 if (size > 1) {ierr = PetscFree(data_w_ghost);CHKERRQ(ierr);} 1050 ierr = PetscFree(flid_fgid);CHKERRQ(ierr); 1051 1052 *a_P_out = Prol; /* out */ 1053 1054 ierr = PetscLogEventEnd(PC_GAMGProlongator_AGG,0,0,0,0);CHKERRQ(ierr); 1055 PetscFunctionReturn(0); 1056 } 1057 1058 /* -------------------------------------------------------------------------- */ 1059 /* 1060 PCGAMGOptProlongator_AGG 1061 1062 Input Parameter: 1063 . pc - this 1064 . Amat - matrix on this fine level 1065 In/Output Parameter: 1066 . a_P - prolongation operator to the next level 1067 */ 1068 static PetscErrorCode PCGAMGOptProlongator_AGG(PC pc,Mat Amat,Mat *a_P) 1069 { 1070 PetscErrorCode ierr; 1071 PC_MG *mg = (PC_MG*)pc->data; 1072 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1073 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 1074 PetscInt jj; 1075 Mat Prol = *a_P; 1076 MPI_Comm comm; 1077 KSP eksp; 1078 Vec bb, xx; 1079 PC epc; 1080 PetscReal alpha, emax, emin; 1081 1082 PetscFunctionBegin; 1083 ierr = PetscObjectGetComm((PetscObject)Amat,&comm);CHKERRQ(ierr); 1084 ierr = PetscLogEventBegin(PC_GAMGOptProlongator_AGG,0,0,0,0);CHKERRQ(ierr); 1085 1086 /* compute maximum singular value of operator to be used in smoother */ 1087 if (0 < pc_gamg_agg->nsmooths) { 1088 /* get eigen estimates */ 1089 if (pc_gamg->emax > 0) { 1090 emin = pc_gamg->emin; 1091 emax = pc_gamg->emax; 1092 } else { 1093 const char *prefix; 1094 1095 ierr = MatCreateVecs(Amat, &bb, NULL);CHKERRQ(ierr); 1096 ierr = MatCreateVecs(Amat, &xx, NULL);CHKERRQ(ierr); 1097 ierr = KSPSetNoisy_Private(bb);CHKERRQ(ierr); 1098 1099 ierr = KSPCreate(comm,&eksp);CHKERRQ(ierr); 1100 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 1101 ierr = KSPSetOptionsPrefix(eksp,prefix);CHKERRQ(ierr); 1102 ierr = KSPAppendOptionsPrefix(eksp,"pc_gamg_esteig_");CHKERRQ(ierr); 1103 { 1104 PetscBool sflg; 1105 ierr = MatGetOption(Amat, MAT_SPD, &sflg);CHKERRQ(ierr); 1106 if (sflg) { 1107 ierr = KSPSetType(eksp, KSPCG);CHKERRQ(ierr); 1108 } 1109 } 1110 ierr = KSPSetErrorIfNotConverged(eksp,pc->erroriffailure);CHKERRQ(ierr); 1111 ierr = KSPSetNormType(eksp, KSP_NORM_NONE);CHKERRQ(ierr); 1112 1113 ierr = KSPSetInitialGuessNonzero(eksp, PETSC_FALSE);CHKERRQ(ierr); 1114 ierr = KSPSetOperators(eksp, Amat, Amat);CHKERRQ(ierr); 1115 1116 ierr = KSPGetPC(eksp, &epc);CHKERRQ(ierr); 1117 ierr = PCSetType(epc, PCJACOBI);CHKERRQ(ierr); /* smoother in smoothed agg. */ 1118 1119 ierr = KSPSetTolerances(eksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT, 10);CHKERRQ(ierr); // 10 is safer, but 5 is often fine, can override with -pc_gamg_esteig_ksp_max_it -mg_levels_ksp_chebyshev_esteig 0,0.25,0,1.2 1120 1121 ierr = KSPSetFromOptions(eksp);CHKERRQ(ierr); 1122 ierr = KSPSetComputeSingularValues(eksp,PETSC_TRUE);CHKERRQ(ierr); 1123 ierr = KSPSolve(eksp, bb, xx);CHKERRQ(ierr); 1124 ierr = KSPCheckSolve(eksp,pc,xx);CHKERRQ(ierr); 1125 1126 ierr = KSPComputeExtremeSingularValues(eksp, &emax, &emin);CHKERRQ(ierr); 1127 ierr = PetscInfo(pc,"%s: Smooth P0: max eigen=%e min=%e PC=%s\n",((PetscObject)pc)->prefix,(double)emax,(double)emin,PCJACOBI);CHKERRQ(ierr); 1128 ierr = VecDestroy(&xx);CHKERRQ(ierr); 1129 ierr = VecDestroy(&bb);CHKERRQ(ierr); 1130 ierr = KSPDestroy(&eksp);CHKERRQ(ierr); 1131 } 1132 if (pc_gamg->use_sa_esteig) { 1133 mg->min_eigen_DinvA[pc_gamg->current_level] = emin; 1134 mg->max_eigen_DinvA[pc_gamg->current_level] = emax; 1135 ierr = PetscInfo(pc,"%s: Smooth P0: level %D, cache spectra %g %g\n",((PetscObject)pc)->prefix,pc_gamg->current_level,(double)emin,(double)emax);CHKERRQ(ierr); 1136 } else { 1137 mg->min_eigen_DinvA[pc_gamg->current_level] = 0; 1138 mg->max_eigen_DinvA[pc_gamg->current_level] = 0; 1139 } 1140 } else { 1141 mg->min_eigen_DinvA[pc_gamg->current_level] = 0; 1142 mg->max_eigen_DinvA[pc_gamg->current_level] = 0; 1143 } 1144 1145 /* smooth P0 */ 1146 for (jj = 0; jj < pc_gamg_agg->nsmooths; jj++) { 1147 Mat tMat; 1148 Vec diag; 1149 1150 ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET9],0,0,0,0);CHKERRQ(ierr); 1151 1152 /* smooth P1 := (I - omega/lam D^{-1}A)P0 */ 1153 ierr = PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2],0,0,0,0);CHKERRQ(ierr); 1154 ierr = MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &tMat);CHKERRQ(ierr); 1155 ierr = PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2],0,0,0,0);CHKERRQ(ierr); 1156 ierr = MatProductClear(tMat);CHKERRQ(ierr); 1157 ierr = MatCreateVecs(Amat, &diag, NULL);CHKERRQ(ierr); 1158 ierr = MatGetDiagonal(Amat, diag);CHKERRQ(ierr); /* effectively PCJACOBI */ 1159 ierr = VecReciprocal(diag);CHKERRQ(ierr); 1160 ierr = MatDiagonalScale(tMat, diag, NULL);CHKERRQ(ierr); 1161 ierr = VecDestroy(&diag);CHKERRQ(ierr); 1162 1163 /* TODO: Set a PCFailedReason and exit the building of the AMG preconditioner */ 1164 PetscCheckFalse(emax == 0.0,PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Computed maximum singular value as zero"); 1165 /* TODO: Document the 1.4 and don't hardwire it in this routine */ 1166 alpha = -1.4/emax; 1167 1168 ierr = MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN);CHKERRQ(ierr); 1169 ierr = MatDestroy(&Prol);CHKERRQ(ierr); 1170 Prol = tMat; 1171 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET9],0,0,0,0);CHKERRQ(ierr); 1172 } 1173 ierr = PetscLogEventEnd(PC_GAMGOptProlongator_AGG,0,0,0,0);CHKERRQ(ierr); 1174 *a_P = Prol; 1175 PetscFunctionReturn(0); 1176 } 1177 1178 /* -------------------------------------------------------------------------- */ 1179 /* 1180 PCCreateGAMG_AGG 1181 1182 Input Parameter: 1183 . pc - 1184 */ 1185 PetscErrorCode PCCreateGAMG_AGG(PC pc) 1186 { 1187 PetscErrorCode ierr; 1188 PC_MG *mg = (PC_MG*)pc->data; 1189 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1190 PC_GAMG_AGG *pc_gamg_agg; 1191 1192 PetscFunctionBegin; 1193 /* create sub context for SA */ 1194 ierr = PetscNewLog(pc,&pc_gamg_agg);CHKERRQ(ierr); 1195 pc_gamg->subctx = pc_gamg_agg; 1196 1197 pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG; 1198 pc_gamg->ops->destroy = PCDestroy_GAMG_AGG; 1199 /* reset does not do anything; setup not virtual */ 1200 1201 /* set internal function pointers */ 1202 pc_gamg->ops->graph = PCGAMGGraph_AGG; 1203 pc_gamg->ops->coarsen = PCGAMGCoarsen_AGG; 1204 pc_gamg->ops->prolongator = PCGAMGProlongator_AGG; 1205 pc_gamg->ops->optprolongator = PCGAMGOptProlongator_AGG; 1206 pc_gamg->ops->createdefaultdata = PCSetData_AGG; 1207 pc_gamg->ops->view = PCView_GAMG_AGG; 1208 1209 pc_gamg_agg->square_graph = 1; 1210 pc_gamg_agg->sym_graph = PETSC_FALSE; 1211 pc_gamg_agg->nsmooths = 1; 1212 1213 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetNSmooths_C",PCGAMGSetNSmooths_AGG);CHKERRQ(ierr); 1214 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetSymGraph_C",PCGAMGSetSymGraph_AGG);CHKERRQ(ierr); 1215 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetSquareGraph_C",PCGAMGSetSquareGraph_AGG);CHKERRQ(ierr); 1216 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_AGG);CHKERRQ(ierr); 1217 PetscFunctionReturn(0); 1218 } 1219