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

[Giacomo Mulas]

Hello.

I am relatively new to petsc and slepc. I am trying to use slepc to solve a
sparse eigenvalue problem, in particular getting all eigenvalues in an
interval (spectral slicing). To run in parallel, the documentation states
that I have to specify an external linear solver, e.g. mumps, and explains
that I can use the command line options

-st_ksp_type preonly -st_pc_type cholesky \
-st_pc_factor_mat_solver_package mumps -mat_mumps_icntl_13 1

to set everything up for it. I was indeed able to get my code to work with
these command line options (even if I will need some additional work to make
it run faster in parallel, but that's another story). However, I would
rather like to be able to put these options in my code itself, as defaults.

So far, I was able to do so for all command line options except for
-mat_mumps_icntl_13 1

How should I do that? Here goes a snippet of the relevant code:

(... omit creation/population of H matrix ...)
       ierr = EPSCreate(mixpars->slepc_comm,&eps);CHKERRQ(ierr);
       ierr = EPSSetOperators(eps,H,PETSC_NULL);CHKERRQ(ierr);
       ierr = EPSSetProblemType(eps,EPS_HEP);CHKERRQ(ierr);
       ierr = EPSSetTolerances(eps, tol, PETSC_DECIDE);
       if (statesinlist>500)
 	ierr = EPSSetDimensions(eps, 500, PETSC_IGNORE, PETSC_IGNORE);
       else
 	ierr = EPSSetDimensions(eps, statesinlist, PETSC_IGNORE,
 				PETSC_IGNORE);
       CHKERRQ(ierr);
       if (statesinlist<50) ierr = EPSSetType(eps, EPSLAPACK);
       else {
 	ierr = EPSGetST(eps,&st);CHKERRQ(ierr);
 	ierr = STSetType(st,STSINVERT);CHKERRQ(ierr);
 	ierr = STGetKSP(st,&ksp);CHKERRQ(ierr);
 	ierr = KSPSetType(ksp,KSPPREONLY);
 	ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
 	ierr = PCSetType(pc, PCCHOLESKY);CHKERRQ(ierr);
 	if (mixpars->slepc_size>1) {
 	  ierr = PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);CHKERRQ(ierr);
 	  /*
 	  Mat F;
 	  ierr = PCFactorGetMatrix(pc,&F);CHKERRQ(ierr);
 	  PetscInt icntl = 13;
 	  PetscInt ival = 1;
 	  ierr = MatMumpsSetIcntl(F,icntl,ival); CHKERRQ(ierr);
 	  */
 	}
       }
       ierr = EPSSetFromOptions(eps); CHKERRQ(ierr);
       ierr = EPSSetInterval(eps,
 			    (PetscScalar) currsubclass->lowtargetlimit[j],
 			    (PetscScalar) currsubclass->uptargetlimit[j]);
       CHKERRQ(ierr);
       ierr = EPSSetWhichEigenpairs(eps, EPS_ALL); CHKERRQ(ierr);
       ierr = EPSSolve(eps);CHKERRQ(ierr);

The part which is commented out is what I tried to use to try to obtain the
same effect as the -mat_mumps_icntl_13 1 command line option. However, if I
decomment it, I get the following error:

[0]PETSC ERROR: --------------------- Error Message
------------------------------------
[0]PETSC ERROR: Operation done in wrong order!
[0]PETSC ERROR: Matrix not yet factored; call after KSPSetUp() or PCSetUp()!
[0]PETSC ERROR:
------------------------------------------------------------------------
[0]PETSC ERROR: Petsc Release Version 3.2.0, Patch 5, Sat Oct 29 13:45:54
CDT 2011 
[0]PETSC ERROR: See docs/changes/index.html for recent updates.
[0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
[0]PETSC ERROR: See docs/index.html for manual pages.

The code works fine if I comment out the lines from 
Mat F;
to
ierr = MatMumpsSetIcntl(F,icntl,ival); CHKERRQ(ierr);
and use the command line option 
-mat_mumps_icntl_13 1

Is it indeed possible to set this as a default in the code, and if yes, is
it documented somewhere? Can someone give me a hint?

A side issue: should I explicitly free the PC, KSP and ST I get using the 
EPSGetST(), STGetKSP() and KSPGetPC() calls respectively? I only use them to
set up default options (instead of giving them on the command line). Or can
I trust that EPSDestroy() will do it?

Thanks in advance, bye
Giacomo Mulas

-- 
_________________________________________________________________

Giacomo Mulas <gmulas at oa-cagliari.inaf.it>
_________________________________________________________________

OSSERVATORIO ASTRONOMICO DI CAGLIARI
Str. 54, Loc. Poggio dei Pini * 09012 Capoterra (CA)

Tel. (OAC): +39 070 71180 248     Fax : +39 070 71180 222
Tel. (UNICA): +39 070 675 4916
_________________________________________________________________

"When the storms are raging around you, stay right where you are"
                          (Freddy Mercury)
_________________________________________________________________