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