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