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