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