<div dir="ltr">Hello,<div><br></div><div>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 <i>common</i> 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:</div><div><br></div><div>[...]</div><div>call PetscInitialize(PETSC_NULL_CHARACTER,perr)<br></div><div><div>ind(1) = 10</div><div>call VecCreate(PETSC_COMM_WORLD,feP,perr)</div><div>call VecSetSizes(feP,PETSC_DECIDE,ind,perr)</div><div>call VecSetFromOptions(feP,perr)</div><div>call VecDuplicate(feP,u1P,perr)</div><div>do jt = 1,10</div><div><span class="gmail-Apple-tab-span" style="white-space:pre">        </span>  ind(1) = jt-1</div><div><span class="gmail-Apple-tab-span" style="white-space:pre">       </span>  fval(1) = jt</div><div>          call VecSetValues(feP,1,ind,fval(1),INSERT_VALUES,perr)<br></div><div>enddo</div><div><div>call VecAssemblyBegin(feP,perr)</div><div>call VecAssemblyEnd(feP,perr)</div></div><div>ind(1) = 10<br></div><div>call MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, ind, ind, PETSC_NULL_INTEGER, TheShellMatrix, perr)</div><div>call MatShellSetOperation(TheShellMatrix, MATOP_MULT, MyMatMult, perr)</div><div>call KSPCreate(PETSC_COMM_WORLD, ksp, perr)<br></div><div>call KSPSetType(ksp,KSPGMRES,perr)</div><div>call KSPSetOperators(ksp,TheShellMatrix,TheShellMatrix,perr)</div><div>call KSPSolve(ksp,feP,u1P,perr)          <br></div><div>call PetscFinalize(PETSC_NULL_CHARACTER,perr)<br></div></div><div>[...]</div><div><br></div><div>subroutine MyMatMult(TheShellMatrix,T,AT,ierr)<br></div><div>[...]</div><div>Vec             T, AT<br></div><div>Mat             TheShellMatrix</div><div>PetscReal   fval(1), u0(1)</div><div>[...]</div><div>call PetscInitialize(PETSC_NULL_CHARACTER,ierr)<br></div><div><div>ind(1) = 10</div><div>call VecCreate(PETSC_COMM_WORLD,AT,ierr)</div><div>call VecSetSizes(AT,PETSC_DECIDE,ind,ierr)</div><div>call VecSetFromOptions(AT,ierr) </div></div><div><div>do i =0,9</div><div><span class="gmail-Apple-tab-span" style="white-space:pre">     </span>ind(1) = i</div><div><span class="gmail-Apple-tab-span" style="white-space:pre">     </span>call VecGetValues(T,1,ind,u0(1),ierr)</div><div><span class="gmail-Apple-tab-span" style="white-space:pre">  </span>fval(1) = u0(1)</div><div><span class="gmail-Apple-tab-span" style="white-space:pre">        </span>call VecSetValues(AT,1,ind,fval(1),INSERT_VALUES,ierr)</div><div>enddo</div><div>call VecAssemblyBegin(AT,ierr)</div><div>call VecAssemblyEnd(AT,ierr)</div><div>return</div><div>end subroutine MyMatMult</div></div><div><br></div><div>The output of this code is something completely invented but in some way related to the actual solution:</div><div><div>5.0964719143762542E-002</div><div>0.10192943828752508     </div><div>0.15289415743128765     </div><div>0.20385887657505017     </div><div>0.25482359571881275     </div><div>0.30578831486257529     </div><div>0.35675303400633784     </div><div>0.40771775315010034     </div><div>0.45868247229386289     </div><div>0.50964719143762549  </div></div><div><br></div><div>Instead, if I use MatMult in MyMatMult I get the right solution. Here's the code</div><div><br></div><div><div>subroutine MyMatMult(TheShellMatrix,T,AT,ierr)</div><div>[...]</div><div>Vec             T, AT</div><div>Mat             TheShellMatrix, IDEN</div><div>PetscReal   fval(1)</div><div>[...]</div><div>call PetscInitialize(PETSC_NULL_CHARACTER,ierr)</div><div>ind(1) = 10</div><div>call MatCreate(PETSC_COMM_WORLD,IDEN,ierr)</div><div>call MatSetSizes(IDEN,PETSC_DECIDE,PETSC_DECIDE,ind,ind,ierr)</div><div>call MatSetUp(IDEN,ierr)</div><div>do i =0,9</div><div><span class="gmail-Apple-tab-span" style="white-space:pre">    </span>ind(1) = i</div><div><span class="gmail-Apple-tab-span" style="white-space:pre">     </span>fval(1) = 1</div><div><span class="gmail-Apple-tab-span" style="white-space:pre">    </span>call MatSetValues(IDEN,1,ind,1,ind,fval(1),INSERT_VALUES,ierr)</div><div>enddo</div><div>call MatAssemblyBegin(IDEN,MAT_FINAL_ASSEMBLY,ierr)</div><div>call MatAssemblyEnd(IDEN,MAT_FINAL_ASSEMBLY,ierr)</div><div>call MatMult(IDEN,T,AT,ierr)</div><div>return</div><div>end subroutine MyMatMult</div></div><div><br></div><div>Thanks in advance for any answer!</div><div>Francesco Migliorini</div><div><br></div></div>