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