1 /* 2 GAMG geometric-algebric multiogrid PC - Mark Adams 2011 3 */ 4 #include <../src/ksp/pc/impls/gamg/gamg.h> /*I "petscpc.h" I*/ 5 #include "private/matimpl.h" 6 7 #if defined PETSC_USE_LOG 8 PetscLogEvent gamg_setup_events[NUM_SET]; 9 #endif 10 #define GAMG_MAXLEVELS 30 11 12 /*#define GAMG_STAGES*/ 13 #if (defined PETSC_USE_LOG && defined GAMG_STAGES) 14 static PetscLogStage gamg_stages[GAMG_MAXLEVELS]; 15 #endif 16 17 /* Private context for the GAMG preconditioner */ 18 typedef struct gamg_TAG { 19 PetscInt m_dim; 20 PetscInt m_Nlevels; 21 PetscInt m_data_sz; 22 PetscInt m_data_rows; 23 PetscInt m_data_cols; 24 PetscInt m_count; 25 PetscInt m_method; /* 0: geomg; 1: plain agg 'pa'; 2: smoothed agg 'sa' */ 26 PetscReal *m_data; /* blocked vector of vertex data on fine grid (coordinates) */ 27 char m_type[64]; 28 PetscBool m_avoid_repart; 29 PetscInt m_min_eq_proc; 30 PetscReal m_threshold; 31 PetscBool m_verbose; 32 } PC_GAMG; 33 34 /* -------------------------------------------------------------------------- */ 35 /* 36 PCSetCoordinates_GAMG 37 38 Input Parameter: 39 . pc - the preconditioner context 40 */ 41 EXTERN_C_BEGIN 42 #undef __FUNCT__ 43 #define __FUNCT__ "PCSetCoordinates_GAMG" 44 PetscErrorCode PCSetCoordinates_GAMG( PC a_pc, PetscInt a_ndm, PetscReal *a_coords ) 45 { 46 PC_MG *mg = (PC_MG*)a_pc->data; 47 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 48 PetscErrorCode ierr; 49 PetscInt arrsz,bs,my0,kk,ii,jj,nloc,Iend; 50 Mat Amat = a_pc->pmat; 51 52 PetscFunctionBegin; 53 PetscValidHeaderSpecific( Amat, MAT_CLASSID, 1 ); 54 ierr = MatGetBlockSize( Amat, &bs ); CHKERRQ( ierr ); 55 ierr = MatGetOwnershipRange( Amat, &my0, &Iend ); CHKERRQ(ierr); 56 nloc = (Iend-my0)/bs; 57 if((Iend-my0)%bs!=0) SETERRQ1(((PetscObject)Amat)->comm,PETSC_ERR_ARG_WRONG, "Bad local size %d.",nloc); 58 59 pc_gamg->m_data_rows = 1; 60 if(a_coords==0 && pc_gamg->m_method==0) { 61 SETERRQ(((PetscObject)Amat)->comm,PETSC_ERR_ARG_WRONG, "Need coordinates for pc_gamg_type 'geo'."); 62 } 63 if( pc_gamg->m_method==0 ) pc_gamg->m_data_cols = a_ndm; /* coordinates */ 64 else{ /* SA: null space vectors */ 65 if(a_coords != 0 && bs==1 ) pc_gamg->m_data_cols = 1; /* scalar w/ coords and SA (not needed) */ 66 else if(a_coords != 0 ) pc_gamg->m_data_cols = (a_ndm==2 ? 3 : 6); /* elasticity */ 67 else pc_gamg->m_data_cols = bs; /* no data, force SA with constant null space vectors */ 68 pc_gamg->m_data_rows = bs; 69 } 70 arrsz = nloc*pc_gamg->m_data_rows*pc_gamg->m_data_cols; 71 72 /* create data - syntactic sugar that should be refactored at some point */ 73 if (pc_gamg->m_data==0 || (pc_gamg->m_data_sz != arrsz)) { 74 ierr = PetscFree( pc_gamg->m_data ); CHKERRQ(ierr); 75 ierr = PetscMalloc((arrsz+1)*sizeof(PetscReal), &pc_gamg->m_data ); CHKERRQ(ierr); 76 } 77 for(kk=0;kk<arrsz;kk++)pc_gamg->m_data[kk] = -999.; 78 pc_gamg->m_data[arrsz] = -99.; 79 /* copy data in - column oriented */ 80 if( pc_gamg->m_method != 0 ) { 81 const PetscInt M = Iend - my0; 82 for(kk=0;kk<nloc;kk++){ 83 PetscReal *data = &pc_gamg->m_data[kk*bs]; 84 if( pc_gamg->m_data_cols==1 ) *data = 1.0; 85 else { 86 for(ii=0;ii<bs;ii++) 87 for(jj=0;jj<bs;jj++) 88 if(ii==jj)data[ii*M + jj] = 1.0; /* translational modes */ 89 else data[ii*M + jj] = 0.0; 90 if( a_coords != 0 ) { 91 if( a_ndm == 2 ){ /* rotational modes */ 92 data += 2*M; 93 data[0] = -a_coords[2*kk+1]; 94 data[1] = a_coords[2*kk]; 95 } 96 else { 97 data += 3*M; 98 data[0] = 0.0; data[M+0] = a_coords[3*kk+2]; data[2*M+0] = -a_coords[3*kk+1]; 99 data[1] = -a_coords[3*kk+2]; data[M+1] = 0.0; data[2*M+1] = a_coords[3*kk]; 100 data[2] = a_coords[3*kk+1]; data[M+2] = -a_coords[3*kk]; data[2*M+2] = 0.0; 101 } 102 } 103 } 104 } 105 } 106 else { 107 for( kk = 0 ; kk < nloc ; kk++ ){ 108 for( ii = 0 ; ii < a_ndm ; ii++ ) { 109 pc_gamg->m_data[ii*nloc + kk] = a_coords[kk*a_ndm + ii]; 110 } 111 } 112 } 113 assert(pc_gamg->m_data[arrsz] == -99.); 114 115 pc_gamg->m_data_sz = arrsz; 116 pc_gamg->m_dim = a_ndm; 117 118 PetscFunctionReturn(0); 119 } 120 EXTERN_C_END 121 122 123 /* ----------------------------------------------------------------------------- */ 124 #undef __FUNCT__ 125 #define __FUNCT__ "PCReset_GAMG" 126 PetscErrorCode PCReset_GAMG(PC pc) 127 { 128 PetscErrorCode ierr; 129 PC_MG *mg = (PC_MG*)pc->data; 130 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 131 132 PetscFunctionBegin; 133 if( pc_gamg->m_data != 0 ) { /* this should not happen, cleaned up in SetUp */ 134 ierr = PetscFree(pc_gamg->m_data); CHKERRQ(ierr); 135 } 136 pc_gamg->m_data = 0; pc_gamg->m_data_sz = 0; 137 PetscFunctionReturn(0); 138 } 139 140 /* -------------------------------------------------------------------------- */ 141 /* 142 PCGAMGPartitionLevel 143 144 Input Parameter: 145 . a_Amat_fine - matrix on this fine (k) level 146 . a_ndata_rows - size of data to move (coarse grid) 147 . a_ndata_cols - size of data to move (coarse grid) 148 In/Output Parameter: 149 . a_P_inout - prolongation operator to the next level (k-1) 150 . a_coarse_data - data that need to be moved 151 . a_nactive_proc - number of active procs 152 Output Parameter: 153 . a_Amat_crs - coarse matrix that is created (k-1) 154 . a_avoid_repart - 155 . a_min_eq_proc - goal for number of eqs. to have per proc. after proc reduction 156 */ 157 158 #undef __FUNCT__ 159 #define __FUNCT__ "PCGAMGPartitionLevel" 160 PetscErrorCode PCGAMGPartitionLevel(PC pc, Mat a_Amat_fine, 161 PetscInt a_ndata_rows, 162 PetscInt a_ndata_cols, 163 PetscInt a_cbs, 164 Mat *a_P_inout, 165 PetscReal **a_coarse_data, 166 PetscMPIInt *a_nactive_proc, 167 Mat *a_Amat_crs, 168 PetscBool a_avoid_repart, 169 PetscInt a_min_eq_proc 170 ) 171 { 172 PC_MG *mg = (PC_MG*)pc->data; 173 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 174 PetscErrorCode ierr; 175 Mat Cmat,Pnew,Pold=*a_P_inout; 176 IS new_indices,isnum; 177 MPI_Comm wcomm = ((PetscObject)a_Amat_fine)->comm; 178 PetscMPIInt mype,npe,new_npe,nactive; 179 PetscInt neq,NN,Istart,Iend,Istart0,Iend0,ncrs_new,ncrs0; 180 181 PetscFunctionBegin; 182 ierr = MPI_Comm_rank( wcomm, &mype ); CHKERRQ(ierr); 183 ierr = MPI_Comm_size( wcomm, &npe ); CHKERRQ(ierr); 184 185 /* RAP */ 186 #ifdef USE_R 187 /* make R wih brute force for now */ 188 ierr = MatTranspose( Pold, Pnew ); 189 ierr = MatDestroy( &Pold ); CHKERRQ(ierr); 190 ierr = MatRARt( a_Amat_fine, Pnew, MAT_INITIAL_MATRIX, 2.0, &Cmat ); CHKERRQ(ierr); 191 Pold = Pnew; 192 #else 193 ierr = MatPtAP( a_Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, &Cmat ); CHKERRQ(ierr); 194 #endif 195 ierr = MatSetBlockSize( Cmat, a_cbs ); CHKERRQ(ierr); 196 ierr = MatGetOwnershipRange( Cmat, &Istart0, &Iend0 ); CHKERRQ(ierr); 197 ncrs0 = (Iend0-Istart0)/a_cbs; assert((Iend0-Istart0)%a_cbs == 0); 198 199 if( a_avoid_repart) { 200 *a_Amat_crs = Cmat; /* output */ 201 } 202 else { 203 /* Repartition Cmat_{k} and move colums of P^{k}_{k-1} and coordinates accordingly */ 204 Mat adj; 205 const PetscInt *idx,data_sz=a_ndata_rows*a_ndata_cols; 206 const PetscInt stride0=ncrs0*a_ndata_rows; 207 PetscInt is_sz,*isnewproc_idx,ii,jj,kk,strideNew,*tidx; 208 /* create sub communicator */ 209 MPI_Comm cm; 210 MPI_Group wg, g2; 211 PetscInt *counts,inpe; 212 PetscMPIInt *ranks; 213 IS isscat; 214 PetscScalar *array; 215 Vec src_crd, dest_crd; 216 PetscReal *data = *a_coarse_data; 217 VecScatter vecscat; 218 IS isnewproc; 219 220 /* get number of PEs to make active, reduce */ 221 ierr = MatGetSize( Cmat, &neq, &NN ); CHKERRQ(ierr); 222 new_npe = neq/a_min_eq_proc; /* hardwire min. number of eq/proc */ 223 if( new_npe == 0 || neq < 2*a_min_eq_proc ) new_npe = 1; 224 else if (new_npe >= *a_nactive_proc ) new_npe = *a_nactive_proc; /* no change, rare */ 225 226 ierr = PetscMalloc( npe*sizeof(PetscMPIInt), &ranks ); CHKERRQ(ierr); 227 ierr = PetscMalloc( npe*sizeof(PetscInt), &counts ); CHKERRQ(ierr); 228 229 ierr = MPI_Allgather( &ncrs0, 1, MPIU_INT, counts, 1, MPIU_INT, wcomm ); CHKERRQ(ierr); 230 assert(counts[mype]==ncrs0); 231 /* count real active pes */ 232 for( nactive = jj = 0 ; jj < npe ; jj++) { 233 if( counts[jj] != 0 ) { 234 ranks[nactive++] = jj; 235 } 236 } 237 238 if (nactive < new_npe) new_npe = nactive; /* this can happen with empty input procs */ 239 240 if (pc_gamg->m_verbose) PetscPrintf(PETSC_COMM_WORLD,"\t[%d]%s npe (active): %d --> %d. new npe = %d, neq = %d\n",mype,__FUNCT__,*a_nactive_proc,nactive,new_npe,neq); 241 242 *a_nactive_proc = new_npe; /* output */ 243 244 ierr = MPI_Comm_group( wcomm, &wg ); CHKERRQ(ierr); 245 ierr = MPI_Group_incl( wg, nactive, ranks, &g2 ); CHKERRQ(ierr); 246 ierr = MPI_Comm_create( wcomm, g2, &cm ); CHKERRQ(ierr); 247 248 ierr = MPI_Group_free( &wg ); CHKERRQ(ierr); 249 ierr = MPI_Group_free( &g2 ); CHKERRQ(ierr); 250 251 /* MatPartitioningApply call MatConvert, which is collective */ 252 #if defined PETSC_USE_LOG 253 ierr = PetscLogEventBegin(gamg_setup_events[SET12],0,0,0,0);CHKERRQ(ierr); 254 #endif 255 if( a_cbs == 1) { 256 ierr = MatConvert( Cmat, MATMPIADJ, MAT_INITIAL_MATRIX, &adj ); CHKERRQ(ierr); 257 } 258 else{ 259 /* make a scalar matrix to partition */ 260 Mat tMat; 261 PetscInt ncols,jj,Ii; 262 const PetscScalar *vals; 263 const PetscInt *idx; 264 PetscInt *d_nnz; 265 static int llev = 0; 266 267 ierr = PetscMalloc( ncrs0*sizeof(PetscInt), &d_nnz ); CHKERRQ(ierr); 268 for ( Ii = Istart0, jj = 0 ; Ii < Iend0 ; Ii += a_cbs, jj++ ) { 269 ierr = MatGetRow(Cmat,Ii,&ncols,0,0); CHKERRQ(ierr); 270 d_nnz[jj] = ncols/a_cbs; 271 if( d_nnz[jj] > ncrs0 ) d_nnz[jj] = ncrs0; 272 ierr = MatRestoreRow(Cmat,Ii,&ncols,0,0); CHKERRQ(ierr); 273 } 274 275 ierr = MatCreateMPIAIJ( wcomm, ncrs0, ncrs0, 276 PETSC_DETERMINE, PETSC_DETERMINE, 277 0, d_nnz, 0, d_nnz, 278 &tMat ); 279 CHKERRQ(ierr); 280 ierr = PetscFree( d_nnz ); CHKERRQ(ierr); 281 282 for ( ii = Istart0; ii < Iend0; ii++ ) { 283 PetscInt dest_row = ii/a_cbs; 284 ierr = MatGetRow(Cmat,ii,&ncols,&idx,&vals); CHKERRQ(ierr); 285 for( jj = 0 ; jj < ncols ; jj++ ){ 286 PetscInt dest_col = idx[jj]/a_cbs; 287 PetscScalar v = 1.0; 288 ierr = MatSetValues(tMat,1,&dest_row,1,&dest_col,&v,ADD_VALUES); CHKERRQ(ierr); 289 } 290 ierr = MatRestoreRow(Cmat,ii,&ncols,&idx,&vals); CHKERRQ(ierr); 291 } 292 ierr = MatAssemblyBegin(tMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 293 ierr = MatAssemblyEnd(tMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 294 295 if( llev++ == -1 ) { 296 PetscViewer viewer; char fname[32]; 297 sprintf(fname,"part_mat_%d.mat",llev); 298 PetscViewerBinaryOpen(wcomm,fname,FILE_MODE_WRITE,&viewer); 299 ierr = MatView( tMat, viewer ); CHKERRQ(ierr); 300 ierr = PetscViewerDestroy( &viewer ); 301 } 302 303 ierr = MatConvert( tMat, MATMPIADJ, MAT_INITIAL_MATRIX, &adj ); CHKERRQ(ierr); 304 305 ierr = MatDestroy( &tMat ); CHKERRQ(ierr); 306 } 307 if( ncrs0 != 0 ){ 308 const PetscInt *is_idx; 309 MatPartitioning mpart; 310 /* hack to fix global data that pmetis.c uses in 'adj' */ 311 for( nactive = jj = 0 ; jj < npe ; jj++ ) { 312 if( counts[jj] != 0 ) { 313 adj->rmap->range[nactive++] = adj->rmap->range[jj]; 314 } 315 } 316 adj->rmap->range[nactive] = adj->rmap->range[npe]; 317 318 ierr = MatPartitioningCreate( cm, &mpart ); CHKERRQ(ierr); 319 ierr = MatPartitioningSetAdjacency( mpart, adj ); CHKERRQ(ierr); 320 ierr = MatPartitioningSetFromOptions( mpart ); CHKERRQ(ierr); 321 ierr = MatPartitioningSetNParts( mpart, new_npe );CHKERRQ(ierr); 322 ierr = MatPartitioningApply( mpart, &isnewproc ); CHKERRQ(ierr); 323 ierr = MatPartitioningDestroy( &mpart ); CHKERRQ(ierr); 324 325 /* collect IS info */ 326 ierr = ISGetLocalSize( isnewproc, &is_sz ); CHKERRQ(ierr); 327 ierr = PetscMalloc( a_cbs*is_sz*sizeof(PetscInt), &isnewproc_idx ); CHKERRQ(ierr); 328 ierr = ISGetIndices( isnewproc, &is_idx ); CHKERRQ(ierr); 329 /* spread partitioning across machine - best way ??? */ 330 NN = 1; /*npe/new_npe;*/ 331 for( kk = jj = 0 ; kk < is_sz ; kk++ ){ 332 for( ii = 0 ; ii < a_cbs ; ii++, jj++ ) { 333 isnewproc_idx[jj] = is_idx[kk] * NN; /* distribution */ 334 } 335 } 336 ierr = ISRestoreIndices( isnewproc, &is_idx ); CHKERRQ(ierr); 337 ierr = ISDestroy( &isnewproc ); CHKERRQ(ierr); 338 ierr = MPI_Comm_free( &cm ); CHKERRQ(ierr); 339 340 is_sz *= a_cbs; 341 } 342 else{ 343 isnewproc_idx = 0; 344 is_sz = 0; 345 } 346 ierr = MatDestroy( &adj ); CHKERRQ(ierr); 347 ierr = ISCreateGeneral( wcomm, is_sz, isnewproc_idx, PETSC_COPY_VALUES, &isnewproc ); 348 if( isnewproc_idx != 0 ) { 349 ierr = PetscFree( isnewproc_idx ); CHKERRQ(ierr); 350 } 351 352 /* 353 Create an index set from the isnewproc index set to indicate the mapping TO 354 */ 355 ierr = ISPartitioningToNumbering( isnewproc, &isnum ); CHKERRQ(ierr); 356 /* 357 Determine how many elements are assigned to each processor 358 */ 359 inpe = npe; 360 ierr = ISPartitioningCount( isnewproc, inpe, counts ); CHKERRQ(ierr); 361 ierr = ISDestroy( &isnewproc ); CHKERRQ(ierr); 362 ncrs_new = counts[mype]/a_cbs; 363 strideNew = ncrs_new*a_ndata_rows; 364 #if defined PETSC_USE_LOG 365 ierr = PetscLogEventEnd(gamg_setup_events[SET12],0,0,0,0); CHKERRQ(ierr); 366 #endif 367 /* Create a vector to contain the newly ordered element information */ 368 ierr = VecCreate( wcomm, &dest_crd ); 369 ierr = VecSetSizes( dest_crd, data_sz*ncrs_new, PETSC_DECIDE ); CHKERRQ(ierr); 370 ierr = VecSetFromOptions( dest_crd ); CHKERRQ(ierr); /* this is needed! */ 371 /* 372 There are 'a_ndata_rows*a_ndata_cols' data items per node, (one can think of the vectors of having 373 a block size of ...). Note, ISs are expanded into equation space by 'a_cbs'. 374 */ 375 ierr = PetscMalloc( (ncrs0*data_sz)*sizeof(PetscInt), &tidx ); CHKERRQ(ierr); 376 ierr = ISGetIndices( isnum, &idx ); CHKERRQ(ierr); 377 for(ii=0,jj=0; ii<ncrs0 ; ii++) { 378 PetscInt id = idx[ii*a_cbs]/a_cbs; /* get node back */ 379 for( kk=0; kk<data_sz ; kk++, jj++) tidx[jj] = id*data_sz + kk; 380 } 381 ierr = ISRestoreIndices( isnum, &idx ); CHKERRQ(ierr); 382 ierr = ISCreateGeneral( wcomm, data_sz*ncrs0, tidx, PETSC_COPY_VALUES, &isscat ); 383 CHKERRQ(ierr); 384 ierr = PetscFree( tidx ); CHKERRQ(ierr); 385 /* 386 Create a vector to contain the original vertex information for each element 387 */ 388 ierr = VecCreateSeq( PETSC_COMM_SELF, data_sz*ncrs0, &src_crd ); CHKERRQ(ierr); 389 for( jj=0; jj<a_ndata_cols ; jj++ ) { 390 for( ii=0 ; ii<ncrs0 ; ii++) { 391 for( kk=0; kk<a_ndata_rows ; kk++ ) { 392 PetscInt ix = ii*a_ndata_rows + kk + jj*stride0, jx = ii*data_sz + kk*a_ndata_cols + jj; 393 PetscScalar tt = (PetscScalar)data[ix]; 394 ierr = VecSetValues( src_crd, 1, &jx, &tt, INSERT_VALUES ); CHKERRQ(ierr); 395 } 396 } 397 } 398 ierr = VecAssemblyBegin(src_crd); CHKERRQ(ierr); 399 ierr = VecAssemblyEnd(src_crd); CHKERRQ(ierr); 400 /* 401 Scatter the element vertex information (still in the original vertex ordering) 402 to the correct processor 403 */ 404 ierr = VecScatterCreate( src_crd, PETSC_NULL, dest_crd, isscat, &vecscat); 405 CHKERRQ(ierr); 406 ierr = ISDestroy( &isscat ); CHKERRQ(ierr); 407 ierr = VecScatterBegin(vecscat,src_crd,dest_crd,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 408 ierr = VecScatterEnd(vecscat,src_crd,dest_crd,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 409 ierr = VecScatterDestroy( &vecscat ); CHKERRQ(ierr); 410 ierr = VecDestroy( &src_crd ); CHKERRQ(ierr); 411 /* 412 Put the element vertex data into a new allocation of the gdata->ele 413 */ 414 ierr = PetscFree( *a_coarse_data ); CHKERRQ(ierr); 415 ierr = PetscMalloc( data_sz*ncrs_new*sizeof(PetscReal), a_coarse_data ); CHKERRQ(ierr); 416 417 ierr = VecGetArray( dest_crd, &array ); CHKERRQ(ierr); 418 data = *a_coarse_data; 419 for( jj=0; jj<a_ndata_cols ; jj++ ) { 420 for( ii=0 ; ii<ncrs_new ; ii++) { 421 for( kk=0; kk<a_ndata_rows ; kk++ ) { 422 PetscInt ix = ii*a_ndata_rows + kk + jj*strideNew, jx = ii*data_sz + kk*a_ndata_cols + jj; 423 data[ix] = PetscRealPart(array[jx]); 424 array[jx] = 1.e300; 425 } 426 } 427 } 428 ierr = VecRestoreArray( dest_crd, &array ); CHKERRQ(ierr); 429 ierr = VecDestroy( &dest_crd ); CHKERRQ(ierr); 430 /* 431 Invert for MatGetSubMatrix 432 */ 433 ierr = ISInvertPermutation( isnum, ncrs_new*a_cbs, &new_indices ); CHKERRQ(ierr); 434 ierr = ISSort( new_indices ); CHKERRQ(ierr); /* is this needed? */ 435 ierr = ISDestroy( &isnum ); CHKERRQ(ierr); 436 /* A_crs output */ 437 ierr = MatGetSubMatrix( Cmat, new_indices, new_indices, MAT_INITIAL_MATRIX, a_Amat_crs ); 438 CHKERRQ(ierr); 439 440 ierr = MatDestroy( &Cmat ); CHKERRQ(ierr); 441 Cmat = *a_Amat_crs; /* output */ 442 ierr = MatSetBlockSize( Cmat, a_cbs ); CHKERRQ(ierr); 443 444 /* prolongator */ 445 ierr = MatGetOwnershipRange( Pold, &Istart, &Iend ); CHKERRQ(ierr); 446 { 447 IS findices; 448 ierr = ISCreateStride(wcomm,Iend-Istart,Istart,1,&findices); CHKERRQ(ierr); 449 #ifdef USE_R 450 ierr = MatGetSubMatrix( Pold, new_indices, findices, MAT_INITIAL_MATRIX, &Pnew ); 451 #else 452 ierr = MatGetSubMatrix( Pold, findices, new_indices, MAT_INITIAL_MATRIX, &Pnew ); 453 #endif 454 CHKERRQ(ierr); 455 456 ierr = ISDestroy( &findices ); CHKERRQ(ierr); 457 } 458 459 ierr = MatDestroy( a_P_inout ); CHKERRQ(ierr); 460 *a_P_inout = Pnew; /* output */ 461 462 ierr = ISDestroy( &new_indices ); CHKERRQ(ierr); 463 ierr = PetscFree( counts ); CHKERRQ(ierr); 464 ierr = PetscFree( ranks ); CHKERRQ(ierr); 465 } 466 467 PetscFunctionReturn(0); 468 } 469 470 /* -------------------------------------------------------------------------- */ 471 /* 472 PCSetUp_GAMG - Prepares for the use of the GAMG preconditioner 473 by setting data structures and options. 474 475 Input Parameter: 476 . pc - the preconditioner context 477 478 Application Interface Routine: PCSetUp() 479 480 Notes: 481 The interface routine PCSetUp() is not usually called directly by 482 the user, but instead is called by PCApply() if necessary. 483 */ 484 #undef __FUNCT__ 485 #define __FUNCT__ "PCSetUp_GAMG" 486 PetscErrorCode PCSetUp_GAMG( PC a_pc ) 487 { 488 PetscErrorCode ierr; 489 PC_MG *mg = (PC_MG*)a_pc->data; 490 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 491 PC_MG_Levels **mglevels = mg->levels; 492 Mat Amat = a_pc->mat, Pmat = a_pc->pmat; 493 PetscInt fine_level, level, level1, M, N, bs, nloc, lidx, Istart, Iend; 494 MPI_Comm wcomm = ((PetscObject)a_pc)->comm; 495 PetscMPIInt mype,npe,nactivepe; 496 PetscBool isOK; 497 Mat Aarr[GAMG_MAXLEVELS], Parr[GAMG_MAXLEVELS]; 498 PetscReal *coarse_data = 0, *data, emaxs[GAMG_MAXLEVELS]; 499 MatInfo info; 500 501 PetscFunctionBegin; 502 pc_gamg->m_count++; 503 if( a_pc->setupcalled > 0 ) { 504 /* just do Galerkin grids */ 505 Mat B,dA,dB; 506 507 /* PCSetUp_MG seems to insists on setting this to GMRES */ 508 ierr = KSPSetType( mglevels[0]->smoothd, KSPPREONLY ); CHKERRQ(ierr); 509 510 /* currently only handle case where mat and pmat are the same on coarser levels */ 511 ierr = KSPGetOperators(mglevels[pc_gamg->m_Nlevels-1]->smoothd,&dA,&dB,PETSC_NULL);CHKERRQ(ierr); 512 /* (re)set to get dirty flag */ 513 ierr = KSPSetOperators(mglevels[pc_gamg->m_Nlevels-1]->smoothd,dA,dB,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 514 ierr = KSPSetUp( mglevels[pc_gamg->m_Nlevels-1]->smoothd ); CHKERRQ(ierr); 515 516 for (level=pc_gamg->m_Nlevels-2; level>-1; level--) { 517 ierr = KSPGetOperators(mglevels[level]->smoothd,PETSC_NULL,&B,PETSC_NULL);CHKERRQ(ierr); 518 /* the first time through the matrix structure has changed from repartitioning */ 519 if( pc_gamg->m_count == 2 ) { 520 ierr = MatDestroy( &B ); CHKERRQ(ierr); 521 ierr = MatPtAP(dB,mglevels[level+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);CHKERRQ(ierr); 522 mglevels[level]->A = B; 523 } 524 else { 525 ierr = MatPtAP(dB,mglevels[level+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);CHKERRQ(ierr); 526 } 527 ierr = KSPSetOperators(mglevels[level]->smoothd,B,B,SAME_NONZERO_PATTERN); CHKERRQ(ierr); 528 dB = B; 529 /* setup KSP/PC */ 530 ierr = KSPSetUp( mglevels[level]->smoothd ); CHKERRQ(ierr); 531 } 532 533 #define PRINT_MATS PETSC_FALSE 534 /* plot levels - A */ 535 if( PRINT_MATS ) { 536 for (lidx=0, level=pc_gamg->m_Nlevels-1; level>0 ; level--,lidx++){ 537 PetscViewer viewer; 538 char fname[32]; KSP smoother; Mat Tmat, TTm; 539 ierr = PCMGGetSmoother( a_pc, lidx, &smoother ); CHKERRQ(ierr); 540 ierr = KSPGetOperators( smoother, &Tmat, &TTm, 0 ); CHKERRQ(ierr); 541 sprintf(fname,"Amat_%d_%d.m",(int)pc_gamg->m_count,(int)level); 542 ierr = PetscViewerASCIIOpen( wcomm, fname, &viewer ); CHKERRQ(ierr); 543 ierr = PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB); CHKERRQ(ierr); 544 ierr = MatView( Tmat, viewer ); CHKERRQ(ierr); 545 ierr = PetscViewerDestroy( &viewer ); 546 } 547 } 548 549 a_pc->setupcalled = 2; 550 551 PetscFunctionReturn(0); 552 } 553 554 ierr = MPI_Comm_rank(wcomm,&mype);CHKERRQ(ierr); 555 ierr = MPI_Comm_size(wcomm,&npe);CHKERRQ(ierr); 556 /* GAMG requires input of fine-grid matrix. It determines nlevels. */ 557 ierr = MatGetBlockSize( Amat, &bs ); CHKERRQ(ierr); 558 ierr = MatGetSize( Amat, &M, &N );CHKERRQ(ierr); 559 ierr = MatGetOwnershipRange( Amat, &Istart, &Iend ); CHKERRQ(ierr); 560 nloc = (Iend-Istart)/bs; assert((Iend-Istart)%bs == 0); 561 562 /* get data of not around */ 563 if( pc_gamg->m_data == 0 && nloc > 0 ) { 564 ierr = PCSetCoordinates_GAMG( a_pc, -1, 0 ); CHKERRQ( ierr ); 565 } 566 data = pc_gamg->m_data; 567 568 /* Get A_i and R_i */ 569 ierr = MatGetInfo(Amat,MAT_GLOBAL_SUM,&info); CHKERRQ(ierr); 570 if (pc_gamg->m_verbose) PetscPrintf(PETSC_COMM_WORLD,"\t[%d]%s level %d N=%d, n data rows=%d, n data cols=%d, nnz/row (ave)=%d, np=%d\n", 571 mype,__FUNCT__,0,N,pc_gamg->m_data_rows,pc_gamg->m_data_cols, 572 (int)(info.nz_used/(PetscReal)N),npe); 573 for ( level=0, Aarr[0] = Pmat, nactivepe = npe; /* hard wired stopping logic */ 574 level < (GAMG_MAXLEVELS-1) && (level==0 || M>2*pc_gamg->m_min_eq_proc); /* && (npe==1 || nactivepe>1); */ 575 level++ ){ 576 level1 = level + 1; 577 #if (defined PETSC_USE_LOG && defined GAMG_STAGES) 578 ierr = PetscLogStagePush(gamg_stages[level]); CHKERRQ( ierr ); 579 #endif 580 #if defined PETSC_USE_LOG 581 ierr = PetscLogEventBegin(gamg_setup_events[SET1],0,0,0,0);CHKERRQ(ierr); 582 #endif 583 ierr = createProlongation(Aarr[level], data, pc_gamg->m_dim, pc_gamg->m_data_cols, pc_gamg->m_method, 584 level, pc_gamg->m_threshold, &bs, &Parr[level1], &coarse_data, &isOK, 585 &emaxs[level], pc_gamg->m_verbose ); 586 CHKERRQ(ierr); 587 ierr = PetscFree( data ); CHKERRQ( ierr ); 588 #if defined PETSC_USE_LOG 589 ierr = PetscLogEventEnd(gamg_setup_events[SET1],0,0,0,0);CHKERRQ(ierr); 590 #endif 591 if(level==0) Aarr[0] = Amat; /* use Pmat for finest level setup, but use mat for solver */ 592 if( isOK ) { 593 #if defined PETSC_USE_LOG 594 ierr = PetscLogEventBegin(gamg_setup_events[SET2],0,0,0,0);CHKERRQ(ierr); 595 #endif 596 ierr = PCGAMGPartitionLevel(a_pc, Aarr[level], (pc_gamg->m_method != 0) ? bs : 1, pc_gamg->m_data_cols, bs, 597 &Parr[level1], &coarse_data, &nactivepe, &Aarr[level1], 598 pc_gamg->m_avoid_repart, pc_gamg->m_min_eq_proc ); 599 CHKERRQ(ierr); 600 #if defined PETSC_USE_LOG 601 ierr = PetscLogEventEnd(gamg_setup_events[SET2],0,0,0,0);CHKERRQ(ierr); 602 #endif 603 ierr = MatGetSize( Aarr[level1], &M, &N );CHKERRQ(ierr); 604 ierr = MatGetInfo(Aarr[level1],MAT_GLOBAL_SUM,&info); CHKERRQ(ierr); 605 if (pc_gamg->m_verbose) PetscPrintf(PETSC_COMM_WORLD,"\t\t[%d]%s %d) N=%d, n data cols=%d, nnz/row (ave)=%d, %d active pes\n", 606 mype,__FUNCT__,(int)level1,N,pc_gamg->m_data_cols, 607 (int)(info.nz_used/(PetscReal)N),nactivepe); 608 /* coarse grids with SA can have zero row/cols from singleton aggregates */ 609 /* aggregation method should gaurrentee this does not happen! */ 610 611 if (pc_gamg->m_verbose) { 612 Vec diag; PetscScalar *data_arr,v; PetscInt Istart,Iend,kk,nloceq,id; 613 v = 1.e-10; /* LU factor has hard wired numbers for small diags so this needs to match (yuk) */ 614 ierr = MatGetOwnershipRange(Aarr[level1], &Istart, &Iend); CHKERRQ(ierr); 615 nloceq = Iend-Istart; 616 ierr = MatGetVecs( Aarr[level1], &diag, 0 ); CHKERRQ(ierr); 617 ierr = MatGetDiagonal( Aarr[level1], diag ); CHKERRQ(ierr); 618 ierr = VecGetArray( diag, &data_arr ); CHKERRQ(ierr); 619 for(kk=0;kk<nloceq;kk++){ 620 if(data_arr[kk]==0.0) { 621 id = kk + Istart; 622 ierr = MatSetValues(Aarr[level1],1,&id,1,&id,&v,INSERT_VALUES); 623 CHKERRQ(ierr); 624 PetscPrintf(PETSC_COMM_SELF,"\t[%d]%s warning: added zero to diag (%d) on level %d \n",mype,__FUNCT__,id,level1); 625 } 626 } 627 ierr = VecRestoreArray( diag, &data_arr ); CHKERRQ(ierr); 628 ierr = VecDestroy( &diag ); CHKERRQ(ierr); 629 ierr = MatAssemblyBegin(Aarr[level1],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 630 ierr = MatAssemblyEnd(Aarr[level1],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 631 } 632 } else { 633 coarse_data = 0; 634 break; 635 } 636 data = coarse_data; 637 638 #if (defined PETSC_USE_LOG && defined GAMG_STAGES) 639 ierr = PetscLogStagePop(); CHKERRQ( ierr ); 640 #endif 641 } 642 if( coarse_data ) { 643 ierr = PetscFree( coarse_data ); CHKERRQ( ierr ); 644 } 645 if (pc_gamg->m_verbose) PetscPrintf(PETSC_COMM_WORLD,"\t[%d]%s %d levels\n",0,__FUNCT__,level + 1); 646 pc_gamg->m_data = 0; /* destroyed coordinate data */ 647 pc_gamg->m_Nlevels = level + 1; 648 fine_level = level; 649 ierr = PCMGSetLevels(a_pc,pc_gamg->m_Nlevels,PETSC_NULL);CHKERRQ(ierr); 650 651 /* set default smoothers */ 652 for ( lidx=1, level = pc_gamg->m_Nlevels-2; 653 lidx <= fine_level; 654 lidx++, level--) { 655 PetscReal emax, emin; 656 KSP smoother; PC subpc; 657 PetscBool isCheb; 658 /* set defaults */ 659 ierr = PCMGGetSmoother( a_pc, lidx, &smoother ); CHKERRQ(ierr); 660 ierr = KSPSetType( smoother, KSPCHEBYCHEV );CHKERRQ(ierr); 661 ierr = KSPGetPC( smoother, &subpc ); CHKERRQ(ierr); 662 /* ierr = KSPSetTolerances(smoother,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,2); CHKERRQ(ierr); */ 663 ierr = PCSetType( subpc, PETSC_GAMG_SMOOTHER ); CHKERRQ(ierr); 664 ierr = KSPSetNormType( smoother, KSP_NORM_NONE ); CHKERRQ(ierr); 665 /* overide defaults with input parameters */ 666 ierr = KSPSetFromOptions( smoother ); CHKERRQ(ierr); 667 668 ierr = KSPSetOperators( smoother, Aarr[level], Aarr[level], SAME_NONZERO_PATTERN ); CHKERRQ(ierr); 669 /* do my own cheby */ 670 ierr = PetscTypeCompare( (PetscObject)smoother, KSPCHEBYCHEV, &isCheb ); CHKERRQ(ierr); 671 if( isCheb ) { 672 ierr = PetscTypeCompare( (PetscObject)subpc, PETSC_GAMG_SMOOTHER, &isCheb ); CHKERRQ(ierr); 673 if( isCheb && emaxs[level] > 0.0 ) emax=emaxs[level]; /* eigen estimate only for diagnal PC */ 674 else{ /* eigen estimate 'emax' */ 675 KSP eksp; Mat Lmat = Aarr[level]; 676 Vec bb, xx; PC pc; 677 const PCType type; 678 679 ierr = PCGetType( subpc, &type ); CHKERRQ(ierr); 680 ierr = MatGetVecs( Lmat, &bb, 0 ); CHKERRQ(ierr); 681 ierr = MatGetVecs( Lmat, &xx, 0 ); CHKERRQ(ierr); 682 { 683 PetscRandom rctx; 684 ierr = PetscRandomCreate(wcomm,&rctx);CHKERRQ(ierr); 685 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 686 ierr = VecSetRandom(bb,rctx);CHKERRQ(ierr); 687 ierr = PetscRandomDestroy( &rctx ); CHKERRQ(ierr); 688 } 689 ierr = KSPCreate(wcomm,&eksp);CHKERRQ(ierr); 690 ierr = KSPSetType( eksp, KSPCG ); CHKERRQ(ierr); 691 ierr = KSPSetTolerances( eksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 10 ); 692 CHKERRQ(ierr); 693 ierr = KSPSetNormType( eksp, KSP_NORM_NONE ); CHKERRQ(ierr); 694 695 ierr = KSPAppendOptionsPrefix( eksp, "est_"); CHKERRQ(ierr); 696 ierr = KSPSetFromOptions( eksp ); CHKERRQ(ierr); 697 698 ierr = KSPSetInitialGuessNonzero( eksp, PETSC_FALSE ); CHKERRQ(ierr); 699 ierr = KSPSetOperators( eksp, Lmat, Lmat, SAME_NONZERO_PATTERN ); CHKERRQ( ierr ); 700 ierr = KSPGetPC( eksp, &pc );CHKERRQ( ierr ); 701 ierr = PCSetType( pc, type ); CHKERRQ(ierr); /* should be same as eigen estimates op. */ 702 703 ierr = KSPSetComputeSingularValues( eksp,PETSC_TRUE ); CHKERRQ(ierr); 704 ierr = KSPSolve( eksp, bb, xx ); CHKERRQ(ierr); 705 ierr = KSPComputeExtremeSingularValues( eksp, &emax, &emin ); CHKERRQ(ierr); 706 ierr = VecDestroy( &xx ); CHKERRQ(ierr); 707 ierr = VecDestroy( &bb ); CHKERRQ(ierr); 708 ierr = KSPDestroy( &eksp ); CHKERRQ(ierr); 709 710 if (pc_gamg->m_verbose) { 711 PetscPrintf(PETSC_COMM_WORLD,"\t\t\t%s PC setup max eigen=%e min=%e PC=%s\n", 712 __FUNCT__,emax,emin,PETSC_GAMG_SMOOTHER); 713 } 714 } 715 { 716 PetscInt N1, N0, tt; 717 ierr = MatGetSize( Aarr[level], &N1, &tt ); CHKERRQ(ierr); 718 ierr = MatGetSize( Aarr[level+1], &N0, &tt ); CHKERRQ(ierr); 719 /* heuristic - is this crap? */ 720 emin = 1.*emax/((PetscReal)N1/(PetscReal)N0); 721 emax *= 1.05; 722 } 723 ierr = KSPChebychevSetEigenvalues( smoother, emax, emin );CHKERRQ(ierr); 724 } 725 } 726 { 727 /* coarse grid */ 728 KSP smoother,*k2; PC subpc,pc2; PetscInt ii,first; 729 Mat Lmat = Aarr[pc_gamg->m_Nlevels-1]; 730 ierr = PCMGGetSmoother( a_pc, 0, &smoother ); CHKERRQ(ierr); 731 ierr = KSPSetOperators( smoother, Lmat, Lmat, SAME_NONZERO_PATTERN ); CHKERRQ(ierr); 732 ierr = KSPSetNormType( smoother, KSP_NORM_NONE ); CHKERRQ(ierr); 733 ierr = KSPGetPC( smoother, &subpc ); CHKERRQ(ierr); 734 ierr = PCSetType( subpc, PCBJACOBI ); CHKERRQ(ierr); 735 ierr = PCSetUp( subpc ); CHKERRQ(ierr); 736 ierr = PCBJacobiGetSubKSP(subpc,&ii,&first,&k2);CHKERRQ(ierr); 737 assert(ii==1); 738 ierr = KSPGetPC(k2[0],&pc2);CHKERRQ(ierr); 739 ierr = PCSetType( pc2, PCLU ); CHKERRQ(ierr); 740 } 741 742 /* should be called in PCSetFromOptions_GAMG(), but cannot be called prior to PCMGSetLevels() */ 743 ierr = PCSetFromOptions_MG(a_pc); CHKERRQ(ierr); 744 { 745 PetscBool galerkin; 746 ierr = PCMGGetGalerkin( a_pc, &galerkin); CHKERRQ(ierr); 747 if(galerkin){ 748 SETERRQ(wcomm,PETSC_ERR_ARG_WRONG, "GAMG does galerkin manually so it must not be used in PC_MG."); 749 } 750 } 751 752 /* plot levels - R/P */ 753 if( PRINT_MATS ) { 754 for (level=pc_gamg->m_Nlevels-1;level>0;level--){ 755 PetscViewer viewer; 756 char fname[32]; 757 sprintf(fname,"Pmat_%d_%d.m",(int)pc_gamg->m_count,(int)level); 758 ierr = PetscViewerASCIIOpen( wcomm, fname, &viewer ); CHKERRQ(ierr); 759 ierr = PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB); CHKERRQ(ierr); 760 ierr = MatView( Parr[level], viewer ); CHKERRQ(ierr); 761 ierr = PetscViewerDestroy( &viewer ); 762 sprintf(fname,"Amat_%d_%d.m",(int)pc_gamg->m_count,(int)level); 763 ierr = PetscViewerASCIIOpen( wcomm, fname, &viewer ); CHKERRQ(ierr); 764 ierr = PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB); CHKERRQ(ierr); 765 ierr = MatView( Aarr[level], viewer ); CHKERRQ(ierr); 766 ierr = PetscViewerDestroy( &viewer ); 767 } 768 } 769 770 /* set interpolation between the levels, clean up */ 771 for (lidx=0,level=pc_gamg->m_Nlevels-1; 772 lidx<fine_level; 773 lidx++, level--){ 774 ierr = PCMGSetInterpolation( a_pc, lidx+1, Parr[level] );CHKERRQ(ierr); 775 ierr = MatDestroy( &Parr[level] ); CHKERRQ(ierr); 776 ierr = MatDestroy( &Aarr[level] ); CHKERRQ(ierr); 777 } 778 779 /* setupcalled is set to 0 so that MG is setup from scratch */ 780 a_pc->setupcalled = 0; 781 ierr = PCSetUp_MG( a_pc );CHKERRQ( ierr ); 782 a_pc->setupcalled = 1; /* use 1 as signal that this has not been re-setup */ 783 784 { 785 KSP smoother; /* PCSetUp_MG seems to insists on setting this to GMRES on coarse grid */ 786 ierr = PCMGGetSmoother( a_pc, 0, &smoother ); CHKERRQ(ierr); 787 ierr = KSPSetType( smoother, KSPPREONLY ); CHKERRQ(ierr); 788 ierr = KSPSetUp( smoother ); CHKERRQ(ierr); 789 } 790 791 PetscFunctionReturn(0); 792 } 793 794 /* ------------------------------------------------------------------------- */ 795 /* 796 PCDestroy_GAMG - Destroys the private context for the GAMG preconditioner 797 that was created with PCCreate_GAMG(). 798 799 Input Parameter: 800 . pc - the preconditioner context 801 802 Application Interface Routine: PCDestroy() 803 */ 804 #undef __FUNCT__ 805 #define __FUNCT__ "PCDestroy_GAMG" 806 PetscErrorCode PCDestroy_GAMG(PC pc) 807 { 808 PetscErrorCode ierr; 809 PC_MG *mg = (PC_MG*)pc->data; 810 PC_GAMG *pc_gamg= (PC_GAMG*)mg->innerctx; 811 812 PetscFunctionBegin; 813 ierr = PCReset_GAMG(pc);CHKERRQ(ierr); 814 ierr = PetscFree(pc_gamg);CHKERRQ(ierr); 815 ierr = PCDestroy_MG(pc);CHKERRQ(ierr); 816 PetscFunctionReturn(0); 817 } 818 819 820 #undef __FUNCT__ 821 #define __FUNCT__ "PCGAMGSetProcEqLim" 822 /*@ 823 PCGAMGSetProcEqLim - Set number of equations to aim for on coarse grids via 824 processor reduction. 825 826 Not Collective on PC 827 828 Input Parameters: 829 . pc - the preconditioner context 830 831 832 Options Database Key: 833 . -pc_gamg_process_eq_limit 834 835 Level: intermediate 836 837 Concepts: Unstructured multrigrid preconditioner 838 839 .seealso: () 840 @*/ 841 PetscErrorCode PCGAMGSetProcEqLim(PC pc, PetscInt n) 842 { 843 PetscErrorCode ierr; 844 845 PetscFunctionBegin; 846 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 847 ierr = PetscTryMethod(pc,"PCGAMGSetProcEqLim_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr); 848 PetscFunctionReturn(0); 849 } 850 851 EXTERN_C_BEGIN 852 #undef __FUNCT__ 853 #define __FUNCT__ "PCGAMGSetProcEqLim_GAMG" 854 PetscErrorCode PCGAMGSetProcEqLim_GAMG(PC pc, PetscInt n) 855 { 856 PC_MG *mg = (PC_MG*)pc->data; 857 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 858 859 PetscFunctionBegin; 860 if(n>0) pc_gamg->m_min_eq_proc = n; 861 PetscFunctionReturn(0); 862 } 863 EXTERN_C_END 864 865 #undef __FUNCT__ 866 #define __FUNCT__ "PCGAMGAvoidRepartitioning" 867 /*@ 868 PCGAMGAvoidRepartitioning - Do not repartition the coarse grids 869 870 Collective on PC 871 872 Input Parameters: 873 . pc - the preconditioner context 874 875 876 Options Database Key: 877 . -pc_gamg_avoid_repartitioning 878 879 Level: intermediate 880 881 Concepts: Unstructured multrigrid preconditioner 882 883 .seealso: () 884 @*/ 885 PetscErrorCode PCGAMGAvoidRepartitioning(PC pc, PetscBool n) 886 { 887 PetscErrorCode ierr; 888 889 PetscFunctionBegin; 890 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 891 ierr = PetscTryMethod(pc,"PCGAMGAvoidRepartitioning_C",(PC,PetscBool),(pc,n));CHKERRQ(ierr); 892 PetscFunctionReturn(0); 893 } 894 895 EXTERN_C_BEGIN 896 #undef __FUNCT__ 897 #define __FUNCT__ "PCGAMGAvoidRepartitioning_GAMG" 898 PetscErrorCode PCGAMGAvoidRepartitioning_GAMG(PC pc, PetscBool n) 899 { 900 PC_MG *mg = (PC_MG*)pc->data; 901 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 902 903 PetscFunctionBegin; 904 pc_gamg->m_avoid_repart = n; 905 PetscFunctionReturn(0); 906 } 907 EXTERN_C_END 908 909 #undef __FUNCT__ 910 #define __FUNCT__ "PCGAMGSetThreshold" 911 /*@ 912 PCGAMGSetThreshold - Relative threshold to use for dropping edges in aggregation graph 913 914 Not collective on PC 915 916 Input Parameters: 917 . pc - the preconditioner context 918 919 920 Options Database Key: 921 . -pc_gamg_threshold 922 923 Level: intermediate 924 925 Concepts: Unstructured multrigrid preconditioner 926 927 .seealso: () 928 @*/ 929 PetscErrorCode PCGAMGSetThreshold(PC pc, PetscReal n) 930 { 931 PetscErrorCode ierr; 932 933 PetscFunctionBegin; 934 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 935 ierr = PetscTryMethod(pc,"PCGAMGSetThreshold_C",(PC,PetscReal),(pc,n));CHKERRQ(ierr); 936 PetscFunctionReturn(0); 937 } 938 939 EXTERN_C_BEGIN 940 #undef __FUNCT__ 941 #define __FUNCT__ "PCGAMGSetThreshold_GAMG" 942 PetscErrorCode PCGAMGSetThreshold_GAMG(PC pc, PetscReal n) 943 { 944 PC_MG *mg = (PC_MG*)pc->data; 945 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 946 947 PetscFunctionBegin; 948 pc_gamg->m_threshold = n; 949 PetscFunctionReturn(0); 950 } 951 EXTERN_C_END 952 953 #undef __FUNCT__ 954 #define __FUNCT__ "PCGAMGSetSolverType" 955 /*@ 956 PCGAMGSetSolverType - Set solution method. 957 958 Collective on PC 959 960 Input Parameters: 961 . pc - the preconditioner context 962 963 964 Options Database Key: 965 . -pc_gamg_type 966 967 Level: intermediate 968 969 Concepts: Unstructured multrigrid preconditioner 970 971 .seealso: () 972 @*/ 973 PetscErrorCode PCGAMGSetSolverType(PC pc, char str[], PetscInt sz ) 974 { 975 PetscErrorCode ierr; 976 977 PetscFunctionBegin; 978 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 979 ierr = PetscTryMethod(pc,"PCGAMGSetSolverType_C",(PC,char[],PetscInt),(pc,str,sz)); 980 CHKERRQ(ierr); 981 PetscFunctionReturn(0); 982 } 983 984 EXTERN_C_BEGIN 985 #undef __FUNCT__ 986 #define __FUNCT__ "PCGAMGSetSolverType_GAMG" 987 PetscErrorCode PCGAMGSetSolverType_GAMG(PC pc, char str[], PetscInt sz ) 988 { 989 PC_MG *mg = (PC_MG*)pc->data; 990 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 991 992 PetscFunctionBegin; 993 if(sz < 64) strcpy(pc_gamg->m_type,str); 994 PetscFunctionReturn(0); 995 } 996 EXTERN_C_END 997 998 #undef __FUNCT__ 999 #define __FUNCT__ "PCSetFromOptions_GAMG" 1000 PetscErrorCode PCSetFromOptions_GAMG(PC pc) 1001 { 1002 PetscErrorCode ierr; 1003 PC_MG *mg = (PC_MG*)pc->data; 1004 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1005 PetscBool flag; 1006 1007 PetscFunctionBegin; 1008 ierr = PetscOptionsHead("GAMG options"); CHKERRQ(ierr); 1009 { 1010 ierr = PetscOptionsString("-pc_gamg_type", 1011 "Solver type: plane aggregation ('pa'), smoothed aggregation ('sa') or geometric multigrid (default)", 1012 "PCGAMGSetSolverType", 1013 pc_gamg->m_type, 1014 pc_gamg->m_type, 1015 64, 1016 &flag ); 1017 CHKERRQ(ierr); 1018 if( flag && pc_gamg->m_data != 0 ) { 1019 if( (strcmp(pc_gamg->m_type,"sa")==0 && pc_gamg->m_method != 2) || 1020 (strcmp(pc_gamg->m_type,"pa")==0 && pc_gamg->m_method != 1) || 1021 (strcmp(pc_gamg->m_type,"geo")==0 && pc_gamg->m_method != 0) ) { 1022 SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "PCSetFromOptions called of PCSetCoordinates (with new method, after data was created)"); 1023 } 1024 } 1025 1026 /* -pc_gamg_verbose */ 1027 ierr = PetscOptionsBool("-pc_gamg_verbose","Verbose (debugging) output for PCGAMG","none",pc_gamg->m_verbose,&pc_gamg->m_verbose,PETSC_NULL);CHKERRQ(ierr); 1028 1029 pc_gamg->m_method = 1; /* default to plane aggregation */ 1030 if (flag ) { 1031 if( strcmp(pc_gamg->m_type,"sa") == 0) pc_gamg->m_method = 2; 1032 else if( strcmp(pc_gamg->m_type,"pa") == 0) pc_gamg->m_method = 1; 1033 else if( strcmp(pc_gamg->m_type,"geo") == 0) pc_gamg->m_method = 0; 1034 else SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "Invalid gamg type: %s",pc_gamg->m_type); 1035 } 1036 /* -pc_gamg_avoid_repartitioning */ 1037 pc_gamg->m_avoid_repart = PETSC_FALSE; 1038 ierr = PetscOptionsBool("-pc_gamg_avoid_repartitioning", 1039 "Do not repartion coarse grids (false)", 1040 "PCGAMGAvoidRepartitioning", 1041 pc_gamg->m_avoid_repart, 1042 &pc_gamg->m_avoid_repart, 1043 &flag); 1044 CHKERRQ(ierr); 1045 1046 /* -pc_gamg_process_eq_limit */ 1047 pc_gamg->m_min_eq_proc = 600; 1048 ierr = PetscOptionsInt("-pc_gamg_process_eq_limit", 1049 "Limit (goal) on number of equations per process on coarse grids", 1050 "PCGAMGSetProcEqLim", 1051 pc_gamg->m_min_eq_proc, 1052 &pc_gamg->m_min_eq_proc, 1053 &flag ); 1054 CHKERRQ(ierr); 1055 1056 /* -pc_gamg_threshold */ 1057 pc_gamg->m_threshold = 0.05; 1058 ierr = PetscOptionsReal("-pc_gamg_threshold", 1059 "Relative threshold to use for dropping edges in aggregation graph", 1060 "PCGAMGSetThreshold", 1061 pc_gamg->m_threshold, 1062 &pc_gamg->m_threshold, 1063 &flag ); 1064 CHKERRQ(ierr); 1065 } 1066 ierr = PetscOptionsTail();CHKERRQ(ierr); 1067 1068 PetscFunctionReturn(0); 1069 } 1070 1071 /* -------------------------------------------------------------------------- */ 1072 /*MC 1073 PCGAMG - Geometric algebraic multigrid (AMG) preconditioning. This preconditioner currently has two 1074 AMG methods: 1) an unstructured geometric method, which requires that you provide coordinates for each 1075 vertex, and 2) smoothed aggregation. Smoothed aggregation (SA) can work without coordinates but it 1076 will generate some common non-trivial null spaces if coordinates are provided. The input fine grid matrix 1077 must have the block size set for 'system' problems (with multiple dof per vertex/cell) to work properly. 1078 SA will generate rotational rigid body mode null space vectors, in addition to the trivial translational 1079 modes, when coordinates are provide in 2D and 3D. 1080 1081 Options Database Keys: 1082 Multigrid options(inherited) 1083 + -pc_mg_cycles <1>: 1 for V cycle, 2 for W-cycle (PCMGSetCycleType) 1084 . -pc_mg_smoothup <1>: Number of post-smoothing steps (PCMGSetNumberSmoothUp) 1085 . -pc_mg_smoothdown <1>: Number of pre-smoothing steps (PCMGSetNumberSmoothDown) 1086 - -pc_mg_type <multiplicative>: (one of) additive multiplicative full cascade kascade 1087 1088 Level: intermediate 1089 1090 Concepts: multigrid 1091 1092 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, 1093 PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(), PCMGSetNumberSmoothDown(), 1094 PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(), 1095 PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(), 1096 PCMGSetCyclesOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR() 1097 M*/ 1098 1099 EXTERN_C_BEGIN 1100 #undef __FUNCT__ 1101 #define __FUNCT__ "PCCreate_GAMG" 1102 PetscErrorCode PCCreate_GAMG(PC pc) 1103 { 1104 PetscErrorCode ierr; 1105 PC_GAMG *pc_gamg; 1106 PC_MG *mg; 1107 PetscClassId cookie; 1108 #if defined PETSC_USE_LOG 1109 static int count = 0; 1110 #endif 1111 1112 PetscFunctionBegin; 1113 /* PCGAMG is an inherited class of PCMG. Initialize pc as PCMG */ 1114 ierr = PCSetType(pc,PCMG);CHKERRQ(ierr); /* calls PCCreate_MG() and MGCreate_Private() */ 1115 ierr = PetscObjectChangeTypeName((PetscObject)pc,PCGAMG);CHKERRQ(ierr); 1116 1117 /* create a supporting struct and attach it to pc */ 1118 ierr = PetscNewLog(pc,PC_GAMG,&pc_gamg);CHKERRQ(ierr); 1119 pc_gamg->m_data_sz = 0; pc_gamg->m_data = 0; pc_gamg->m_count = 0; 1120 mg = (PC_MG*)pc->data; 1121 mg->innerctx = pc_gamg; 1122 1123 pc_gamg->m_Nlevels = -1; 1124 1125 /* overwrite the pointers of PCMG by the functions of PCGAMG */ 1126 pc->ops->setfromoptions = PCSetFromOptions_GAMG; 1127 pc->ops->setup = PCSetUp_GAMG; 1128 pc->ops->reset = PCReset_GAMG; 1129 pc->ops->destroy = PCDestroy_GAMG; 1130 1131 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1132 "PCSetCoordinates_C", 1133 "PCSetCoordinates_GAMG", 1134 PCSetCoordinates_GAMG); 1135 CHKERRQ(ierr); 1136 1137 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1138 "PCGAMGSetProcEqLim_C", 1139 "PCGAMGSetProcEqLim_GAMG", 1140 PCGAMGSetProcEqLim_GAMG); 1141 CHKERRQ(ierr); 1142 1143 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1144 "PCGAMGAvoidRepartitioning_C", 1145 "PCGAMGAvoidRepartitioning_GAMG", 1146 PCGAMGAvoidRepartitioning_GAMG); 1147 CHKERRQ(ierr); 1148 1149 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1150 "PCGAMGSetThreshold_C", 1151 "PCGAMGSetThreshold_GAMG", 1152 PCGAMGSetThreshold_GAMG); 1153 CHKERRQ(ierr); 1154 1155 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1156 "PCGAMGSetSolverType_C", 1157 "PCGAMGSetSolverType_GAMG", 1158 PCGAMGSetSolverType_GAMG); 1159 CHKERRQ(ierr); 1160 1161 #if defined PETSC_USE_LOG 1162 if( count++ == 0 ) { 1163 PetscClassIdRegister("GAMG Setup",&cookie); 1164 PetscLogEventRegister("GAMG: createProl", cookie, &gamg_setup_events[SET1]); 1165 PetscLogEventRegister(" Graph", cookie, &gamg_setup_events[GRAPH]); 1166 PetscLogEventRegister(" G.Mat", cookie, &gamg_setup_events[GRAPH_MAT]); 1167 PetscLogEventRegister(" G.Filter", cookie, &gamg_setup_events[GRAPH_FILTER]); 1168 PetscLogEventRegister(" G.Square", cookie, &gamg_setup_events[GRAPH_SQR]); 1169 PetscLogEventRegister(" MIS/Agg", cookie, &gamg_setup_events[SET4]); 1170 PetscLogEventRegister(" geo: growSupp", cookie, &gamg_setup_events[SET5]); 1171 PetscLogEventRegister(" geo: triangle", cookie, &gamg_setup_events[SET6]); 1172 PetscLogEventRegister(" search&set", cookie, &gamg_setup_events[FIND_V]); 1173 PetscLogEventRegister(" SA: init", cookie, &gamg_setup_events[SET7]); 1174 /* PetscLogEventRegister(" SA: frmProl0", cookie, &gamg_setup_events[SET8]); */ 1175 PetscLogEventRegister(" SA: smooth", cookie, &gamg_setup_events[SET9]); 1176 PetscLogEventRegister("GAMG: partLevel", cookie, &gamg_setup_events[SET2]); 1177 PetscLogEventRegister(" PL repartition", cookie, &gamg_setup_events[SET12]); 1178 /* PetscLogEventRegister(" PL move data", cookie, &gamg_setup_events[SET13]); */ 1179 /* PetscLogEventRegister("GAMG: fix", cookie, &gamg_setup_events[SET10]); */ 1180 /* PetscLogEventRegister("GAMG: set levels", cookie, &gamg_setup_events[SET11]); */ 1181 1182 /* create timer stages */ 1183 #if defined GAMG_STAGES 1184 { 1185 char str[32]; 1186 sprintf(str,"MG Level %d (finest)",0); 1187 PetscLogStageRegister(str, &gamg_stages[0]); 1188 PetscInt lidx; 1189 for (lidx=1;lidx<9;lidx++){ 1190 sprintf(str,"MG Level %d",lidx); 1191 PetscLogStageRegister(str, &gamg_stages[lidx]); 1192 } 1193 } 1194 #endif 1195 } 1196 #endif 1197 PetscFunctionReturn(0); 1198 } 1199 EXTERN_C_END 1200