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) 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->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); 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->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", 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], pc_gamg->m_verbose ); 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->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", 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->m_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->m_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 = KSPAppendOptionsPrefix( eksp, "est_"); CHKERRQ(ierr); 671 ierr = KSPSetType( eksp, KSPCG ); CHKERRQ(ierr); 672 ierr = KSPSetFromOptions( eksp ); CHKERRQ(ierr); 673 ierr = KSPSetInitialGuessNonzero( eksp, PETSC_FALSE ); CHKERRQ(ierr); 674 ierr = KSPSetOperators( eksp, Lmat, Lmat, SAME_NONZERO_PATTERN ); CHKERRQ( ierr ); 675 ierr = KSPGetPC( eksp, &pc );CHKERRQ( ierr ); 676 ierr = PCSetType( pc, PETSC_GAMG_SMOOTHER ); CHKERRQ(ierr); /* should be same as eigen estimates op. */ 677 ierr = KSPSetTolerances( eksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 10 ); 678 CHKERRQ(ierr); 679 ierr = KSPSetNormType( eksp, KSP_NORM_NONE ); CHKERRQ(ierr); 680 681 ierr = KSPSetComputeSingularValues( eksp,PETSC_TRUE ); CHKERRQ(ierr); 682 ierr = KSPSolve( eksp, bb, xx ); CHKERRQ(ierr); 683 ierr = KSPComputeExtremeSingularValues( eksp, &emax, &emin ); CHKERRQ(ierr); 684 ierr = VecDestroy( &xx ); CHKERRQ(ierr); 685 ierr = VecDestroy( &bb ); CHKERRQ(ierr); 686 ierr = KSPDestroy( &eksp ); CHKERRQ(ierr); 687 if (pc_gamg->m_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); 688 } 689 { 690 PetscInt N1, N0, tt; 691 ierr = MatGetSize( Aarr[level], &N1, &tt ); CHKERRQ(ierr); 692 ierr = MatGetSize( Aarr[level+1], &N0, &tt ); CHKERRQ(ierr); 693 emin = 1.*emax/((PetscReal)N1/(PetscReal)N0); /* this should be about the coarsening rate */ 694 emax *= 1.05; 695 } 696 697 ierr = KSPSetOperators( smoother, Aarr[level], Aarr[level], SAME_NONZERO_PATTERN ); 698 ierr = KSPChebychevSetEigenvalues( smoother, emax, emin );CHKERRQ(ierr); 699 /* ierr = KSPSetTolerances(smoother,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,2); CHKERRQ(ierr); */ 700 ierr = KSPGetPC( smoother, &subpc ); CHKERRQ(ierr); 701 ierr = PCSetType( subpc, PETSC_GAMG_SMOOTHER ); CHKERRQ(ierr); 702 ierr = KSPSetNormType( smoother, KSP_NORM_NONE ); CHKERRQ(ierr); 703 } 704 { 705 /* coarse grid */ 706 KSP smoother,*k2; PC subpc,pc2; PetscInt ii,first; 707 Mat Lmat = Aarr[pc_gamg->m_Nlevels-1]; 708 ierr = PCMGGetSmoother( a_pc, 0, &smoother ); CHKERRQ(ierr); 709 ierr = KSPSetOperators( smoother, Lmat, Lmat, SAME_NONZERO_PATTERN ); CHKERRQ(ierr); 710 ierr = KSPSetNormType( smoother, KSP_NORM_NONE ); CHKERRQ(ierr); 711 ierr = KSPGetPC( smoother, &subpc ); CHKERRQ(ierr); 712 ierr = PCSetType( subpc, PCBJACOBI ); CHKERRQ(ierr); 713 ierr = PCSetUp( subpc ); CHKERRQ(ierr); 714 ierr = PCBJacobiGetSubKSP(subpc,&ii,&first,&k2);CHKERRQ(ierr); 715 assert(ii==1); 716 ierr = KSPGetPC(k2[0],&pc2);CHKERRQ(ierr); 717 ierr = PCSetType( pc2, PCLU ); CHKERRQ(ierr); 718 } 719 720 /* should be called in PCSetFromOptions_GAMG(), but cannot be called prior to PCMGSetLevels() */ 721 ierr = PCSetFromOptions_MG(a_pc); CHKERRQ(ierr); 722 { 723 PetscBool galerkin; 724 ierr = PCMGGetGalerkin( a_pc, &galerkin); CHKERRQ(ierr); 725 if(galerkin){ 726 SETERRQ(wcomm,PETSC_ERR_ARG_WRONG, "GAMG does galerkin manually so it must not be used in PC_MG."); 727 } 728 } 729 730 /* plot levels - R/P */ 731 if( PRINT_MATS ) { 732 for (level=pc_gamg->m_Nlevels-1;level>0;level--){ 733 PetscViewer viewer; 734 char fname[32]; 735 sprintf(fname,"Pmat_%d_%d.m",(int)pc_gamg->m_count,(int)level); 736 ierr = PetscViewerASCIIOpen( wcomm, fname, &viewer ); CHKERRQ(ierr); 737 ierr = PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB); CHKERRQ(ierr); 738 ierr = MatView( Parr[level], viewer ); CHKERRQ(ierr); 739 ierr = PetscViewerDestroy( &viewer ); 740 sprintf(fname,"Amat_%d_%d.m",(int)pc_gamg->m_count,(int)level); 741 ierr = PetscViewerASCIIOpen( wcomm, fname, &viewer ); CHKERRQ(ierr); 742 ierr = PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB); CHKERRQ(ierr); 743 ierr = MatView( Aarr[level], viewer ); CHKERRQ(ierr); 744 ierr = PetscViewerDestroy( &viewer ); 745 } 746 } 747 748 /* set interpolation between the levels, clean up */ 749 for (lidx=0,level=pc_gamg->m_Nlevels-1; 750 lidx<fine_level; 751 lidx++, level--){ 752 ierr = PCMGSetInterpolation( a_pc, lidx+1, Parr[level] );CHKERRQ(ierr); 753 ierr = MatDestroy( &Parr[level] ); CHKERRQ(ierr); 754 ierr = MatDestroy( &Aarr[level] ); CHKERRQ(ierr); 755 } 756 757 /* setupcalled is set to 0 so that MG is setup from scratch */ 758 a_pc->setupcalled = 0; 759 ierr = PCSetUp_MG( a_pc );CHKERRQ( ierr ); 760 a_pc->setupcalled = 1; /* use 1 as signal that this has not been re-setup */ 761 762 { 763 KSP smoother; /* PCSetUp_MG seems to insists on setting this to GMRES on coarse grid */ 764 ierr = PCMGGetSmoother( a_pc, 0, &smoother ); CHKERRQ(ierr); 765 ierr = KSPSetType( smoother, KSPPREONLY ); CHKERRQ(ierr); 766 ierr = KSPSetUp( smoother ); CHKERRQ(ierr); 767 } 768 769 PetscFunctionReturn(0); 770 } 771 772 /* ------------------------------------------------------------------------- */ 773 /* 774 PCDestroy_GAMG - Destroys the private context for the GAMG preconditioner 775 that was created with PCCreate_GAMG(). 776 777 Input Parameter: 778 . pc - the preconditioner context 779 780 Application Interface Routine: PCDestroy() 781 */ 782 #undef __FUNCT__ 783 #define __FUNCT__ "PCDestroy_GAMG" 784 PetscErrorCode PCDestroy_GAMG(PC pc) 785 { 786 PetscErrorCode ierr; 787 PC_MG *mg = (PC_MG*)pc->data; 788 PC_GAMG *pc_gamg= (PC_GAMG*)mg->innerctx; 789 790 PetscFunctionBegin; 791 ierr = PCReset_GAMG(pc);CHKERRQ(ierr); 792 ierr = PetscFree(pc_gamg);CHKERRQ(ierr); 793 ierr = PCDestroy_MG(pc);CHKERRQ(ierr); 794 PetscFunctionReturn(0); 795 } 796 797 798 #undef __FUNCT__ 799 #define __FUNCT__ "PCGAMGSetProcEqLim" 800 /*@ 801 PCGAMGSetProcEqLim - Set number of equations to aim for on coarse grids via 802 processor reduction. 803 804 Not Collective on PC 805 806 Input Parameters: 807 . pc - the preconditioner context 808 809 810 Options Database Key: 811 . -pc_gamg_process_eq_limit 812 813 Level: intermediate 814 815 Concepts: Unstructured multrigrid preconditioner 816 817 .seealso: () 818 @*/ 819 PetscErrorCode PCGAMGSetProcEqLim(PC pc, PetscInt n) 820 { 821 PetscErrorCode ierr; 822 823 PetscFunctionBegin; 824 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 825 ierr = PetscTryMethod(pc,"PCGAMGSetProcEqLim_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr); 826 PetscFunctionReturn(0); 827 } 828 829 EXTERN_C_BEGIN 830 #undef __FUNCT__ 831 #define __FUNCT__ "PCGAMGSetProcEqLim_GAMG" 832 PetscErrorCode PCGAMGSetProcEqLim_GAMG(PC pc, PetscInt n) 833 { 834 PC_MG *mg = (PC_MG*)pc->data; 835 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 836 837 PetscFunctionBegin; 838 if(n>0) pc_gamg->m_min_eq_proc = n; 839 PetscFunctionReturn(0); 840 } 841 EXTERN_C_END 842 843 #undef __FUNCT__ 844 #define __FUNCT__ "PCGAMGAvoidRepartitioning" 845 /*@ 846 PCGAMGAvoidRepartitioning - Do not repartition the coarse grids 847 848 Collective on PC 849 850 Input Parameters: 851 . pc - the preconditioner context 852 853 854 Options Database Key: 855 . -pc_gamg_avoid_repartitioning 856 857 Level: intermediate 858 859 Concepts: Unstructured multrigrid preconditioner 860 861 .seealso: () 862 @*/ 863 PetscErrorCode PCGAMGAvoidRepartitioning(PC pc, PetscBool n) 864 { 865 PetscErrorCode ierr; 866 867 PetscFunctionBegin; 868 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 869 ierr = PetscTryMethod(pc,"PCGAMGAvoidRepartitioning_C",(PC,PetscBool),(pc,n));CHKERRQ(ierr); 870 PetscFunctionReturn(0); 871 } 872 873 EXTERN_C_BEGIN 874 #undef __FUNCT__ 875 #define __FUNCT__ "PCGAMGAvoidRepartitioning_GAMG" 876 PetscErrorCode PCGAMGAvoidRepartitioning_GAMG(PC pc, PetscBool n) 877 { 878 PC_MG *mg = (PC_MG*)pc->data; 879 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 880 881 PetscFunctionBegin; 882 pc_gamg->m_avoid_repart = n; 883 PetscFunctionReturn(0); 884 } 885 EXTERN_C_END 886 887 #undef __FUNCT__ 888 #define __FUNCT__ "PCGAMGSetThreshold" 889 /*@ 890 PCGAMGSetThreshold - Relative threshold to use for dropping edges in aggregation graph 891 892 Not collective on PC 893 894 Input Parameters: 895 . pc - the preconditioner context 896 897 898 Options Database Key: 899 . -pc_gamg_threshold 900 901 Level: intermediate 902 903 Concepts: Unstructured multrigrid preconditioner 904 905 .seealso: () 906 @*/ 907 PetscErrorCode PCGAMGSetThreshold(PC pc, PetscReal n) 908 { 909 PetscErrorCode ierr; 910 911 PetscFunctionBegin; 912 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 913 ierr = PetscTryMethod(pc,"PCGAMGSetThreshold_C",(PC,PetscReal),(pc,n));CHKERRQ(ierr); 914 PetscFunctionReturn(0); 915 } 916 917 EXTERN_C_BEGIN 918 #undef __FUNCT__ 919 #define __FUNCT__ "PCGAMGSetThreshold_GAMG" 920 PetscErrorCode PCGAMGSetThreshold_GAMG(PC pc, PetscReal n) 921 { 922 PC_MG *mg = (PC_MG*)pc->data; 923 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 924 925 PetscFunctionBegin; 926 pc_gamg->m_threshold = n; 927 PetscFunctionReturn(0); 928 } 929 EXTERN_C_END 930 931 #undef __FUNCT__ 932 #define __FUNCT__ "PCGAMGSetSolverType" 933 /*@ 934 PCGAMGSetSolverType - Set solution method. 935 936 Collective on PC 937 938 Input Parameters: 939 . pc - the preconditioner context 940 941 942 Options Database Key: 943 . -pc_gamg_type 944 945 Level: intermediate 946 947 Concepts: Unstructured multrigrid preconditioner 948 949 .seealso: () 950 @*/ 951 PetscErrorCode PCGAMGSetSolverType(PC pc, char str[], PetscInt sz ) 952 { 953 PetscErrorCode ierr; 954 955 PetscFunctionBegin; 956 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 957 ierr = PetscTryMethod(pc,"PCGAMGSetSolverType_C",(PC,char[],PetscInt),(pc,str,sz)); 958 CHKERRQ(ierr); 959 PetscFunctionReturn(0); 960 } 961 962 EXTERN_C_BEGIN 963 #undef __FUNCT__ 964 #define __FUNCT__ "PCGAMGSetSolverType_GAMG" 965 PetscErrorCode PCGAMGSetSolverType_GAMG(PC pc, char str[], PetscInt sz ) 966 { 967 PC_MG *mg = (PC_MG*)pc->data; 968 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 969 970 PetscFunctionBegin; 971 if(sz < 64) strcpy(pc_gamg->m_type,str); 972 PetscFunctionReturn(0); 973 } 974 EXTERN_C_END 975 976 #undef __FUNCT__ 977 #define __FUNCT__ "PCSetFromOptions_GAMG" 978 PetscErrorCode PCSetFromOptions_GAMG(PC pc) 979 { 980 PetscErrorCode ierr; 981 PC_MG *mg = (PC_MG*)pc->data; 982 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 983 PetscBool flag; 984 985 PetscFunctionBegin; 986 ierr = PetscOptionsHead("GAMG options"); CHKERRQ(ierr); 987 { 988 ierr = PetscOptionsString("-pc_gamg_type", 989 "Solver type: plane aggregation ('pa'), smoothed aggregation ('sa') or geometric multigrid (default)", 990 "PCGAMGSetSolverType", 991 pc_gamg->m_type, 992 pc_gamg->m_type, 993 64, 994 &flag ); 995 CHKERRQ(ierr); 996 if( flag && pc_gamg->m_data != 0 ) { 997 if( (strcmp(pc_gamg->m_type,"sa")==0 && pc_gamg->m_method != 2) || 998 (strcmp(pc_gamg->m_type,"pa")==0 && pc_gamg->m_method != 1) || 999 (strcmp(pc_gamg->m_type,"geo")==0 && pc_gamg->m_method != 0) ) { 1000 SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "PCSetFromOptions called of PCSetCoordinates (with new method, after data was created)"); 1001 } 1002 } 1003 1004 /* -pc_gamg_verbose */ 1005 ierr = PetscOptionsBool("-pc_gamg_verbose","Verbose (debugging) output for PCGAMG","none",pc_gamg->m_verbose,&pc_gamg->m_verbose,PETSC_NULL);CHKERRQ(ierr); 1006 1007 if (flag && strcmp(pc_gamg->m_type,"sa") == 0) pc_gamg->m_method = 2; 1008 else if (flag && strcmp(pc_gamg->m_type,"pa") == 0) pc_gamg->m_method = 1; 1009 else pc_gamg->m_method = 0; 1010 1011 /* -pc_gamg_avoid_repartitioning */ 1012 pc_gamg->m_avoid_repart = PETSC_FALSE; 1013 ierr = PetscOptionsBool("-pc_gamg_avoid_repartitioning", 1014 "Do not repartion coarse grids (false)", 1015 "PCGAMGAvoidRepartitioning", 1016 pc_gamg->m_avoid_repart, 1017 &pc_gamg->m_avoid_repart, 1018 &flag); 1019 CHKERRQ(ierr); 1020 1021 /* -pc_gamg_process_eq_limit */ 1022 pc_gamg->m_min_eq_proc = 600; 1023 ierr = PetscOptionsInt("-pc_gamg_process_eq_limit", 1024 "Limit (goal) on number of equations per process on coarse grids", 1025 "PCGAMGSetProcEqLim", 1026 pc_gamg->m_min_eq_proc, 1027 &pc_gamg->m_min_eq_proc, 1028 &flag ); 1029 CHKERRQ(ierr); 1030 1031 /* -pc_gamg_threshold */ 1032 pc_gamg->m_threshold = 0.05; 1033 ierr = PetscOptionsReal("-pc_gamg_threshold", 1034 "Relative threshold to use for dropping edges in aggregation graph", 1035 "PCGAMGSetThreshold", 1036 pc_gamg->m_threshold, 1037 &pc_gamg->m_threshold, 1038 &flag ); 1039 CHKERRQ(ierr); 1040 } 1041 ierr = PetscOptionsTail();CHKERRQ(ierr); 1042 1043 PetscFunctionReturn(0); 1044 } 1045 1046 /* -------------------------------------------------------------------------- */ 1047 /*MC 1048 PCGAMG - Geometric algebraic multigrid (AMG) preconditioning. This preconditioner currently has two 1049 AMG methods: 1) an unstructured geometric method, which requires that you provide coordinates for each 1050 vertex, and 2) smoothed aggregation. Smoothed aggregation (SA) can work without coordinates but it 1051 will generate some common non-trivial null spaces if coordinates are provided. The input fine grid matrix 1052 must have the block size set for 'system' problems (with multiple dof per vertex/cell) to work properly. 1053 SA will generate rotational rigid body mode null space vectors, in addition to the trivial translational 1054 modes, when coordinates are provide in 2D and 3D. 1055 1056 Options Database Keys: 1057 Multigrid options(inherited) 1058 + -pc_mg_cycles <1>: 1 for V cycle, 2 for W-cycle (PCMGSetCycleType) 1059 . -pc_mg_smoothup <1>: Number of post-smoothing steps (PCMGSetNumberSmoothUp) 1060 . -pc_mg_smoothdown <1>: Number of pre-smoothing steps (PCMGSetNumberSmoothDown) 1061 - -pc_mg_type <multiplicative>: (one of) additive multiplicative full cascade kascade 1062 1063 Level: intermediate 1064 1065 Concepts: multigrid 1066 1067 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, 1068 PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(), PCMGSetNumberSmoothDown(), 1069 PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(), 1070 PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(), 1071 PCMGSetCyclesOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR() 1072 M*/ 1073 1074 EXTERN_C_BEGIN 1075 #undef __FUNCT__ 1076 #define __FUNCT__ "PCCreate_GAMG" 1077 PetscErrorCode PCCreate_GAMG(PC pc) 1078 { 1079 PetscErrorCode ierr; 1080 PC_GAMG *pc_gamg; 1081 PC_MG *mg; 1082 PetscClassId cookie; 1083 #if defined PETSC_USE_LOG 1084 static int count = 0; 1085 #endif 1086 1087 PetscFunctionBegin; 1088 /* PCGAMG is an inherited class of PCMG. Initialize pc as PCMG */ 1089 ierr = PCSetType(pc,PCMG);CHKERRQ(ierr); /* calls PCCreate_MG() and MGCreate_Private() */ 1090 ierr = PetscObjectChangeTypeName((PetscObject)pc,PCGAMG);CHKERRQ(ierr); 1091 1092 /* create a supporting struct and attach it to pc */ 1093 ierr = PetscNewLog(pc,PC_GAMG,&pc_gamg);CHKERRQ(ierr); 1094 pc_gamg->m_data_sz = 0; pc_gamg->m_data = 0; pc_gamg->m_count = 0; 1095 mg = (PC_MG*)pc->data; 1096 mg->innerctx = pc_gamg; 1097 1098 pc_gamg->m_Nlevels = -1; 1099 1100 /* overwrite the pointers of PCMG by the functions of PCGAMG */ 1101 pc->ops->setfromoptions = PCSetFromOptions_GAMG; 1102 pc->ops->setup = PCSetUp_GAMG; 1103 pc->ops->reset = PCReset_GAMG; 1104 pc->ops->destroy = PCDestroy_GAMG; 1105 1106 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1107 "PCSetCoordinates_C", 1108 "PCSetCoordinates_GAMG", 1109 PCSetCoordinates_GAMG); 1110 CHKERRQ(ierr); 1111 1112 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1113 "PCGAMGSetProcEqLim_C", 1114 "PCGAMGSetProcEqLim_GAMG", 1115 PCGAMGSetProcEqLim_GAMG); 1116 CHKERRQ(ierr); 1117 1118 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1119 "PCGAMGAvoidRepartitioning_C", 1120 "PCGAMGAvoidRepartitioning_GAMG", 1121 PCGAMGAvoidRepartitioning_GAMG); 1122 CHKERRQ(ierr); 1123 1124 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1125 "PCGAMGSetThreshold_C", 1126 "PCGAMGSetThreshold_GAMG", 1127 PCGAMGSetThreshold_GAMG); 1128 CHKERRQ(ierr); 1129 1130 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1131 "PCGAMGSetSolverType_C", 1132 "PCGAMGSetSolverType_GAMG", 1133 PCGAMGSetSolverType_GAMG); 1134 CHKERRQ(ierr); 1135 1136 #if defined PETSC_USE_LOG 1137 if( count++ == 0 ) { 1138 PetscClassIdRegister("GAMG Setup",&cookie); 1139 PetscLogEventRegister("GAMG: createProl", cookie, &gamg_setup_events[SET1]); 1140 PetscLogEventRegister(" Graph", cookie, &gamg_setup_events[GRAPH]); 1141 PetscLogEventRegister(" G.Mat", cookie, &gamg_setup_events[GRAPH_MAT]); 1142 PetscLogEventRegister(" G.Filter", cookie, &gamg_setup_events[GRAPH_FILTER]); 1143 PetscLogEventRegister(" G.Square", cookie, &gamg_setup_events[GRAPH_SQR]); 1144 PetscLogEventRegister(" MIS/Agg", cookie, &gamg_setup_events[SET4]); 1145 PetscLogEventRegister(" geo: growSupp", cookie, &gamg_setup_events[SET5]); 1146 PetscLogEventRegister(" geo: triangle", cookie, &gamg_setup_events[SET6]); 1147 PetscLogEventRegister(" search&set", cookie, &gamg_setup_events[FIND_V]); 1148 PetscLogEventRegister(" SA: init", cookie, &gamg_setup_events[SET7]); 1149 /* PetscLogEventRegister(" SA: frmProl0", cookie, &gamg_setup_events[SET8]); */ 1150 PetscLogEventRegister(" SA: smooth", cookie, &gamg_setup_events[SET9]); 1151 PetscLogEventRegister("GAMG: partLevel", cookie, &gamg_setup_events[SET2]); 1152 PetscLogEventRegister(" PL repartition", cookie, &gamg_setup_events[SET12]); 1153 /* PetscLogEventRegister(" PL move data", cookie, &gamg_setup_events[SET13]); */ 1154 /* PetscLogEventRegister("GAMG: fix", cookie, &gamg_setup_events[SET10]); */ 1155 /* PetscLogEventRegister("GAMG: set levels", cookie, &gamg_setup_events[SET11]); */ 1156 1157 /* create timer stages */ 1158 #if defined GAMG_STAGES 1159 { 1160 char str[32]; 1161 sprintf(str,"MG Level %d (finest)",0); 1162 PetscLogStageRegister(str, &gamg_stages[0]); 1163 PetscInt lidx; 1164 for (lidx=1;lidx<9;lidx++){ 1165 sprintf(str,"MG Level %d",lidx); 1166 PetscLogStageRegister(str, &gamg_stages[lidx]); 1167 } 1168 } 1169 #endif 1170 } 1171 #endif 1172 PetscFunctionReturn(0); 1173 } 1174 EXTERN_C_END 1175