On Fri, Aug 26, 2011 at 3:50 PM, John Fettig <john.fettig at gmail.com> wrote: > On Fri, Aug 26, 2011 at 2:27 PM, Jed Brown <jedbrown at mcs.anl.gov> wrote: >> With ML, there is a technical problem that we don't get the interpolation >> back in a PETSc Mat, so we can't use PETSc MatPtAP(). I don't know if they >> have a similar function that is callable by users, or another way to only >> update the Galerkin operators without redoing the smoothed aggregation. > > There must be a way for them to do it.
I think the way to do it is similar to the way it is done in the function ML_Gen_MultiLevelHierarchy_UsingSmoothedAggr_ReuseExistingAgg, although that function rebuilds R, A, and P. So we just have to rebuild A. Here's some code that does that, it seems to work as well as the previous code I sent that only worked in serial (this works in parallel too). Could this be made to be the default if you set the flag "SAME_NONZERO_PATTERN"? Then if you want to rebuild the preconditioner, call with "DIFFERENT_NONZERO_PATTERN". Regards, John Nlevels = pc_ml->Nlevels; fine_level = Nlevels - 1; ml_object = pc_ml->ml_object; gridctx = pc_ml->gridctx; gridctx[fine_level].A = A; ierr = PetscTypeCompare((PetscObject) A, MATSEQAIJ, &isSeq);CHKERRQ(ierr); ierr = PetscTypeCompare((PetscObject) A, MATMPIAIJ, &isMPI);CHKERRQ(ierr); if (isMPI){ ierr = MatConvert_MPIAIJ_ML(A,PETSC_NULL,MAT_INITIAL_MATRIX,&Aloc);CHKERRQ(ierr); } else if (isSeq) { Aloc = A; } else SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "Matrix type '%s' cannot be used with ML. ML can only handle AIJ matrices.",((PetscObject)A)->type_name); ierr = MatGetSize(Aloc,&m,&nlocal_allcols);CHKERRQ(ierr); PetscMLdata = pc_ml->PetscMLdata; PetscMLdata->A = A; PetscMLdata->Aloc = Aloc; ML_Init_Amatrix(ml_object,0,m,m,PetscMLdata); ML_Set_Amatrix_Matvec(ml_object,0,PetscML_matvec); mesh_level = ml_object->ML_finest_level; while( ml_object->SingleLevel[mesh_level].Rmat->to != NULL) { old_mesh_level = mesh_level; mesh_level = ml_object->SingleLevel[mesh_level].Rmat->to->levelnum; // clean and regenerate A mlmat = &(ml_object->Amat[mesh_level]); ML_Operator_Clean(mlmat); ML_Operator_Init(mlmat,ml_object->comm); ML_Gen_AmatrixRAP(ml_object, old_mesh_level, mesh_level); } level = fine_level - 1; if (size == 1){ // convert ML P, R and A into seqaij format for (mllevel=1; mllevel<Nlevels; mllevel++){ mlmat = &(ml_object->Amat[mllevel]); ierr = MatWrapML_SeqAIJ(mlmat,MAT_INITIAL_MATRIX,&gridctx[level].A);CHKERRQ(ierr); level--; } } else { /* convert ML P and R into shell format, ML A into mpiaij format */ for (mllevel=1; mllevel<Nlevels; mllevel++){ mlmat = &(ml_object->Amat[mllevel]); ierr = MatWrapML_MPIAIJ(mlmat,&gridctx[level].A);CHKERRQ(ierr); level--; } } for (level=0; level<fine_level; level++){ if (level > 0){ ierr = PCMGSetResidual(pc,level,PCMGDefaultResidual,gridctx[level].A);CHKERRQ(ierr); } ierr = KSPSetOperators(gridctx[level].ksp,gridctx[level].A,gridctx[level].A,SAME_NONZERO_PATTERN);CHKERRQ(ierr); } ierr = PCMGSetResidual(pc,fine_level,PCMGDefaultResidual,gridctx[fine_level].A);CHKERRQ(ierr); ierr = KSPSetOperators(gridctx[fine_level].ksp,gridctx[level].A,gridctx[fine_level].A,SAME_NONZERO_PATTERN);CHKERRQ(ierr); // setupcalled is set to 0 so that MG is setup from scratch pc->setupcalled = 0; ierr = PCSetUp_MG(pc);CHKERRQ(ierr); PetscFunctionReturn(0);