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