1 /* 2 GAMG geometric-algebric multiogrid PC - Mark Adams 2011 3 */ 4 5 #include <../src/ksp/pc/impls/gamg/gamg.h> /*I "petscpc.h" I*/ 6 #include <petsc-private/kspimpl.h> 7 8 #include <assert.h> 9 #include <petscblaslapack.h> 10 11 typedef struct { 12 PetscInt nsmooths; 13 PetscBool sym_graph; 14 PetscBool square_graph; 15 }PC_GAMG_AGG; 16 17 #undef __FUNCT__ 18 #define __FUNCT__ "PCGAMGSetNSmooths" 19 /*@ 20 PCGAMGSetNSmooths - Set number of smoothing steps (1 is typical) 21 22 Not Collective on PC 23 24 Input Parameters: 25 . pc - the preconditioner context 26 27 Options Database Key: 28 . -pc_gamg_agg_nsmooths 29 30 Level: intermediate 31 32 Concepts: Aggregation AMG preconditioner 33 34 .seealso: () 35 @*/ 36 PetscErrorCode PCGAMGSetNSmooths(PC pc, PetscInt n) 37 { 38 PetscErrorCode ierr; 39 40 PetscFunctionBegin; 41 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 42 ierr = PetscTryMethod(pc,"PCGAMGSetNSmooths_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr); 43 PetscFunctionReturn(0); 44 } 45 46 EXTERN_C_BEGIN 47 #undef __FUNCT__ 48 #define __FUNCT__ "PCGAMGSetNSmooths_GAMG" 49 PetscErrorCode PCGAMGSetNSmooths_GAMG(PC pc, PetscInt n) 50 { 51 PC_MG *mg = (PC_MG*)pc->data; 52 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 53 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 54 55 PetscFunctionBegin; 56 pc_gamg_agg->nsmooths = n; 57 PetscFunctionReturn(0); 58 } 59 EXTERN_C_END 60 61 #undef __FUNCT__ 62 #define __FUNCT__ "PCGAMGSetSymGraph" 63 /*@ 64 PCGAMGSetSymGraph - 65 66 Not Collective on PC 67 68 Input Parameters: 69 . pc - the preconditioner context 70 71 Options Database Key: 72 . -pc_gamg_sym_graph 73 74 Level: intermediate 75 76 Concepts: Aggregation AMG preconditioner 77 78 .seealso: () 79 @*/ 80 PetscErrorCode PCGAMGSetSymGraph(PC pc, PetscBool n) 81 { 82 PetscErrorCode ierr; 83 84 PetscFunctionBegin; 85 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 86 ierr = PetscTryMethod(pc,"PCGAMGSetSymGraph_C",(PC,PetscBool),(pc,n));CHKERRQ(ierr); 87 PetscFunctionReturn(0); 88 } 89 90 EXTERN_C_BEGIN 91 #undef __FUNCT__ 92 #define __FUNCT__ "PCGAMGSetSymGraph_GAMG" 93 PetscErrorCode PCGAMGSetSymGraph_GAMG(PC pc, PetscBool n) 94 { 95 PC_MG *mg = (PC_MG*)pc->data; 96 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 97 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 98 99 PetscFunctionBegin; 100 pc_gamg_agg->sym_graph = n; 101 PetscFunctionReturn(0); 102 } 103 EXTERN_C_END 104 105 #undef __FUNCT__ 106 #define __FUNCT__ "PCGAMGSetSquareGraph" 107 /*@ 108 PCGAMGSetSquareGraph - 109 110 Not Collective on PC 111 112 Input Parameters: 113 . pc - the preconditioner context 114 115 Options Database Key: 116 . -pc_gamg_square_graph 117 118 Level: intermediate 119 120 Concepts: Aggregation AMG preconditioner 121 122 .seealso: () 123 @*/ 124 PetscErrorCode PCGAMGSetSquareGraph(PC pc, PetscBool n) 125 { 126 PetscErrorCode ierr; 127 128 PetscFunctionBegin; 129 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 130 ierr = PetscTryMethod(pc,"PCGAMGSetSquareGraph_C",(PC,PetscBool),(pc,n));CHKERRQ(ierr); 131 PetscFunctionReturn(0); 132 } 133 134 EXTERN_C_BEGIN 135 #undef __FUNCT__ 136 #define __FUNCT__ "PCGAMGSetSquareGraph_GAMG" 137 PetscErrorCode PCGAMGSetSquareGraph_GAMG(PC pc, PetscBool n) 138 { 139 PC_MG *mg = (PC_MG*)pc->data; 140 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 141 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 142 143 PetscFunctionBegin; 144 pc_gamg_agg->square_graph = n; 145 PetscFunctionReturn(0); 146 } 147 EXTERN_C_END 148 149 /* -------------------------------------------------------------------------- */ 150 /* 151 PCSetFromOptions_GAMG_AGG 152 153 Input Parameter: 154 . pc - 155 */ 156 #undef __FUNCT__ 157 #define __FUNCT__ "PCSetFromOptions_GAMG_AGG" 158 PetscErrorCode PCSetFromOptions_GAMG_AGG(PC pc) 159 { 160 PetscErrorCode ierr; 161 PC_MG *mg = (PC_MG*)pc->data; 162 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 163 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 164 PetscBool flag; 165 166 PetscFunctionBegin; 167 /* call base class */ 168 ierr = PCSetFromOptions_GAMG(pc);CHKERRQ(ierr); 169 170 ierr = PetscOptionsHead("GAMG-AGG options");CHKERRQ(ierr); 171 { 172 /* -pc_gamg_agg_nsmooths */ 173 pc_gamg_agg->nsmooths = 0; 174 ierr = PetscOptionsInt("-pc_gamg_agg_nsmooths", 175 "smoothing steps for smoothed aggregation, usually 1 (0)", 176 "PCGAMGSetNSmooths", 177 pc_gamg_agg->nsmooths, 178 &pc_gamg_agg->nsmooths, 179 &flag);CHKERRQ(ierr); 180 181 /* -pc_gamg_sym_graph */ 182 pc_gamg_agg->sym_graph = PETSC_FALSE; 183 ierr = PetscOptionsBool("-pc_gamg_sym_graph", 184 "Set for asymmetric matrices", 185 "PCGAMGSetSymGraph", 186 pc_gamg_agg->sym_graph, 187 &pc_gamg_agg->sym_graph, 188 &flag);CHKERRQ(ierr); 189 190 /* -pc_gamg_square_graph */ 191 pc_gamg_agg->square_graph = PETSC_TRUE; 192 ierr = PetscOptionsBool("-pc_gamg_square_graph", 193 "For faster coarsening and lower coarse grid complexity", 194 "PCGAMGSetSquareGraph", 195 pc_gamg_agg->square_graph, 196 &pc_gamg_agg->square_graph, 197 &flag);CHKERRQ(ierr); 198 } 199 ierr = PetscOptionsTail();CHKERRQ(ierr); 200 201 PetscFunctionReturn(0); 202 } 203 204 /* -------------------------------------------------------------------------- */ 205 /* 206 PCDestroy_AGG 207 208 Input Parameter: 209 . pc - 210 */ 211 #undef __FUNCT__ 212 #define __FUNCT__ "PCDestroy_AGG" 213 PetscErrorCode PCDestroy_AGG(PC pc) 214 { 215 PetscErrorCode ierr; 216 PC_MG *mg = (PC_MG*)pc->data; 217 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 218 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 219 220 PetscFunctionBegin; 221 if (pc_gamg_agg) { 222 ierr = PetscFree(pc_gamg_agg);CHKERRQ(ierr); 223 pc_gamg_agg = 0; 224 } 225 226 /* call base class */ 227 ierr = PCDestroy_GAMG(pc);CHKERRQ(ierr); 228 229 PetscFunctionReturn(0); 230 } 231 232 /* -------------------------------------------------------------------------- */ 233 /* 234 PCSetCoordinates_AGG 235 - collective 236 237 Input Parameter: 238 . pc - the preconditioner context 239 . ndm - dimesion of data (used for dof/vertex for Stokes) 240 . a_nloc - number of vertices local 241 . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...} 242 */ 243 EXTERN_C_BEGIN 244 #undef __FUNCT__ 245 #define __FUNCT__ "PCSetCoordinates_AGG" 246 PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords) 247 { 248 PC_MG *mg = (PC_MG*)pc->data; 249 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 250 PetscErrorCode ierr; 251 PetscInt arrsz,kk,ii,jj,nloc,ndatarows,ndf; 252 Mat mat = pc->pmat; 253 /* MPI_Comm wcomm = ((PetscObject)pc)->comm; */ 254 255 PetscFunctionBegin; 256 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 257 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 258 nloc = a_nloc; 259 260 /* SA: null space vectors */ 261 ierr = MatGetBlockSize(mat, &ndf);CHKERRQ(ierr); /* this does not work for Stokes */ 262 if (coords && ndf==1) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */ 263 else if (coords) { 264 if (ndm > ndf) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"degrees of motion %d > block size %d",ndm,ndf); 265 pc_gamg->data_cell_cols = (ndm==2 ? 3 : 6); /* displacement elasticity */ 266 if (ndm != ndf) { 267 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); 268 } 269 } else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */ 270 pc_gamg->data_cell_rows = ndatarows = ndf; assert(pc_gamg->data_cell_cols>0); 271 arrsz = nloc*pc_gamg->data_cell_rows*pc_gamg->data_cell_cols; 272 273 /* create data - syntactic sugar that should be refactored at some point */ 274 if (pc_gamg->data==0 || (pc_gamg->data_sz != arrsz)) { 275 ierr = PetscFree(pc_gamg->data);CHKERRQ(ierr); 276 ierr = PetscMalloc((arrsz+1)*sizeof(PetscReal), &pc_gamg->data);CHKERRQ(ierr); 277 } 278 /* copy data in - column oriented */ 279 for (kk=0;kk<nloc;kk++) { 280 const PetscInt M = nloc*pc_gamg->data_cell_rows; /* stride into data */ 281 PetscReal *data = &pc_gamg->data[kk*ndatarows]; /* start of cell */ 282 if (pc_gamg->data_cell_cols==1) *data = 1.0; 283 else { 284 /* translational modes */ 285 for (ii=0;ii<ndatarows;ii++) 286 for (jj=0;jj<ndatarows;jj++) 287 if (ii==jj)data[ii*M + jj] = 1.0; 288 else data[ii*M + jj] = 0.0; 289 290 /* rotational modes */ 291 if (coords) { 292 if (ndm == 2) { 293 data += 2*M; 294 data[0] = -coords[2*kk+1]; 295 data[1] = coords[2*kk]; 296 } else { 297 data += 3*M; 298 data[0] = 0.0; data[M+0] = coords[3*kk+2]; data[2*M+0] = -coords[3*kk+1]; 299 data[1] = -coords[3*kk+2]; data[M+1] = 0.0; data[2*M+1] = coords[3*kk]; 300 data[2] = coords[3*kk+1]; data[M+2] = -coords[3*kk]; data[2*M+2] = 0.0; 301 } 302 } 303 } 304 } 305 306 pc_gamg->data_sz = arrsz; 307 308 PetscFunctionReturn(0); 309 } 310 EXTERN_C_END 311 312 typedef PetscInt NState; 313 static const NState NOT_DONE=-2; 314 static const NState DELETED=-1; 315 static const NState REMOVED=-3; 316 #define IS_SELECTED(s) (s!=DELETED && s!=NOT_DONE && s!=REMOVED) 317 318 /* -------------------------------------------------------------------------- */ 319 /* 320 smoothAggs - greedy grab of with G1 (unsquared graph) -- AIJ specific 321 - AGG-MG specific: clears singletons out of 'selected_2' 322 323 Input Parameter: 324 . Gmat_2 - glabal matrix of graph (data not defined) 325 . Gmat_1 - base graph to grab with 326 Input/Output Parameter: 327 . aggs_2 - linked list of aggs with gids) 328 */ 329 #undef __FUNCT__ 330 #define __FUNCT__ "smoothAggs" 331 static PetscErrorCode smoothAggs(const Mat Gmat_2, /* base (squared) graph */ 332 const Mat Gmat_1, /* base graph */ 333 /* const IS selected_2, [nselected local] selected vertices */ 334 PetscCoarsenData *aggs_2) /* [nselected local] global ID of aggregate */ 335 { 336 PetscErrorCode ierr; 337 PetscBool isMPI; 338 Mat_SeqAIJ *matA_1, *matB_1=0, *matA_2, *matB_2=0; 339 MPI_Comm wcomm = ((PetscObject)Gmat_2)->comm; 340 PetscMPIInt mype,npe; 341 PetscInt lid,*ii,*idx,ix,Iend,my0,kk,n,j; 342 Mat_MPIAIJ *mpimat_2 = 0, *mpimat_1=0; 343 const PetscInt nloc = Gmat_2->rmap->n; 344 PetscScalar *cpcol_1_state,*cpcol_2_state,*cpcol_2_par_orig,*lid_parent_gid; 345 PetscInt *lid_cprowID_1; 346 NState *lid_state; 347 Vec ghost_par_orig2; 348 349 PetscFunctionBegin; 350 ierr = MPI_Comm_rank(wcomm, &mype);CHKERRQ(ierr); 351 ierr = MPI_Comm_size(wcomm, &npe);CHKERRQ(ierr); 352 ierr = MatGetOwnershipRange(Gmat_1,&my0,&Iend);CHKERRQ(ierr); 353 354 if (PETSC_FALSE) { 355 PetscViewer viewer; char fname[32]; static int llev=0; 356 sprintf(fname,"Gmat2_%d.m",llev++); 357 PetscViewerASCIIOpen(wcomm,fname,&viewer); 358 ierr = PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr); 359 ierr = MatView(Gmat_2, viewer);CHKERRQ(ierr); 360 ierr = PetscViewerDestroy(&viewer); 361 } 362 363 /* get submatrices */ 364 ierr = PetscObjectTypeCompare((PetscObject)Gmat_1, MATMPIAIJ, &isMPI);CHKERRQ(ierr); 365 if (isMPI) { 366 /* grab matrix objects */ 367 mpimat_2 = (Mat_MPIAIJ*)Gmat_2->data; 368 mpimat_1 = (Mat_MPIAIJ*)Gmat_1->data; 369 matA_1 = (Mat_SeqAIJ*)mpimat_1->A->data; 370 matB_1 = (Mat_SeqAIJ*)mpimat_1->B->data; 371 matA_2 = (Mat_SeqAIJ*)mpimat_2->A->data; 372 matB_2 = (Mat_SeqAIJ*)mpimat_2->B->data; 373 374 /* force compressed row storage for B matrix in AuxMat */ 375 matB_1->compressedrow.check = PETSC_TRUE; 376 ierr = MatCheckCompressedRow(mpimat_1->B,&matB_1->compressedrow,matB_1->i,Gmat_1->rmap->n,-1.0);CHKERRQ(ierr); 377 378 ierr = PetscMalloc(nloc*sizeof(PetscInt), &lid_cprowID_1);CHKERRQ(ierr); 379 for (lid = 0 ; lid < nloc ; lid++) lid_cprowID_1[lid] = -1; 380 for (ix=0; ix<matB_1->compressedrow.nrows; ix++) { 381 PetscInt lid = matB_1->compressedrow.rindex[ix]; 382 lid_cprowID_1[lid] = ix; 383 } 384 } else { 385 matA_1 = (Mat_SeqAIJ*)Gmat_1->data; 386 matA_2 = (Mat_SeqAIJ*)Gmat_2->data; 387 lid_cprowID_1 = PETSC_NULL; 388 } 389 assert(matA_1 && !matA_1->compressedrow.use); 390 assert(matB_1==0 || matB_1->compressedrow.use); 391 assert(matA_2 && !matA_2->compressedrow.use); 392 assert(matB_2==0 || matB_2->compressedrow.use); 393 394 /* get state of locals and selected gid for deleted */ 395 ierr = PetscMalloc(nloc*sizeof(NState), &lid_state);CHKERRQ(ierr); 396 ierr = PetscMalloc(nloc*sizeof(PetscScalar), &lid_parent_gid);CHKERRQ(ierr); 397 for (lid = 0 ; lid < nloc ; lid++) { 398 lid_parent_gid[lid] = -1.0; 399 lid_state[lid] = DELETED; 400 } 401 402 /* set lid_state */ 403 for (lid = 0 ; lid < nloc ; lid++) { 404 PetscCDPos pos; 405 ierr = PetscCDGetHeadPos(aggs_2,lid,&pos);CHKERRQ(ierr); 406 if (pos) { 407 PetscInt gid1; 408 ierr = PetscLLNGetID(pos, &gid1);CHKERRQ(ierr); assert(gid1==lid+my0); 409 lid_state[lid] = gid1; 410 } 411 } 412 413 /* map local to selected local, DELETED means a ghost owns it */ 414 for (lid=kk=0;lid<nloc;lid++) { 415 NState state = lid_state[lid]; 416 if (IS_SELECTED(state)) { 417 PetscCDPos pos; 418 ierr = PetscCDGetHeadPos(aggs_2,lid,&pos);CHKERRQ(ierr); 419 while (pos) { 420 PetscInt gid1; 421 ierr = PetscLLNGetID(pos, &gid1);CHKERRQ(ierr); 422 ierr = PetscCDGetNextPos(aggs_2,lid,&pos);CHKERRQ(ierr); 423 424 if (gid1 >= my0 && gid1 < Iend) { 425 lid_parent_gid[gid1-my0] = (PetscScalar)(lid + my0); 426 } 427 } 428 } 429 } 430 /* get 'cpcol_1/2_state' & cpcol_2_par_orig - uses mpimat_1/2->lvec for temp space */ 431 if (isMPI) { 432 Vec tempVec; 433 /* get 'cpcol_1_state' */ 434 ierr = MatGetVecs(Gmat_1, &tempVec, 0);CHKERRQ(ierr); 435 for (kk=0,j=my0;kk<nloc;kk++,j++) { 436 PetscScalar v = (PetscScalar)lid_state[kk]; 437 ierr = VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES);CHKERRQ(ierr); 438 } 439 ierr = VecAssemblyBegin(tempVec);CHKERRQ(ierr); 440 ierr = VecAssemblyEnd(tempVec);CHKERRQ(ierr); 441 ierr = VecScatterBegin(mpimat_1->Mvctx,tempVec, mpimat_1->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 442 ierr = VecScatterEnd(mpimat_1->Mvctx,tempVec, mpimat_1->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 443 ierr = VecGetArray(mpimat_1->lvec, &cpcol_1_state);CHKERRQ(ierr); 444 /* get 'cpcol_2_state' */ 445 ierr = VecScatterBegin(mpimat_2->Mvctx,tempVec, mpimat_2->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 446 ierr = VecScatterEnd(mpimat_2->Mvctx,tempVec, mpimat_2->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 447 ierr = VecGetArray(mpimat_2->lvec, &cpcol_2_state);CHKERRQ(ierr); 448 /* get 'cpcol_2_par_orig' */ 449 for (kk=0,j=my0;kk<nloc;kk++,j++) { 450 PetscScalar v = (PetscScalar)lid_parent_gid[kk]; 451 ierr = VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES);CHKERRQ(ierr); 452 } 453 ierr = VecAssemblyBegin(tempVec);CHKERRQ(ierr); 454 ierr = VecAssemblyEnd(tempVec);CHKERRQ(ierr); 455 ierr = VecDuplicate(mpimat_2->lvec, &ghost_par_orig2);CHKERRQ(ierr); 456 ierr = VecScatterBegin(mpimat_2->Mvctx,tempVec, ghost_par_orig2,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 457 ierr = VecScatterEnd(mpimat_2->Mvctx,tempVec, ghost_par_orig2,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 458 ierr = VecGetArray(ghost_par_orig2, &cpcol_2_par_orig);CHKERRQ(ierr); 459 460 ierr = VecDestroy(&tempVec);CHKERRQ(ierr); 461 } /* ismpi */ 462 463 /* doit */ 464 for (lid=0;lid<nloc;lid++) { 465 NState state = lid_state[lid]; 466 if (IS_SELECTED(state)) { 467 /* steal locals */ 468 ii = matA_1->i; n = ii[lid+1] - ii[lid]; 469 idx = matA_1->j + ii[lid]; 470 for (j=0; j<n; j++) { 471 PetscInt lidj = idx[j], sgid; 472 NState statej = lid_state[lidj]; 473 if (statej==DELETED && (sgid=(PetscInt)PetscRealPart(lid_parent_gid[lidj])) != lid+my0) { /* steal local */ 474 lid_parent_gid[lidj] = (PetscScalar)(lid+my0); /* send this if sgid is not local */ 475 if (sgid >= my0 && sgid < Iend) { /* I'm stealing this local from a local sgid */ 476 PetscInt hav=0,slid=sgid-my0,gidj=lidj+my0; 477 PetscCDPos pos,last=PETSC_NULL; 478 /* looking for local from local so id_llist_2 works */ 479 ierr = PetscCDGetHeadPos(aggs_2,slid,&pos);CHKERRQ(ierr); 480 while (pos) { 481 PetscInt gid; 482 ierr = PetscLLNGetID(pos, &gid);CHKERRQ(ierr); 483 if (gid == gidj) { 484 assert(last); 485 ierr = PetscCDRemoveNextNode(aggs_2, slid, last);CHKERRQ(ierr); 486 ierr = PetscCDAppendNode(aggs_2, lid, pos);CHKERRQ(ierr); 487 hav = 1; 488 break; 489 } else last = pos; 490 491 ierr = PetscCDGetNextPos(aggs_2,slid,&pos);CHKERRQ(ierr); 492 } 493 if (hav!=1) { 494 if (hav==0)SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"failed to find adj in 'selected' lists - structurally unsymmetric matrix"); 495 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"found node %d times???",hav); 496 } 497 } else { /* I'm stealing this local, owned by a ghost */ 498 if (sgid != -1) { 499 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Have un-symmetric graph (apparently). Use '-pc_gamg_sym_graph true' to symetrize the graph or '-pc_gamg_threshold 0.0' if the matrix is structurally symmetric."); 500 } 501 ierr = PetscCDAppendID(aggs_2, lid, lidj+my0);CHKERRQ(ierr); 502 } 503 } 504 } /* local neighbors */ 505 } else if (state == DELETED && lid_cprowID_1) { 506 PetscInt sgidold = (PetscInt)PetscRealPart(lid_parent_gid[lid]); 507 /* see if I have a selected ghost neighbor that will steal me */ 508 if ((ix=lid_cprowID_1[lid]) != -1) { 509 ii = matB_1->compressedrow.i; n = ii[ix+1] - ii[ix]; 510 idx = matB_1->j + ii[ix]; 511 for (j=0 ; j<n ; j++) { 512 PetscInt cpid = idx[j]; 513 NState statej = (NState)PetscRealPart(cpcol_1_state[cpid]); 514 if (IS_SELECTED(statej) && sgidold != (PetscInt)statej) { /* ghost will steal this, remove from my list */ 515 lid_parent_gid[lid] = (PetscScalar)statej; /* send who selected */ 516 if (sgidold>=my0 && sgidold<Iend) { /* this was mine */ 517 PetscInt hav=0,oldslidj=sgidold-my0; 518 PetscCDPos pos,last=PETSC_NULL; 519 /* remove from 'oldslidj' list */ 520 ierr = PetscCDGetHeadPos(aggs_2,oldslidj,&pos);CHKERRQ(ierr); 521 while (pos) { 522 PetscInt gid; 523 ierr = PetscLLNGetID(pos, &gid);CHKERRQ(ierr); 524 if (lid+my0 == gid) { 525 /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove lid from oldslidj list *\/ */ 526 assert(last); 527 ierr = PetscCDRemoveNextNode(aggs_2, oldslidj, last);CHKERRQ(ierr); 528 /* ghost (PetscScalar)statej will add this later */ 529 hav = 1; 530 break; 531 } else last = pos; 532 533 ierr = PetscCDGetNextPos(aggs_2,oldslidj,&pos);CHKERRQ(ierr); 534 } 535 if (hav!=1) { 536 if (hav==0)SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"failed to find adj in 'selected' lists - structurally unsymmetric matrix"); 537 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"found node %d times???",hav); 538 } 539 } else { 540 /* ghosts remove this later */ 541 } 542 } 543 } 544 } 545 } /* selected/deleted */ 546 } /* node loop */ 547 548 if (isMPI) { 549 PetscScalar *cpcol_2_parent,*cpcol_2_gid; 550 Vec tempVec,ghostgids2,ghostparents2; 551 PetscInt cpid,nghost_2; 552 GAMGHashTable gid_cpid; 553 554 ierr = VecGetSize(mpimat_2->lvec, &nghost_2);CHKERRQ(ierr); 555 ierr = MatGetVecs(Gmat_2, &tempVec, 0);CHKERRQ(ierr); 556 557 /* get 'cpcol_2_parent' */ 558 for (kk=0,j=my0;kk<nloc;kk++,j++) { 559 ierr = VecSetValues(tempVec, 1, &j, &lid_parent_gid[kk], INSERT_VALUES);CHKERRQ(ierr); 560 } 561 ierr = VecAssemblyBegin(tempVec);CHKERRQ(ierr); 562 ierr = VecAssemblyEnd(tempVec);CHKERRQ(ierr); 563 ierr = VecDuplicate(mpimat_2->lvec, &ghostparents2);CHKERRQ(ierr); 564 ierr = VecScatterBegin(mpimat_2->Mvctx,tempVec, ghostparents2,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 565 ierr = VecScatterEnd(mpimat_2->Mvctx,tempVec, ghostparents2,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 566 ierr = VecGetArray(ghostparents2, &cpcol_2_parent);CHKERRQ(ierr); 567 568 /* get 'cpcol_2_gid' */ 569 for (kk=0,j=my0;kk<nloc;kk++,j++) { 570 PetscScalar v = (PetscScalar)j; 571 ierr = VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES);CHKERRQ(ierr); 572 } 573 ierr = VecAssemblyBegin(tempVec);CHKERRQ(ierr); 574 ierr = VecAssemblyEnd(tempVec);CHKERRQ(ierr); 575 ierr = VecDuplicate(mpimat_2->lvec, &ghostgids2);CHKERRQ(ierr); 576 ierr = VecScatterBegin(mpimat_2->Mvctx,tempVec, ghostgids2,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 577 ierr = VecScatterEnd(mpimat_2->Mvctx,tempVec, ghostgids2,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 578 ierr = VecGetArray(ghostgids2, &cpcol_2_gid);CHKERRQ(ierr); 579 580 ierr = VecDestroy(&tempVec);CHKERRQ(ierr); 581 582 /* look for deleted ghosts and add to table */ 583 ierr = GAMGTableCreate(2*nghost_2+1, &gid_cpid);CHKERRQ(ierr); 584 for (cpid = 0 ; cpid < nghost_2 ; cpid++) { 585 NState state = (NState)PetscRealPart(cpcol_2_state[cpid]); 586 if (state==DELETED) { 587 PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]); 588 PetscInt sgid_old = (PetscInt)PetscRealPart(cpcol_2_par_orig[cpid]); 589 if (sgid_old == -1 && sgid_new != -1) { 590 PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]); 591 ierr = GAMGTableAdd(&gid_cpid, gid, cpid);CHKERRQ(ierr); 592 } 593 } 594 } 595 596 /* look for deleted ghosts and see if they moved - remove it */ 597 for (lid=0;lid<nloc;lid++) { 598 NState state = lid_state[lid]; 599 if (IS_SELECTED(state)) { 600 PetscCDPos pos,last=PETSC_NULL; 601 /* look for deleted ghosts and see if they moved */ 602 ierr = PetscCDGetHeadPos(aggs_2,lid,&pos);CHKERRQ(ierr); 603 while (pos) { 604 PetscInt gid; 605 ierr = PetscLLNGetID(pos, &gid);CHKERRQ(ierr); 606 607 if (gid < my0 || gid >= Iend) { 608 ierr = GAMGTableFind(&gid_cpid, gid, &cpid);CHKERRQ(ierr); 609 if (cpid != -1) { 610 /* a moved ghost - */ 611 /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove 'flid' from list *\/ */ 612 ierr = PetscCDRemoveNextNode(aggs_2, lid, last);CHKERRQ(ierr); 613 } else last = pos; 614 } else last = pos; 615 616 ierr = PetscCDGetNextPos(aggs_2,lid,&pos);CHKERRQ(ierr); 617 } /* loop over list of deleted */ 618 } /* selected */ 619 } 620 ierr = GAMGTableDestroy(&gid_cpid);CHKERRQ(ierr); 621 622 /* look at ghosts, see if they changed - and it */ 623 for (cpid = 0 ; cpid < nghost_2 ; cpid++) { 624 PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]); 625 if (sgid_new >= my0 && sgid_new < Iend) { /* this is mine */ 626 PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]); 627 PetscInt slid_new=sgid_new-my0,hav=0; 628 PetscCDPos pos; 629 /* search for this gid to see if I have it */ 630 ierr = PetscCDGetHeadPos(aggs_2,slid_new,&pos);CHKERRQ(ierr); 631 while (pos) { 632 PetscInt gidj; 633 ierr = PetscLLNGetID(pos, &gidj);CHKERRQ(ierr); 634 ierr = PetscCDGetNextPos(aggs_2,slid_new,&pos);CHKERRQ(ierr); 635 636 if (gidj == gid) { hav = 1; break; } 637 } 638 if (hav != 1) { 639 /* insert 'flidj' into head of llist */ 640 ierr = PetscCDAppendID(aggs_2, slid_new, gid);CHKERRQ(ierr); 641 } 642 } 643 } 644 645 ierr = VecRestoreArray(mpimat_1->lvec, &cpcol_1_state);CHKERRQ(ierr); 646 ierr = VecRestoreArray(mpimat_2->lvec, &cpcol_2_state);CHKERRQ(ierr); 647 ierr = VecRestoreArray(ghostparents2, &cpcol_2_parent);CHKERRQ(ierr); 648 ierr = VecRestoreArray(ghostgids2, &cpcol_2_gid);CHKERRQ(ierr); 649 ierr = PetscFree(lid_cprowID_1);CHKERRQ(ierr); 650 ierr = VecDestroy(&ghostgids2);CHKERRQ(ierr); 651 ierr = VecDestroy(&ghostparents2);CHKERRQ(ierr); 652 ierr = VecDestroy(&ghost_par_orig2);CHKERRQ(ierr); 653 } 654 655 ierr = PetscFree(lid_parent_gid);CHKERRQ(ierr); 656 ierr = PetscFree(lid_state);CHKERRQ(ierr); 657 658 PetscFunctionReturn(0); 659 } 660 661 /* -------------------------------------------------------------------------- */ 662 /* 663 PCSetData_AGG - called if data is not set with PCSetCoordinates. 664 Looks in Mat for near null space. 665 Does not work for Stokes 666 667 Input Parameter: 668 . pc - 669 . a_A - matrix to get (near) null space out of. 670 */ 671 #undef __FUNCT__ 672 #define __FUNCT__ "PCSetData_AGG" 673 PetscErrorCode PCSetData_AGG(PC pc, Mat a_A) 674 { 675 PetscErrorCode ierr; 676 PC_MG *mg = (PC_MG*)pc->data; 677 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 678 MatNullSpace mnull; 679 680 PetscFunctionBegin; 681 ierr = MatGetNearNullSpace(a_A, &mnull);CHKERRQ(ierr); 682 if (!mnull) { 683 PetscInt bs,NN,MM; 684 ierr = MatGetBlockSize(a_A, &bs);CHKERRQ(ierr); 685 ierr = MatGetLocalSize(a_A, &MM, &NN);CHKERRQ(ierr); assert(MM%bs==0); 686 ierr = PCSetCoordinates_AGG(pc, bs, MM/bs, PETSC_NULL);CHKERRQ(ierr); 687 } else { 688 PetscReal *nullvec; 689 PetscBool has_const; 690 PetscInt i,j,mlocal,nvec,bs; 691 const Vec *vecs; const PetscScalar *v; 692 ierr = MatGetLocalSize(a_A,&mlocal,PETSC_NULL);CHKERRQ(ierr); 693 ierr = MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);CHKERRQ(ierr); 694 ierr = PetscMalloc((nvec+!!has_const)*mlocal*sizeof(*nullvec),&nullvec);CHKERRQ(ierr); 695 if (has_const) for (i=0; i<mlocal; i++) nullvec[i] = 1.0; 696 for (i=0; i<nvec; i++) { 697 ierr = VecGetArrayRead(vecs[i],&v);CHKERRQ(ierr); 698 for (j=0; j<mlocal; j++) nullvec[(i+!!has_const)*mlocal + j] = PetscRealPart(v[j]); 699 ierr = VecRestoreArrayRead(vecs[i],&v);CHKERRQ(ierr); 700 } 701 pc_gamg->data = nullvec; 702 pc_gamg->data_cell_cols = (nvec+!!has_const); 703 ierr = MatGetBlockSize(a_A, &bs);CHKERRQ(ierr); /* this does not work for Stokes */ 704 pc_gamg->data_cell_rows = bs; 705 } 706 PetscFunctionReturn(0); 707 } 708 709 /* -------------------------------------------------------------------------- */ 710 /* 711 formProl0 712 713 Input Parameter: 714 . agg_llists - list of arrays with aggregates 715 . bs - block size 716 . nSAvec - column bs of new P 717 . my0crs - global index of start of locals 718 . data_stride - bs*(nloc nodes + ghost nodes) 719 . data_in[data_stride*nSAvec] - local data on fine grid 720 . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist' 721 Output Parameter: 722 . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data 723 . a_Prol - prolongation operator 724 */ 725 #undef __FUNCT__ 726 #define __FUNCT__ "formProl0" 727 static PetscErrorCode formProl0(const PetscCoarsenData *agg_llists,/* list from selected vertices of aggregate unselected vertices */ 728 const PetscInt bs, /* (row) block size */ 729 const PetscInt nSAvec, /* column bs */ 730 const PetscInt my0crs, /* global index of start of locals */ 731 const PetscInt data_stride, /* (nloc+nghost)*bs */ 732 PetscReal data_in[], /* [data_stride][nSAvec] */ 733 const PetscInt flid_fgid[], /* [data_stride/bs] */ 734 PetscReal **a_data_out, 735 Mat a_Prol) /* prolongation operator (output)*/ 736 { 737 PetscErrorCode ierr; 738 PetscInt Istart,my0,Iend,nloc,clid,flid,aggID,kk,jj,ii,mm,ndone,nSelected,minsz,nghosts,out_data_stride; 739 MPI_Comm wcomm = ((PetscObject)a_Prol)->comm; 740 PetscMPIInt mype, npe; 741 PetscReal *out_data; 742 PetscCDPos pos; 743 GAMGHashTable fgid_flid; 744 745 /* #define OUT_AGGS */ 746 #if defined(OUT_AGGS) 747 static PetscInt llev = 0; char fname[32]; FILE *file = PETSC_NULL; PetscInt pM; 748 #endif 749 750 PetscFunctionBegin; 751 ierr = MPI_Comm_rank(wcomm,&mype);CHKERRQ(ierr); 752 ierr = MPI_Comm_size(wcomm,&npe);CHKERRQ(ierr); 753 ierr = MatGetOwnershipRange(a_Prol, &Istart, &Iend);CHKERRQ(ierr); 754 nloc = (Iend-Istart)/bs; my0 = Istart/bs; assert((Iend-Istart)%bs==0); 755 Iend /= bs; 756 nghosts = data_stride/bs - nloc; 757 758 ierr = GAMGTableCreate(2*nghosts+1, &fgid_flid);CHKERRQ(ierr); 759 for (kk=0;kk<nghosts;kk++) { 760 ierr = GAMGTableAdd(&fgid_flid, flid_fgid[nloc+kk], nloc+kk);CHKERRQ(ierr); 761 } 762 763 #if defined(OUT_AGGS) 764 sprintf(fname,"aggs_%d_%d.m",llev++,mype); 765 if (llev==1) { 766 file = fopen(fname,"w"); 767 } 768 MatGetSize(a_Prol, &pM, &jj); 769 #endif 770 771 /* count selected -- same as number of cols of P */ 772 for (nSelected=mm=0;mm<nloc;mm++) { 773 PetscBool ise; 774 ierr = PetscCDEmptyAt(agg_llists, mm, &ise);CHKERRQ(ierr); 775 if (!ise) nSelected++; 776 } 777 ierr = MatGetOwnershipRangeColumn(a_Prol, &ii, &jj);CHKERRQ(ierr); 778 assert((ii/nSAvec)==my0crs); assert(nSelected==(jj-ii)/nSAvec); 779 780 /* aloc space for coarse point data (output) */ 781 out_data_stride = nSelected*nSAvec; 782 ierr = PetscMalloc(out_data_stride*nSAvec*sizeof(PetscReal), &out_data);CHKERRQ(ierr); 783 for (ii=0;ii<out_data_stride*nSAvec;ii++) { 784 out_data[ii]=1.e300; 785 } 786 *a_data_out = out_data; /* output - stride nSelected*nSAvec */ 787 788 /* find points and set prolongation */ 789 minsz = 100; 790 ndone = 0; 791 for (mm = clid = 0 ; mm < nloc ; mm++) { 792 ierr = PetscCDSizeAt(agg_llists, mm, &jj);CHKERRQ(ierr); 793 if (jj > 0) { 794 const PetscInt lid = mm, cgid = my0crs + clid; 795 PetscInt cids[100]; /* max bs */ 796 PetscBLASInt asz=jj,M=asz*bs,N=nSAvec,INFO; 797 PetscBLASInt Mdata=M+((N-M>0)?N-M:0),LDA=Mdata,LWORK=N*bs; 798 PetscScalar *qqc,*qqr,*TAU,*WORK; 799 PetscInt *fids; 800 PetscReal *data; 801 /* count agg */ 802 if (asz<minsz) minsz = asz; 803 804 /* get block */ 805 ierr = PetscMalloc((Mdata*N)*sizeof(PetscScalar), &qqc);CHKERRQ(ierr); 806 ierr = PetscMalloc((M*N)*sizeof(PetscScalar), &qqr);CHKERRQ(ierr); 807 ierr = PetscMalloc(N*sizeof(PetscScalar), &TAU);CHKERRQ(ierr); 808 ierr = PetscMalloc(LWORK*sizeof(PetscScalar), &WORK);CHKERRQ(ierr); 809 ierr = PetscMalloc(M*sizeof(PetscInt), &fids);CHKERRQ(ierr); 810 811 aggID = 0; 812 ierr = PetscCDGetHeadPos(agg_llists,lid,&pos);CHKERRQ(ierr); 813 while (pos) { 814 PetscInt gid1; 815 ierr = PetscLLNGetID(pos, &gid1);CHKERRQ(ierr); 816 ierr = PetscCDGetNextPos(agg_llists,lid,&pos);CHKERRQ(ierr); 817 818 if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0; 819 else { 820 ierr = GAMGTableFind(&fgid_flid, gid1, &flid);CHKERRQ(ierr); 821 assert(flid>=0); 822 } 823 /* copy in B_i matrix - column oriented */ 824 data = &data_in[flid*bs]; 825 for (kk = ii = 0; ii < bs ; ii++) { 826 for (jj = 0; jj < N ; jj++) { 827 PetscReal d = data[jj*data_stride + ii]; 828 qqc[jj*Mdata + aggID*bs + ii] = d; 829 } 830 } 831 #if defined(OUT_AGGS) 832 if (llev==1) { 833 char str[] = "plot(%e,%e,'r*'), hold on,\n", col[] = "rgbkmc", sim[] = "*os+h>d<vx^"; 834 PetscInt MM,pi,pj; 835 str[12] = col[(clid+17*mype)%6]; str[13] = sim[(clid+17*mype)%11]; 836 MM = (PetscInt)(PetscSqrtReal((PetscReal)pM)); 837 pj = gid1/MM; pi = gid1%MM; 838 fprintf(file,str,(double)pi,(double)pj); 839 /* fprintf(file,str,data[2*data_stride+1],-data[2*data_stride]); */ 840 } 841 #endif 842 /* set fine IDs */ 843 for (kk=0;kk<bs;kk++) fids[aggID*bs + kk] = flid_fgid[flid]*bs + kk; 844 845 aggID++; 846 } 847 848 /* pad with zeros */ 849 for (ii = asz*bs; ii < Mdata ; ii++) { 850 for (jj = 0; jj < N ; jj++, kk++) { 851 qqc[jj*Mdata + ii] = .0; 852 } 853 } 854 855 ndone += aggID; 856 /* QR */ 857 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 858 LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO); 859 ierr = PetscFPTrapPop();CHKERRQ(ierr); 860 if (INFO != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xGEQRF error"); 861 /* get R - column oriented - output B_{i+1} */ 862 { 863 PetscReal *data = &out_data[clid*nSAvec]; 864 for (jj = 0; jj < nSAvec ; jj++) { 865 for (ii = 0; ii < nSAvec ; ii++) { 866 assert(data[jj*out_data_stride + ii] == 1.e300); 867 if (ii <= jj) data[jj*out_data_stride + ii] = PetscRealPart(qqc[jj*Mdata + ii]); 868 else data[jj*out_data_stride + ii] = 0.; 869 } 870 } 871 } 872 873 /* get Q - row oriented */ 874 LAPACKungqr_(&Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO); 875 if (INFO != 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xORGQR error arg %d",-INFO); 876 877 for (ii = 0 ; ii < M ; ii++) { 878 for (jj = 0 ; jj < N ; jj++) { 879 qqr[N*ii + jj] = qqc[jj*Mdata + ii]; 880 } 881 } 882 883 /* add diagonal block of P0 */ 884 for (kk=0;kk<N;kk++) { 885 cids[kk] = N*cgid + kk; /* global col IDs in P0 */ 886 } 887 ierr = MatSetValues(a_Prol,M,fids,N,cids,qqr,INSERT_VALUES);CHKERRQ(ierr); 888 889 ierr = PetscFree(qqc);CHKERRQ(ierr); 890 ierr = PetscFree(qqr);CHKERRQ(ierr); 891 ierr = PetscFree(TAU);CHKERRQ(ierr); 892 ierr = PetscFree(WORK);CHKERRQ(ierr); 893 ierr = PetscFree(fids);CHKERRQ(ierr); 894 clid++; 895 } /* coarse agg */ 896 } /* for all fine nodes */ 897 ierr = MatAssemblyBegin(a_Prol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 898 ierr = MatAssemblyEnd(a_Prol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 899 900 /* ierr = MPI_Allreduce(&ndone, &ii, 1, MPIU_INT, MPIU_SUM, wcomm); */ 901 /* MatGetSize(a_Prol, &kk, &jj); */ 902 /* ierr = MPI_Allreduce(&minsz, &jj, 1, MPIU_INT, MPIU_MIN, wcomm); */ 903 /* PetscPrintf(wcomm," **** [%d]%s %d total done, %d nodes (%d local done), min agg. size = %d\n",mype,__FUNCT__,ii,kk/bs,ndone,jj); */ 904 905 #if defined(OUT_AGGS) 906 if (llev==1) fclose(file); 907 #endif 908 ierr = GAMGTableDestroy(&fgid_flid);CHKERRQ(ierr); 909 910 PetscFunctionReturn(0); 911 } 912 913 /* -------------------------------------------------------------------------- */ 914 /* 915 PCGAMGgraph_AGG 916 917 Input Parameter: 918 . pc - this 919 . Amat - matrix on this fine level 920 Output Parameter: 921 . a_Gmat - 922 */ 923 #undef __FUNCT__ 924 #define __FUNCT__ "PCGAMGgraph_AGG" 925 PetscErrorCode PCGAMGgraph_AGG(PC pc, 926 const Mat Amat, 927 Mat *a_Gmat) 928 { 929 PetscErrorCode ierr; 930 PC_MG *mg = (PC_MG*)pc->data; 931 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 932 const PetscInt verbose = pc_gamg->verbose; 933 const PetscReal vfilter = pc_gamg->threshold; 934 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 935 PetscMPIInt mype,npe; 936 Mat Gmat; 937 MPI_Comm wcomm = ((PetscObject)Amat)->comm; 938 PetscBool /* set,flg , */symm; 939 940 PetscFunctionBegin; 941 #if defined PETSC_USE_LOG 942 ierr = PetscLogEventBegin(PC_GAMGGgraph_AGG,0,0,0,0);CHKERRQ(ierr); 943 #endif 944 ierr = MPI_Comm_rank(wcomm, &mype);CHKERRQ(ierr); 945 ierr = MPI_Comm_size(wcomm, &npe);CHKERRQ(ierr); 946 947 /* ierr = MatIsSymmetricKnown(Amat, &set, &flg);CHKERRQ(ierr); || !(set && flg) -- this causes lot of symm calls */ 948 symm = (PetscBool)(pc_gamg_agg->sym_graph); /* && !pc_gamg_agg->square_graph; */ 949 950 ierr = PCGAMGCreateGraph(Amat, &Gmat);CHKERRQ(ierr); 951 ierr = PCGAMGFilterGraph(&Gmat, vfilter, symm, verbose);CHKERRQ(ierr); 952 953 *a_Gmat = Gmat; 954 955 #if defined PETSC_USE_LOG 956 ierr = PetscLogEventEnd(PC_GAMGGgraph_AGG,0,0,0,0);CHKERRQ(ierr); 957 #endif 958 PetscFunctionReturn(0); 959 } 960 961 /* -------------------------------------------------------------------------- */ 962 /* 963 PCGAMGCoarsen_AGG 964 965 Input Parameter: 966 . a_pc - this 967 Input/Output Parameter: 968 . a_Gmat1 - graph on this fine level - coarsening can change this (squares it) 969 Output Parameter: 970 . agg_lists - list of aggregates 971 */ 972 #undef __FUNCT__ 973 #define __FUNCT__ "PCGAMGCoarsen_AGG" 974 PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc, 975 Mat *a_Gmat1, 976 PetscCoarsenData **agg_lists) 977 { 978 PetscErrorCode ierr; 979 PC_MG *mg = (PC_MG*)a_pc->data; 980 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 981 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 982 Mat mat,Gmat2, Gmat1 = *a_Gmat1; /* squared graph */ 983 const PetscInt verbose = pc_gamg->verbose; 984 IS perm; 985 PetscInt Ii,nloc,bs,n,m; 986 PetscInt *permute; 987 PetscBool *bIndexSet; 988 MatCoarsen crs; 989 MPI_Comm wcomm = ((PetscObject)Gmat1)->comm; 990 PetscMPIInt mype,npe; 991 992 PetscFunctionBegin; 993 #if defined PETSC_USE_LOG 994 ierr = PetscLogEventBegin(PC_GAMGCoarsen_AGG,0,0,0,0);CHKERRQ(ierr); 995 #endif 996 ierr = MPI_Comm_rank(wcomm, &mype);CHKERRQ(ierr); 997 ierr = MPI_Comm_size(wcomm, &npe);CHKERRQ(ierr); 998 ierr = MatGetLocalSize(Gmat1, &n, &m);CHKERRQ(ierr); 999 ierr = MatGetBlockSize(Gmat1, &bs);CHKERRQ(ierr); assert(bs==1); 1000 nloc = n/bs; 1001 1002 if (pc_gamg_agg->square_graph) { 1003 if (verbose > 1) { 1004 PetscPrintf(wcomm,"[%d]%s square graph\n",mype,__FUNCT__); 1005 } 1006 /* ierr = MatMatTransposeMult(Gmat1, Gmat1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat2); */ 1007 ierr = MatTransposeMatMult(Gmat1, Gmat1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat2);CHKERRQ(ierr); 1008 if (verbose > 2) { 1009 ierr = PetscPrintf(wcomm,"[%d]%s square graph done\n",mype,__FUNCT__);CHKERRQ(ierr); 1010 /* check for symetry */ 1011 if (verbose > 4) { 1012 1013 } 1014 } 1015 } else Gmat2 = Gmat1; 1016 1017 /* get MIS aggs */ 1018 /* randomize */ 1019 ierr = PetscMalloc(nloc*sizeof(PetscInt), &permute);CHKERRQ(ierr); 1020 ierr = PetscMalloc(nloc*sizeof(PetscBool), &bIndexSet);CHKERRQ(ierr); 1021 for (Ii = 0; Ii < nloc ; Ii++) { 1022 bIndexSet[Ii] = PETSC_FALSE; 1023 permute[Ii] = Ii; 1024 } 1025 srand(1); /* make deterministic */ 1026 for (Ii = 0; Ii < nloc ; Ii++) { 1027 PetscInt iSwapIndex = rand()%nloc; 1028 if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) { 1029 PetscInt iTemp = permute[iSwapIndex]; 1030 permute[iSwapIndex] = permute[Ii]; 1031 permute[Ii] = iTemp; 1032 bIndexSet[iSwapIndex] = PETSC_TRUE; 1033 } 1034 } 1035 ierr = PetscFree(bIndexSet);CHKERRQ(ierr); 1036 1037 if (verbose > 1) { 1038 PetscPrintf(wcomm,"[%d]%s coarsen graph\n",mype,__FUNCT__); 1039 } 1040 1041 ierr = ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm);CHKERRQ(ierr); 1042 #if defined PETSC_GAMG_USE_LOG 1043 ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET4],0,0,0,0);CHKERRQ(ierr); 1044 #endif 1045 ierr = MatCoarsenCreate(wcomm, &crs);CHKERRQ(ierr); 1046 /* ierr = MatCoarsenSetType(crs, MATCOARSENMIS);CHKERRQ(ierr); */ 1047 ierr = MatCoarsenSetFromOptions(crs);CHKERRQ(ierr); 1048 ierr = MatCoarsenSetGreedyOrdering(crs, perm);CHKERRQ(ierr); 1049 ierr = MatCoarsenSetAdjacency(crs, Gmat2);CHKERRQ(ierr); 1050 ierr = MatCoarsenSetVerbose(crs, pc_gamg->verbose);CHKERRQ(ierr); 1051 ierr = MatCoarsenSetStrictAggs(crs, PETSC_TRUE);CHKERRQ(ierr); 1052 ierr = MatCoarsenApply(crs);CHKERRQ(ierr); 1053 ierr = MatCoarsenGetData(crs, agg_lists);CHKERRQ(ierr); /* output */ 1054 ierr = MatCoarsenDestroy(&crs);CHKERRQ(ierr); 1055 1056 ierr = ISDestroy(&perm);CHKERRQ(ierr); 1057 ierr = PetscFree(permute);CHKERRQ(ierr); 1058 #if defined PETSC_GAMG_USE_LOG 1059 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET4],0,0,0,0);CHKERRQ(ierr); 1060 #endif 1061 if (verbose > 2) { 1062 PetscPrintf(wcomm,"[%d]%s coarsen graph done\n",mype,__FUNCT__); 1063 } 1064 1065 /* smooth aggs */ 1066 if (Gmat2 != Gmat1) { 1067 const PetscCoarsenData *llist = *agg_lists; 1068 ierr = smoothAggs(Gmat2, Gmat1, *agg_lists);CHKERRQ(ierr); 1069 ierr = MatDestroy(&Gmat1);CHKERRQ(ierr); 1070 *a_Gmat1 = Gmat2; /* output */ 1071 ierr = PetscCDGetMat(llist, &mat);CHKERRQ(ierr); 1072 if (mat) SETERRQ(wcomm,PETSC_ERR_ARG_WRONG, "Auxilary matrix with squared graph????"); 1073 } else { 1074 const PetscCoarsenData *llist = *agg_lists; 1075 /* see if we have a matrix that takes precedence (returned from MatCoarsenAppply) */ 1076 ierr = PetscCDGetMat(llist, &mat);CHKERRQ(ierr); 1077 if (mat) { 1078 ierr = MatDestroy(&Gmat1);CHKERRQ(ierr); 1079 *a_Gmat1 = mat; /* output */ 1080 } 1081 } 1082 #if defined PETSC_USE_LOG 1083 ierr = PetscLogEventEnd(PC_GAMGCoarsen_AGG,0,0,0,0);CHKERRQ(ierr); 1084 #endif 1085 if (verbose > 2) { 1086 PetscPrintf(wcomm,"[%d]%s PCGAMGCoarsen_AGG done\n",mype,__FUNCT__); 1087 } 1088 1089 PetscFunctionReturn(0); 1090 } 1091 1092 /* -------------------------------------------------------------------------- */ 1093 /* 1094 PCGAMGProlongator_AGG 1095 1096 Input Parameter: 1097 . pc - this 1098 . Amat - matrix on this fine level 1099 . Graph - used to get ghost data for nodes in 1100 . agg_lists - list of aggregates 1101 Output Parameter: 1102 . a_P_out - prolongation operator to the next level 1103 */ 1104 #undef __FUNCT__ 1105 #define __FUNCT__ "PCGAMGProlongator_AGG" 1106 PetscErrorCode PCGAMGProlongator_AGG(PC pc, 1107 const Mat Amat, 1108 const Mat Gmat, 1109 PetscCoarsenData *agg_lists, 1110 Mat *a_P_out) 1111 { 1112 PC_MG *mg = (PC_MG*)pc->data; 1113 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1114 const PetscInt verbose = pc_gamg->verbose; 1115 const PetscInt data_cols = pc_gamg->data_cell_cols; 1116 PetscErrorCode ierr; 1117 PetscInt Istart,Iend,nloc,ii,jj,kk,my0,nLocalSelected,bs; 1118 Mat Prol; 1119 PetscMPIInt mype, npe; 1120 MPI_Comm wcomm = ((PetscObject)Amat)->comm; 1121 const PetscInt col_bs = data_cols; 1122 PetscReal *data_w_ghost; 1123 PetscInt myCrs0, nbnodes=0, *flid_fgid; 1124 1125 PetscFunctionBegin; 1126 #if defined PETSC_USE_LOG 1127 ierr = PetscLogEventBegin(PC_GAMGProlongator_AGG,0,0,0,0);CHKERRQ(ierr); 1128 #endif 1129 ierr = MPI_Comm_rank(wcomm, &mype);CHKERRQ(ierr); 1130 ierr = MPI_Comm_size(wcomm, &npe);CHKERRQ(ierr); 1131 ierr = MatGetOwnershipRange(Amat, &Istart, &Iend);CHKERRQ(ierr); 1132 ierr = MatGetBlockSize(Amat, &bs);CHKERRQ(ierr); 1133 nloc = (Iend-Istart)/bs; my0 = Istart/bs; assert((Iend-Istart)%bs==0); 1134 1135 /* get 'nLocalSelected' */ 1136 for (ii=0, nLocalSelected = 0 ; ii < nloc ; ii++) { 1137 PetscBool ise; 1138 /* filter out singletons 0 or 1? */ 1139 ierr = PetscCDEmptyAt(agg_lists, ii, &ise);CHKERRQ(ierr); 1140 if (!ise) nLocalSelected++; 1141 } 1142 1143 /* create prolongator, create P matrix */ 1144 ierr = MatCreate(wcomm, &Prol);CHKERRQ(ierr); 1145 ierr = MatSetSizes(Prol,nloc*bs,nLocalSelected*col_bs,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1146 ierr = MatSetBlockSizes(Prol, bs, col_bs);CHKERRQ(ierr); 1147 ierr = MatSetType(Prol, MATAIJ);CHKERRQ(ierr); 1148 ierr = MatSeqAIJSetPreallocation(Prol, data_cols, PETSC_NULL);CHKERRQ(ierr); 1149 ierr = MatMPIAIJSetPreallocation(Prol,data_cols, PETSC_NULL,data_cols, PETSC_NULL);CHKERRQ(ierr); 1150 /* nloc*bs, nLocalSelected*col_bs, */ 1151 /* PETSC_DETERMINE, PETSC_DETERMINE, */ 1152 /* data_cols, PETSC_NULL, data_cols, PETSC_NULL, */ 1153 /* &Prol); */ 1154 1155 /* can get all points "removed" */ 1156 ierr = MatGetSize(Prol, &kk, &ii);CHKERRQ(ierr); 1157 if (ii==0) { 1158 if (verbose > 0) { 1159 PetscPrintf(wcomm,"[%d]%s no selected points on coarse grid\n",mype,__FUNCT__); 1160 } 1161 ierr = MatDestroy(&Prol);CHKERRQ(ierr); 1162 *a_P_out = PETSC_NULL; /* out */ 1163 PetscFunctionReturn(0); 1164 } 1165 if (verbose > 0) { 1166 PetscPrintf(wcomm,"\t\t[%d]%s New grid %d nodes\n",mype,__FUNCT__,ii/col_bs); 1167 } 1168 ierr = MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk);CHKERRQ(ierr); 1169 1170 assert((kk-myCrs0)%col_bs==0); 1171 myCrs0 = myCrs0/col_bs; 1172 assert((kk/col_bs-myCrs0)==nLocalSelected); 1173 1174 /* create global vector of data in 'data_w_ghost' */ 1175 #if defined PETSC_GAMG_USE_LOG 1176 ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET7],0,0,0,0);CHKERRQ(ierr); 1177 #endif 1178 if (npe > 1) { /* */ 1179 PetscReal *tmp_gdata,*tmp_ldata,*tp2; 1180 ierr = PetscMalloc(nloc*sizeof(PetscReal), &tmp_ldata);CHKERRQ(ierr); 1181 for (jj = 0 ; jj < data_cols ; jj++) { 1182 for (kk = 0 ; kk < bs ; kk++) { 1183 PetscInt ii,stride; 1184 const PetscReal *tp = pc_gamg->data + jj*bs*nloc + kk; 1185 for (ii = 0 ; ii < nloc ; ii++, tp += bs) { 1186 tmp_ldata[ii] = *tp; 1187 } 1188 ierr = PCGAMGGetDataWithGhosts(Gmat, 1, tmp_ldata, &stride, &tmp_gdata);CHKERRQ(ierr); 1189 1190 if (jj==0 && kk==0) { /* now I know how many todal nodes - allocate */ 1191 ierr = PetscMalloc(stride*bs*data_cols*sizeof(PetscReal), &data_w_ghost);CHKERRQ(ierr); 1192 nbnodes = bs*stride; 1193 } 1194 tp2 = data_w_ghost + jj*bs*stride + kk; 1195 for (ii = 0 ; ii < stride ; ii++, tp2 += bs) { 1196 *tp2 = tmp_gdata[ii]; 1197 } 1198 ierr = PetscFree(tmp_gdata);CHKERRQ(ierr); 1199 } 1200 } 1201 ierr = PetscFree(tmp_ldata);CHKERRQ(ierr); 1202 } else { 1203 nbnodes = bs*nloc; 1204 data_w_ghost = (PetscReal*)pc_gamg->data; 1205 } 1206 1207 /* get P0 */ 1208 if (npe > 1) { 1209 PetscReal *fid_glid_loc,*fiddata; 1210 PetscInt stride; 1211 1212 ierr = PetscMalloc(nloc*sizeof(PetscReal), &fid_glid_loc);CHKERRQ(ierr); 1213 for (kk=0;kk<nloc;kk++) fid_glid_loc[kk] = (PetscReal)(my0+kk); 1214 ierr = PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata);CHKERRQ(ierr); 1215 ierr = PetscMalloc(stride*sizeof(PetscInt), &flid_fgid);CHKERRQ(ierr); 1216 for (kk=0;kk<stride;kk++) flid_fgid[kk] = (PetscInt)fiddata[kk]; 1217 ierr = PetscFree(fiddata);CHKERRQ(ierr); 1218 1219 assert(stride==nbnodes/bs); 1220 ierr = PetscFree(fid_glid_loc);CHKERRQ(ierr); 1221 } else { 1222 ierr = PetscMalloc(nloc*sizeof(PetscInt), &flid_fgid);CHKERRQ(ierr); 1223 for (kk=0;kk<nloc;kk++) flid_fgid[kk] = my0 + kk; 1224 } 1225 #if defined PETSC_GAMG_USE_LOG 1226 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET7],0,0,0,0);CHKERRQ(ierr); 1227 ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET8],0,0,0,0);CHKERRQ(ierr); 1228 #endif 1229 { 1230 PetscReal *data_out = PETSC_NULL; 1231 ierr = formProl0(agg_lists, bs, data_cols, myCrs0, nbnodes,data_w_ghost, flid_fgid, &data_out, Prol);CHKERRQ(ierr); 1232 ierr = PetscFree(pc_gamg->data);CHKERRQ(ierr); 1233 pc_gamg->data = data_out; 1234 pc_gamg->data_cell_rows = data_cols; 1235 pc_gamg->data_sz = data_cols*data_cols*nLocalSelected; 1236 } 1237 #if defined PETSC_GAMG_USE_LOG 1238 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET8],0,0,0,0);CHKERRQ(ierr); 1239 #endif 1240 if (npe > 1) ierr = PetscFree(data_w_ghost);CHKERRQ(ierr); 1241 ierr = PetscFree(flid_fgid);CHKERRQ(ierr); 1242 1243 *a_P_out = Prol; /* out */ 1244 #if defined PETSC_USE_LOG 1245 ierr = PetscLogEventEnd(PC_GAMGProlongator_AGG,0,0,0,0);CHKERRQ(ierr); 1246 #endif 1247 PetscFunctionReturn(0); 1248 } 1249 1250 /* -------------------------------------------------------------------------- */ 1251 /* 1252 PCGAMGOptprol_AGG 1253 1254 Input Parameter: 1255 . pc - this 1256 . Amat - matrix on this fine level 1257 In/Output Parameter: 1258 . a_P_out - prolongation operator to the next level 1259 */ 1260 #undef __FUNCT__ 1261 #define __FUNCT__ "PCGAMGOptprol_AGG" 1262 PetscErrorCode PCGAMGOptprol_AGG(PC pc, 1263 const Mat Amat, 1264 Mat *a_P) 1265 { 1266 PetscErrorCode ierr; 1267 PC_MG *mg = (PC_MG*)pc->data; 1268 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1269 const PetscInt verbose = pc_gamg->verbose; 1270 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 1271 PetscInt jj; 1272 PetscMPIInt mype,npe; 1273 Mat Prol = *a_P; 1274 MPI_Comm wcomm = ((PetscObject)Amat)->comm; 1275 1276 PetscFunctionBegin; 1277 #if defined PETSC_USE_LOG 1278 ierr = PetscLogEventBegin(PC_GAMGOptprol_AGG,0,0,0,0);CHKERRQ(ierr); 1279 #endif 1280 ierr = MPI_Comm_rank(wcomm, &mype);CHKERRQ(ierr); 1281 ierr = MPI_Comm_size(wcomm, &npe);CHKERRQ(ierr); 1282 1283 /* smooth P0 */ 1284 for (jj = 0 ; jj < pc_gamg_agg->nsmooths ; jj++) { 1285 Mat tMat; 1286 Vec diag; 1287 PetscReal alpha, emax, emin; 1288 #if defined PETSC_GAMG_USE_LOG 1289 ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET9],0,0,0,0);CHKERRQ(ierr); 1290 #endif 1291 if (jj == 0) { 1292 KSP eksp; 1293 Vec bb, xx; 1294 PC pc; 1295 ierr = MatGetVecs(Amat, &bb, 0);CHKERRQ(ierr); 1296 ierr = MatGetVecs(Amat, &xx, 0);CHKERRQ(ierr); 1297 { 1298 PetscRandom rctx; 1299 ierr = PetscRandomCreate(wcomm,&rctx);CHKERRQ(ierr); 1300 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 1301 ierr = VecSetRandom(bb,rctx);CHKERRQ(ierr); 1302 ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr); 1303 } 1304 1305 /* zeroing out BC rows -- needed for crazy matrices */ 1306 { 1307 PetscInt Istart,Iend,ncols,jj,Ii; 1308 PetscScalar zero = 0.0; 1309 ierr = MatGetOwnershipRange(Amat, &Istart, &Iend);CHKERRQ(ierr); 1310 for (Ii = Istart, jj = 0 ; Ii < Iend ; Ii++, jj++) { 1311 ierr = MatGetRow(Amat,Ii,&ncols,0,0);CHKERRQ(ierr); 1312 if (ncols <= 1) { 1313 ierr = VecSetValues(bb, 1, &Ii, &zero, INSERT_VALUES);CHKERRQ(ierr); 1314 } 1315 ierr = MatRestoreRow(Amat,Ii,&ncols,0,0);CHKERRQ(ierr); 1316 } 1317 ierr = VecAssemblyBegin(bb);CHKERRQ(ierr); 1318 ierr = VecAssemblyEnd(bb);CHKERRQ(ierr); 1319 } 1320 1321 ierr = KSPCreate(wcomm,&eksp);CHKERRQ(ierr); 1322 ierr = KSPSetTolerances(eksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,10);CHKERRQ(ierr); 1323 ierr = KSPSetNormType(eksp, KSP_NORM_NONE);CHKERRQ(ierr); 1324 ierr = KSPAppendOptionsPrefix(eksp, "gamg_est_");CHKERRQ(ierr); 1325 ierr = KSPSetFromOptions(eksp);CHKERRQ(ierr); 1326 1327 ierr = KSPSetInitialGuessNonzero(eksp, PETSC_FALSE);CHKERRQ(ierr); 1328 ierr = KSPSetOperators(eksp, Amat, Amat, SAME_NONZERO_PATTERN);CHKERRQ(ierr); 1329 ierr = KSPSetComputeSingularValues(eksp,PETSC_TRUE);CHKERRQ(ierr); 1330 1331 ierr = KSPGetPC(eksp, &pc);CHKERRQ(ierr); 1332 ierr = PCSetType(pc, PCJACOBI);CHKERRQ(ierr); /* smoother in smoothed agg. */ 1333 1334 /* solve - keep stuff out of logging */ 1335 ierr = PetscLogEventDeactivate(KSP_Solve);CHKERRQ(ierr); 1336 ierr = PetscLogEventDeactivate(PC_Apply);CHKERRQ(ierr); 1337 ierr = KSPSolve(eksp, bb, xx);CHKERRQ(ierr); 1338 ierr = PetscLogEventActivate(KSP_Solve);CHKERRQ(ierr); 1339 ierr = PetscLogEventActivate(PC_Apply);CHKERRQ(ierr); 1340 1341 ierr = KSPComputeExtremeSingularValues(eksp, &emax, &emin);CHKERRQ(ierr); 1342 if (verbose > 0) { 1343 PetscPrintf(wcomm,"\t\t\t%s smooth P0: max eigen=%e min=%e PC=%s\n", 1344 __FUNCT__,emax,emin,PCJACOBI); 1345 } 1346 ierr = VecDestroy(&xx);CHKERRQ(ierr); 1347 ierr = VecDestroy(&bb);CHKERRQ(ierr); 1348 ierr = KSPDestroy(&eksp);CHKERRQ(ierr); 1349 1350 if (pc_gamg->emax_id == -1) { 1351 ierr = PetscObjectComposedDataRegister(&pc_gamg->emax_id);CHKERRQ(ierr); 1352 assert(pc_gamg->emax_id != -1); 1353 } 1354 ierr = PetscObjectComposedDataSetScalar((PetscObject)Amat, pc_gamg->emax_id, emax);CHKERRQ(ierr); 1355 } 1356 1357 /* smooth P1 := (I - omega/lam D^{-1}A)P0 */ 1358 ierr = MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &tMat);CHKERRQ(ierr); 1359 ierr = MatGetVecs(Amat, &diag, 0);CHKERRQ(ierr); 1360 ierr = MatGetDiagonal(Amat, diag);CHKERRQ(ierr); /* effectively PCJACOBI */ 1361 ierr = VecReciprocal(diag);CHKERRQ(ierr); 1362 ierr = MatDiagonalScale(tMat, diag, 0);CHKERRQ(ierr); 1363 ierr = VecDestroy(&diag);CHKERRQ(ierr); 1364 alpha = -1.4/emax; 1365 ierr = MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN);CHKERRQ(ierr); 1366 ierr = MatDestroy(&Prol);CHKERRQ(ierr); 1367 Prol = tMat; 1368 #if defined PETSC_GAMG_USE_LOG 1369 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET9],0,0,0,0);CHKERRQ(ierr); 1370 #endif 1371 } 1372 #if defined PETSC_USE_LOG 1373 ierr = PetscLogEventEnd(PC_GAMGOptprol_AGG,0,0,0,0);CHKERRQ(ierr); 1374 #endif 1375 *a_P = Prol; 1376 1377 PetscFunctionReturn(0); 1378 } 1379 1380 /* -------------------------------------------------------------------------- */ 1381 /* 1382 PCGAMGKKTProl_AGG 1383 1384 Input Parameter: 1385 . pc - this 1386 . Prol11 - matrix on this fine level 1387 . A21 - matrix on this fine level 1388 In/Output Parameter: 1389 . a_P22 - prolongation operator to the next level 1390 */ 1391 #undef __FUNCT__ 1392 #define __FUNCT__ "PCGAMGKKTProl_AGG" 1393 PetscErrorCode PCGAMGKKTProl_AGG(PC pc, 1394 const Mat Prol11, 1395 const Mat A21, 1396 Mat *a_P22) 1397 { 1398 PetscErrorCode ierr; 1399 PC_MG *mg = (PC_MG*)pc->data; 1400 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1401 const PetscInt verbose = pc_gamg->verbose; 1402 /* PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; */ 1403 PetscMPIInt mype,npe; 1404 Mat Prol22,Tmat,Gmat; 1405 MPI_Comm wcomm = ((PetscObject)pc)->comm; 1406 PetscCoarsenData *agg_lists; 1407 1408 PetscFunctionBegin; 1409 #if defined PETSC_USE_LOG 1410 ierr = PetscLogEventBegin(PC_GAMGKKTProl_AGG,0,0,0,0);CHKERRQ(ierr); 1411 #endif 1412 ierr = MPI_Comm_rank(wcomm, &mype);CHKERRQ(ierr); 1413 ierr = MPI_Comm_size(wcomm, &npe);CHKERRQ(ierr); 1414 1415 /* form C graph */ 1416 ierr = MatMatMult(A21, Prol11, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Tmat);CHKERRQ(ierr); 1417 ierr = MatMatTransposeMult(Tmat, Tmat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat);CHKERRQ(ierr); 1418 ierr = MatDestroy(&Tmat);CHKERRQ(ierr); 1419 ierr = PCGAMGFilterGraph(&Gmat, 0.0, PETSC_FALSE, verbose);CHKERRQ(ierr); 1420 1421 /* coarsen constraints */ 1422 { 1423 MatCoarsen crs; 1424 ierr = MatCoarsenCreate(wcomm, &crs);CHKERRQ(ierr); 1425 ierr = MatCoarsenSetType(crs, MATCOARSENMIS);CHKERRQ(ierr); 1426 ierr = MatCoarsenSetAdjacency(crs, Gmat);CHKERRQ(ierr); 1427 ierr = MatCoarsenSetVerbose(crs, verbose);CHKERRQ(ierr); 1428 ierr = MatCoarsenSetStrictAggs(crs, PETSC_TRUE);CHKERRQ(ierr); 1429 ierr = MatCoarsenApply(crs);CHKERRQ(ierr); 1430 ierr = MatCoarsenGetData(crs, &agg_lists);CHKERRQ(ierr); 1431 ierr = MatCoarsenDestroy(&crs);CHKERRQ(ierr); 1432 } 1433 1434 /* form simple prolongation 'Prol22' */ 1435 { 1436 PetscInt ii,mm,clid,my0,nloc,nLocalSelected; 1437 PetscScalar val = 1.0; 1438 /* get 'nLocalSelected' */ 1439 ierr = MatGetLocalSize(Gmat, &nloc, &ii);CHKERRQ(ierr); 1440 for (ii=0, nLocalSelected = 0 ; ii < nloc ; ii++) { 1441 PetscBool ise; 1442 /* filter out singletons 0 or 1? */ 1443 ierr = PetscCDEmptyAt(agg_lists, ii, &ise);CHKERRQ(ierr); 1444 if (!ise) nLocalSelected++; 1445 } 1446 1447 ierr = MatCreate(wcomm,&Prol22);CHKERRQ(ierr); 1448 ierr = MatSetSizes(Prol22,nloc, nLocalSelected,PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1449 ierr = MatSetType(Prol22, MATAIJ);CHKERRQ(ierr); 1450 ierr = MatSeqAIJSetPreallocation(Prol22,1,PETSC_NULL);CHKERRQ(ierr); 1451 ierr = MatMPIAIJSetPreallocation(Prol22,1,PETSC_NULL,1,PETSC_NULL);CHKERRQ(ierr); 1452 /* ierr = MatCreateAIJ(wcomm, */ 1453 /* nloc, nLocalSelected, */ 1454 /* PETSC_DETERMINE, PETSC_DETERMINE, */ 1455 /* 1, PETSC_NULL, 1, PETSC_NULL, */ 1456 /* &Prol22); */ 1457 1458 ierr = MatGetOwnershipRange(Prol22, &my0, &ii);CHKERRQ(ierr); 1459 nloc = ii - my0; 1460 1461 /* make aggregates */ 1462 for (mm = clid = 0 ; mm < nloc ; mm++) { 1463 ierr = PetscCDSizeAt(agg_lists, mm, &ii);CHKERRQ(ierr); 1464 if (ii > 0) { 1465 PetscInt asz=ii,cgid=my0+clid,rids[1000]; 1466 PetscCDPos pos; 1467 if (asz>1000)SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Very large aggregate: %d",asz); 1468 ii = 0; 1469 ierr = PetscCDGetHeadPos(agg_lists,mm,&pos);CHKERRQ(ierr); 1470 while (pos) { 1471 PetscInt gid1; 1472 ierr = PetscLLNGetID(pos, &gid1);CHKERRQ(ierr); 1473 ierr = PetscCDGetNextPos(agg_lists,mm,&pos);CHKERRQ(ierr); 1474 1475 rids[ii++] = gid1; 1476 } 1477 assert(ii==asz); 1478 /* add diagonal block of P0 */ 1479 ierr = MatSetValues(Prol22,asz,rids,1,&cgid,&val,INSERT_VALUES);CHKERRQ(ierr); 1480 1481 clid++; 1482 } /* coarse agg */ 1483 } /* for all fine nodes */ 1484 ierr = MatAssemblyBegin(Prol22,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1485 ierr = MatAssemblyEnd(Prol22,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1486 } 1487 1488 /* clean up */ 1489 ierr = MatDestroy(&Gmat);CHKERRQ(ierr); 1490 ierr = PetscCDDestroy(agg_lists);CHKERRQ(ierr); 1491 #if defined PETSC_USE_LOG 1492 ierr = PetscLogEventEnd(PC_GAMGKKTProl_AGG,0,0,0,0);CHKERRQ(ierr); 1493 #endif 1494 *a_P22 = Prol22; 1495 1496 PetscFunctionReturn(0); 1497 } 1498 1499 /* -------------------------------------------------------------------------- */ 1500 /* 1501 PCCreateGAMG_AGG 1502 1503 Input Parameter: 1504 . pc - 1505 */ 1506 #undef __FUNCT__ 1507 #define __FUNCT__ "PCCreateGAMG_AGG" 1508 PetscErrorCode PCCreateGAMG_AGG(PC pc) 1509 { 1510 PetscErrorCode ierr; 1511 PC_MG *mg = (PC_MG*)pc->data; 1512 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1513 PC_GAMG_AGG *pc_gamg_agg; 1514 1515 PetscFunctionBegin; 1516 if (pc_gamg->subctx) { 1517 /* call base class */ 1518 ierr = PCDestroy_GAMG(pc);CHKERRQ(ierr); 1519 } 1520 1521 /* create sub context for SA */ 1522 ierr = PetscNewLog(pc, PC_GAMG_AGG, &pc_gamg_agg);CHKERRQ(ierr); 1523 pc_gamg->subctx = pc_gamg_agg; 1524 1525 pc->ops->setfromoptions = PCSetFromOptions_GAMG_AGG; 1526 pc->ops->destroy = PCDestroy_AGG; 1527 /* reset does not do anything; setup not virtual */ 1528 1529 /* set internal function pointers */ 1530 pc_gamg->graph = PCGAMGgraph_AGG; 1531 pc_gamg->coarsen = PCGAMGCoarsen_AGG; 1532 pc_gamg->prolongator = PCGAMGProlongator_AGG; 1533 pc_gamg->optprol = PCGAMGOptprol_AGG; 1534 pc_gamg->formkktprol = PCGAMGKKTProl_AGG; 1535 1536 pc_gamg->createdefaultdata = PCSetData_AGG; 1537 1538 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCSetCoordinates_C","PCSetCoordinates_AGG",PCSetCoordinates_AGG);CHKERRQ(ierr); 1539 PetscFunctionReturn(0); 1540 } 1541