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