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