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