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