pc.c

Go to the documentation of this file.

/*! \page pc The preconditioner

  Choosing a preconditioner changes a linear system into a preconditioned
  system. However, this is not a transformation of any coefficient matrix,
  so this preprocessor is handled a bit differently from the previous ones.
 */

#include <stdlib.h>
#include <stdio.h>
#include "syspro.h"
#include "syspro_impl.h"
#include "sysprolinear.h"
#include "sysprotransform.h"
#include "linear_impl.h"
#include "linpc.h"
#include "linksp.h"
#include "petsc.h"
#include "petscmat.h"
#include "petscpc.h"
#include "petscksp.h"
#include "anamod.h"

#define PetscObjectStateDecrease(obj) ((obj)->state--,0)

#define PREPROCESSOR "pc"

#undef __FUNCT__
#define __FUNCT__ "setup_pc_choices"
static PetscErrorCode setup_pc_choices()
{
  SalsaTransform pc; SalsaTransformObject cur;
  PetscErrorCode ierr;

  PetscFunctionBegin;

  ierr = TransformGetByName(PREPROCESSOR,&pc); CHKERRQ(ierr);

  /*
   * Set up various properties of PCs; these can be used
   * later to disable KSP choices
   */
  /* symmetric substitutions */
  ierr = SysProDefineIntAnnotation(PREPROCESSOR,"notrans"); CHKERRQ(ierr);
  /* requires transpose */
  ierr = SysProDefineIntAnnotation(PREPROCESSOR,"noparallel"); CHKERRQ(ierr);
  /* is direct solve */
  ierr = SysProDefineIntAnnotation(PREPROCESSOR,"directsolve"); CHKERRQ(ierr);
  /* marked as breaking down */
  ierr = SysProDefineIntAnnotation(PREPROCESSOR,"marked"); CHKERRQ(ierr);
  
#if defined(PETSC_HAVE_BLOCKSOLVE)
  /*
   * BlockSolve 95
   */
  ierr = NewTransformObject(PREPROCESSOR,PCBS95,&cur); CHKERRQ(ierr);
  ierr = TransformObjectSetExplanation(cur,"BlockSolve 95"); CHKERRQ(ierr);
  ierr = TransformObjectIntAnnotate(cur,"notrans",1); CHKERRQ(ierr);
#endif
  
#if defined(PETSC_HAVE_HYPRE)
  /*
   * Hypre methods: boomeramg and pilut
   */
  ierr = NewTransformObject(PREPROCESSOR,PCBOOMERAMG,&cur); CHKERRQ(ierr);
  ierr = TransformObjectSetExplanation(cur,"Hypre BoomerAMG"); CHKERRQ(ierr);
  ierr = TransformObjectIntAnnotate(cur,"notrans",1); CHKERRQ(ierr);

  ierr = NewTransformObject(PREPROCESSOR,PCPILUT,&cur); CHKERRQ(ierr);
  ierr = TransformObjectSetExplanation(cur,"Hypre pILUt"); CHKERRQ(ierr);
  ierr = TransformObjectIntAnnotate(cur,"notrans",1); CHKERRQ(ierr);

  ierr = NewTransformObject(PREPROCESSOR,PCEUCLID,&cur); CHKERRQ(ierr);
  ierr = TransformObjectSetExplanation(cur,"Hypre Euclid"); CHKERRQ(ierr);
  ierr = TransformObjectIntAnnotate(cur,"notrans",1); CHKERRQ(ierr);
  ierr = TransformObjectDefineOption(cur,"fill"); CHKERRQ(ierr);
  ierr = TransformObjectAddOption(cur,1); CHKERRQ(ierr);

  ierr = NewTransformObject(PREPROCESSOR,PCPARASAILS,&cur); CHKERRQ(ierr);
  ierr = TransformObjectSetExplanation(cur,"Hypre ParaSails"); CHKERRQ(ierr);
  ierr = TransformObjectIntAnnotate(cur,"notrans",1); CHKERRQ(ierr);
  ierr = TransformObjectDefineOption(cur,"fill"); CHKERRQ(ierr);
  ierr = TransformObjectAddOption(cur,1); CHKERRQ(ierr);
  ierr = TransformObjectAddOption(cur,2); CHKERRQ(ierr);
  ierr = TransformObjectAddOption(cur,3); CHKERRQ(ierr);
#endif

  /*
   * ILU sequential
   */
  ierr = NewTransformObject(PREPROCESSOR,PCILU,&cur); CHKERRQ(ierr);
  ierr = TransformObjectSetExplanation(cur,"Incomplete LU"); CHKERRQ(ierr);
  ierr = TransformObjectIntAnnotate(cur,"noparallel",1); CHKERRQ(ierr);
  ierr = TransformObjectDefineOption(cur,"fill"); CHKERRQ(ierr);
  ierr = TransformObjectAddOption(cur,1); CHKERRQ(ierr);
  ierr = TransformObjectAddOption(cur,2); CHKERRQ(ierr);
  ierr = TransformObjectAddOption(cur,3); CHKERRQ(ierr);
  /* shifted version */
  ierr = NewTransformObject(PREPROCESSOR,PCSILU,&cur); CHKERRQ(ierr);
  ierr = TransformObjectSetExplanation(cur,"shifted Incomplete LU"); CHKERRQ(ierr);
  ierr = TransformObjectIntAnnotate(cur,"noparallel",1); CHKERRQ(ierr);
  ierr = TransformObjectDefineOption(cur,"fill"); CHKERRQ(ierr);
  ierr = TransformObjectAddOption(cur,1); CHKERRQ(ierr);
  ierr = TransformObjectAddOption(cur,2); CHKERRQ(ierr);
  ierr = TransformObjectAddOption(cur,3); CHKERRQ(ierr);
  
#if 0
#if defined(PETSC_HAVE_SPAI)
  /*
   * SPAI
   */
  ierr = NewTransformObject(PREPROCESSOR,PCSPAI,&cur); CHKERRQ(ierr);
  ierr = TransformObjectSetExplanation(cur,"Sparse Approximate Inverse"); CHKERRQ(ierr);
  ierr = TransformObjectIntAnnotate(cur,"notrans",1); CHKERRQ(ierr);
  ierr = TransformObjectIntAnnotate(cur,"noparallel",1); CHKERRQ(ierr);
  ierr = TransformObjectDefineOption(cur,"blocksize"); CHKERRQ(ierr);
  ierr = TransformObjectAddOption(cur,1); CHKERRQ(ierr);
  ierr = TransformObjectAddOption(cur,2); CHKERRQ(ierr);
  ierr = TransformObjectAddOption(cur,3); CHKERRQ(ierr);
#endif
#endif

  /*
   * Additive Schwarz
   */
  ierr = NewTransformObject
    (PREPROCESSOR,PCASM,&cur); CHKERRQ(ierr);
  ierr = TransformObjectSetExplanation(cur,"Additive Schwarz method"); CHKERRQ(ierr);
  ierr = TransformObjectDefineOption(cur,"local-method"); CHKERRQ(ierr);
  /*
  ierr = TransformCurrentItemDefineOption(cur,"local method: -1: LU; $\\geq1$ fill levels$+1$; $\\leq-1$: negative of local CG iterations$+1$","local method"); CHKERRQ(ierr);
  */
  ierr = TransformObjectAddOption(cur,1); CHKERRQ(ierr);
  ierr = TransformObjectAddOptionExplanation(cur,1,"ILU(0)"); CHKERRQ(ierr); 
  ierr = TransformObjectAddOption(cur,2); CHKERRQ(ierr);
  ierr = TransformObjectAddOptionExplanation(cur,2,"ILU(1)"); CHKERRQ(ierr); 
  ierr = TransformObjectAddOption(cur,3); CHKERRQ(ierr);
  ierr = TransformObjectAddOptionExplanation(cur,3,"ILU(2)"); CHKERRQ(ierr); 
  ierr = TransformObjectAddOption(cur,-1); CHKERRQ(ierr);
  ierr = TransformObjectAddOptionExplanation(cur,-1,"LU"); CHKERRQ(ierr); 
  /*
  ierr = TransformObjectAddOption(cur,-7); CHKERRQ(ierr);
  ierr = TransformObjectAddOptionExplanation
    (cur,"6 iterations GMRES"); CHKERRQ(ierr); 
  */

  /*
   * Restricted Additive Schwarz
   */
  ierr = NewTransformObject(PREPROCESSOR,PCRASM,&cur); CHKERRQ(ierr);
  ierr = TransformObjectSetExplanation(cur,"Restrictied Additive Schwarz method"); CHKERRQ(ierr);
  ierr = TransformObjectDefineOption(cur,"local-method"); CHKERRQ(ierr);
  ierr = TransformObjectAddOption(cur,1); CHKERRQ(ierr);
  ierr = TransformObjectAddOptionExplanation(cur,1,"ILU(0)"); CHKERRQ(ierr); 
  ierr = TransformObjectAddOption(cur,2); CHKERRQ(ierr);
  ierr = TransformObjectAddOptionExplanation(cur,2,"ILU(1)"); CHKERRQ(ierr); 
  ierr = TransformObjectAddOption(cur,3); CHKERRQ(ierr);
  ierr = TransformObjectAddOptionExplanation(cur,3,"ILU(2)"); CHKERRQ(ierr); 
  ierr = TransformObjectAddOption(cur,-1); CHKERRQ(ierr);
  ierr = TransformObjectAddOptionExplanation(cur,-1,"LU"); CHKERRQ(ierr); 
  /*
  ierr = TransformObjectAddOption(cur,-7); CHKERRQ(ierr);
  ierr = TransformObjectAddOptionExplanation
    (cur,"6 iterations GMRES"); CHKERRQ(ierr); 
  */

  /*
   * Block Jacobi
   */
  ierr = NewTransformObject(PREPROCESSOR,PCBJACOBI,&cur); CHKERRQ(ierr);
  ierr = TransformObjectSetExplanation(cur,"Block Jacobi"); CHKERRQ(ierr);
  ierr = TransformObjectDefineOption(cur,"local-method"); CHKERRQ(ierr);
  ierr = TransformObjectIntAnnotate(cur,"notrans",1); CHKERRQ(ierr);
  ierr = TransformObjectAddOption(cur,1); CHKERRQ(ierr);
  ierr = TransformObjectAddOption(cur,2); CHKERRQ(ierr);
  ierr = TransformObjectAddOption(cur,3); CHKERRQ(ierr);
  ierr = TransformObjectAddOption(cur,-1); CHKERRQ(ierr);
  /*
  ierr = TransformObjectAddOption(cur,-7); CHKERRQ(ierr);
  */

  ierr = NewTransformObject(PREPROCESSOR,PCJACOBI,&cur); CHKERRQ(ierr);
  ierr = TransformObjectSetExplanation(cur,"Point Jacobi"); CHKERRQ(ierr);
  ierr = NewTransformObject(PREPROCESSOR,PCNONE,&cur); CHKERRQ(ierr);
  ierr = TransformObjectSetExplanation(cur,"None"); CHKERRQ(ierr);
  
  /*
   * Direct methods
   */
  ierr = NewTransformObject(PREPROCESSOR,PCLU,&cur); CHKERRQ(ierr);
  ierr = TransformObjectSetExplanation(cur,"Petsc LU"); CHKERRQ(ierr);
  ierr = TransformObjectIntAnnotate(cur,"directsolve",1); CHKERRQ(ierr);
#if defined(PETSC_HAVE_MUMPS)
  ierr = NewTransformObject(PREPROCESSOR,PCMUMPS,&cur); CHKERRQ(ierr);
  ierr = TransformObjectSetExplanation(cur,"Mumps LU"); CHKERRQ(ierr);
  ierr = TransformObjectIntAnnotate(cur,"directsolve",1); CHKERRQ(ierr);
#endif
#if defined(PETSC_HAVE_SPOOLES)
  ierr = NewTransformObject(PREPROCESSOR,PCSPOOLES,&cur); CHKERRQ(ierr);
  ierr = TransformObjectSetExplanation(cur,"Spooles LU"); CHKERRQ(ierr);
  ierr = TransformObjectIntAnnotate(cur,"directsolve",1); CHKERRQ(ierr);
#endif
#if defined(PETSC_HAVE_SUPERLU_DIST)
  ierr = NewTransformObject(PREPROCESSOR,PCSUPERLU,&cur); CHKERRQ(ierr);
  ierr = TransformObjectSetExplanation(cur,"SuperLU LU"); CHKERRQ(ierr);
  ierr = TransformObjectIntAnnotate(cur,"directsolve",1); CHKERRQ(ierr);
#endif
#if defined(PETSC_HAVE_UMFPACK)
  ierr = NewTransformObject(PREPROCESSOR,PCUMFPACK,&cur); CHKERRQ(ierr);
  ierr = TransformObjectSetExplanation(cur,"Umfpack LU"); CHKERRQ(ierr);
  ierr = TransformObjectIntAnnotate(cur,"directsolve",1); CHKERRQ(ierr);
#endif

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "disable_pcs"
static PetscErrorCode disable_pcs
(NumericalProblem theproblem,SalsaTransform pc)
{
  SalsaTransformObject curpc; int ds;
  PetscBool flg; PetscErrorCode ierr;
  PetscFunctionBegin;

  ierr = SysProRetrieveQuantity
    (theproblem,"variance","diagonal-sign",(void*)&ds,PETSC_NULL,
     &flg); CHKERRQ(ierr);
  ierr = TransformObjectGetByName("pc","euclid",&curpc); CHKERRQ(ierr);
  if (0) {
    ierr = TransformObjectMark(curpc);  CHKERRQ(ierr);}

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "pcoptionshandling"
/*! Disable certain preconditioners based on commandline options.
  
  At the moment this is only disabling of direct solvers if the 
  user asks for iterative only.
*/
static PetscErrorCode pcoptionshandling()
{
  char **names; int nnames,iname; PetscBool onlyiterative,onlydirect;
  PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = PetscOptionsHasName
    (PETSC_NULL,"-syspro_pc_iterative",&onlyiterative); CHKERRQ(ierr);
  ierr = PetscOptionsHasName
    (PETSC_NULL,"-syspro_pc_direct",&onlydirect); CHKERRQ(ierr);
  if (onlyiterative || onlydirect) { 
    ierr = RetrieveAllPreprocessorValues
      (PREPROCESSOR,(const char***)&names,&nnames); CHKERRQ(ierr);
    for (iname=0; iname<nnames; iname++) {
      SalsaTransformObject to; int isdirect; PetscBool flg;
      ierr = TransformObjectGetByName
        (PREPROCESSOR,names[iname],&to); CHKERRQ(ierr);
      ierr = TransformObjectGetIntAnnotation
        (to,"directsolve",&isdirect,&flg); CHKERRQ(ierr);
      if ( (flg && isdirect && onlyiterative) ||
           ( onlydirect && ( !flg || !isdirect) ) 
           ){
        ierr = TransformObjectMark(to); CHKERRQ(ierr);
      }
    }
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "create_solver"
/*! Create a solver and 
  install a monitor that dynamically increases the maximum
  number of iterations.
*/
static PetscErrorCode create_solver(NumericalProblem prob,void **ctx)
{
  KSP solver; PetscErrorCode ierr;
  PetscFunctionBegin;

  ierr = KSPCreate(prob->comm,&solver); CHKERRQ(ierr);
  /*ierr = SysProLinearInstallCustomKSPMonitor(solver); CHKERRQ(ierr);*/
  *(KSP*)ctx = solver;
#if 0
  { 
    PC pc;
    ierr = PCCreate(prob->comm,&pc); CHKERRQ(ierr);
    *(PC*)ctx = pc;
  }
#endif
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "destroy_solver"
static PetscErrorCode destroy_solver(void *ctx)
{
  KSP solver; PetscErrorCode ierr;
  PetscFunctionBegin;

  solver = (KSP)ctx;
  ierr = KSPDestroy(&solver); CHKERRQ(ierr);
#if 0
  {
    PC pc;
    pc = (PC)ctx;
    ierr = PCDestroy(pc); CHKERRQ(ierr);
  }
#endif
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "setup_pc"
static PetscErrorCode setup_pc
(const char *type,int pcv,PetscBool overwrite,
 NumericalProblem inproblem,NumericalProblem *outproblem,
 void *gctx,void **ctx,PetscBool *success)
{
  LinearSystem newproblem,problem = (LinearSystem)inproblem;
  KSP solver; PC pc;  Mat A,B,Buse;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  /*
   * handling the problem data is trivial
   */
  ierr = LinearSystemDuplicatePointers(problem,&newproblem); CHKERRQ(ierr);

  ierr = LinearSystemGetParts
    (problem,&A,&B,PETSC_NULL,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
  ierr = set_preconditioner_base_matrix((PCType)type,B,&Buse); CHKERRQ(ierr);
  ierr = LinearSystemSetParts
    (newproblem,PETSC_NULL,Buse,PETSC_NULL,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
  /*  ierr = PurgeAttachedArrays(A); CHKERRQ(ierr);*/
  ierr = PetscObjectStateIncrease((PetscObject)A); CHKERRQ(ierr);

  *outproblem = (NumericalProblem)newproblem;
  
  /*
   * the fun part is setting up the Petsc PC
   */
  solver = (KSP)gctx;
  ierr = KSPSetOperators(solver,A,Buse); CHKERRQ(ierr);
  ierr = KSPGetPC(solver,&pc); CHKERRQ(ierr);
#if 0
  pc = (PC)gctx;
  ierr = PCSetOperators(pc,A,Buse,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
#endif
  ierr = SetPetscOptionsForPC(pc,(PCType)type,pcv,6); CHKERRQ(ierr);
  {
    PetscLogDouble t1,t2,*tstore; PetscErrorCode ierr_save1,ierr_save2;
    ierr = PetscTime(&t1); CHKERRQ(ierr);
    ierr = PCSetFromOptions(pc); CHKERRQ(ierr);
    ierr = PetscPushErrorHandler(&PetscReturnErrorHandler,NULL); CHKERRQ(ierr);
    ierr = PCSetUp(pc); ierr_save1 = ierr;
    ierr = PCSetUpOnBlocks(pc); ierr_save2 = ierr;
    ierr = PetscTime(&t2); CHKERRQ(ierr);
    ierr = PetscPopErrorHandler(); CHKERRQ(ierr);
    *success = TRUTH(!ierr_save1 && !ierr_save2);
    ierr = PreprocessorGetContext("solution",(void**)&tstore); CHKERRQ(ierr);
    *tstore = t2-t1;
  }
  
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "unset_pc"
static PetscErrorCode unset_pc
(const char *type,PetscBool overwrite,void *gctx,void *ctx,
 NumericalProblem thisproblem,NumericalProblem upproblem,
 NumericalSolution old,NumericalSolution nnew)
{
  LinearSystem problem = (LinearSystem)thisproblem;
  Mat A;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = LinearSystemGetParts
    (problem,&A,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
  /*  ierr = PurgeAttachedArrays(A); CHKERRQ(ierr);*/
  ierr = PetscObjectStateDecrease((PetscObject)A); CHKERRQ(ierr);

  ierr = LinearSolutionCopy
    ((LinearSolution)old,(LinearSolution)nnew); CHKERRQ(ierr);
  /*ierr = LinearSolutionDelete((LinearSolution)old); CHKERRQ(ierr);*/
  //ierr = DeleteLinearSystem((LinearSystem)thisproblem); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DeclarePCPreprocessor"
PetscErrorCode DeclarePCPreprocessor(void)
{
  SystemPreprocessor pp; PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = DeclarePreprocessor
    (PREPROCESSOR,
     &setup_pc_choices,&disable_pcs,0,0,
     &create_solver,&destroy_solver,
     &setup_pc,&unset_pc); CHKERRQ(ierr);
  ierr = PreprocessorSetPreservedCategories
    (PREPROCESSOR,"laplace"); CHKERRQ(ierr);
  ierr = SystemPreprocessorGetByName(PREPROCESSOR,&pp); CHKERRQ(ierr);
  pp->optionshandling = &pcoptionshandling;
  PetscFunctionReturn(0);
}