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