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