[petsc-users] Using MatShell without MatMult

[Francesco Migliorini]

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?
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)
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)
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)
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170407/7f5c74c7/attachment-0001.html>