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