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