[petsc-users] Using MatShell without MatMult

[Matthew Knepley]

On Fri, Apr 7, 2017 at 5:11 AM, Francesco Migliorini <
francescomigliorini93 at gmail.com> wrote:

> Hello,
>
> I need to solve a linear system using GMRES without creating explicitly
> the matrix because very large. So, I am trying to use the MatShell strategy
> but I am stucked. The problem is that it seems to me that inside the
> user-defined MyMatMult it is required to use MatMult and this would
> honestly vanish all the gain from using this strategy. Indeed, I would need
> to access directly to the entries of the input vector, multiply them by
> some parameters imported in MyMatMult with *common* and finally compose
> the output vector without creating any matrix. First of all, is it
> possible?
>

Yes.


> Secondly, if so, where is my mistake? Here's an example of my code with a
> very simple 10x10 system with the identity matrix:
>
> [...]
> call PetscInitialize(PETSC_NULL_CHARACTER,perr)
> ind(1) = 10
> call VecCreate(PETSC_COMM_WORLD,feP,perr)
> call VecSetSizes(feP,PETSC_DECIDE,ind,perr)
> call VecSetFromOptions(feP,perr)
> call VecDuplicate(feP,u1P,perr)
> do jt = 1,10
>  ind(1) = jt-1
>  fval(1) = jt
>           call VecSetValues(feP,1,ind,fval(1),INSERT_VALUES,perr)
> enddo
> call VecAssemblyBegin(feP,perr)
> call VecAssemblyEnd(feP,perr)
> ind(1) = 10
> call MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, ind,
> ind, PETSC_NULL_INTEGER, TheShellMatrix, perr)
> call MatShellSetOperation(TheShellMatrix, MATOP_MULT, MyMatMult, perr)
>

Here I would probably use


http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatShellSetContext.html

instead of a common block, but that works too.


> call KSPCreate(PETSC_COMM_WORLD, ksp, perr)
> call KSPSetType(ksp,KSPGMRES,perr)
> call KSPSetOperators(ksp,TheShellMatrix,TheShellMatrix,perr)
> call KSPSolve(ksp,feP,u1P,perr)
> call PetscFinalize(PETSC_NULL_CHARACTER,perr)
> [...]
>
> subroutine MyMatMult(TheShellMatrix,T,AT,ierr)
> [...]
> Vec             T, AT
> Mat             TheShellMatrix
> PetscReal   fval(1), u0(1)
> [...]
> call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
> ind(1) = 10
> call VecCreate(PETSC_COMM_WORLD,AT,ierr)
> call VecSetSizes(AT,PETSC_DECIDE,ind,ierr)
> call VecSetFromOptions(AT,ierr)
>

Its not your job to create AT. We are passing it in, so just use it.


> do i =0,9
> ind(1) = i
> call VecGetValues(T,1,ind,u0(1),ierr)
> fval(1) = u0(1)
> call VecSetValues(AT,1,ind,fval(1),INSERT_VALUES,ierr)
>

You can do it this way, but its easier to use


http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecGetArray.html

outside the loop for both vectors.

   Matt


> enddo
> call VecAssemblyBegin(AT,ierr)
> call VecAssemblyEnd(AT,ierr)
> return
> end subroutine MyMatMult
>
> The output of this code is something completely invented but in some way
> related to the actual solution:
> 5.0964719143762542E-002
> 0.10192943828752508
> 0.15289415743128765
> 0.20385887657505017
> 0.25482359571881275
> 0.30578831486257529
> 0.35675303400633784
> 0.40771775315010034
> 0.45868247229386289
> 0.50964719143762549
>
> Instead, if I use MatMult in MyMatMult I get the right solution. Here's
> the code
>
> subroutine MyMatMult(TheShellMatrix,T,AT,ierr)
> [...]
> Vec             T, AT
> Mat             TheShellMatrix, IDEN
> PetscReal   fval(1)
> [...]
> call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
> ind(1) = 10
> call MatCreate(PETSC_COMM_WORLD,IDEN,ierr)
> call MatSetSizes(IDEN,PETSC_DECIDE,PETSC_DECIDE,ind,ind,ierr)
> call MatSetUp(IDEN,ierr)
> do i =0,9
> ind(1) = i
> fval(1) = 1
> call MatSetValues(IDEN,1,ind,1,ind,fval(1),INSERT_VALUES,ierr)
> enddo
> call MatAssemblyBegin(IDEN,MAT_FINAL_ASSEMBLY,ierr)
> call MatAssemblyEnd(IDEN,MAT_FINAL_ASSEMBLY,ierr)
> call MatMult(IDEN,T,AT,ierr)
> return
> end subroutine MyMatMult
>
> Thanks in advance for any answer!
> Francesco Migliorini
>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170407/4c866233/attachment.html>