You should also not call PetscInitialize() from within your user MatMult function.
On Fri, 7 Apr 2017 at 13:24, Matthew Knepley <knep...@gmail.com> wrote: > On Fri, Apr 7, 2017 at 5:11 AM, Francesco Migliorini < > francescomigliorin...@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 >
