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