1 /* 2 GAMG geometric-algebric multiogrid 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 = neq/min_eq_proc; /* 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 /* aggregation method should gaurrentee this does not happen! */ 616 617 if (pc_gamg->m_verbose) { 618 Vec diag; PetscScalar *data_arr,v; PetscInt Istart,Iend,kk,nloceq,id; 619 v = 1.e-10; /* LU factor has hard wired numbers for small diags so this needs to match (yuk) */ 620 ierr = MatGetOwnershipRange(Aarr[level1], &Istart, &Iend); CHKERRQ(ierr); 621 nloceq = Iend-Istart; 622 ierr = MatGetVecs( Aarr[level1], &diag, 0 ); CHKERRQ(ierr); 623 ierr = MatGetDiagonal( Aarr[level1], diag ); CHKERRQ(ierr); 624 ierr = VecGetArray( diag, &data_arr ); CHKERRQ(ierr); 625 for(kk=0;kk<nloceq;kk++){ 626 if(data_arr[kk]==0.0) { 627 id = kk + Istart; 628 ierr = MatSetValues(Aarr[level1],1,&id,1,&id,&v,INSERT_VALUES); 629 CHKERRQ(ierr); 630 PetscPrintf(PETSC_COMM_SELF,"\t[%d]%s warning: added zero to diag (%d) on level %d \n",mype,__FUNCT__,id,level1); 631 } 632 } 633 ierr = VecRestoreArray( diag, &data_arr ); CHKERRQ(ierr); 634 ierr = VecDestroy( &diag ); CHKERRQ(ierr); 635 ierr = MatAssemblyBegin(Aarr[level1],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 636 ierr = MatAssemblyEnd(Aarr[level1],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 637 } 638 } else { 639 coarse_data = 0; 640 break; 641 } 642 data = coarse_data; 643 644 #if (defined PETSC_USE_LOG && defined GAMG_STAGES) 645 ierr = PetscLogStagePop(); CHKERRQ( ierr ); 646 #endif 647 } 648 if( coarse_data ) { 649 ierr = PetscFree( coarse_data ); CHKERRQ( ierr ); 650 } 651 if (pc_gamg->m_verbose) PetscPrintf(PETSC_COMM_WORLD,"\t[%d]%s %d levels\n",0,__FUNCT__,level + 1); 652 pc_gamg->m_data = 0; /* destroyed coordinate data */ 653 pc_gamg->m_Nlevels = level + 1; 654 fine_level = level; 655 ierr = PCMGSetLevels(a_pc,pc_gamg->m_Nlevels,PETSC_NULL);CHKERRQ(ierr); 656 657 /* set default smoothers */ 658 for ( lidx=1, level = pc_gamg->m_Nlevels-2; 659 lidx <= fine_level; 660 lidx++, level--) { 661 PetscReal emax, emin; 662 KSP smoother; PC subpc; 663 PetscBool isCheb; 664 /* set defaults */ 665 ierr = PCMGGetSmoother( a_pc, lidx, &smoother ); CHKERRQ(ierr); 666 ierr = KSPSetType( smoother, KSPCHEBYCHEV );CHKERRQ(ierr); 667 ierr = KSPGetPC( smoother, &subpc ); CHKERRQ(ierr); 668 /* ierr = KSPSetTolerances(smoother,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,2); CHKERRQ(ierr); */ 669 ierr = PCSetType( subpc, PETSC_GAMG_SMOOTHER ); CHKERRQ(ierr); 670 ierr = KSPSetNormType( smoother, KSP_NORM_NONE ); CHKERRQ(ierr); 671 /* overide defaults with input parameters */ 672 ierr = KSPSetFromOptions( smoother ); CHKERRQ(ierr); 673 674 ierr = KSPSetOperators( smoother, Aarr[level], Aarr[level], SAME_NONZERO_PATTERN ); CHKERRQ(ierr); 675 /* do my own cheby */ 676 ierr = PetscTypeCompare( (PetscObject)smoother, KSPCHEBYCHEV, &isCheb ); CHKERRQ(ierr); 677 678 if( isCheb ) { 679 ierr = PetscTypeCompare( (PetscObject)subpc, PETSC_GAMG_SMOOTHER, &isCheb ); CHKERRQ(ierr); 680 if( isCheb && emaxs[level] > 0.0 ) emax=emaxs[level]; /* eigen estimate only for diagnal PC */ 681 else{ /* eigen estimate 'emax' */ 682 KSP eksp; Mat Lmat = Aarr[level]; 683 Vec bb, xx; PC pc; 684 const PCType type; 685 686 ierr = PCGetType( subpc, &type ); CHKERRQ(ierr); 687 ierr = MatGetVecs( Lmat, &bb, 0 ); CHKERRQ(ierr); 688 ierr = MatGetVecs( Lmat, &xx, 0 ); CHKERRQ(ierr); 689 { 690 PetscRandom rctx; 691 ierr = PetscRandomCreate(wcomm,&rctx);CHKERRQ(ierr); 692 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 693 ierr = VecSetRandom(bb,rctx);CHKERRQ(ierr); 694 ierr = PetscRandomDestroy( &rctx ); CHKERRQ(ierr); 695 } 696 ierr = KSPCreate(wcomm,&eksp);CHKERRQ(ierr); 697 ierr = KSPSetType( eksp, KSPCG ); CHKERRQ(ierr); 698 ierr = KSPSetTolerances( eksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 10 ); 699 CHKERRQ(ierr); 700 ierr = KSPSetNormType( eksp, KSP_NORM_NONE ); CHKERRQ(ierr); 701 702 ierr = KSPAppendOptionsPrefix( eksp, "est_"); CHKERRQ(ierr); 703 ierr = KSPSetFromOptions( eksp ); CHKERRQ(ierr); 704 705 ierr = KSPSetInitialGuessNonzero( eksp, PETSC_FALSE ); CHKERRQ(ierr); 706 ierr = KSPSetOperators( eksp, Lmat, Lmat, SAME_NONZERO_PATTERN ); CHKERRQ( ierr ); 707 ierr = KSPGetPC( eksp, &pc );CHKERRQ( ierr ); 708 ierr = PCSetType( pc, type ); CHKERRQ(ierr); /* should be same as eigen estimates op. */ 709 710 ierr = KSPSetComputeSingularValues( eksp,PETSC_TRUE ); CHKERRQ(ierr); 711 ierr = KSPSolve( eksp, bb, xx ); CHKERRQ(ierr); 712 ierr = KSPComputeExtremeSingularValues( eksp, &emax, &emin ); CHKERRQ(ierr); 713 ierr = VecDestroy( &xx ); CHKERRQ(ierr); 714 ierr = VecDestroy( &bb ); CHKERRQ(ierr); 715 ierr = KSPDestroy( &eksp ); CHKERRQ(ierr); 716 717 if (pc_gamg->m_verbose) { 718 PetscPrintf(PETSC_COMM_WORLD,"\t\t\t%s PC setup max eigen=%e min=%e\n", 719 __FUNCT__,emax,emin); 720 } 721 } 722 { 723 PetscInt N1, N0, tt; 724 ierr = MatGetSize( Aarr[level], &N1, &tt ); CHKERRQ(ierr); 725 ierr = MatGetSize( Aarr[level+1], &N0, &tt ); CHKERRQ(ierr); 726 /* heuristic - is this crap? */ 727 emin = 1.*emax/((PetscReal)N1/(PetscReal)N0); 728 emax *= 1.05; 729 } 730 ierr = KSPChebychevSetEigenvalues( smoother, emax, emin );CHKERRQ(ierr); 731 } 732 } 733 { 734 /* coarse grid */ 735 KSP smoother,*k2; PC subpc,pc2; PetscInt ii,first; 736 Mat Lmat = Aarr[pc_gamg->m_Nlevels-1]; 737 ierr = PCMGGetSmoother( a_pc, 0, &smoother ); CHKERRQ(ierr); 738 ierr = KSPSetOperators( smoother, Lmat, Lmat, SAME_NONZERO_PATTERN ); CHKERRQ(ierr); 739 ierr = KSPSetNormType( smoother, KSP_NORM_NONE ); CHKERRQ(ierr); 740 ierr = KSPGetPC( smoother, &subpc ); CHKERRQ(ierr); 741 ierr = PCSetType( subpc, PCBJACOBI ); CHKERRQ(ierr); 742 ierr = PCSetUp( subpc ); CHKERRQ(ierr); 743 ierr = PCBJacobiGetSubKSP(subpc,&ii,&first,&k2);CHKERRQ(ierr); 744 assert(ii==1); 745 ierr = KSPGetPC(k2[0],&pc2);CHKERRQ(ierr); 746 ierr = PCSetType( pc2, PCLU ); CHKERRQ(ierr); 747 } 748 749 /* should be called in PCSetFromOptions_GAMG(), but cannot be called prior to PCMGSetLevels() */ 750 ierr = PCSetFromOptions_MG(a_pc); CHKERRQ(ierr); 751 { 752 PetscBool galerkin; 753 ierr = PCMGGetGalerkin( a_pc, &galerkin); CHKERRQ(ierr); 754 if(galerkin){ 755 SETERRQ(wcomm,PETSC_ERR_ARG_WRONG, "GAMG does galerkin manually so it must not be used in PC_MG."); 756 } 757 } 758 759 /* plot levels - R/P */ 760 if( PRINT_MATS ) { 761 for (level=pc_gamg->m_Nlevels-1;level>0;level--){ 762 PetscViewer viewer; 763 char fname[32]; 764 sprintf(fname,"Pmat_%d_%d.m",(int)pc_gamg->m_count,(int)level); 765 ierr = PetscViewerASCIIOpen( wcomm, fname, &viewer ); CHKERRQ(ierr); 766 ierr = PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB); CHKERRQ(ierr); 767 ierr = MatView( Parr[level], viewer ); CHKERRQ(ierr); 768 ierr = PetscViewerDestroy( &viewer ); 769 sprintf(fname,"Amat_%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( Aarr[level], viewer ); CHKERRQ(ierr); 773 ierr = PetscViewerDestroy( &viewer ); 774 } 775 } 776 777 /* set interpolation between the levels, clean up */ 778 for (lidx=0,level=pc_gamg->m_Nlevels-1; 779 lidx<fine_level; 780 lidx++, level--){ 781 ierr = PCMGSetInterpolation( a_pc, lidx+1, Parr[level] );CHKERRQ(ierr); 782 ierr = MatDestroy( &Parr[level] ); CHKERRQ(ierr); 783 ierr = MatDestroy( &Aarr[level] ); CHKERRQ(ierr); 784 } 785 786 /* setupcalled is set to 0 so that MG is setup from scratch */ 787 a_pc->setupcalled = 0; 788 ierr = PCSetUp_MG( a_pc );CHKERRQ( ierr ); 789 a_pc->setupcalled = 1; /* use 1 as signal that this has not been re-setup */ 790 791 { 792 KSP smoother; /* PCSetUp_MG seems to insists on setting this to GMRES on coarse grid */ 793 ierr = PCMGGetSmoother( a_pc, 0, &smoother ); CHKERRQ(ierr); 794 ierr = KSPSetType( smoother, KSPPREONLY ); CHKERRQ(ierr); 795 ierr = KSPSetUp( smoother ); CHKERRQ(ierr); 796 } 797 798 PetscFunctionReturn(0); 799 } 800 801 /* ------------------------------------------------------------------------- */ 802 /* 803 PCDestroy_GAMG - Destroys the private context for the GAMG preconditioner 804 that was created with PCCreate_GAMG(). 805 806 Input Parameter: 807 . pc - the preconditioner context 808 809 Application Interface Routine: PCDestroy() 810 */ 811 #undef __FUNCT__ 812 #define __FUNCT__ "PCDestroy_GAMG" 813 PetscErrorCode PCDestroy_GAMG(PC pc) 814 { 815 PetscErrorCode ierr; 816 PC_MG *mg = (PC_MG*)pc->data; 817 PC_GAMG *pc_gamg= (PC_GAMG*)mg->innerctx; 818 819 PetscFunctionBegin; 820 ierr = PCReset_GAMG(pc);CHKERRQ(ierr); 821 ierr = PetscFree(pc_gamg);CHKERRQ(ierr); 822 ierr = PCDestroy_MG(pc);CHKERRQ(ierr); 823 PetscFunctionReturn(0); 824 } 825 826 827 #undef __FUNCT__ 828 #define __FUNCT__ "PCGAMGSetProcEqLim" 829 /*@ 830 PCGAMGSetProcEqLim - Set number of equations to aim for on coarse grids via 831 processor reduction. 832 833 Not Collective on PC 834 835 Input Parameters: 836 . pc - the preconditioner context 837 838 839 Options Database Key: 840 . -pc_gamg_process_eq_limit 841 842 Level: intermediate 843 844 Concepts: Unstructured multrigrid preconditioner 845 846 .seealso: () 847 @*/ 848 PetscErrorCode PCGAMGSetProcEqLim(PC pc, PetscInt n) 849 { 850 PetscErrorCode ierr; 851 852 PetscFunctionBegin; 853 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 854 ierr = PetscTryMethod(pc,"PCGAMGSetProcEqLim_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr); 855 PetscFunctionReturn(0); 856 } 857 858 EXTERN_C_BEGIN 859 #undef __FUNCT__ 860 #define __FUNCT__ "PCGAMGSetProcEqLim_GAMG" 861 PetscErrorCode PCGAMGSetProcEqLim_GAMG(PC pc, PetscInt n) 862 { 863 PC_MG *mg = (PC_MG*)pc->data; 864 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 865 866 PetscFunctionBegin; 867 if(n>0) pc_gamg->m_min_eq_proc = n; 868 PetscFunctionReturn(0); 869 } 870 EXTERN_C_END 871 872 #undef __FUNCT__ 873 #define __FUNCT__ "PCGAMGSetCoarseEqLim" 874 /*@ 875 PCGAMGSetCoarseEqLim - Set max number of equations on coarse grids. 876 877 Collective on PC 878 879 Input Parameters: 880 . pc - the preconditioner context 881 882 883 Options Database Key: 884 . -pc_gamg_coarse_eq_limit 885 886 Level: intermediate 887 888 Concepts: Unstructured multrigrid preconditioner 889 890 .seealso: () 891 @*/ 892 PetscErrorCode PCGAMGSetCoarseEqLim(PC pc, PetscInt n) 893 { 894 PetscErrorCode ierr; 895 896 PetscFunctionBegin; 897 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 898 ierr = PetscTryMethod(pc,"PCGAMGSetCoarseEqLim_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr); 899 PetscFunctionReturn(0); 900 } 901 902 EXTERN_C_BEGIN 903 #undef __FUNCT__ 904 #define __FUNCT__ "PCGAMGSetCoarseEqLim_GAMG" 905 PetscErrorCode PCGAMGSetCoarseEqLim_GAMG(PC pc, PetscInt n) 906 { 907 PC_MG *mg = (PC_MG*)pc->data; 908 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 909 910 PetscFunctionBegin; 911 if(n>0) pc_gamg->m_coarse_eq_limit = n; 912 PetscFunctionReturn(0); 913 } 914 EXTERN_C_END 915 916 #undef __FUNCT__ 917 #define __FUNCT__ "PCGAMGSetRepartitioning" 918 /*@ 919 PCGAMGSetRepartitioning - Repartition the coarse grids 920 921 Collective on PC 922 923 Input Parameters: 924 . pc - the preconditioner context 925 926 927 Options Database Key: 928 . -pc_gamg_repartition 929 930 Level: intermediate 931 932 Concepts: Unstructured multrigrid preconditioner 933 934 .seealso: () 935 @*/ 936 PetscErrorCode PCGAMGSetRepartitioning(PC pc, PetscBool n) 937 { 938 PetscErrorCode ierr; 939 940 PetscFunctionBegin; 941 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 942 ierr = PetscTryMethod(pc,"PCGAMGSetRepartitioning_C",(PC,PetscBool),(pc,n));CHKERRQ(ierr); 943 PetscFunctionReturn(0); 944 } 945 946 EXTERN_C_BEGIN 947 #undef __FUNCT__ 948 #define __FUNCT__ "PCGAMGSetRepartitioning_GAMG" 949 PetscErrorCode PCGAMGSetRepartitioning_GAMG(PC pc, PetscBool n) 950 { 951 PC_MG *mg = (PC_MG*)pc->data; 952 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 953 954 PetscFunctionBegin; 955 pc_gamg->m_repart = n; 956 PetscFunctionReturn(0); 957 } 958 EXTERN_C_END 959 960 #undef __FUNCT__ 961 #define __FUNCT__ "PCGAMGSetNlevels" 962 /*@ 963 PCGAMGSetNlevels - 964 965 Not collective on PC 966 967 Input Parameters: 968 . pc - the preconditioner context 969 970 971 Options Database Key: 972 . -pc_mg_levels 973 974 Level: intermediate 975 976 Concepts: Unstructured multrigrid preconditioner 977 978 .seealso: () 979 @*/ 980 PetscErrorCode PCGAMGSetNlevels(PC pc, PetscInt n) 981 { 982 PetscErrorCode ierr; 983 984 PetscFunctionBegin; 985 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 986 ierr = PetscTryMethod(pc,"PCGAMGSetNlevels_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr); 987 PetscFunctionReturn(0); 988 } 989 990 EXTERN_C_BEGIN 991 #undef __FUNCT__ 992 #define __FUNCT__ "PCGAMGSetNlevels_GAMG" 993 PetscErrorCode PCGAMGSetNlevels_GAMG(PC pc, PetscInt n) 994 { 995 PC_MG *mg = (PC_MG*)pc->data; 996 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 997 998 PetscFunctionBegin; 999 pc_gamg->m_Nlevels = n; 1000 PetscFunctionReturn(0); 1001 } 1002 EXTERN_C_END 1003 1004 #undef __FUNCT__ 1005 #define __FUNCT__ "PCGAMGSetThreshold" 1006 /*@ 1007 PCGAMGSetThreshold - Relative threshold to use for dropping edges in aggregation graph 1008 1009 Not collective on PC 1010 1011 Input Parameters: 1012 . pc - the preconditioner context 1013 1014 1015 Options Database Key: 1016 . -pc_gamg_threshold 1017 1018 Level: intermediate 1019 1020 Concepts: Unstructured multrigrid preconditioner 1021 1022 .seealso: () 1023 @*/ 1024 PetscErrorCode PCGAMGSetThreshold(PC pc, PetscReal n) 1025 { 1026 PetscErrorCode ierr; 1027 1028 PetscFunctionBegin; 1029 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1030 ierr = PetscTryMethod(pc,"PCGAMGSetThreshold_C",(PC,PetscReal),(pc,n));CHKERRQ(ierr); 1031 PetscFunctionReturn(0); 1032 } 1033 1034 EXTERN_C_BEGIN 1035 #undef __FUNCT__ 1036 #define __FUNCT__ "PCGAMGSetThreshold_GAMG" 1037 PetscErrorCode PCGAMGSetThreshold_GAMG(PC pc, PetscReal n) 1038 { 1039 PC_MG *mg = (PC_MG*)pc->data; 1040 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1041 1042 PetscFunctionBegin; 1043 pc_gamg->m_threshold = n; 1044 PetscFunctionReturn(0); 1045 } 1046 EXTERN_C_END 1047 1048 #undef __FUNCT__ 1049 #define __FUNCT__ "PCGAMGSetSolverType" 1050 /*@ 1051 PCGAMGSetSolverType - Set solution method. 1052 1053 Collective on PC 1054 1055 Input Parameters: 1056 . pc - the preconditioner context 1057 1058 1059 Options Database Key: 1060 . -pc_gamg_type 1061 1062 Level: intermediate 1063 1064 Concepts: Unstructured multrigrid preconditioner 1065 1066 .seealso: () 1067 @*/ 1068 PetscErrorCode PCGAMGSetSolverType(PC pc, char str[], PetscInt sz ) 1069 { 1070 PetscErrorCode ierr; 1071 1072 PetscFunctionBegin; 1073 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1074 ierr = PetscTryMethod(pc,"PCGAMGSetSolverType_C",(PC,char[],PetscInt),(pc,str,sz)); 1075 CHKERRQ(ierr); 1076 PetscFunctionReturn(0); 1077 } 1078 1079 EXTERN_C_BEGIN 1080 #undef __FUNCT__ 1081 #define __FUNCT__ "PCGAMGSetSolverType_GAMG" 1082 PetscErrorCode PCGAMGSetSolverType_GAMG(PC pc, char str[], PetscInt sz ) 1083 { 1084 PC_MG *mg = (PC_MG*)pc->data; 1085 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1086 1087 PetscFunctionBegin; 1088 if(sz < 64) strcpy(pc_gamg->m_type,str); 1089 PetscFunctionReturn(0); 1090 } 1091 EXTERN_C_END 1092 1093 #undef __FUNCT__ 1094 #define __FUNCT__ "PCSetFromOptions_GAMG" 1095 PetscErrorCode PCSetFromOptions_GAMG(PC pc) 1096 { 1097 PetscErrorCode ierr; 1098 PC_MG *mg = (PC_MG*)pc->data; 1099 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1100 PetscBool flag; 1101 1102 PetscFunctionBegin; 1103 ierr = PetscOptionsHead("GAMG options"); CHKERRQ(ierr); 1104 { 1105 ierr = PetscOptionsString("-pc_gamg_type", 1106 "Solver type: plane aggregation ('pa'), smoothed aggregation ('sa') or geometric multigrid (default)", 1107 "PCGAMGSetSolverType", 1108 pc_gamg->m_type, 1109 pc_gamg->m_type, 1110 64, 1111 &flag ); 1112 CHKERRQ(ierr); 1113 if( flag && pc_gamg->m_data != 0 ) { 1114 if( (strcmp(pc_gamg->m_type,"sa")==0 && pc_gamg->m_method != 2) || 1115 (strcmp(pc_gamg->m_type,"pa")==0 && pc_gamg->m_method != 1) || 1116 (strcmp(pc_gamg->m_type,"geo")==0 && pc_gamg->m_method != 0) ) { 1117 SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "PCSetFromOptions called of PCSetCoordinates (with new method, after data was created)"); 1118 } 1119 } 1120 1121 /* -pc_gamg_verbose */ 1122 ierr = PetscOptionsBool("-pc_gamg_verbose","Verbose (debugging) output for PCGAMG","none",pc_gamg->m_verbose,&pc_gamg->m_verbose,PETSC_NULL);CHKERRQ(ierr); 1123 1124 pc_gamg->m_method = 1; /* default to plane aggregation */ 1125 if (flag ) { 1126 if( strcmp(pc_gamg->m_type,"sa") == 0) pc_gamg->m_method = 2; 1127 else if( strcmp(pc_gamg->m_type,"pa") == 0) pc_gamg->m_method = 1; 1128 else if( strcmp(pc_gamg->m_type,"geo") == 0) pc_gamg->m_method = 0; 1129 else SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "Invalid gamg type: %s",pc_gamg->m_type); 1130 } 1131 /* -pc_gamg_repartition */ 1132 pc_gamg->m_repart = PETSC_FALSE; 1133 ierr = PetscOptionsBool("-pc_gamg_repartition", 1134 "Repartion coarse grids (false)", 1135 "PCGAMGRepartitioning", 1136 pc_gamg->m_repart, 1137 &pc_gamg->m_repart, 1138 &flag); 1139 CHKERRQ(ierr); 1140 1141 /* -pc_gamg_process_eq_limit */ 1142 pc_gamg->m_min_eq_proc = 600; 1143 ierr = PetscOptionsInt("-pc_gamg_process_eq_limit", 1144 "Limit (goal) on number of equations per process on coarse grids", 1145 "PCGAMGSetProcEqLim", 1146 pc_gamg->m_min_eq_proc, 1147 &pc_gamg->m_min_eq_proc, 1148 &flag ); 1149 CHKERRQ(ierr); 1150 1151 /* -pc_gamg_coarse_eq_limit */ 1152 pc_gamg->m_coarse_eq_limit = 1500; 1153 ierr = PetscOptionsInt("-pc_gamg_coarse_eq_limit", 1154 "Limit on number of equations for the coarse grid", 1155 "PCGAMGSetCoarseEqLim", 1156 pc_gamg->m_coarse_eq_limit, 1157 &pc_gamg->m_coarse_eq_limit, 1158 &flag ); 1159 CHKERRQ(ierr); 1160 1161 /* -pc_gamg_threshold */ 1162 pc_gamg->m_threshold = 0.05; 1163 ierr = PetscOptionsReal("-pc_gamg_threshold", 1164 "Relative threshold to use for dropping edges in aggregation graph", 1165 "PCGAMGSetThreshold", 1166 pc_gamg->m_threshold, 1167 &pc_gamg->m_threshold, 1168 &flag ); 1169 CHKERRQ(ierr); 1170 1171 pc_gamg->m_Nlevels = GAMG_MAXLEVELS; 1172 ierr = PetscOptionsInt("-pc_mg_levels", 1173 "Set number of MG levels", 1174 "PCGAMGSetNlevels", 1175 pc_gamg->m_Nlevels, 1176 &pc_gamg->m_Nlevels, 1177 &flag ); 1178 } 1179 ierr = PetscOptionsTail();CHKERRQ(ierr); 1180 1181 PetscFunctionReturn(0); 1182 } 1183 1184 /* -------------------------------------------------------------------------- */ 1185 /*MC 1186 PCGAMG - Geometric algebraic multigrid (AMG) preconditioning. This preconditioner currently has two 1187 AMG methods: 1) an unstructured geometric method, which requires that you provide coordinates for each 1188 vertex, and 2) smoothed aggregation. Smoothed aggregation (SA) can work without coordinates but it 1189 will generate some common non-trivial null spaces if coordinates are provided. The input fine grid matrix 1190 must have the block size set for 'system' problems (with multiple dof per vertex/cell) to work properly. 1191 SA will generate rotational rigid body mode null space vectors, in addition to the trivial translational 1192 modes, when coordinates are provide in 2D and 3D. 1193 1194 Options Database Keys: 1195 Multigrid options(inherited) 1196 + -pc_mg_cycles <1>: 1 for V cycle, 2 for W-cycle (PCMGSetCycleType) 1197 . -pc_mg_smoothup <1>: Number of post-smoothing steps (PCMGSetNumberSmoothUp) 1198 . -pc_mg_smoothdown <1>: Number of pre-smoothing steps (PCMGSetNumberSmoothDown) 1199 - -pc_mg_type <multiplicative>: (one of) additive multiplicative full cascade kascade 1200 1201 Level: intermediate 1202 1203 Concepts: multigrid 1204 1205 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, 1206 PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(), PCMGSetNumberSmoothDown(), 1207 PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(), 1208 PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(), 1209 PCMGSetCyclesOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR() 1210 M*/ 1211 1212 EXTERN_C_BEGIN 1213 #undef __FUNCT__ 1214 #define __FUNCT__ "PCCreate_GAMG" 1215 PetscErrorCode PCCreate_GAMG(PC pc) 1216 { 1217 PetscErrorCode ierr; 1218 PC_GAMG *pc_gamg; 1219 PC_MG *mg; 1220 PetscClassId cookie; 1221 #if defined PETSC_USE_LOG 1222 static int count = 0; 1223 #endif 1224 1225 PetscFunctionBegin; 1226 /* PCGAMG is an inherited class of PCMG. Initialize pc as PCMG */ 1227 ierr = PCSetType(pc,PCMG);CHKERRQ(ierr); /* calls PCCreate_MG() and MGCreate_Private() */ 1228 ierr = PetscObjectChangeTypeName((PetscObject)pc,PCGAMG);CHKERRQ(ierr); 1229 1230 /* create a supporting struct and attach it to pc */ 1231 ierr = PetscNewLog(pc,PC_GAMG,&pc_gamg);CHKERRQ(ierr); 1232 pc_gamg->m_data_sz = 0; pc_gamg->m_data = 0; pc_gamg->m_count = 0; 1233 mg = (PC_MG*)pc->data; 1234 mg->innerctx = pc_gamg; 1235 1236 pc_gamg->m_Nlevels = -1; 1237 1238 /* overwrite the pointers of PCMG by the functions of PCGAMG */ 1239 pc->ops->setfromoptions = PCSetFromOptions_GAMG; 1240 pc->ops->setup = PCSetUp_GAMG; 1241 pc->ops->reset = PCReset_GAMG; 1242 pc->ops->destroy = PCDestroy_GAMG; 1243 1244 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1245 "PCSetCoordinates_C", 1246 "PCSetCoordinates_GAMG", 1247 PCSetCoordinates_GAMG); 1248 CHKERRQ(ierr); 1249 1250 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1251 "PCGAMGSetProcEqLim_C", 1252 "PCGAMGSetProcEqLim_GAMG", 1253 PCGAMGSetProcEqLim_GAMG); 1254 CHKERRQ(ierr); 1255 1256 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1257 "PCGAMGSetCoarseEqLim_C", 1258 "PCGAMGSetCoarseEqLim_GAMG", 1259 PCGAMGSetCoarseEqLim_GAMG); 1260 CHKERRQ(ierr); 1261 1262 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1263 "PCGAMGSetRepartitioning_C", 1264 "PCGAMGSetRepartitioning_GAMG", 1265 PCGAMGSetRepartitioning_GAMG); 1266 CHKERRQ(ierr); 1267 1268 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1269 "PCGAMGSetThreshold_C", 1270 "PCGAMGSetThreshold_GAMG", 1271 PCGAMGSetThreshold_GAMG); 1272 CHKERRQ(ierr); 1273 1274 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1275 "PCGAMGSetSolverType_C", 1276 "PCGAMGSetSolverType_GAMG", 1277 PCGAMGSetSolverType_GAMG); 1278 CHKERRQ(ierr); 1279 1280 #if defined PETSC_USE_LOG 1281 if( count++ == 0 ) { 1282 PetscClassIdRegister("GAMG Setup",&cookie); 1283 PetscLogEventRegister("GAMG: createProl", cookie, &gamg_setup_events[SET1]); 1284 PetscLogEventRegister(" Graph", cookie, &gamg_setup_events[GRAPH]); 1285 PetscLogEventRegister(" G.Mat", cookie, &gamg_setup_events[GRAPH_MAT]); 1286 PetscLogEventRegister(" G.Filter", cookie, &gamg_setup_events[GRAPH_FILTER]); 1287 PetscLogEventRegister(" G.Square", cookie, &gamg_setup_events[GRAPH_SQR]); 1288 PetscLogEventRegister(" MIS/Agg", cookie, &gamg_setup_events[SET4]); 1289 PetscLogEventRegister(" geo: growSupp", cookie, &gamg_setup_events[SET5]); 1290 PetscLogEventRegister(" geo: triangle", cookie, &gamg_setup_events[SET6]); 1291 PetscLogEventRegister(" search&set", cookie, &gamg_setup_events[FIND_V]); 1292 PetscLogEventRegister(" SA: init", cookie, &gamg_setup_events[SET7]); 1293 /* PetscLogEventRegister(" SA: frmProl0", cookie, &gamg_setup_events[SET8]); */ 1294 PetscLogEventRegister(" SA: smooth", cookie, &gamg_setup_events[SET9]); 1295 PetscLogEventRegister("GAMG: partLevel", cookie, &gamg_setup_events[SET2]); 1296 PetscLogEventRegister(" PL repartition", cookie, &gamg_setup_events[SET12]); 1297 1298 /* PetscLogEventRegister(" PL 1", cookie, &gamg_setup_events[SET13]); */ 1299 /* PetscLogEventRegister(" PL 2", cookie, &gamg_setup_events[SET14]); */ 1300 /* PetscLogEventRegister(" PL 3", cookie, &gamg_setup_events[SET15]); */ 1301 1302 /* PetscLogEventRegister(" PL move data", cookie, &gamg_setup_events[SET13]); */ 1303 /* PetscLogEventRegister("GAMG: fix", cookie, &gamg_setup_events[SET10]); */ 1304 /* PetscLogEventRegister("GAMG: set levels", cookie, &gamg_setup_events[SET11]); */ 1305 1306 /* create timer stages */ 1307 #if defined GAMG_STAGES 1308 { 1309 char str[32]; 1310 sprintf(str,"MG Level %d (finest)",0); 1311 PetscLogStageRegister(str, &gamg_stages[0]); 1312 PetscInt lidx; 1313 for (lidx=1;lidx<9;lidx++){ 1314 sprintf(str,"MG Level %d",lidx); 1315 PetscLogStageRegister(str, &gamg_stages[lidx]); 1316 } 1317 } 1318 #endif 1319 } 1320 #endif 1321 PetscFunctionReturn(0); 1322 } 1323 EXTERN_C_END 1324