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