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