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