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