[petsc-users] How to set up -mat_mumps_icntl13 1 in the code?

[Jose E. Roman]

El 14/05/2012, a las 19:45, Giacomo Mulas escribió:

> On Mon, 14 May 2012, Matthew Knepley wrote:
> 
>> The problems below are exactly why we do not advocate this strategy.
>> However, you
>> can get this effectusing http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/Sys/PetscOpti
>> onsSetValue.html
> 
> Thanks Matt!
> 
> Hmm, this would be more or less the same as adding the options to argc and
> argv before initialising slepc and petsc, I gather. Functional, even if not
> quite as clean as I would have liked. My idea was to set some working
> default in the code, but still honor command line options if different, whereas what you suggest overrides any conflicting command line options. Indeed, in my code I call EPSSetFromOptions() _after_ my explicit function
> calls, to make sure that command line options override my defaults (at
> least, that's my intended logic...).
> 
> Maybe I should first check for appropriate command line options, and then
> set them in the code _only if_ none were given on the command line?
> What about:
> 	  PetscInt ival;
> 	  PetscBool isset;
> 	  ierr = PetscOptionsGetInt(PETSC_NULL, "-mat_mumps_icntl_13",
> 	  &ival,
> 				    &isset);CHKERRQ(ierr);
> 	  if (isset==PETSC_FALSE) {
> 	    ierr = PetscOptionsSetValue("-mat_mumps_icntl_13","1");
> 	    CHKERRQ(ierr);
> 	  }
> 
> Does it look sensible?
> 
> Also, could you please specify a bit more what strategy you do not advocate,
> and for what problems? Do you mean setting up options by direct function
> calls, as opposed to setting them via PetscOptionsSetValue() ?
> Or do you mean it is deprecated to set options in the code in any way, be it
> via direct function calls or via PetscOptionsSetValue() ?
> 
> And the problems you are referring to are essentially the runtime error I
> get, i.e. calling functions in the wrong order? Or something else?

In order to set MUMPS options procedurally with MatMumpsSetIcntl (as in PETSc's ex52) one must first call KSPSetOperators() and PCSetType(). The problem is that in your example KSPSetOperators() will be called in STSetUp(). You could try to call STSetUp() explicitly but you may end up factoring the matrix twice. The design of ST is not prepared for a solution to your problem.

> 
> Finally, if I call EPSGetST(), should I then call STDestroy() on the ST I
> obtained, or will it be done by EPSDestroy()? And the same question iterates
> for the KSP I obtained via STGetKSP(), and the PC I obtained via KSPGetPC().

EPSGetST() does not increase the reference count, so EPSDestroy() will free it for you. When in doubt, run with -malloc_dump to see if there are uncleaned objects.

Jose

> 
> Thanks, bye
> Giacomo