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