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=0; 261 MPI_Comm comm; 262 PetscInt lid,*ii,*idx,ix,Iend,my0,kk,n,j; 263 Mat_MPIAIJ *mpimat_2 = 0, *mpimat_1=0; 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, 0);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, 0);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 #if defined(PETSC_MISSING_LAPACK_ORGQR) 763 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"ORGQR - Lapack routine is unavailable."); 764 #else 765 PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO)); 766 if (INFO != 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xORGQR error arg %d",-INFO); 767 #endif 768 769 for (ii = 0; ii < M; ii++) { 770 for (jj = 0; jj < N; jj++) { 771 qqr[N*ii + jj] = qqc[jj*Mdata + ii]; 772 } 773 } 774 775 /* add diagonal block of P0 */ 776 for (kk=0; kk<N; kk++) { 777 cids[kk] = N*cgid + kk; /* global col IDs in P0 */ 778 } 779 ierr = MatSetValues(a_Prol,M,fids,N,cids,qqr,INSERT_VALUES);CHKERRQ(ierr); 780 ierr = PetscFree5(qqc,qqr,TAU,WORK,fids);CHKERRQ(ierr); 781 clid++; 782 } /* coarse agg */ 783 } /* for all fine nodes */ 784 ierr = MatAssemblyBegin(a_Prol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 785 ierr = MatAssemblyEnd(a_Prol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 786 ierr = PCGAMGHashTableDestroy(&fgid_flid);CHKERRQ(ierr); 787 PetscFunctionReturn(0); 788 } 789 790 static PetscErrorCode PCView_GAMG_AGG(PC pc,PetscViewer viewer) 791 { 792 PetscErrorCode ierr; 793 PC_MG *mg = (PC_MG*)pc->data; 794 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 795 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 796 797 PetscFunctionBegin; 798 ierr = PetscViewerASCIIPrintf(viewer," AGG specific options\n");CHKERRQ(ierr); 799 ierr = PetscViewerASCIIPrintf(viewer," Symmetric graph %s\n",pc_gamg_agg->sym_graph ? "true" : "false");CHKERRQ(ierr); 800 ierr = PetscViewerASCIIPrintf(viewer," Number of levels to square graph %D\n",pc_gamg_agg->square_graph);CHKERRQ(ierr); 801 ierr = PetscViewerASCIIPrintf(viewer," Number smoothing steps %D\n",pc_gamg_agg->nsmooths);CHKERRQ(ierr); 802 PetscFunctionReturn(0); 803 } 804 805 /* -------------------------------------------------------------------------- */ 806 /* 807 PCGAMGGraph_AGG 808 809 Input Parameter: 810 . pc - this 811 . Amat - matrix on this fine level 812 Output Parameter: 813 . a_Gmat - 814 */ 815 static PetscErrorCode PCGAMGGraph_AGG(PC pc,Mat Amat,Mat *a_Gmat) 816 { 817 PetscErrorCode ierr; 818 PC_MG *mg = (PC_MG*)pc->data; 819 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 820 const PetscReal vfilter = pc_gamg->threshold[pc_gamg->current_level]; 821 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 822 Mat Gmat; 823 MPI_Comm comm; 824 PetscBool /* set,flg , */ symm; 825 826 PetscFunctionBegin; 827 ierr = PetscObjectGetComm((PetscObject)Amat,&comm);CHKERRQ(ierr); 828 ierr = PetscLogEventBegin(PC_GAMGGraph_AGG,0,0,0,0);CHKERRQ(ierr); 829 830 /* ierr = MatIsSymmetricKnown(Amat, &set, &flg);CHKERRQ(ierr); || !(set && flg) -- this causes lot of symm calls */ 831 symm = (PetscBool)(pc_gamg_agg->sym_graph); /* && !pc_gamg_agg->square_graph; */ 832 833 ierr = PCGAMGCreateGraph(Amat, &Gmat);CHKERRQ(ierr); 834 ierr = PCGAMGFilterGraph(&Gmat, vfilter, symm);CHKERRQ(ierr); 835 *a_Gmat = Gmat; 836 ierr = PetscLogEventEnd(PC_GAMGGraph_AGG,0,0,0,0);CHKERRQ(ierr); 837 PetscFunctionReturn(0); 838 } 839 840 /* -------------------------------------------------------------------------- */ 841 /* 842 PCGAMGCoarsen_AGG 843 844 Input Parameter: 845 . a_pc - this 846 Input/Output Parameter: 847 . a_Gmat1 - graph on this fine level - coarsening can change this (squares it) 848 Output Parameter: 849 . agg_lists - list of aggregates 850 */ 851 static PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc,Mat *a_Gmat1,PetscCoarsenData **agg_lists) 852 { 853 PetscErrorCode ierr; 854 PC_MG *mg = (PC_MG*)a_pc->data; 855 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 856 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 857 Mat mat,Gmat2, Gmat1 = *a_Gmat1; /* squared graph */ 858 IS perm; 859 PetscInt Istart,Iend,Ii,nloc,bs,n,m; 860 PetscInt *permute; 861 PetscBool *bIndexSet; 862 MatCoarsen crs; 863 MPI_Comm comm; 864 PetscMPIInt rank; 865 PetscReal hashfact; 866 PetscInt iSwapIndex; 867 PetscRandom random; 868 869 PetscFunctionBegin; 870 ierr = PetscLogEventBegin(PC_GAMGCoarsen_AGG,0,0,0,0);CHKERRQ(ierr); 871 ierr = PetscObjectGetComm((PetscObject)Gmat1,&comm);CHKERRQ(ierr); 872 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 873 ierr = MatGetLocalSize(Gmat1, &n, &m);CHKERRQ(ierr); 874 ierr = MatGetBlockSize(Gmat1, &bs);CHKERRQ(ierr); 875 if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"bs %D must be 1",bs); 876 nloc = n/bs; 877 878 if (pc_gamg->current_level < pc_gamg_agg->square_graph) { 879 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); 880 ierr = MatTransposeMatMult(Gmat1, Gmat1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat2);CHKERRQ(ierr); 881 } else Gmat2 = Gmat1; 882 883 /* get MIS aggs - randomize */ 884 ierr = PetscMalloc1(nloc, &permute);CHKERRQ(ierr); 885 ierr = PetscCalloc1(nloc, &bIndexSet);CHKERRQ(ierr); 886 for (Ii = 0; Ii < nloc; Ii++) { 887 permute[Ii] = Ii; 888 } 889 ierr = PetscRandomCreate(PETSC_COMM_SELF,&random);CHKERRQ(ierr); 890 ierr = MatGetOwnershipRange(Gmat1, &Istart, &Iend);CHKERRQ(ierr); 891 for (Ii = 0; Ii < nloc; Ii++) { 892 ierr = PetscRandomGetValueReal(random,&hashfact);CHKERRQ(ierr); 893 iSwapIndex = (PetscInt) (hashfact*nloc)%nloc; 894 if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) { 895 PetscInt iTemp = permute[iSwapIndex]; 896 permute[iSwapIndex] = permute[Ii]; 897 permute[Ii] = iTemp; 898 bIndexSet[iSwapIndex] = PETSC_TRUE; 899 } 900 } 901 ierr = PetscFree(bIndexSet);CHKERRQ(ierr); 902 ierr = PetscRandomDestroy(&random);CHKERRQ(ierr); 903 ierr = ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm);CHKERRQ(ierr); 904 #if defined PETSC_GAMG_USE_LOG 905 ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET4],0,0,0,0);CHKERRQ(ierr); 906 #endif 907 ierr = MatCoarsenCreate(comm, &crs);CHKERRQ(ierr); 908 ierr = MatCoarsenSetFromOptions(crs);CHKERRQ(ierr); 909 ierr = MatCoarsenSetGreedyOrdering(crs, perm);CHKERRQ(ierr); 910 ierr = MatCoarsenSetAdjacency(crs, Gmat2);CHKERRQ(ierr); 911 ierr = MatCoarsenSetStrictAggs(crs, PETSC_TRUE);CHKERRQ(ierr); 912 ierr = MatCoarsenApply(crs);CHKERRQ(ierr); 913 ierr = MatCoarsenGetData(crs, agg_lists);CHKERRQ(ierr); /* output */ 914 ierr = MatCoarsenDestroy(&crs);CHKERRQ(ierr); 915 916 ierr = ISDestroy(&perm);CHKERRQ(ierr); 917 ierr = PetscFree(permute);CHKERRQ(ierr); 918 #if defined PETSC_GAMG_USE_LOG 919 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET4],0,0,0,0);CHKERRQ(ierr); 920 #endif 921 922 /* smooth aggs */ 923 if (Gmat2 != Gmat1) { 924 const PetscCoarsenData *llist = *agg_lists; 925 ierr = smoothAggs(a_pc,Gmat2, Gmat1, *agg_lists);CHKERRQ(ierr); 926 ierr = MatDestroy(&Gmat1);CHKERRQ(ierr); 927 *a_Gmat1 = Gmat2; /* output */ 928 ierr = PetscCDGetMat(llist, &mat);CHKERRQ(ierr); 929 if (mat) SETERRQ(comm,PETSC_ERR_ARG_WRONG, "Auxilary matrix with squared graph????"); 930 } else { 931 const PetscCoarsenData *llist = *agg_lists; 932 /* see if we have a matrix that takes precedence (returned from MatCoarsenApply) */ 933 ierr = PetscCDGetMat(llist, &mat);CHKERRQ(ierr); 934 if (mat) { 935 ierr = MatDestroy(&Gmat1);CHKERRQ(ierr); 936 *a_Gmat1 = mat; /* output */ 937 } 938 } 939 ierr = PetscLogEventEnd(PC_GAMGCoarsen_AGG,0,0,0,0);CHKERRQ(ierr); 940 PetscFunctionReturn(0); 941 } 942 943 /* -------------------------------------------------------------------------- */ 944 /* 945 PCGAMGProlongator_AGG 946 947 Input Parameter: 948 . pc - this 949 . Amat - matrix on this fine level 950 . Graph - used to get ghost data for nodes in 951 . agg_lists - list of aggregates 952 Output Parameter: 953 . a_P_out - prolongation operator to the next level 954 */ 955 static PetscErrorCode PCGAMGProlongator_AGG(PC pc,Mat Amat,Mat Gmat,PetscCoarsenData *agg_lists,Mat *a_P_out) 956 { 957 PC_MG *mg = (PC_MG*)pc->data; 958 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 959 const PetscInt col_bs = pc_gamg->data_cell_cols; 960 PetscErrorCode ierr; 961 PetscInt Istart,Iend,nloc,ii,jj,kk,my0,nLocalSelected,bs; 962 Mat Prol; 963 PetscMPIInt rank, size; 964 MPI_Comm comm; 965 PetscReal *data_w_ghost; 966 PetscInt myCrs0, nbnodes=0, *flid_fgid; 967 MatType mtype; 968 969 PetscFunctionBegin; 970 ierr = PetscObjectGetComm((PetscObject)Amat,&comm);CHKERRQ(ierr); 971 if (col_bs < 1) SETERRQ(comm,PETSC_ERR_PLIB,"Column bs cannot be less than 1"); 972 ierr = PetscLogEventBegin(PC_GAMGProlongator_AGG,0,0,0,0);CHKERRQ(ierr); 973 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 974 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 975 ierr = MatGetOwnershipRange(Amat, &Istart, &Iend);CHKERRQ(ierr); 976 ierr = MatGetBlockSize(Amat, &bs);CHKERRQ(ierr); 977 nloc = (Iend-Istart)/bs; my0 = Istart/bs; 978 if ((Iend-Istart) % bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"(Iend %D - Istart %D) not divisible by bs %D",Iend,Istart,bs); 979 980 /* get 'nLocalSelected' */ 981 for (ii=0, nLocalSelected = 0; ii < nloc; ii++) { 982 PetscBool ise; 983 /* filter out singletons 0 or 1? */ 984 ierr = PetscCDEmptyAt(agg_lists, ii, &ise);CHKERRQ(ierr); 985 if (!ise) nLocalSelected++; 986 } 987 988 /* create prolongator, create P matrix */ 989 ierr = MatGetType(Amat,&mtype);CHKERRQ(ierr); 990 ierr = MatCreate(comm, &Prol);CHKERRQ(ierr); 991 ierr = MatSetSizes(Prol,nloc*bs,nLocalSelected*col_bs,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 992 ierr = MatSetBlockSizes(Prol, bs, col_bs);CHKERRQ(ierr); 993 ierr = MatSetType(Prol, mtype);CHKERRQ(ierr); 994 ierr = MatSeqAIJSetPreallocation(Prol,col_bs, NULL);CHKERRQ(ierr); 995 ierr = MatMPIAIJSetPreallocation(Prol,col_bs, NULL, col_bs, NULL);CHKERRQ(ierr); 996 997 /* can get all points "removed" */ 998 ierr = MatGetSize(Prol, &kk, &ii);CHKERRQ(ierr); 999 if (!ii) { 1000 ierr = PetscInfo(pc,"No selected points on coarse grid\n");CHKERRQ(ierr); 1001 ierr = MatDestroy(&Prol);CHKERRQ(ierr); 1002 *a_P_out = NULL; /* out */ 1003 ierr = PetscLogEventEnd(PC_GAMGProlongator_AGG,0,0,0,0);CHKERRQ(ierr); 1004 PetscFunctionReturn(0); 1005 } 1006 ierr = PetscInfo1(pc,"New grid %D nodes\n",ii/col_bs);CHKERRQ(ierr); 1007 ierr = MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk);CHKERRQ(ierr); 1008 1009 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); 1010 myCrs0 = myCrs0/col_bs; 1011 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); 1012 1013 /* create global vector of data in 'data_w_ghost' */ 1014 #if defined PETSC_GAMG_USE_LOG 1015 ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET7],0,0,0,0);CHKERRQ(ierr); 1016 #endif 1017 if (size > 1) { /* */ 1018 PetscReal *tmp_gdata,*tmp_ldata,*tp2; 1019 ierr = PetscMalloc1(nloc, &tmp_ldata);CHKERRQ(ierr); 1020 for (jj = 0; jj < col_bs; jj++) { 1021 for (kk = 0; kk < bs; kk++) { 1022 PetscInt ii,stride; 1023 const PetscReal *tp = pc_gamg->data + jj*bs*nloc + kk; 1024 for (ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp; 1025 1026 ierr = PCGAMGGetDataWithGhosts(Gmat, 1, tmp_ldata, &stride, &tmp_gdata);CHKERRQ(ierr); 1027 1028 if (!jj && !kk) { /* now I know how many todal nodes - allocate */ 1029 ierr = PetscMalloc1(stride*bs*col_bs, &data_w_ghost);CHKERRQ(ierr); 1030 nbnodes = bs*stride; 1031 } 1032 tp2 = data_w_ghost + jj*bs*stride + kk; 1033 for (ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii]; 1034 ierr = PetscFree(tmp_gdata);CHKERRQ(ierr); 1035 } 1036 } 1037 ierr = PetscFree(tmp_ldata);CHKERRQ(ierr); 1038 } else { 1039 nbnodes = bs*nloc; 1040 data_w_ghost = (PetscReal*)pc_gamg->data; 1041 } 1042 1043 /* get P0 */ 1044 if (size > 1) { 1045 PetscReal *fid_glid_loc,*fiddata; 1046 PetscInt stride; 1047 1048 ierr = PetscMalloc1(nloc, &fid_glid_loc);CHKERRQ(ierr); 1049 for (kk=0; kk<nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0+kk); 1050 ierr = PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata);CHKERRQ(ierr); 1051 ierr = PetscMalloc1(stride, &flid_fgid);CHKERRQ(ierr); 1052 for (kk=0; kk<stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk]; 1053 ierr = PetscFree(fiddata);CHKERRQ(ierr); 1054 1055 if (stride != nbnodes/bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"stride %D != nbnodes %D/bs %D",stride,nbnodes,bs); 1056 ierr = PetscFree(fid_glid_loc);CHKERRQ(ierr); 1057 } else { 1058 ierr = PetscMalloc1(nloc, &flid_fgid);CHKERRQ(ierr); 1059 for (kk=0; kk<nloc; kk++) flid_fgid[kk] = my0 + kk; 1060 } 1061 #if defined PETSC_GAMG_USE_LOG 1062 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET7],0,0,0,0);CHKERRQ(ierr); 1063 ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET8],0,0,0,0);CHKERRQ(ierr); 1064 #endif 1065 { 1066 PetscReal *data_out = NULL; 1067 ierr = formProl0(agg_lists, bs, col_bs, myCrs0, nbnodes,data_w_ghost, flid_fgid, &data_out, Prol);CHKERRQ(ierr); 1068 ierr = PetscFree(pc_gamg->data);CHKERRQ(ierr); 1069 1070 pc_gamg->data = data_out; 1071 pc_gamg->data_cell_rows = col_bs; 1072 pc_gamg->data_sz = col_bs*col_bs*nLocalSelected; 1073 } 1074 #if defined PETSC_GAMG_USE_LOG 1075 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET8],0,0,0,0);CHKERRQ(ierr); 1076 #endif 1077 if (size > 1) {ierr = PetscFree(data_w_ghost);CHKERRQ(ierr);} 1078 ierr = PetscFree(flid_fgid);CHKERRQ(ierr); 1079 1080 *a_P_out = Prol; /* out */ 1081 1082 ierr = PetscLogEventEnd(PC_GAMGProlongator_AGG,0,0,0,0);CHKERRQ(ierr); 1083 PetscFunctionReturn(0); 1084 } 1085 1086 /* -------------------------------------------------------------------------- */ 1087 /* 1088 PCGAMGOptProlongator_AGG 1089 1090 Input Parameter: 1091 . pc - this 1092 . Amat - matrix on this fine level 1093 In/Output Parameter: 1094 . a_P - prolongation operator to the next level 1095 */ 1096 static PetscErrorCode PCGAMGOptProlongator_AGG(PC pc,Mat Amat,Mat *a_P) 1097 { 1098 PetscErrorCode ierr; 1099 PC_MG *mg = (PC_MG*)pc->data; 1100 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1101 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 1102 PetscInt jj; 1103 Mat Prol = *a_P; 1104 MPI_Comm comm; 1105 KSP eksp; 1106 Vec bb, xx; 1107 PC epc; 1108 PetscReal alpha, emax, emin; 1109 PetscRandom random; 1110 1111 PetscFunctionBegin; 1112 ierr = PetscObjectGetComm((PetscObject)Amat,&comm);CHKERRQ(ierr); 1113 ierr = PetscLogEventBegin(PC_GAMGOptProlongator_AGG,0,0,0,0);CHKERRQ(ierr); 1114 1115 /* compute maximum value of operator to be used in smoother */ 1116 if (0 < pc_gamg_agg->nsmooths) { 1117 ierr = MatCreateVecs(Amat, &bb, 0);CHKERRQ(ierr); 1118 ierr = MatCreateVecs(Amat, &xx, 0);CHKERRQ(ierr); 1119 ierr = PetscRandomCreate(PETSC_COMM_SELF,&random);CHKERRQ(ierr); 1120 ierr = VecSetRandom(bb,random);CHKERRQ(ierr); 1121 ierr = PetscRandomDestroy(&random);CHKERRQ(ierr); 1122 1123 ierr = KSPCreate(comm,&eksp);CHKERRQ(ierr); 1124 ierr = KSPSetErrorIfNotConverged(eksp,pc->erroriffailure);CHKERRQ(ierr); 1125 ierr = KSPSetTolerances(eksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,10);CHKERRQ(ierr); 1126 ierr = KSPSetNormType(eksp, KSP_NORM_NONE);CHKERRQ(ierr); 1127 1128 ierr = KSPSetInitialGuessNonzero(eksp, PETSC_FALSE);CHKERRQ(ierr); 1129 ierr = KSPSetOperators(eksp, Amat, Amat);CHKERRQ(ierr); 1130 ierr = KSPSetComputeSingularValues(eksp,PETSC_TRUE);CHKERRQ(ierr); 1131 1132 ierr = KSPGetPC(eksp, &epc);CHKERRQ(ierr); 1133 ierr = PCSetType(epc, PCJACOBI);CHKERRQ(ierr); /* smoother in smoothed agg. */ 1134 1135 ierr = KSPSetOptionsPrefix(eksp,((PetscObject)pc)->prefix);CHKERRQ(ierr); 1136 ierr = KSPAppendOptionsPrefix(eksp, "gamg_est_");CHKERRQ(ierr); 1137 ierr = KSPSetFromOptions(eksp);CHKERRQ(ierr); 1138 1139 /* solve - keep stuff out of logging */ 1140 ierr = PetscLogEventDeactivate(KSP_Solve);CHKERRQ(ierr); 1141 ierr = PetscLogEventDeactivate(PC_Apply);CHKERRQ(ierr); 1142 ierr = KSPSolve(eksp, bb, xx);CHKERRQ(ierr); 1143 ierr = KSPCheckSolve(eksp,pc,xx);CHKERRQ(ierr); 1144 ierr = PetscLogEventActivate(KSP_Solve);CHKERRQ(ierr); 1145 ierr = PetscLogEventActivate(PC_Apply);CHKERRQ(ierr); 1146 1147 ierr = KSPComputeExtremeSingularValues(eksp, &emax, &emin);CHKERRQ(ierr); 1148 ierr = PetscInfo3(pc,"Smooth P0: max eigen=%e min=%e PC=%s\n",emax,emin,PCJACOBI);CHKERRQ(ierr); 1149 ierr = VecDestroy(&xx);CHKERRQ(ierr); 1150 ierr = VecDestroy(&bb);CHKERRQ(ierr); 1151 ierr = KSPDestroy(&eksp);CHKERRQ(ierr); 1152 } 1153 1154 /* smooth P0 */ 1155 for (jj = 0; jj < pc_gamg_agg->nsmooths; jj++) { 1156 Mat tMat; 1157 Vec diag; 1158 1159 #if defined PETSC_GAMG_USE_LOG 1160 ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET9],0,0,0,0);CHKERRQ(ierr); 1161 #endif 1162 1163 /* smooth P1 := (I - omega/lam D^{-1}A)P0 */ 1164 ierr = MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &tMat);CHKERRQ(ierr); 1165 ierr = MatCreateVecs(Amat, &diag, 0);CHKERRQ(ierr); 1166 ierr = MatGetDiagonal(Amat, diag);CHKERRQ(ierr); /* effectively PCJACOBI */ 1167 ierr = VecReciprocal(diag);CHKERRQ(ierr); 1168 ierr = MatDiagonalScale(tMat, diag, 0);CHKERRQ(ierr); 1169 ierr = VecDestroy(&diag);CHKERRQ(ierr); 1170 alpha = -1.4/emax; 1171 ierr = MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN);CHKERRQ(ierr); 1172 ierr = MatDestroy(&Prol);CHKERRQ(ierr); 1173 Prol = tMat; 1174 #if defined PETSC_GAMG_USE_LOG 1175 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET9],0,0,0,0);CHKERRQ(ierr); 1176 #endif 1177 } 1178 ierr = PetscLogEventEnd(PC_GAMGOptProlongator_AGG,0,0,0,0);CHKERRQ(ierr); 1179 *a_P = Prol; 1180 PetscFunctionReturn(0); 1181 } 1182 1183 /* -------------------------------------------------------------------------- */ 1184 /* 1185 PCCreateGAMG_AGG 1186 1187 Input Parameter: 1188 . pc - 1189 */ 1190 PetscErrorCode PCCreateGAMG_AGG(PC pc) 1191 { 1192 PetscErrorCode ierr; 1193 PC_MG *mg = (PC_MG*)pc->data; 1194 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1195 PC_GAMG_AGG *pc_gamg_agg; 1196 1197 PetscFunctionBegin; 1198 /* create sub context for SA */ 1199 ierr = PetscNewLog(pc,&pc_gamg_agg);CHKERRQ(ierr); 1200 pc_gamg->subctx = pc_gamg_agg; 1201 1202 pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG; 1203 pc_gamg->ops->destroy = PCDestroy_GAMG_AGG; 1204 /* reset does not do anything; setup not virtual */ 1205 1206 /* set internal function pointers */ 1207 pc_gamg->ops->graph = PCGAMGGraph_AGG; 1208 pc_gamg->ops->coarsen = PCGAMGCoarsen_AGG; 1209 pc_gamg->ops->prolongator = PCGAMGProlongator_AGG; 1210 pc_gamg->ops->optprolongator = PCGAMGOptProlongator_AGG; 1211 pc_gamg->ops->createdefaultdata = PCSetData_AGG; 1212 pc_gamg->ops->view = PCView_GAMG_AGG; 1213 1214 pc_gamg_agg->square_graph = 1; 1215 pc_gamg_agg->sym_graph = PETSC_FALSE; 1216 pc_gamg_agg->nsmooths = 1; 1217 1218 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetNSmooths_C",PCGAMGSetNSmooths_AGG);CHKERRQ(ierr); 1219 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetSymGraph_C",PCGAMGSetSymGraph_AGG);CHKERRQ(ierr); 1220 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetSquareGraph_C",PCGAMGSetSquareGraph_AGG);CHKERRQ(ierr); 1221 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_AGG);CHKERRQ(ierr); 1222 PetscFunctionReturn(0); 1223 } 1224