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 #endif 584 if(level==0) Aarr[0] = Amat; /* use Pmat for finest level setup, but use mat for solver */ 585 if( isOK ) { 586 #if defined PETSC_USE_LOG 587 ierr = PetscLogEventBegin(gamg_setup_events[SET2],0,0,0,0);CHKERRQ(ierr); 588 #endif 589 ierr = createLevel( a_pc, Aarr[level], (pc_gamg->m_method != 0) ? bs : 1, pc_gamg->m_data_cols, bs, 590 (PetscBool)(level == pc_gamg->m_Nlevels-2), 591 &Parr[level1], &coarse_data, &nactivepe, &Aarr[level1] ); 592 CHKERRQ(ierr); 593 #if defined PETSC_USE_LOG 594 ierr = PetscLogEventEnd(gamg_setup_events[SET2],0,0,0,0);CHKERRQ(ierr); 595 #endif 596 ierr = MatGetSize( Aarr[level1], &M, &N );CHKERRQ(ierr); 597 ierr = MatGetInfo(Aarr[level1],MAT_GLOBAL_SUM,&info); CHKERRQ(ierr); 598 if (pc_gamg->m_verbose) PetscPrintf(PETSC_COMM_WORLD,"\t\t[%d]%s %d) N=%d, n data cols=%d, nnz/row (ave)=%d, %d active pes\n", 599 mype,__FUNCT__,(int)level1,N,pc_gamg->m_data_cols, 600 (int)(info.nz_used/(PetscReal)N),nactivepe); 601 /* coarse grids with SA can have zero row/cols from singleton aggregates */ 602 /* aggregation method should gaurrentee this does not happen! */ 603 604 if (pc_gamg->m_verbose) { 605 Vec diag; PetscScalar *data_arr,v; PetscInt Istart,Iend,kk,nloceq,id; 606 v = 1.e-10; /* LU factor has hard wired numbers for small diags so this needs to match (yuk) */ 607 ierr = MatGetOwnershipRange(Aarr[level1], &Istart, &Iend); CHKERRQ(ierr); 608 nloceq = Iend-Istart; 609 ierr = MatGetVecs( Aarr[level1], &diag, 0 ); CHKERRQ(ierr); 610 ierr = MatGetDiagonal( Aarr[level1], diag ); CHKERRQ(ierr); 611 ierr = VecGetArray( diag, &data_arr ); CHKERRQ(ierr); 612 for(kk=0;kk<nloceq;kk++){ 613 if(data_arr[kk]==0.0) { 614 id = kk + Istart; 615 ierr = MatSetValues(Aarr[level1],1,&id,1,&id,&v,INSERT_VALUES); 616 CHKERRQ(ierr); 617 PetscPrintf(PETSC_COMM_SELF,"\t[%d]%s warning: added zero to diag (%d) on level %d \n",mype,__FUNCT__,id,level1); 618 } 619 } 620 ierr = VecRestoreArray( diag, &data_arr ); CHKERRQ(ierr); 621 ierr = VecDestroy( &diag ); CHKERRQ(ierr); 622 ierr = MatAssemblyBegin(Aarr[level1],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 623 ierr = MatAssemblyEnd(Aarr[level1],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 624 } 625 } else { 626 coarse_data = 0; 627 break; 628 } 629 data = coarse_data; 630 631 #if (defined PETSC_USE_LOG && defined GAMG_STAGES) 632 ierr = PetscLogStagePop(); CHKERRQ( ierr ); 633 #endif 634 } 635 if( coarse_data ) { 636 ierr = PetscFree( coarse_data ); CHKERRQ( ierr ); 637 } 638 if (pc_gamg->m_verbose) PetscPrintf(PETSC_COMM_WORLD,"\t[%d]%s %d levels\n",0,__FUNCT__,level + 1); 639 pc_gamg->m_data = 0; /* destroyed coordinate data */ 640 pc_gamg->m_Nlevels = level + 1; 641 fine_level = level; 642 ierr = PCMGSetLevels(a_pc,pc_gamg->m_Nlevels,PETSC_NULL);CHKERRQ(ierr); 643 644 /* set default smoothers */ 645 for ( lidx=1, level = pc_gamg->m_Nlevels-2; 646 lidx <= fine_level; 647 lidx++, level--) { 648 PetscReal emax, emin; 649 KSP smoother; PC subpc; 650 PetscBool isCheb; 651 /* set defaults */ 652 ierr = PCMGGetSmoother( a_pc, lidx, &smoother ); CHKERRQ(ierr); 653 ierr = KSPSetType( smoother, KSPCHEBYCHEV );CHKERRQ(ierr); 654 ierr = KSPGetPC( smoother, &subpc ); CHKERRQ(ierr); 655 /* ierr = KSPSetTolerances(smoother,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,2); CHKERRQ(ierr); */ 656 ierr = PCSetType( subpc, PETSC_GAMG_SMOOTHER ); CHKERRQ(ierr); 657 ierr = KSPSetNormType( smoother, KSP_NORM_NONE ); CHKERRQ(ierr); 658 /* overide defaults with input parameters */ 659 ierr = KSPSetFromOptions( smoother ); CHKERRQ(ierr); 660 661 ierr = KSPSetOperators( smoother, Aarr[level], Aarr[level], SAME_NONZERO_PATTERN ); CHKERRQ(ierr); 662 /* do my own cheby */ 663 ierr = PetscTypeCompare( (PetscObject)smoother, KSPCHEBYCHEV, &isCheb ); CHKERRQ(ierr); 664 if( isCheb ) { 665 ierr = PetscTypeCompare( (PetscObject)subpc, PETSC_GAMG_SMOOTHER, &isCheb ); CHKERRQ(ierr); 666 if( isCheb && emaxs[level] > 0.0 ) emax=emaxs[level]; /* eigen estimate only for diagnal PC */ 667 else{ /* eigen estimate 'emax' */ 668 KSP eksp; Mat Lmat = Aarr[level]; 669 Vec bb, xx; PC pc; 670 const PCType type; 671 672 ierr = PCGetType( subpc, &type ); CHKERRQ(ierr); 673 ierr = MatGetVecs( Lmat, &bb, 0 ); CHKERRQ(ierr); 674 ierr = MatGetVecs( Lmat, &xx, 0 ); CHKERRQ(ierr); 675 { 676 PetscRandom rctx; 677 ierr = PetscRandomCreate(wcomm,&rctx);CHKERRQ(ierr); 678 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 679 ierr = VecSetRandom(bb,rctx);CHKERRQ(ierr); 680 ierr = PetscRandomDestroy( &rctx ); CHKERRQ(ierr); 681 } 682 ierr = KSPCreate(wcomm,&eksp);CHKERRQ(ierr); 683 ierr = KSPSetType( eksp, KSPCG ); CHKERRQ(ierr); 684 ierr = KSPSetTolerances( eksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 10 ); 685 CHKERRQ(ierr); 686 ierr = KSPSetNormType( eksp, KSP_NORM_NONE ); CHKERRQ(ierr); 687 688 ierr = KSPAppendOptionsPrefix( eksp, "est_"); CHKERRQ(ierr); 689 ierr = KSPSetFromOptions( eksp ); CHKERRQ(ierr); 690 691 ierr = KSPSetInitialGuessNonzero( eksp, PETSC_FALSE ); CHKERRQ(ierr); 692 ierr = KSPSetOperators( eksp, Lmat, Lmat, SAME_NONZERO_PATTERN ); CHKERRQ( ierr ); 693 ierr = KSPGetPC( eksp, &pc );CHKERRQ( ierr ); 694 ierr = PCSetType( pc, type ); CHKERRQ(ierr); /* should be same as eigen estimates op. */ 695 696 ierr = KSPSetComputeSingularValues( eksp,PETSC_TRUE ); CHKERRQ(ierr); 697 ierr = KSPSolve( eksp, bb, xx ); CHKERRQ(ierr); 698 ierr = KSPComputeExtremeSingularValues( eksp, &emax, &emin ); CHKERRQ(ierr); 699 ierr = VecDestroy( &xx ); CHKERRQ(ierr); 700 ierr = VecDestroy( &bb ); CHKERRQ(ierr); 701 ierr = KSPDestroy( &eksp ); CHKERRQ(ierr); 702 703 if (pc_gamg->m_verbose) { 704 PetscPrintf(PETSC_COMM_WORLD,"\t\t\t%s PC setup max eigen=%e min=%e PC=%s\n", 705 __FUNCT__,emax,emin,PETSC_GAMG_SMOOTHER); 706 } 707 } 708 { 709 PetscInt N1, N0, tt; 710 ierr = MatGetSize( Aarr[level], &N1, &tt ); CHKERRQ(ierr); 711 ierr = MatGetSize( Aarr[level+1], &N0, &tt ); CHKERRQ(ierr); 712 /* heuristic - is this crap? */ 713 emin = 1.*emax/((PetscReal)N1/(PetscReal)N0); 714 emax *= 1.05; 715 } 716 ierr = KSPChebychevSetEigenvalues( smoother, emax, emin );CHKERRQ(ierr); 717 } 718 } 719 { 720 /* coarse grid */ 721 KSP smoother,*k2; PC subpc,pc2; PetscInt ii,first; 722 Mat Lmat = Aarr[pc_gamg->m_Nlevels-1]; 723 ierr = PCMGGetSmoother( a_pc, 0, &smoother ); CHKERRQ(ierr); 724 ierr = KSPSetOperators( smoother, Lmat, Lmat, SAME_NONZERO_PATTERN ); CHKERRQ(ierr); 725 ierr = KSPSetNormType( smoother, KSP_NORM_NONE ); CHKERRQ(ierr); 726 ierr = KSPGetPC( smoother, &subpc ); CHKERRQ(ierr); 727 ierr = PCSetType( subpc, PCBJACOBI ); CHKERRQ(ierr); 728 ierr = PCSetUp( subpc ); CHKERRQ(ierr); 729 ierr = PCBJacobiGetSubKSP(subpc,&ii,&first,&k2);CHKERRQ(ierr); 730 assert(ii==1); 731 ierr = KSPGetPC(k2[0],&pc2);CHKERRQ(ierr); 732 ierr = PCSetType( pc2, PCLU ); CHKERRQ(ierr); 733 } 734 735 /* should be called in PCSetFromOptions_GAMG(), but cannot be called prior to PCMGSetLevels() */ 736 ierr = PCSetFromOptions_MG(a_pc); CHKERRQ(ierr); 737 { 738 PetscBool galerkin; 739 ierr = PCMGGetGalerkin( a_pc, &galerkin); CHKERRQ(ierr); 740 if(galerkin){ 741 SETERRQ(wcomm,PETSC_ERR_ARG_WRONG, "GAMG does galerkin manually so it must not be used in PC_MG."); 742 } 743 } 744 745 /* plot levels - R/P */ 746 if( PRINT_MATS ) { 747 for (level=pc_gamg->m_Nlevels-1;level>0;level--){ 748 PetscViewer viewer; 749 char fname[32]; 750 sprintf(fname,"Pmat_%d_%d.m",(int)pc_gamg->m_count,(int)level); 751 ierr = PetscViewerASCIIOpen( wcomm, fname, &viewer ); CHKERRQ(ierr); 752 ierr = PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB); CHKERRQ(ierr); 753 ierr = MatView( Parr[level], viewer ); CHKERRQ(ierr); 754 ierr = PetscViewerDestroy( &viewer ); 755 sprintf(fname,"Amat_%d_%d.m",(int)pc_gamg->m_count,(int)level); 756 ierr = PetscViewerASCIIOpen( wcomm, fname, &viewer ); CHKERRQ(ierr); 757 ierr = PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB); CHKERRQ(ierr); 758 ierr = MatView( Aarr[level], viewer ); CHKERRQ(ierr); 759 ierr = PetscViewerDestroy( &viewer ); 760 } 761 } 762 763 /* set interpolation between the levels, clean up */ 764 for (lidx=0,level=pc_gamg->m_Nlevels-1; 765 lidx<fine_level; 766 lidx++, level--){ 767 ierr = PCMGSetInterpolation( a_pc, lidx+1, Parr[level] );CHKERRQ(ierr); 768 ierr = MatDestroy( &Parr[level] ); CHKERRQ(ierr); 769 ierr = MatDestroy( &Aarr[level] ); CHKERRQ(ierr); 770 } 771 772 /* setupcalled is set to 0 so that MG is setup from scratch */ 773 a_pc->setupcalled = 0; 774 ierr = PCSetUp_MG( a_pc );CHKERRQ( ierr ); 775 a_pc->setupcalled = 1; /* use 1 as signal that this has not been re-setup */ 776 777 { 778 KSP smoother; /* PCSetUp_MG seems to insists on setting this to GMRES on coarse grid */ 779 ierr = PCMGGetSmoother( a_pc, 0, &smoother ); CHKERRQ(ierr); 780 ierr = KSPSetType( smoother, KSPPREONLY ); CHKERRQ(ierr); 781 ierr = KSPSetUp( smoother ); CHKERRQ(ierr); 782 } 783 784 PetscFunctionReturn(0); 785 } 786 787 /* ------------------------------------------------------------------------- */ 788 /* 789 PCDestroy_GAMG - Destroys the private context for the GAMG preconditioner 790 that was created with PCCreate_GAMG(). 791 792 Input Parameter: 793 . pc - the preconditioner context 794 795 Application Interface Routine: PCDestroy() 796 */ 797 #undef __FUNCT__ 798 #define __FUNCT__ "PCDestroy_GAMG" 799 PetscErrorCode PCDestroy_GAMG(PC pc) 800 { 801 PetscErrorCode ierr; 802 PC_MG *mg = (PC_MG*)pc->data; 803 PC_GAMG *pc_gamg= (PC_GAMG*)mg->innerctx; 804 805 PetscFunctionBegin; 806 ierr = PCReset_GAMG(pc);CHKERRQ(ierr); 807 ierr = PetscFree(pc_gamg);CHKERRQ(ierr); 808 ierr = PCDestroy_MG(pc);CHKERRQ(ierr); 809 PetscFunctionReturn(0); 810 } 811 812 813 #undef __FUNCT__ 814 #define __FUNCT__ "PCGAMGSetProcEqLim" 815 /*@ 816 PCGAMGSetProcEqLim - Set number of equations to aim for on coarse grids via 817 processor reduction. 818 819 Not Collective on PC 820 821 Input Parameters: 822 . pc - the preconditioner context 823 824 825 Options Database Key: 826 . -pc_gamg_process_eq_limit 827 828 Level: intermediate 829 830 Concepts: Unstructured multrigrid preconditioner 831 832 .seealso: () 833 @*/ 834 PetscErrorCode PCGAMGSetProcEqLim(PC pc, PetscInt n) 835 { 836 PetscErrorCode ierr; 837 838 PetscFunctionBegin; 839 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 840 ierr = PetscTryMethod(pc,"PCGAMGSetProcEqLim_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr); 841 PetscFunctionReturn(0); 842 } 843 844 EXTERN_C_BEGIN 845 #undef __FUNCT__ 846 #define __FUNCT__ "PCGAMGSetProcEqLim_GAMG" 847 PetscErrorCode PCGAMGSetProcEqLim_GAMG(PC pc, PetscInt n) 848 { 849 PC_MG *mg = (PC_MG*)pc->data; 850 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 851 852 PetscFunctionBegin; 853 if(n>0) pc_gamg->m_min_eq_proc = n; 854 PetscFunctionReturn(0); 855 } 856 EXTERN_C_END 857 858 #undef __FUNCT__ 859 #define __FUNCT__ "PCGAMGSetCoarseEqLim" 860 /*@ 861 PCGAMGSetCoarseEqLim - Set max number of equations on coarse grids. 862 863 Collective on PC 864 865 Input Parameters: 866 . pc - the preconditioner context 867 868 869 Options Database Key: 870 . -pc_gamg_coarse_eq_limit 871 872 Level: intermediate 873 874 Concepts: Unstructured multrigrid preconditioner 875 876 .seealso: () 877 @*/ 878 PetscErrorCode PCGAMGSetCoarseEqLim(PC pc, PetscInt n) 879 { 880 PetscErrorCode ierr; 881 882 PetscFunctionBegin; 883 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 884 ierr = PetscTryMethod(pc,"PCGAMGSetCoarseEqLim_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr); 885 PetscFunctionReturn(0); 886 } 887 888 EXTERN_C_BEGIN 889 #undef __FUNCT__ 890 #define __FUNCT__ "PCGAMGSetCoarseEqLim_GAMG" 891 PetscErrorCode PCGAMGSetCoarseEqLim_GAMG(PC pc, PetscInt n) 892 { 893 PC_MG *mg = (PC_MG*)pc->data; 894 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 895 896 PetscFunctionBegin; 897 if(n>0) pc_gamg->m_coarse_eq_limit = n; 898 PetscFunctionReturn(0); 899 } 900 EXTERN_C_END 901 902 #undef __FUNCT__ 903 #define __FUNCT__ "PCGAMGSetRepartitioning" 904 /*@ 905 PCGAMGSetRepartitioning - Repartition the coarse grids 906 907 Collective on PC 908 909 Input Parameters: 910 . pc - the preconditioner context 911 912 913 Options Database Key: 914 . -pc_gamg_repartition 915 916 Level: intermediate 917 918 Concepts: Unstructured multrigrid preconditioner 919 920 .seealso: () 921 @*/ 922 PetscErrorCode PCGAMGSetRepartitioning(PC pc, PetscBool n) 923 { 924 PetscErrorCode ierr; 925 926 PetscFunctionBegin; 927 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 928 ierr = PetscTryMethod(pc,"PCGAMGSetRepartitioning_C",(PC,PetscBool),(pc,n));CHKERRQ(ierr); 929 PetscFunctionReturn(0); 930 } 931 932 EXTERN_C_BEGIN 933 #undef __FUNCT__ 934 #define __FUNCT__ "PCGAMGSetRepartitioning_GAMG" 935 PetscErrorCode PCGAMGSetRepartitioning_GAMG(PC pc, PetscBool n) 936 { 937 PC_MG *mg = (PC_MG*)pc->data; 938 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 939 940 PetscFunctionBegin; 941 pc_gamg->m_repart = n; 942 PetscFunctionReturn(0); 943 } 944 EXTERN_C_END 945 946 #undef __FUNCT__ 947 #define __FUNCT__ "PCGAMGSetNlevels" 948 /*@ 949 PCGAMGSetNlevels - 950 951 Not collective on PC 952 953 Input Parameters: 954 . pc - the preconditioner context 955 956 957 Options Database Key: 958 . -pc_mg_levels 959 960 Level: intermediate 961 962 Concepts: Unstructured multrigrid preconditioner 963 964 .seealso: () 965 @*/ 966 PetscErrorCode PCGAMGSetNlevels(PC pc, PetscInt n) 967 { 968 PetscErrorCode ierr; 969 970 PetscFunctionBegin; 971 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 972 ierr = PetscTryMethod(pc,"PCGAMGSetNlevels_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr); 973 PetscFunctionReturn(0); 974 } 975 976 EXTERN_C_BEGIN 977 #undef __FUNCT__ 978 #define __FUNCT__ "PCGAMGSetNlevels_GAMG" 979 PetscErrorCode PCGAMGSetNlevels_GAMG(PC pc, PetscInt n) 980 { 981 PC_MG *mg = (PC_MG*)pc->data; 982 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 983 984 PetscFunctionBegin; 985 pc_gamg->m_Nlevels = n; 986 PetscFunctionReturn(0); 987 } 988 EXTERN_C_END 989 990 #undef __FUNCT__ 991 #define __FUNCT__ "PCGAMGSetThreshold" 992 /*@ 993 PCGAMGSetThreshold - Relative threshold to use for dropping edges in aggregation graph 994 995 Not collective on PC 996 997 Input Parameters: 998 . pc - the preconditioner context 999 1000 1001 Options Database Key: 1002 . -pc_gamg_threshold 1003 1004 Level: intermediate 1005 1006 Concepts: Unstructured multrigrid preconditioner 1007 1008 .seealso: () 1009 @*/ 1010 PetscErrorCode PCGAMGSetThreshold(PC pc, PetscReal n) 1011 { 1012 PetscErrorCode ierr; 1013 1014 PetscFunctionBegin; 1015 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1016 ierr = PetscTryMethod(pc,"PCGAMGSetThreshold_C",(PC,PetscReal),(pc,n));CHKERRQ(ierr); 1017 PetscFunctionReturn(0); 1018 } 1019 1020 EXTERN_C_BEGIN 1021 #undef __FUNCT__ 1022 #define __FUNCT__ "PCGAMGSetThreshold_GAMG" 1023 PetscErrorCode PCGAMGSetThreshold_GAMG(PC pc, PetscReal n) 1024 { 1025 PC_MG *mg = (PC_MG*)pc->data; 1026 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1027 1028 PetscFunctionBegin; 1029 pc_gamg->m_threshold = n; 1030 PetscFunctionReturn(0); 1031 } 1032 EXTERN_C_END 1033 1034 #undef __FUNCT__ 1035 #define __FUNCT__ "PCGAMGSetSolverType" 1036 /*@ 1037 PCGAMGSetSolverType - Set solution method. 1038 1039 Collective on PC 1040 1041 Input Parameters: 1042 . pc - the preconditioner context 1043 1044 1045 Options Database Key: 1046 . -pc_gamg_type 1047 1048 Level: intermediate 1049 1050 Concepts: Unstructured multrigrid preconditioner 1051 1052 .seealso: () 1053 @*/ 1054 PetscErrorCode PCGAMGSetSolverType(PC pc, char str[], PetscInt sz ) 1055 { 1056 PetscErrorCode ierr; 1057 1058 PetscFunctionBegin; 1059 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1060 ierr = PetscTryMethod(pc,"PCGAMGSetSolverType_C",(PC,char[],PetscInt),(pc,str,sz)); 1061 CHKERRQ(ierr); 1062 PetscFunctionReturn(0); 1063 } 1064 1065 EXTERN_C_BEGIN 1066 #undef __FUNCT__ 1067 #define __FUNCT__ "PCGAMGSetSolverType_GAMG" 1068 PetscErrorCode PCGAMGSetSolverType_GAMG(PC pc, char str[], PetscInt sz ) 1069 { 1070 PC_MG *mg = (PC_MG*)pc->data; 1071 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1072 1073 PetscFunctionBegin; 1074 if(sz < 64) strcpy(pc_gamg->m_type,str); 1075 PetscFunctionReturn(0); 1076 } 1077 EXTERN_C_END 1078 1079 #undef __FUNCT__ 1080 #define __FUNCT__ "PCSetFromOptions_GAMG" 1081 PetscErrorCode PCSetFromOptions_GAMG(PC pc) 1082 { 1083 PetscErrorCode ierr; 1084 PC_MG *mg = (PC_MG*)pc->data; 1085 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1086 PetscBool flag; 1087 1088 PetscFunctionBegin; 1089 ierr = PetscOptionsHead("GAMG options"); CHKERRQ(ierr); 1090 { 1091 ierr = PetscOptionsString("-pc_gamg_type", 1092 "Solver type: plane aggregation ('pa'), smoothed aggregation ('sa') or geometric multigrid (default)", 1093 "PCGAMGSetSolverType", 1094 pc_gamg->m_type, 1095 pc_gamg->m_type, 1096 64, 1097 &flag ); 1098 CHKERRQ(ierr); 1099 if( flag && pc_gamg->m_data != 0 ) { 1100 if( (strcmp(pc_gamg->m_type,"sa")==0 && pc_gamg->m_method != 2) || 1101 (strcmp(pc_gamg->m_type,"pa")==0 && pc_gamg->m_method != 1) || 1102 (strcmp(pc_gamg->m_type,"geo")==0 && pc_gamg->m_method != 0) ) { 1103 SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "PCSetFromOptions called of PCSetCoordinates (with new method, after data was created)"); 1104 } 1105 } 1106 1107 /* -pc_gamg_verbose */ 1108 ierr = PetscOptionsBool("-pc_gamg_verbose","Verbose (debugging) output for PCGAMG","none",pc_gamg->m_verbose,&pc_gamg->m_verbose,PETSC_NULL);CHKERRQ(ierr); 1109 1110 pc_gamg->m_method = 1; /* default to plane aggregation */ 1111 if (flag ) { 1112 if( strcmp(pc_gamg->m_type,"sa") == 0) pc_gamg->m_method = 2; 1113 else if( strcmp(pc_gamg->m_type,"pa") == 0) pc_gamg->m_method = 1; 1114 else if( strcmp(pc_gamg->m_type,"geo") == 0) pc_gamg->m_method = 0; 1115 else SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "Invalid gamg type: %s",pc_gamg->m_type); 1116 } 1117 /* -pc_gamg_repartition */ 1118 pc_gamg->m_repart = PETSC_FALSE; 1119 ierr = PetscOptionsBool("-pc_gamg_repartition", 1120 "Repartion coarse grids (false)", 1121 "PCGAMGRepartitioning", 1122 pc_gamg->m_repart, 1123 &pc_gamg->m_repart, 1124 &flag); 1125 CHKERRQ(ierr); 1126 1127 /* -pc_gamg_process_eq_limit */ 1128 pc_gamg->m_min_eq_proc = 600; 1129 ierr = PetscOptionsInt("-pc_gamg_process_eq_limit", 1130 "Limit (goal) on number of equations per process on coarse grids", 1131 "PCGAMGSetProcEqLim", 1132 pc_gamg->m_min_eq_proc, 1133 &pc_gamg->m_min_eq_proc, 1134 &flag ); 1135 CHKERRQ(ierr); 1136 1137 /* -pc_gamg_coarse_eq_limit */ 1138 pc_gamg->m_coarse_eq_limit = 1500; 1139 ierr = PetscOptionsInt("-pc_gamg_coarse_eq_limit", 1140 "Limit on number of equations for the coarse grid", 1141 "PCGAMGSetCoarseEqLim", 1142 pc_gamg->m_coarse_eq_limit, 1143 &pc_gamg->m_coarse_eq_limit, 1144 &flag ); 1145 CHKERRQ(ierr); 1146 1147 /* -pc_gamg_threshold */ 1148 pc_gamg->m_threshold = 0.05; 1149 ierr = PetscOptionsReal("-pc_gamg_threshold", 1150 "Relative threshold to use for dropping edges in aggregation graph", 1151 "PCGAMGSetThreshold", 1152 pc_gamg->m_threshold, 1153 &pc_gamg->m_threshold, 1154 &flag ); 1155 CHKERRQ(ierr); 1156 1157 pc_gamg->m_Nlevels = GAMG_MAXLEVELS; 1158 ierr = PetscOptionsInt("-pc_mg_levels", 1159 "Set number of MG levels", 1160 "PCGAMGSetNlevels", 1161 pc_gamg->m_Nlevels, 1162 &pc_gamg->m_Nlevels, 1163 &flag ); 1164 } 1165 ierr = PetscOptionsTail();CHKERRQ(ierr); 1166 1167 PetscFunctionReturn(0); 1168 } 1169 1170 /* -------------------------------------------------------------------------- */ 1171 /*MC 1172 PCGAMG - Geometric algebraic multigrid (AMG) preconditioning. This preconditioner currently has two 1173 AMG methods: 1) an unstructured geometric method, which requires that you provide coordinates for each 1174 vertex, and 2) smoothed aggregation. Smoothed aggregation (SA) can work without coordinates but it 1175 will generate some common non-trivial null spaces if coordinates are provided. The input fine grid matrix 1176 must have the block size set for 'system' problems (with multiple dof per vertex/cell) to work properly. 1177 SA will generate rotational rigid body mode null space vectors, in addition to the trivial translational 1178 modes, when coordinates are provide in 2D and 3D. 1179 1180 Options Database Keys: 1181 Multigrid options(inherited) 1182 + -pc_mg_cycles <1>: 1 for V cycle, 2 for W-cycle (PCMGSetCycleType) 1183 . -pc_mg_smoothup <1>: Number of post-smoothing steps (PCMGSetNumberSmoothUp) 1184 . -pc_mg_smoothdown <1>: Number of pre-smoothing steps (PCMGSetNumberSmoothDown) 1185 - -pc_mg_type <multiplicative>: (one of) additive multiplicative full cascade kascade 1186 1187 Level: intermediate 1188 1189 Concepts: multigrid 1190 1191 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, 1192 PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(), PCMGSetNumberSmoothDown(), 1193 PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(), 1194 PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(), 1195 PCMGSetCyclesOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR() 1196 M*/ 1197 1198 EXTERN_C_BEGIN 1199 #undef __FUNCT__ 1200 #define __FUNCT__ "PCCreate_GAMG" 1201 PetscErrorCode PCCreate_GAMG(PC pc) 1202 { 1203 PetscErrorCode ierr; 1204 PC_GAMG *pc_gamg; 1205 PC_MG *mg; 1206 PetscClassId cookie; 1207 #if defined PETSC_USE_LOG 1208 static int count = 0; 1209 #endif 1210 1211 PetscFunctionBegin; 1212 /* PCGAMG is an inherited class of PCMG. Initialize pc as PCMG */ 1213 ierr = PCSetType(pc,PCMG);CHKERRQ(ierr); /* calls PCCreate_MG() and MGCreate_Private() */ 1214 ierr = PetscObjectChangeTypeName((PetscObject)pc,PCGAMG);CHKERRQ(ierr); 1215 1216 /* create a supporting struct and attach it to pc */ 1217 ierr = PetscNewLog(pc,PC_GAMG,&pc_gamg);CHKERRQ(ierr); 1218 pc_gamg->m_data_sz = 0; pc_gamg->m_data = 0; pc_gamg->m_count = 0; 1219 mg = (PC_MG*)pc->data; 1220 mg->innerctx = pc_gamg; 1221 1222 pc_gamg->m_Nlevels = -1; 1223 1224 /* overwrite the pointers of PCMG by the functions of PCGAMG */ 1225 pc->ops->setfromoptions = PCSetFromOptions_GAMG; 1226 pc->ops->setup = PCSetUp_GAMG; 1227 pc->ops->reset = PCReset_GAMG; 1228 pc->ops->destroy = PCDestroy_GAMG; 1229 1230 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1231 "PCSetCoordinates_C", 1232 "PCSetCoordinates_GAMG", 1233 PCSetCoordinates_GAMG); 1234 CHKERRQ(ierr); 1235 1236 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1237 "PCGAMGSetProcEqLim_C", 1238 "PCGAMGSetProcEqLim_GAMG", 1239 PCGAMGSetProcEqLim_GAMG); 1240 CHKERRQ(ierr); 1241 1242 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1243 "PCGAMGSetCoarseEqLim_C", 1244 "PCGAMGSetCoarseEqLim_GAMG", 1245 PCGAMGSetCoarseEqLim_GAMG); 1246 CHKERRQ(ierr); 1247 1248 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1249 "PCGAMGSetRepartitioning_C", 1250 "PCGAMGSetRepartitioning_GAMG", 1251 PCGAMGSetRepartitioning_GAMG); 1252 CHKERRQ(ierr); 1253 1254 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1255 "PCGAMGSetThreshold_C", 1256 "PCGAMGSetThreshold_GAMG", 1257 PCGAMGSetThreshold_GAMG); 1258 CHKERRQ(ierr); 1259 1260 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1261 "PCGAMGSetSolverType_C", 1262 "PCGAMGSetSolverType_GAMG", 1263 PCGAMGSetSolverType_GAMG); 1264 CHKERRQ(ierr); 1265 1266 #if defined PETSC_USE_LOG 1267 if( count++ == 0 ) { 1268 PetscClassIdRegister("GAMG Setup",&cookie); 1269 PetscLogEventRegister("GAMG: createProl", cookie, &gamg_setup_events[SET1]); 1270 PetscLogEventRegister(" Graph", cookie, &gamg_setup_events[GRAPH]); 1271 PetscLogEventRegister(" G.Mat", cookie, &gamg_setup_events[GRAPH_MAT]); 1272 PetscLogEventRegister(" G.Filter", cookie, &gamg_setup_events[GRAPH_FILTER]); 1273 PetscLogEventRegister(" G.Square", cookie, &gamg_setup_events[GRAPH_SQR]); 1274 PetscLogEventRegister(" MIS/Agg", cookie, &gamg_setup_events[SET4]); 1275 PetscLogEventRegister(" geo: growSupp", cookie, &gamg_setup_events[SET5]); 1276 PetscLogEventRegister(" geo: triangle", cookie, &gamg_setup_events[SET6]); 1277 PetscLogEventRegister(" search&set", cookie, &gamg_setup_events[FIND_V]); 1278 PetscLogEventRegister(" SA: init", cookie, &gamg_setup_events[SET7]); 1279 /* PetscLogEventRegister(" SA: frmProl0", cookie, &gamg_setup_events[SET8]); */ 1280 PetscLogEventRegister(" SA: smooth", cookie, &gamg_setup_events[SET9]); 1281 PetscLogEventRegister("GAMG: partLevel", cookie, &gamg_setup_events[SET2]); 1282 PetscLogEventRegister(" PL repartition", cookie, &gamg_setup_events[SET12]); 1283 1284 /* PetscLogEventRegister(" PL 1", cookie, &gamg_setup_events[SET13]); */ 1285 /* PetscLogEventRegister(" PL 2", cookie, &gamg_setup_events[SET14]); */ 1286 /* PetscLogEventRegister(" PL 3", cookie, &gamg_setup_events[SET15]); */ 1287 1288 /* PetscLogEventRegister(" PL move data", cookie, &gamg_setup_events[SET13]); */ 1289 /* PetscLogEventRegister("GAMG: fix", cookie, &gamg_setup_events[SET10]); */ 1290 /* PetscLogEventRegister("GAMG: set levels", cookie, &gamg_setup_events[SET11]); */ 1291 1292 /* create timer stages */ 1293 #if defined GAMG_STAGES 1294 { 1295 char str[32]; 1296 sprintf(str,"MG Level %d (finest)",0); 1297 PetscLogStageRegister(str, &gamg_stages[0]); 1298 PetscInt lidx; 1299 for (lidx=1;lidx<9;lidx++){ 1300 sprintf(str,"MG Level %d",lidx); 1301 PetscLogStageRegister(str, &gamg_stages[lidx]); 1302 } 1303 } 1304 #endif 1305 } 1306 #endif 1307 PetscFunctionReturn(0); 1308 } 1309 EXTERN_C_END 1310