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