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