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