[petsc-dev] ML Reuse redux

[John Fettig]

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);