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