[petsc-users] wrap MATMPIAIJ in MatShellSetOperation and use in SLEPc

[Jose E. Roman]

El 13/10/2014, a las 19:16, Denis Davydov escribió:

> Dear all, 
> 
> I am solving a generalised eigenvalue problem 
> 
>    S \cdot v = a M \cdot v  (1)
> 
> using Krylov-Schur eigensolver with the shift-and-invert transformation.  
> Now I would like to expand the problem to the case   
> 
>    [S + c^\ast \otimes c ] \cdot v = a M \cdot v      (2),
> 
> where “c” is some vector. 
> 
> As a first step, I was trying to “wrap" MATMPIAIJ matrix in the MatCreateShell() matrix 
> with MatShellSetOperation() which then does multiplication. 
> So far I do not add any action due to “c” (e.g. c = 0), but only aim to workout how to approach the problem.  
> Unfortunately, when I use shift-and-invert, i hit some errors (below). 
> In another test problem without shift-and-invert, it seems to be working fine (I am still testing). 
> Supposedly, the user-defined matrices are not supported in full in SLEPc and the error is triggered. 
> I am quite sure it is related to the usage shift-and-invert, but I don’t know implementation details to tell why is it the case. 
> 
> Is there any other way I could consider the problem (2) within SLEPc with shift-and-invert or 
> do some additional tricks to make it work within MatCreateShell()? 
> 
> p.s. this problem arises from the density functional theory with non-local pseudo-potentials, S is hermitian.
> 
> Kind regards,
> Denis 
> 

This case is discussed in section 3.4.2 of SLEPc's users guide. Basically, the default behaviour in shift-and-invert is that SLEPc builds matrix C=A-sigma*B (where sigma is the shift). In this case, the "A" matrix (the shell matrix in your case) needs two additional operations: duplicate and AXPY. You could implement these two operations in your shell matrix. You can also try with "-st_matmode inplace" in which case only AXPY is required. But probably what you want is to set "-st_matmode shell" so that you do not need any additional operation; the downside in all cases is that direct linear solves cannot be used and you must be cautious when using shift-and-invert with iterative linear solvers (see discussion at section 3.4.1).

Jose