<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Fri, Apr 7, 2017 at 5:11 AM, Francesco Migliorini <span dir="ltr"><<a href="mailto:francescomigliorini93@gmail.com" target="_blank">francescomigliorini93@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><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? </div></div></blockquote><div><br></div><div>Yes.</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>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_<wbr>CHARACTER,perr)<br></div><div><div>ind(1) = 10</div><div>call VecCreate(PETSC_COMM_WORLD,<wbr>feP,perr)</div><div>call VecSetSizes(feP,PETSC_DECIDE,<wbr>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-m_-4755789698870337980gmail-Apple-tab-span" style="white-space:pre-wrap"> </span> ind(1) = jt-1</div><div><span class="gmail-m_-4755789698870337980gmail-Apple-tab-span" style="white-space:pre-wrap"> </span> fval(1) = jt</div><div> call VecSetValues(feP,1,ind,fval(1)<wbr>,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_<wbr>WORLD, PETSC_DECIDE, PETSC_DECIDE, ind, ind, PETSC_NULL_INTEGER, TheShellMatrix, perr)</div><div>call MatShellSetOperation(<wbr>TheShellMatrix, MATOP_MULT, MyMatMult, perr)</div></div></div></blockquote><div><br></div><div>Here I would probably use</div><div><br></div><div> <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatShellSetContext.html">http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatShellSetContext.html</a></div><div><br></div><div>instead of a common block, but that works too.</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><div>call KSPCreate(PETSC_COMM_WORLD, ksp, perr)<br></div><div>call KSPSetType(ksp,KSPGMRES,perr)</div><div>call KSPSetOperators(ksp,<wbr>TheShellMatrix,TheShellMatrix,<wbr>perr)</div><div>call KSPSolve(ksp,feP,u1P,perr) <br></div><div>call PetscFinalize(PETSC_NULL_<wbr>CHARACTER,perr)<br></div></div><div>[...]</div><div><br></div><div>subroutine MyMatMult(TheShellMatrix,T,AT,<wbr>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_<wbr>CHARACTER,ierr)<br></div><div><div>ind(1) = 10</div><div>call VecCreate(PETSC_COMM_WORLD,AT,<wbr>ierr)</div><div>call VecSetSizes(AT,PETSC_DECIDE,<wbr>ind,ierr)</div><div>call VecSetFromOptions(AT,ierr) </div></div></div></blockquote><div><br></div><div>Its not your job to create AT. We are passing it in, so just use it.</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><div>do i =0,9</div><div><span class="gmail-m_-4755789698870337980gmail-Apple-tab-span" style="white-space:pre-wrap"> </span>ind(1) = i</div><div><span class="gmail-m_-4755789698870337980gmail-Apple-tab-span" style="white-space:pre-wrap"> </span>call VecGetValues(T,1,ind,u0(1),<wbr>ierr)</div><div><span class="gmail-m_-4755789698870337980gmail-Apple-tab-span" style="white-space:pre-wrap"> </span>fval(1) = u0(1)</div><div><span class="gmail-m_-4755789698870337980gmail-Apple-tab-span" style="white-space:pre-wrap"> </span>call VecSetValues(AT,1,ind,fval(1),<wbr>INSERT_VALUES,ierr)</div></div></div></blockquote><div><br></div><div>You can do it this way, but its easier to use</div><div><br></div><div> <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecGetArray.html">http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecGetArray.html</a></div><div><br></div><div>outside the loop for both vectors.</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><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,<wbr>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_<wbr>CHARACTER,ierr)</div><div>ind(1) = 10</div><div>call MatCreate(PETSC_COMM_WORLD,<wbr>IDEN,ierr)</div><div>call MatSetSizes(IDEN,PETSC_DECIDE,<wbr>PETSC_DECIDE,ind,ind,ierr)</div><div>call MatSetUp(IDEN,ierr)</div><div>do i =0,9</div><div><span class="gmail-m_-4755789698870337980gmail-Apple-tab-span" style="white-space:pre-wrap"> </span>ind(1) = i</div><div><span class="gmail-m_-4755789698870337980gmail-Apple-tab-span" style="white-space:pre-wrap"> </span>fval(1) = 1</div><div><span class="gmail-m_-4755789698870337980gmail-Apple-tab-span" style="white-space:pre-wrap"> </span>call MatSetValues(IDEN,1,ind,1,ind,<wbr>fval(1),INSERT_VALUES,ierr)</div><div>enddo</div><div>call MatAssemblyBegin(IDEN,MAT_<wbr>FINAL_ASSEMBLY,ierr)</div><div>call MatAssemblyEnd(IDEN,MAT_FINAL_<wbr>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><span class="gmail-HOEnZb"><font color="#888888"><div>Francesco Migliorini</div><div><br></div></font></span></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</div></div>