1 /* 2 GAMG geometric-algebric multiogrid PC - Mark Adams 2011 3 */ 4 #include <../src/ksp/pc/impls/gamg/gamg.h> /*I "petscpc.h" I*/ 5 #include "private/matimpl.h" 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 /* Private context for the GAMG preconditioner */ 18 static PetscBool s_avoid_repart = PETSC_FALSE; 19 static PetscInt s_min_eq_proc = 600; 20 static PetscReal s_threshold = 0.05; 21 typedef struct gamg_TAG { 22 PetscInt m_dim; 23 PetscInt m_Nlevels; 24 PetscInt m_data_sz; 25 PetscInt m_data_rows; 26 PetscInt m_data_cols; 27 PetscInt m_count; 28 PetscInt m_method; /* 0: geomg; 1: plain agg 'pa'; 2: smoothed agg 'sa' */ 29 PetscReal *m_data; /* blocked vector of vertex data on fine grid (coordinates) */ 30 char m_type[64]; 31 } PC_GAMG; 32 33 /* -------------------------------------------------------------------------- */ 34 /* 35 PCSetCoordinates_GAMG 36 37 Input Parameter: 38 . pc - the preconditioner context 39 */ 40 EXTERN_C_BEGIN 41 #undef __FUNCT__ 42 #define __FUNCT__ "PCSetCoordinates_GAMG" 43 PetscErrorCode PCSetCoordinates_GAMG( PC a_pc, PetscInt a_ndm, PetscReal *a_coords ) 44 { 45 PC_MG *mg = (PC_MG*)a_pc->data; 46 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 47 PetscErrorCode ierr; 48 PetscInt arrsz,bs,my0,kk,ii,jj,nloc,Iend; 49 Mat Amat = a_pc->pmat; 50 51 PetscFunctionBegin; 52 PetscValidHeaderSpecific( Amat, MAT_CLASSID, 1 ); 53 ierr = MatGetBlockSize( Amat, &bs ); CHKERRQ( ierr ); 54 ierr = MatGetOwnershipRange( Amat, &my0, &Iend ); CHKERRQ(ierr); 55 nloc = (Iend-my0)/bs; 56 if((Iend-my0)%bs!=0) SETERRQ1(((PetscObject)Amat)->comm,PETSC_ERR_ARG_WRONG, "Bad local size %d.",nloc); 57 58 pc_gamg->m_data_rows = 1; 59 if(a_coords==0 && pc_gamg->m_method==0) pc_gamg->m_method = 2; /* use SA if no coords */ 60 if( pc_gamg->m_method==0 ) pc_gamg->m_data_cols = a_ndm; /* coordinates */ 61 else{ /* SA: null space vectors */ 62 if(a_coords != 0 && bs==1 ) pc_gamg->m_data_cols = 1; /* scalar w/ coords and SA (not needed) */ 63 else if(a_coords != 0 ) pc_gamg->m_data_cols = (a_ndm==2 ? 3 : 6); /* elasticity */ 64 else pc_gamg->m_data_cols = bs; /* no data, force SA with constant null space vectors */ 65 pc_gamg->m_data_rows = bs; 66 } 67 arrsz = nloc*pc_gamg->m_data_rows*pc_gamg->m_data_cols; 68 69 /* create data - syntactic sugar that should be refactored at some point */ 70 if (pc_gamg->m_data==0 || (pc_gamg->m_data_sz != arrsz)) { 71 ierr = PetscFree( pc_gamg->m_data ); CHKERRQ(ierr); 72 ierr = PetscMalloc((arrsz+1)*sizeof(double), &pc_gamg->m_data ); CHKERRQ(ierr); 73 } 74 for(kk=0;kk<arrsz;kk++)pc_gamg->m_data[kk] = -999.; 75 pc_gamg->m_data[arrsz] = -99.; 76 /* copy data in - column oriented */ 77 if( pc_gamg->m_method != 0 ) { 78 const PetscInt M = Iend - my0; 79 for(kk=0;kk<nloc;kk++){ 80 PetscReal *data = &pc_gamg->m_data[kk*bs]; 81 if( pc_gamg->m_data_cols==1 ) *data = 1.0; 82 else { 83 for(ii=0;ii<bs;ii++) 84 for(jj=0;jj<bs;jj++) 85 if(ii==jj)data[ii*M + jj] = 1.0; /* translational modes */ 86 else data[ii*M + jj] = 0.0; 87 if( a_coords != 0 ) { 88 if( a_ndm == 2 ){ /* rotational modes */ 89 data += 2*M; 90 data[0] = -a_coords[2*kk+1]; 91 data[1] = a_coords[2*kk]; 92 } 93 else { 94 data += 3*M; 95 data[0] = 0.0; data[M+0] = a_coords[3*kk+2]; data[2*M+0] = -a_coords[3*kk+1]; 96 data[1] = -a_coords[3*kk+2]; data[M+1] = 0.0; data[2*M+1] = a_coords[3*kk]; 97 data[2] = a_coords[3*kk+1]; data[M+2] = -a_coords[3*kk]; data[2*M+2] = 0.0; 98 } 99 } 100 } 101 } 102 } 103 else { 104 for( kk = 0 ; kk < nloc ; kk++ ){ 105 for( ii = 0 ; ii < a_ndm ; ii++ ) { 106 pc_gamg->m_data[ii*nloc + kk] = a_coords[kk*a_ndm + ii]; 107 } 108 } 109 } 110 assert(pc_gamg->m_data[arrsz] == -99.); 111 112 pc_gamg->m_data_sz = arrsz; 113 pc_gamg->m_dim = a_ndm; 114 115 PetscFunctionReturn(0); 116 } 117 EXTERN_C_END 118 119 120 /* ----------------------------------------------------------------------------- */ 121 #undef __FUNCT__ 122 #define __FUNCT__ "PCReset_GAMG" 123 PetscErrorCode PCReset_GAMG(PC pc) 124 { 125 PetscErrorCode ierr; 126 PC_MG *mg = (PC_MG*)pc->data; 127 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 128 129 PetscFunctionBegin; 130 if( pc_gamg->m_data != 0 ) { /* this should not happen, cleaned up in SetUp */ 131 ierr = PetscFree(pc_gamg->m_data); CHKERRQ(ierr); 132 } 133 pc_gamg->m_data = 0; pc_gamg->m_data_sz = 0; 134 PetscFunctionReturn(0); 135 } 136 137 /* -------------------------------------------------------------------------- */ 138 /* 139 partitionLevel 140 141 Input Parameter: 142 . a_Amat_fine - matrix on this fine (k) level 143 . a_ndata_rows - size of data to move (coarse grid) 144 . a_ndata_cols - size of data to move (coarse grid) 145 In/Output Parameter: 146 . a_P_inout - prolongation operator to the next level (k-1) 147 . a_coarse_data - data that need to be moved 148 . a_nactive_proc - number of active procs 149 Output Parameter: 150 . a_Amat_crs - coarse matrix that is created (k-1) 151 */ 152 153 #define TOP_GRID_LIM 2*s_min_eq_proc /* this will happen anyway */ 154 155 #undef __FUNCT__ 156 #define __FUNCT__ "partitionLevel" 157 PetscErrorCode partitionLevel( Mat a_Amat_fine, 158 PetscInt a_ndata_rows, 159 PetscInt a_ndata_cols, 160 PetscInt a_cbs, 161 Mat *a_P_inout, 162 PetscReal **a_coarse_data, 163 PetscMPIInt *a_nactive_proc, 164 Mat *a_Amat_crs 165 ) 166 { 167 PetscErrorCode ierr; 168 Mat Cmat,Pnew,Pold=*a_P_inout; 169 IS new_indices,isnum; 170 MPI_Comm wcomm = ((PetscObject)a_Amat_fine)->comm; 171 PetscMPIInt mype,npe,new_npe,nactive; 172 PetscInt neq,NN,Istart,Iend,Istart0,Iend0,ncrs_new,ncrs0; 173 174 PetscFunctionBegin; 175 ierr = MPI_Comm_rank( wcomm, &mype ); CHKERRQ(ierr); 176 ierr = MPI_Comm_size( wcomm, &npe ); CHKERRQ(ierr); 177 /* RAP */ 178 #ifdef USE_R 179 /* make R wih brute force for now */ 180 ierr = MatTranspose( Pold, Pnew ); 181 ierr = MatDestroy( &Pold ); CHKERRQ(ierr); 182 ierr = MatRARt( a_Amat_fine, Pnew, MAT_INITIAL_MATRIX, 2.0, &Cmat ); CHKERRQ(ierr); 183 Pold = Pnew; 184 #else 185 ierr = MatPtAP( a_Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, &Cmat ); CHKERRQ(ierr); 186 #endif 187 ierr = MatSetBlockSize( Cmat, a_cbs ); CHKERRQ(ierr); 188 ierr = MatGetOwnershipRange( Cmat, &Istart0, &Iend0 ); CHKERRQ(ierr); 189 ncrs0 = (Iend0-Istart0)/a_cbs; assert((Iend0-Istart0)%a_cbs == 0); 190 191 if( s_avoid_repart) { 192 *a_Amat_crs = Cmat; /* output */ 193 } 194 else { 195 /* Repartition Cmat_{k} and move colums of P^{k}_{k-1} and coordinates accordingly */ 196 Mat adj; 197 const PetscInt *idx,data_sz=a_ndata_rows*a_ndata_cols; 198 const PetscInt stride0=ncrs0*a_ndata_rows; 199 PetscInt is_sz,*isnewproc_idx,ii,jj,kk,strideNew,*tidx; 200 /* create sub communicator */ 201 MPI_Comm cm; 202 MPI_Group wg, g2; 203 PetscInt *counts,inpe; 204 PetscMPIInt *ranks; 205 IS isscat; 206 PetscScalar *array; 207 Vec src_crd, dest_crd; 208 PetscReal *data = *a_coarse_data; 209 VecScatter vecscat; 210 IS isnewproc; 211 212 /* get number of PEs to make active, reduce */ 213 ierr = MatGetSize( Cmat, &neq, &NN ); CHKERRQ(ierr); 214 new_npe = neq/s_min_eq_proc; /* hardwire min. number of eq/proc */ 215 if( new_npe == 0 || neq < TOP_GRID_LIM ) new_npe = 1; 216 else if (new_npe >= *a_nactive_proc ) new_npe = *a_nactive_proc; /* no change, rare */ 217 218 ierr = PetscMalloc( npe*sizeof(PetscMPIInt), &ranks ); CHKERRQ(ierr); 219 ierr = PetscMalloc( npe*sizeof(PetscInt), &counts ); CHKERRQ(ierr); 220 221 ierr = MPI_Allgather( &ncrs0, 1, MPIU_INT, counts, 1, MPIU_INT, wcomm ); CHKERRQ(ierr); 222 assert(counts[mype]==ncrs0); 223 /* count real active pes */ 224 for( nactive = jj = 0 ; jj < npe ; jj++) { 225 if( counts[jj] != 0 ) { 226 ranks[nactive++] = jj; 227 } 228 } 229 230 if (nactive < new_npe) new_npe = nactive; /* this can happen with empty input procs */ 231 232 #ifdef VERBOSE 233 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); 234 #endif 235 236 *a_nactive_proc = new_npe; /* output */ 237 238 ierr = MPI_Comm_group( wcomm, &wg ); CHKERRQ(ierr); 239 ierr = MPI_Group_incl( wg, nactive, ranks, &g2 ); CHKERRQ(ierr); 240 ierr = MPI_Comm_create( wcomm, g2, &cm ); CHKERRQ(ierr); 241 242 ierr = MPI_Group_free( &wg ); CHKERRQ(ierr); 243 ierr = MPI_Group_free( &g2 ); CHKERRQ(ierr); 244 245 /* MatPartitioningApply call MatConvert, which is collective */ 246 #if defined PETSC_USE_LOG 247 ierr = PetscLogEventBegin(gamg_setup_events[SET12],0,0,0,0);CHKERRQ(ierr); 248 #endif 249 if( a_cbs == 1) { 250 ierr = MatConvert( Cmat, MATMPIADJ, MAT_INITIAL_MATRIX, &adj ); CHKERRQ(ierr); 251 } 252 else{ 253 /* make a scalar matrix to partition */ 254 Mat tMat; 255 PetscInt ncols,jj,Ii; 256 const PetscScalar *vals; 257 const PetscInt *idx; 258 PetscInt *d_nnz; 259 static int llev = 0; 260 261 ierr = PetscMalloc( ncrs0*sizeof(PetscInt), &d_nnz ); CHKERRQ(ierr); 262 for ( Ii = Istart0, jj = 0 ; Ii < Iend0 ; Ii += a_cbs, jj++ ) { 263 ierr = MatGetRow(Cmat,Ii,&ncols,0,0); CHKERRQ(ierr); 264 d_nnz[jj] = ncols/a_cbs; 265 if( d_nnz[jj] > ncrs0 ) d_nnz[jj] = ncrs0; 266 ierr = MatRestoreRow(Cmat,Ii,&ncols,0,0); CHKERRQ(ierr); 267 } 268 269 ierr = MatCreateMPIAIJ( wcomm, ncrs0, ncrs0, 270 PETSC_DETERMINE, PETSC_DETERMINE, 271 0, d_nnz, 0, d_nnz, 272 &tMat ); 273 CHKERRQ(ierr); 274 ierr = PetscFree( d_nnz ); CHKERRQ(ierr); 275 276 for ( ii = Istart0; ii < Iend0; ii++ ) { 277 PetscInt dest_row = ii/a_cbs; 278 ierr = MatGetRow(Cmat,ii,&ncols,&idx,&vals); CHKERRQ(ierr); 279 for( jj = 0 ; jj < ncols ; jj++ ){ 280 PetscInt dest_col = idx[jj]/a_cbs; 281 PetscScalar v = 1.0; 282 ierr = MatSetValues(tMat,1,&dest_row,1,&dest_col,&v,ADD_VALUES); CHKERRQ(ierr); 283 } 284 ierr = MatRestoreRow(Cmat,ii,&ncols,&idx,&vals); CHKERRQ(ierr); 285 } 286 ierr = MatAssemblyBegin(tMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 287 ierr = MatAssemblyEnd(tMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 288 289 if( llev++ == -1 ) { 290 PetscViewer viewer; char fname[32]; 291 sprintf(fname,"part_mat_%d.mat",llev); 292 PetscViewerBinaryOpen(wcomm,fname,FILE_MODE_WRITE,&viewer); 293 ierr = MatView( tMat, viewer ); CHKERRQ(ierr); 294 ierr = PetscViewerDestroy( &viewer ); 295 } 296 297 ierr = MatConvert( tMat, MATMPIADJ, MAT_INITIAL_MATRIX, &adj ); CHKERRQ(ierr); 298 299 ierr = MatDestroy( &tMat ); CHKERRQ(ierr); 300 } 301 if( ncrs0 != 0 ){ 302 const PetscInt *is_idx; 303 MatPartitioning mpart; 304 /* hack to fix global data that pmetis.c uses in 'adj' */ 305 for( nactive = jj = 0 ; jj < npe ; jj++ ) { 306 if( counts[jj] != 0 ) { 307 adj->rmap->range[nactive++] = adj->rmap->range[jj]; 308 } 309 } 310 adj->rmap->range[nactive] = adj->rmap->range[npe]; 311 312 ierr = MatPartitioningCreate( cm, &mpart ); CHKERRQ(ierr); 313 ierr = MatPartitioningSetAdjacency( mpart, adj ); CHKERRQ(ierr); 314 ierr = MatPartitioningSetFromOptions( mpart ); CHKERRQ(ierr); 315 ierr = MatPartitioningSetNParts( mpart, new_npe );CHKERRQ(ierr); 316 ierr = MatPartitioningApply( mpart, &isnewproc ); CHKERRQ(ierr); 317 ierr = MatPartitioningDestroy( &mpart ); CHKERRQ(ierr); 318 319 /* collect IS info */ 320 ierr = ISGetLocalSize( isnewproc, &is_sz ); CHKERRQ(ierr); 321 ierr = PetscMalloc( a_cbs*is_sz*sizeof(PetscInt), &isnewproc_idx ); CHKERRQ(ierr); 322 ierr = ISGetIndices( isnewproc, &is_idx ); CHKERRQ(ierr); 323 /* spread partitioning across machine - best way ??? */ 324 NN = 1; /*npe/new_npe;*/ 325 for( kk = jj = 0 ; kk < is_sz ; kk++ ){ 326 for( ii = 0 ; ii < a_cbs ; ii++, jj++ ) { 327 isnewproc_idx[jj] = is_idx[kk] * NN; /* distribution */ 328 } 329 } 330 ierr = ISRestoreIndices( isnewproc, &is_idx ); CHKERRQ(ierr); 331 ierr = ISDestroy( &isnewproc ); CHKERRQ(ierr); 332 ierr = MPI_Comm_free( &cm ); CHKERRQ(ierr); 333 334 is_sz *= a_cbs; 335 } 336 else{ 337 isnewproc_idx = 0; 338 is_sz = 0; 339 } 340 ierr = MatDestroy( &adj ); CHKERRQ(ierr); 341 ierr = ISCreateGeneral( wcomm, is_sz, isnewproc_idx, PETSC_COPY_VALUES, &isnewproc ); 342 if( isnewproc_idx != 0 ) { 343 ierr = PetscFree( isnewproc_idx ); CHKERRQ(ierr); 344 } 345 346 /* 347 Create an index set from the isnewproc index set to indicate the mapping TO 348 */ 349 ierr = ISPartitioningToNumbering( isnewproc, &isnum ); CHKERRQ(ierr); 350 /* 351 Determine how many elements are assigned to each processor 352 */ 353 inpe = npe; 354 ierr = ISPartitioningCount( isnewproc, inpe, counts ); CHKERRQ(ierr); 355 ierr = ISDestroy( &isnewproc ); CHKERRQ(ierr); 356 ncrs_new = counts[mype]/a_cbs; 357 strideNew = ncrs_new*a_ndata_rows; 358 #if defined PETSC_USE_LOG 359 ierr = PetscLogEventEnd(gamg_setup_events[SET12],0,0,0,0); CHKERRQ(ierr); 360 #endif 361 /* Create a vector to contain the newly ordered element information */ 362 ierr = VecCreate( wcomm, &dest_crd ); 363 ierr = VecSetSizes( dest_crd, data_sz*ncrs_new, PETSC_DECIDE ); CHKERRQ(ierr); 364 ierr = VecSetFromOptions( dest_crd ); CHKERRQ(ierr); /* this is needed! */ 365 /* 366 There are 'a_ndata_rows*a_ndata_cols' data items per node, (one can think of the vectors of having 367 a block size of ...). Note, ISs are expanded into equation space by 'a_cbs'. 368 */ 369 ierr = PetscMalloc( (ncrs0*data_sz)*sizeof(PetscInt), &tidx ); CHKERRQ(ierr); 370 ierr = ISGetIndices( isnum, &idx ); CHKERRQ(ierr); 371 for(ii=0,jj=0; ii<ncrs0 ; ii++) { 372 PetscInt id = idx[ii*a_cbs]/a_cbs; /* get node back */ 373 for( kk=0; kk<data_sz ; kk++, jj++) tidx[jj] = id*data_sz + kk; 374 } 375 ierr = ISRestoreIndices( isnum, &idx ); CHKERRQ(ierr); 376 ierr = ISCreateGeneral( wcomm, data_sz*ncrs0, tidx, PETSC_COPY_VALUES, &isscat ); 377 CHKERRQ(ierr); 378 ierr = PetscFree( tidx ); CHKERRQ(ierr); 379 /* 380 Create a vector to contain the original vertex information for each element 381 */ 382 ierr = VecCreateSeq( PETSC_COMM_SELF, data_sz*ncrs0, &src_crd ); CHKERRQ(ierr); 383 for( jj=0; jj<a_ndata_cols ; jj++ ) { 384 for( ii=0 ; ii<ncrs0 ; ii++) { 385 for( kk=0; kk<a_ndata_rows ; kk++ ) { 386 PetscInt ix = ii*a_ndata_rows + kk + jj*stride0, jx = ii*data_sz + kk*a_ndata_cols + jj; 387 PetscScalar tt = (PetscScalar)data[ix]; 388 ierr = VecSetValues( src_crd, 1, &jx, &tt, INSERT_VALUES ); CHKERRQ(ierr); 389 } 390 } 391 } 392 ierr = VecAssemblyBegin(src_crd); CHKERRQ(ierr); 393 ierr = VecAssemblyEnd(src_crd); CHKERRQ(ierr); 394 /* 395 Scatter the element vertex information (still in the original vertex ordering) 396 to the correct processor 397 */ 398 ierr = VecScatterCreate( src_crd, PETSC_NULL, dest_crd, isscat, &vecscat); 399 CHKERRQ(ierr); 400 ierr = ISDestroy( &isscat ); CHKERRQ(ierr); 401 ierr = VecScatterBegin(vecscat,src_crd,dest_crd,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 402 ierr = VecScatterEnd(vecscat,src_crd,dest_crd,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 403 ierr = VecScatterDestroy( &vecscat ); CHKERRQ(ierr); 404 ierr = VecDestroy( &src_crd ); CHKERRQ(ierr); 405 /* 406 Put the element vertex data into a new allocation of the gdata->ele 407 */ 408 ierr = PetscFree( *a_coarse_data ); CHKERRQ(ierr); 409 ierr = PetscMalloc( data_sz*ncrs_new*sizeof(PetscReal), a_coarse_data ); CHKERRQ(ierr); 410 411 ierr = VecGetArray( dest_crd, &array ); CHKERRQ(ierr); 412 data = *a_coarse_data; 413 for( jj=0; jj<a_ndata_cols ; jj++ ) { 414 for( ii=0 ; ii<ncrs_new ; ii++) { 415 for( kk=0; kk<a_ndata_rows ; kk++ ) { 416 PetscInt ix = ii*a_ndata_rows + kk + jj*strideNew, jx = ii*data_sz + kk*a_ndata_cols + jj; 417 data[ix] = PetscRealPart(array[jx]); 418 array[jx] = 1.e300; 419 } 420 } 421 } 422 ierr = VecRestoreArray( dest_crd, &array ); CHKERRQ(ierr); 423 ierr = VecDestroy( &dest_crd ); CHKERRQ(ierr); 424 /* 425 Invert for MatGetSubMatrix 426 */ 427 ierr = ISInvertPermutation( isnum, ncrs_new*a_cbs, &new_indices ); CHKERRQ(ierr); 428 ierr = ISSort( new_indices ); CHKERRQ(ierr); /* is this needed? */ 429 ierr = ISDestroy( &isnum ); CHKERRQ(ierr); 430 /* A_crs output */ 431 ierr = MatGetSubMatrix( Cmat, new_indices, new_indices, MAT_INITIAL_MATRIX, a_Amat_crs ); 432 CHKERRQ(ierr); 433 434 ierr = MatDestroy( &Cmat ); CHKERRQ(ierr); 435 Cmat = *a_Amat_crs; /* output */ 436 ierr = MatSetBlockSize( Cmat, a_cbs ); CHKERRQ(ierr); 437 438 /* prolongator */ 439 ierr = MatGetOwnershipRange( Pold, &Istart, &Iend ); CHKERRQ(ierr); 440 { 441 IS findices; 442 ierr = ISCreateStride(wcomm,Iend-Istart,Istart,1,&findices); CHKERRQ(ierr); 443 #ifdef USE_R 444 ierr = MatGetSubMatrix( Pold, new_indices, findices, MAT_INITIAL_MATRIX, &Pnew ); 445 #else 446 ierr = MatGetSubMatrix( Pold, findices, new_indices, MAT_INITIAL_MATRIX, &Pnew ); 447 #endif 448 CHKERRQ(ierr); 449 450 ierr = ISDestroy( &findices ); CHKERRQ(ierr); 451 } 452 453 ierr = MatDestroy( a_P_inout ); CHKERRQ(ierr); 454 *a_P_inout = Pnew; /* output */ 455 456 ierr = ISDestroy( &new_indices ); CHKERRQ(ierr); 457 ierr = PetscFree( counts ); CHKERRQ(ierr); 458 ierr = PetscFree( ranks ); CHKERRQ(ierr); 459 } 460 461 PetscFunctionReturn(0); 462 } 463 464 /* -------------------------------------------------------------------------- */ 465 /* 466 PCSetUp_GAMG - Prepares for the use of the GAMG preconditioner 467 by setting data structures and options. 468 469 Input Parameter: 470 . pc - the preconditioner context 471 472 Application Interface Routine: PCSetUp() 473 474 Notes: 475 The interface routine PCSetUp() is not usually called directly by 476 the user, but instead is called by PCApply() if necessary. 477 */ 478 #undef __FUNCT__ 479 #define __FUNCT__ "PCSetUp_GAMG" 480 PetscErrorCode PCSetUp_GAMG( PC a_pc ) 481 { 482 PetscErrorCode ierr; 483 PC_MG *mg = (PC_MG*)a_pc->data; 484 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 485 PC_MG_Levels **mglevels = mg->levels; 486 Mat Amat = a_pc->mat, Pmat = a_pc->pmat; 487 PetscInt fine_level, level, level1, M, N, bs, nloc, lidx, Istart, Iend; 488 MPI_Comm wcomm = ((PetscObject)a_pc)->comm; 489 PetscMPIInt mype,npe,nactivepe; 490 PetscBool isOK; 491 Mat Aarr[GAMG_MAXLEVELS], Parr[GAMG_MAXLEVELS]; 492 PetscReal *coarse_data = 0, *data, emaxs[GAMG_MAXLEVELS]; 493 MatInfo info; 494 495 PetscFunctionBegin; 496 pc_gamg->m_count++; 497 if( a_pc->setupcalled > 0 ) { 498 /* just do Galerkin grids */ 499 Mat B,dA,dB; 500 501 /* PCSetUp_MG seems to insists on setting this to GMRES */ 502 ierr = KSPSetType( mglevels[0]->smoothd, KSPPREONLY ); CHKERRQ(ierr); 503 504 /* currently only handle case where mat and pmat are the same on coarser levels */ 505 ierr = KSPGetOperators(mglevels[pc_gamg->m_Nlevels-1]->smoothd,&dA,&dB,PETSC_NULL);CHKERRQ(ierr); 506 /* (re)set to get dirty flag */ 507 ierr = KSPSetOperators(mglevels[pc_gamg->m_Nlevels-1]->smoothd,dA,dB,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 508 ierr = KSPSetUp( mglevels[pc_gamg->m_Nlevels-1]->smoothd ); CHKERRQ(ierr); 509 510 for (level=pc_gamg->m_Nlevels-2; level>-1; level--) { 511 ierr = KSPGetOperators(mglevels[level]->smoothd,PETSC_NULL,&B,PETSC_NULL);CHKERRQ(ierr); 512 /* the first time through the matrix structure has changed from repartitioning */ 513 if( pc_gamg->m_count == 2 ) { 514 ierr = MatDestroy( &B ); CHKERRQ(ierr); 515 ierr = MatPtAP(dB,mglevels[level+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);CHKERRQ(ierr); 516 mglevels[level]->A = B; 517 } 518 else { 519 ierr = MatPtAP(dB,mglevels[level+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);CHKERRQ(ierr); 520 } 521 ierr = KSPSetOperators(mglevels[level]->smoothd,B,B,SAME_NONZERO_PATTERN); CHKERRQ(ierr); 522 dB = B; 523 /* setup KSP/PC */ 524 ierr = KSPSetUp( mglevels[level]->smoothd ); CHKERRQ(ierr); 525 } 526 527 #define PRINT_MATS !PETSC_TRUE 528 /* plot levels - A */ 529 if( PRINT_MATS ) { 530 for (lidx=0, level=pc_gamg->m_Nlevels-1; level>0 ; level--,lidx++){ 531 PetscViewer viewer; 532 char fname[32]; KSP smoother; Mat Tmat, TTm; 533 ierr = PCMGGetSmoother( a_pc, lidx, &smoother ); CHKERRQ(ierr); 534 ierr = KSPGetOperators( smoother, &Tmat, &TTm, 0 ); CHKERRQ(ierr); 535 sprintf(fname,"Amat_%d_%d.m",(int)pc_gamg->m_count,(int)level); 536 ierr = PetscViewerASCIIOpen( wcomm, fname, &viewer ); CHKERRQ(ierr); 537 ierr = PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB); CHKERRQ(ierr); 538 ierr = MatView( Tmat, viewer ); CHKERRQ(ierr); 539 ierr = PetscViewerDestroy( &viewer ); 540 } 541 } 542 543 a_pc->setupcalled = 2; 544 545 PetscFunctionReturn(0); 546 } 547 ierr = MPI_Comm_rank(wcomm,&mype);CHKERRQ(ierr); 548 ierr = MPI_Comm_size(wcomm,&npe);CHKERRQ(ierr); 549 /* GAMG requires input of fine-grid matrix. It determines nlevels. */ 550 ierr = MatGetBlockSize( Amat, &bs ); CHKERRQ(ierr); 551 ierr = MatGetSize( Amat, &M, &N );CHKERRQ(ierr); 552 ierr = MatGetOwnershipRange( Amat, &Istart, &Iend ); CHKERRQ(ierr); 553 nloc = (Iend-Istart)/bs; assert((Iend-Istart)%bs == 0); 554 555 /* get data of not around */ 556 if( pc_gamg->m_data == 0 && nloc > 0 ) { 557 ierr = PCSetCoordinates_GAMG( a_pc, -1, 0 ); CHKERRQ( ierr ); 558 } 559 data = pc_gamg->m_data; 560 561 /* Get A_i and R_i */ 562 ierr = MatGetInfo(Amat,MAT_GLOBAL_SUM,&info); CHKERRQ(ierr); 563 #ifdef VERBOSE 564 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", 565 mype,__FUNCT__,0,N,pc_gamg->m_data_rows,pc_gamg->m_data_cols, 566 (int)(info.nz_used/(PetscReal)N),npe); 567 #endif 568 for ( level=0, Aarr[0] = Pmat, nactivepe = npe; /* hard wired stopping logic */ 569 level < (GAMG_MAXLEVELS-1) && (level==0 || M>TOP_GRID_LIM); /* && (npe==1 || nactivepe>1); */ 570 level++ ){ 571 level1 = level + 1; 572 #if (defined PETSC_USE_LOG && defined GAMG_STAGES) 573 ierr = PetscLogStagePush(gamg_stages[level]); CHKERRQ( ierr ); 574 #endif 575 #if defined PETSC_USE_LOG 576 ierr = PetscLogEventBegin(gamg_setup_events[SET1],0,0,0,0);CHKERRQ(ierr); 577 #endif 578 ierr = createProlongation(Aarr[level], data, pc_gamg->m_dim, pc_gamg->m_data_cols, pc_gamg->m_method, 579 level, s_threshold, &bs, &Parr[level1], &coarse_data, &isOK, &emaxs[level] ); 580 CHKERRQ(ierr); 581 ierr = PetscFree( data ); CHKERRQ( ierr ); 582 #if defined PETSC_USE_LOG 583 ierr = PetscLogEventEnd(gamg_setup_events[SET1],0,0,0,0);CHKERRQ(ierr); 584 #endif 585 if(level==0) Aarr[0] = Amat; /* use Pmat for finest level setup, but use mat for solver */ 586 if( isOK ) { 587 #if defined PETSC_USE_LOG 588 ierr = PetscLogEventBegin(gamg_setup_events[SET2],0,0,0,0);CHKERRQ(ierr); 589 #endif 590 ierr = partitionLevel( Aarr[level], (pc_gamg->m_method != 0) ? bs : 1, pc_gamg->m_data_cols, bs, 591 &Parr[level1], &coarse_data, &nactivepe, &Aarr[level1] ); 592 CHKERRQ(ierr); 593 #if defined PETSC_USE_LOG 594 ierr = PetscLogEventEnd(gamg_setup_events[SET2],0,0,0,0);CHKERRQ(ierr); 595 #endif 596 ierr = MatGetSize( Aarr[level1], &M, &N );CHKERRQ(ierr); 597 ierr = MatGetInfo(Aarr[level1],MAT_GLOBAL_SUM,&info); CHKERRQ(ierr); 598 #ifdef VERBOSE 599 PetscPrintf(PETSC_COMM_WORLD,"\t\t[%d]%s %d) N=%d, n data cols=%d, nnz/row (ave)=%d, %d active pes\n", 600 mype,__FUNCT__,(int)level1,N,pc_gamg->m_data_cols, 601 (int)(info.nz_used/(PetscReal)N),nactivepe); 602 #endif 603 /* coarse grids with SA can have zero row/cols from singleton aggregates */ 604 /* aggregation method should gaurrentee this does not happen! */ 605 606 #ifdef VERBOSE 607 if( PETSC_TRUE ){ 608 Vec diag; PetscScalar *data_arr,v; PetscInt Istart,Iend,kk,nloceq,id; 609 v = 1.e-10; /* LU factor has hard wired numbers for small diags so this needs to match (yuk) */ 610 ierr = MatGetOwnershipRange(Aarr[level1], &Istart, &Iend); CHKERRQ(ierr); 611 nloceq = Iend-Istart; 612 ierr = MatGetVecs( Aarr[level1], &diag, 0 ); CHKERRQ(ierr); 613 ierr = MatGetDiagonal( Aarr[level1], diag ); CHKERRQ(ierr); 614 ierr = VecGetArray( diag, &data_arr ); CHKERRQ(ierr); 615 for(kk=0;kk<nloceq;kk++){ 616 if(data_arr[kk]==0.0) { 617 id = kk + Istart; 618 ierr = MatSetValues(Aarr[level1],1,&id,1,&id,&v,INSERT_VALUES); 619 CHKERRQ(ierr); 620 PetscPrintf(PETSC_COMM_SELF,"\t[%d]%s warning: added zero to diag (%d) on level %d \n",mype,__FUNCT__,id,level1); 621 } 622 } 623 ierr = VecRestoreArray( diag, &data_arr ); CHKERRQ(ierr); 624 ierr = VecDestroy( &diag ); CHKERRQ(ierr); 625 ierr = MatAssemblyBegin(Aarr[level1],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 626 ierr = MatAssemblyEnd(Aarr[level1],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 627 } 628 #endif 629 } 630 else{ 631 coarse_data = 0; 632 break; 633 } 634 data = coarse_data; 635 636 #if (defined PETSC_USE_LOG && defined GAMG_STAGES) 637 ierr = PetscLogStagePop(); CHKERRQ( ierr ); 638 #endif 639 } 640 if( coarse_data ) { 641 ierr = PetscFree( coarse_data ); CHKERRQ( ierr ); 642 } 643 #ifdef VERBOSE 644 PetscPrintf(PETSC_COMM_WORLD,"\t[%d]%s %d levels\n",0,__FUNCT__,level + 1); 645 #endif 646 pc_gamg->m_data = 0; /* destroyed coordinate data */ 647 pc_gamg->m_Nlevels = level + 1; 648 fine_level = level; 649 ierr = PCMGSetLevels(a_pc,pc_gamg->m_Nlevels,PETSC_NULL);CHKERRQ(ierr); 650 651 /* set default smoothers */ 652 for ( lidx=1, level = pc_gamg->m_Nlevels-2; 653 lidx <= fine_level; 654 lidx++, level--) { 655 PetscReal emax, emin; 656 KSP smoother; PC subpc; 657 ierr = PCMGGetSmoother( a_pc, lidx, &smoother ); CHKERRQ(ierr); 658 ierr = KSPSetType( smoother, KSPCHEBYCHEV );CHKERRQ(ierr); 659 if( emaxs[level] > 0.0 ) emax=emaxs[level]; 660 else{ /* eigen estimate 'emax' */ 661 KSP eksp; Mat Lmat = Aarr[level]; 662 Vec bb, xx; PC pc; 663 664 ierr = MatGetVecs( Lmat, &bb, 0 ); CHKERRQ(ierr); 665 ierr = MatGetVecs( Lmat, &xx, 0 ); CHKERRQ(ierr); 666 { 667 PetscRandom rctx; 668 ierr = PetscRandomCreate(wcomm,&rctx);CHKERRQ(ierr); 669 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 670 ierr = VecSetRandom(bb,rctx);CHKERRQ(ierr); 671 ierr = PetscRandomDestroy( &rctx ); CHKERRQ(ierr); 672 } 673 ierr = KSPCreate(wcomm,&eksp);CHKERRQ(ierr); 674 ierr = KSPSetType( eksp, KSPCG ); CHKERRQ(ierr); 675 ierr = KSPSetInitialGuessNonzero( eksp, PETSC_FALSE ); CHKERRQ(ierr); 676 ierr = KSPSetOperators( eksp, Lmat, Lmat, SAME_NONZERO_PATTERN ); CHKERRQ( ierr ); 677 ierr = KSPGetPC( eksp, &pc );CHKERRQ( ierr ); 678 ierr = PCSetType( pc, PETSC_GAMG_SMOOTHER ); CHKERRQ(ierr); /* should be same as eigen estimates op. */ 679 ierr = KSPSetTolerances( eksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 10 ); 680 CHKERRQ(ierr); 681 ierr = KSPSetNormType( eksp, KSP_NORM_NONE ); CHKERRQ(ierr); 682 683 ierr = KSPSetComputeSingularValues( eksp,PETSC_TRUE ); CHKERRQ(ierr); 684 ierr = KSPSolve( eksp, bb, xx ); CHKERRQ(ierr); 685 ierr = KSPComputeExtremeSingularValues( eksp, &emax, &emin ); CHKERRQ(ierr); 686 ierr = VecDestroy( &xx ); CHKERRQ(ierr); 687 ierr = VecDestroy( &bb ); CHKERRQ(ierr); 688 ierr = KSPDestroy( &eksp ); CHKERRQ(ierr); 689 #ifdef VERBOSE 690 PetscPrintf(PETSC_COMM_WORLD,"\t\t\t%s PC setup max eigen=%e min=%e PC=%s\n",__FUNCT__,emax,emin,PETSC_GAMG_SMOOTHER); 691 #endif 692 } 693 { 694 PetscInt N1, N0, tt; 695 ierr = MatGetSize( Aarr[level], &N1, &tt ); CHKERRQ(ierr); 696 ierr = MatGetSize( Aarr[level+1], &N0, &tt ); CHKERRQ(ierr); 697 emin = 1.*emax/((PetscReal)N1/(PetscReal)N0); /* this should be about the coarsening rate */ 698 emax *= 1.05; 699 } 700 701 ierr = KSPSetOperators( smoother, Aarr[level], Aarr[level], SAME_NONZERO_PATTERN ); 702 ierr = KSPChebychevSetEigenvalues( smoother, emax, emin );CHKERRQ(ierr); 703 /* ierr = KSPSetTolerances(smoother,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,2); CHKERRQ(ierr); */ 704 ierr = KSPGetPC( smoother, &subpc ); CHKERRQ(ierr); 705 ierr = PCSetType( subpc, PETSC_GAMG_SMOOTHER ); CHKERRQ(ierr); 706 ierr = KSPSetNormType( smoother, KSP_NORM_NONE ); CHKERRQ(ierr); 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) s_min_eq_proc = n; 843 PetscFunctionReturn(0); 844 } 845 EXTERN_C_END 846 847 #undef __FUNCT__ 848 #define __FUNCT__ "PCGAMGAvoidRepartitioning" 849 /*@ 850 PCGAMGAvoidRepartitioning - Do not repartition the 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_avoid_repartitioning 860 861 Level: intermediate 862 863 Concepts: Unstructured multrigrid preconditioner 864 865 .seealso: () 866 @*/ 867 PetscErrorCode PCGAMGAvoidRepartitioning(PC pc, PetscBool n) 868 { 869 PetscErrorCode ierr; 870 871 PetscFunctionBegin; 872 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 873 ierr = PetscTryMethod(pc,"PCGAMGAvoidRepartitioning_C",(PC,PetscBool),(pc,n));CHKERRQ(ierr); 874 PetscFunctionReturn(0); 875 } 876 877 EXTERN_C_BEGIN 878 #undef __FUNCT__ 879 #define __FUNCT__ "PCGAMGAvoidRepartitioning_GAMG" 880 PetscErrorCode PCGAMGAvoidRepartitioning_GAMG(PC pc, PetscBool n) 881 { 882 /* PC_MG *mg = (PC_MG*)pc->data; */ 883 /* PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; */ 884 885 PetscFunctionBegin; 886 s_avoid_repart = n; 887 PetscFunctionReturn(0); 888 } 889 EXTERN_C_END 890 891 #undef __FUNCT__ 892 #define __FUNCT__ "PCGAMGSetThreshold" 893 /*@ 894 PCGAMGSetThreshold - Relative threshold to use for dropping edges in aggregation graph 895 896 Not collective on PC 897 898 Input Parameters: 899 . pc - the preconditioner context 900 901 902 Options Database Key: 903 . -pc_gamg_threshold 904 905 Level: intermediate 906 907 Concepts: Unstructured multrigrid preconditioner 908 909 .seealso: () 910 @*/ 911 PetscErrorCode PCGAMGSetThreshold(PC pc, PetscReal n) 912 { 913 PetscErrorCode ierr; 914 915 PetscFunctionBegin; 916 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 917 ierr = PetscTryMethod(pc,"PCGAMGSetThreshold_C",(PC,PetscReal),(pc,n));CHKERRQ(ierr); 918 PetscFunctionReturn(0); 919 } 920 921 EXTERN_C_BEGIN 922 #undef __FUNCT__ 923 #define __FUNCT__ "PCGAMGSetThreshold_GAMG" 924 PetscErrorCode PCGAMGSetThreshold_GAMG(PC pc, PetscReal n) 925 { 926 /* PC_MG *mg = (PC_MG*)pc->data; */ 927 /* PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; */ 928 929 PetscFunctionBegin; 930 s_threshold = n; 931 PetscFunctionReturn(0); 932 } 933 EXTERN_C_END 934 935 #undef __FUNCT__ 936 #define __FUNCT__ "PCGAMGSetSolverType" 937 /*@ 938 PCGAMGSetSolverType - Set solution method. 939 940 Collective on PC 941 942 Input Parameters: 943 . pc - the preconditioner context 944 945 946 Options Database Key: 947 . -pc_gamg_type 948 949 Level: intermediate 950 951 Concepts: Unstructured multrigrid preconditioner 952 953 .seealso: () 954 @*/ 955 PetscErrorCode PCGAMGSetSolverType(PC pc, char str[], PetscInt sz ) 956 { 957 PetscErrorCode ierr; 958 959 PetscFunctionBegin; 960 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 961 ierr = PetscTryMethod(pc,"PCGAMGSetSolverType_C",(PC,char[],PetscInt),(pc,str,sz)); 962 CHKERRQ(ierr); 963 PetscFunctionReturn(0); 964 } 965 966 EXTERN_C_BEGIN 967 #undef __FUNCT__ 968 #define __FUNCT__ "PCGAMGSetSolverType_GAMG" 969 PetscErrorCode PCGAMGSetSolverType_GAMG(PC pc, char str[], PetscInt sz ) 970 { 971 PC_MG *mg = (PC_MG*)pc->data; 972 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 973 974 PetscFunctionBegin; 975 if(sz < 64) strcpy(pc_gamg->m_type,str); 976 PetscFunctionReturn(0); 977 } 978 EXTERN_C_END 979 980 #undef __FUNCT__ 981 #define __FUNCT__ "PCSetFromOptions_GAMG" 982 PetscErrorCode PCSetFromOptions_GAMG(PC pc) 983 { 984 PetscErrorCode ierr; 985 PC_MG *mg = (PC_MG*)pc->data; 986 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 987 PetscBool flag; 988 989 PetscFunctionBegin; 990 ierr = PetscOptionsHead("GAMG options"); CHKERRQ(ierr); 991 { 992 ierr = PetscOptionsString("-pc_gamg_type", 993 "Solver type: plane aggregation ('pa'), smoothed aggregation ('sa') or geometric multigrid (default)", 994 "PCGAMGSetSolverType", 995 pc_gamg->m_type, 996 pc_gamg->m_type, 997 64, 998 &flag ); 999 CHKERRQ(ierr); 1000 if( flag && pc_gamg->m_data != 0 ) { 1001 if( (strcmp(pc_gamg->m_type,"sa")==0 && pc_gamg->m_method != 2) || 1002 (strcmp(pc_gamg->m_type,"pa")==0 && pc_gamg->m_method != 1) || 1003 (strcmp(pc_gamg->m_type,"geo")==0 && pc_gamg->m_method != 0) ) { 1004 SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "PCSetFromOptions called of PCSetCoordinates (with new method, after data was created)"); 1005 } 1006 } 1007 1008 if (flag && strcmp(pc_gamg->m_type,"sa") == 0) pc_gamg->m_method = 2; 1009 else if (flag && strcmp(pc_gamg->m_type,"pa") == 0) pc_gamg->m_method = 1; 1010 else pc_gamg->m_method = 0; 1011 1012 /* common (static) variable */ 1013 ierr = PetscOptionsBool("-pc_gamg_avoid_repartitioning", 1014 "Do not repartion coarse grids (false)", 1015 "PCGAMGAvoidRepartitioning", 1016 s_avoid_repart, 1017 &s_avoid_repart, 1018 &flag); 1019 CHKERRQ(ierr); 1020 1021 /* common (static) variable */ 1022 ierr = PetscOptionsInt("-pc_gamg_process_eq_limit", 1023 "Limit (goal) on number of equations per process on coarse grids", 1024 "PCGAMGSetProcEqLim", 1025 s_min_eq_proc, 1026 &s_min_eq_proc, 1027 &flag ); 1028 CHKERRQ(ierr); 1029 1030 /* common (static) variable */ 1031 ierr = PetscOptionsReal("-pc_gamg_threshold", 1032 "Relative threshold to use for dropping edges in aggregation graph", 1033 "PCGAMGSetThreshold", 1034 s_threshold, 1035 &s_threshold, 1036 &flag ); 1037 CHKERRQ(ierr); 1038 } 1039 ierr = PetscOptionsTail();CHKERRQ(ierr); 1040 1041 PetscFunctionReturn(0); 1042 } 1043 1044 /* -------------------------------------------------------------------------- */ 1045 /* 1046 PCCreate_GAMG - Creates a GAMG preconditioner context, PC_GAMG 1047 1048 Input Parameter: 1049 . pc - the preconditioner context 1050 1051 Application Interface Routine: PCCreate() 1052 1053 */ 1054 /* MC 1055 PCGAMG - Geometric algebraic multigrid (AMG) preconditioning. This preconditioner currently has two 1056 AMG methods: 1) an unstructured geometric method, which requires that you provide coordinates for each 1057 vertex, and 2) smoothed aggregation. Smoothed aggregation (SA) can work without coordinates but it 1058 will generate some common non-trivial null spaces if coordinates are provided. The input fine grid matrix 1059 must have the block size set for 'system' problems (with multiple dof per vertex/cell) to work properly. 1060 SA will generate rotational rigid body mode null space vectors, in addition to the trivial translational 1061 modes, when coordinates are provide in 2D and 3D. 1062 1063 Options Database Key: 1064 Multigrid options(inherited) 1065 + -pc_mg_cycles <1>: 1 for V cycle, 2 for W-cycle (MGSetCycles) 1066 . -pc_mg_smoothup <1>: Number of post-smoothing steps (MGSetNumberSmoothUp) 1067 . -pc_mg_smoothdown <1>: Number of pre-smoothing steps (MGSetNumberSmoothDown) 1068 -pc_mg_type <multiplicative>: (one of) additive multiplicative full cascade kascade 1069 GAMG options: 1070 1071 Level: intermediate 1072 Concepts: multigrid 1073 1074 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, 1075 PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), MPSetCycles(), PCMGSetNumberSmoothDown(), 1076 PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(), 1077 PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(), 1078 PCMGSetCyclesOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR() 1079 M */ 1080 1081 EXTERN_C_BEGIN 1082 #undef __FUNCT__ 1083 #define __FUNCT__ "PCCreate_GAMG" 1084 PetscErrorCode PCCreate_GAMG(PC pc) 1085 { 1086 PetscErrorCode ierr; 1087 PC_GAMG *pc_gamg; 1088 PC_MG *mg; 1089 PetscClassId cookie; 1090 #if defined PETSC_USE_LOG 1091 static int count = 0; 1092 #endif 1093 1094 PetscFunctionBegin; 1095 /* PCGAMG is an inherited class of PCMG. Initialize pc as PCMG */ 1096 ierr = PCSetType(pc,PCMG);CHKERRQ(ierr); /* calls PCCreate_MG() and MGCreate_Private() */ 1097 ierr = PetscObjectChangeTypeName((PetscObject)pc,PCGAMG);CHKERRQ(ierr); 1098 1099 /* create a supporting struct and attach it to pc */ 1100 ierr = PetscNewLog(pc,PC_GAMG,&pc_gamg);CHKERRQ(ierr); 1101 pc_gamg->m_data_sz = 0; pc_gamg->m_data = 0; pc_gamg->m_count = 0; 1102 mg = (PC_MG*)pc->data; 1103 mg->innerctx = pc_gamg; 1104 1105 pc_gamg->m_Nlevels = -1; 1106 1107 /* overwrite the pointers of PCMG by the functions of PCGAMG */ 1108 pc->ops->setfromoptions = PCSetFromOptions_GAMG; 1109 pc->ops->setup = PCSetUp_GAMG; 1110 pc->ops->reset = PCReset_GAMG; 1111 pc->ops->destroy = PCDestroy_GAMG; 1112 1113 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1114 "PCSetCoordinates_C", 1115 "PCSetCoordinates_GAMG", 1116 PCSetCoordinates_GAMG); 1117 CHKERRQ(ierr); 1118 1119 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1120 "PCGAMGSetProcEqLim_C", 1121 "PCGAMGSetProcEqLim_GAMG", 1122 PCGAMGSetProcEqLim_GAMG); 1123 CHKERRQ(ierr); 1124 1125 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1126 "PCGAMGAvoidRepartitioning_C", 1127 "PCGAMGAvoidRepartitioning_GAMG", 1128 PCGAMGAvoidRepartitioning_GAMG); 1129 CHKERRQ(ierr); 1130 1131 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1132 "PCGAMGSetThreshold_C", 1133 "PCGAMGSetThreshold_GAMG", 1134 PCGAMGSetThreshold_GAMG); 1135 CHKERRQ(ierr); 1136 1137 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1138 "PCGAMGSetSolverType_C", 1139 "PCGAMGSetSolverType_GAMG", 1140 PCGAMGSetSolverType_GAMG); 1141 CHKERRQ(ierr); 1142 1143 #if defined PETSC_USE_LOG 1144 if( count++ == 0 ) { 1145 PetscClassIdRegister("GAMG Setup",&cookie); 1146 PetscLogEventRegister("GAMG: createProl", cookie, &gamg_setup_events[SET1]); 1147 PetscLogEventRegister(" Graph", cookie, &gamg_setup_events[GRAPH]); 1148 PetscLogEventRegister(" G.Mat", cookie, &gamg_setup_events[GRAPH_MAT]); 1149 PetscLogEventRegister(" G.Filter", cookie, &gamg_setup_events[GRAPH_FILTER]); 1150 PetscLogEventRegister(" G.Square", cookie, &gamg_setup_events[GRAPH_SQR]); 1151 PetscLogEventRegister(" MIS/Agg", cookie, &gamg_setup_events[SET4]); 1152 PetscLogEventRegister(" geo: growSupp", cookie, &gamg_setup_events[SET5]); 1153 PetscLogEventRegister(" geo: triangle", cookie, &gamg_setup_events[SET6]); 1154 PetscLogEventRegister(" search&set", cookie, &gamg_setup_events[FIND_V]); 1155 PetscLogEventRegister(" SA: init", cookie, &gamg_setup_events[SET7]); 1156 /* PetscLogEventRegister(" SA: frmProl0", cookie, &gamg_setup_events[SET8]); */ 1157 PetscLogEventRegister(" SA: smooth", cookie, &gamg_setup_events[SET9]); 1158 PetscLogEventRegister("GAMG: partLevel", cookie, &gamg_setup_events[SET2]); 1159 PetscLogEventRegister(" PL repartition", cookie, &gamg_setup_events[SET12]); 1160 /* PetscLogEventRegister(" PL move data", cookie, &gamg_setup_events[SET13]); */ 1161 /* PetscLogEventRegister("GAMG: fix", cookie, &gamg_setup_events[SET10]); */ 1162 /* PetscLogEventRegister("GAMG: set levels", cookie, &gamg_setup_events[SET11]); */ 1163 1164 /* create timer stages */ 1165 #if defined GAMG_STAGES 1166 { 1167 char str[32]; 1168 sprintf(str,"MG Level %d (finest)",0); 1169 PetscLogStageRegister(str, &gamg_stages[0]); 1170 PetscInt lidx; 1171 for (lidx=1;lidx<9;lidx++){ 1172 sprintf(str,"MG Level %d",lidx); 1173 PetscLogStageRegister(str, &gamg_stages[lidx]); 1174 } 1175 } 1176 #endif 1177 } 1178 #endif 1179 PetscFunctionReturn(0); 1180 } 1181 EXTERN_C_END 1182