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