<div>You should also not call PetscInitialize() from within your user MatMult function.</div><div><br></div><div><br></div><div><br></div><div><br><div class="gmail_quote"><div>On Fri, 7 Apr 2017 at 13:24, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg">On Fri, Apr 7, 2017 at 5:11 AM, Francesco Migliorini <span class="gmail_msg"><<a href="mailto:francescomigliorini93@gmail.com" class="gmail_msg" target="_blank">francescomigliorini93@gmail.com</a>></span> wrote:<br class="gmail_msg"><blockquote class="gmail_quote gmail_msg" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="gmail_msg">Hello,<div class="gmail_msg"><br class="gmail_msg"></div><div class="gmail_msg">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 class="gmail_msg">common</i> and finally compose the output vector without creating any matrix. First of all, is it possible? </div></div></blockquote><div class="gmail_msg"><br class="gmail_msg"></div></div></div></div><div class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg"><div class="gmail_msg">Yes.</div></div></div></div><div class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg"><div class="gmail_msg"> </div><blockquote class="gmail_quote gmail_msg" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="gmail_msg"><div class="gmail_msg">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 class="gmail_msg"><br class="gmail_msg"></div><div class="gmail_msg">[...]</div><div class="gmail_msg">call PetscInitialize(PETSC_NULL_CHARACTER,perr)<br class="gmail_msg"></div><div class="gmail_msg"><div class="gmail_msg">ind(1) = 10</div><div class="gmail_msg">call VecCreate(PETSC_COMM_WORLD,feP,perr)</div><div class="gmail_msg">call VecSetSizes(feP,PETSC_DECIDE,ind,perr)</div><div class="gmail_msg">call VecSetFromOptions(feP,perr)</div><div class="gmail_msg">call VecDuplicate(feP,u1P,perr)</div><div class="gmail_msg">do jt = 1,10</div><div class="gmail_msg"><span class="m_6294238755145690588gmail-m_-4755789698870337980gmail-Apple-tab-span gmail_msg" style="white-space:pre-wrap"> </span> ind(1) = jt-1</div><div class="gmail_msg"><span class="m_6294238755145690588gmail-m_-4755789698870337980gmail-Apple-tab-span gmail_msg" style="white-space:pre-wrap"> </span> fval(1) = jt</div><div class="gmail_msg"> call VecSetValues(feP,1,ind,fval(1),INSERT_VALUES,perr)<br class="gmail_msg"></div><div class="gmail_msg">enddo</div><div class="gmail_msg"><div class="gmail_msg">call VecAssemblyBegin(feP,perr)</div><div class="gmail_msg">call VecAssemblyEnd(feP,perr)</div></div><div class="gmail_msg">ind(1) = 10<br class="gmail_msg"></div><div class="gmail_msg">call MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, ind, ind, PETSC_NULL_INTEGER, TheShellMatrix, perr)</div><div class="gmail_msg">call MatShellSetOperation(TheShellMatrix, MATOP_MULT, MyMatMult, perr)</div></div></div></blockquote><div class="gmail_msg"><br class="gmail_msg"></div></div></div></div><div class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg"><div class="gmail_msg">Here I would probably use</div><div class="gmail_msg"><br class="gmail_msg"></div><div class="gmail_msg"> <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatShellSetContext.html" class="gmail_msg" target="_blank">http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatShellSetContext.html</a></div><div class="gmail_msg"><br class="gmail_msg"></div><div class="gmail_msg">instead of a common block, but that works too.</div></div></div></div><div class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg"><div class="gmail_msg"> </div><blockquote class="gmail_quote gmail_msg" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="gmail_msg"><div class="gmail_msg"><div class="gmail_msg">call KSPCreate(PETSC_COMM_WORLD, ksp, perr)<br class="gmail_msg"></div><div class="gmail_msg">call KSPSetType(ksp,KSPGMRES,perr)</div><div class="gmail_msg">call KSPSetOperators(ksp,TheShellMatrix,TheShellMatrix,perr)</div><div class="gmail_msg">call KSPSolve(ksp,feP,u1P,perr) <br class="gmail_msg"></div><div class="gmail_msg">call PetscFinalize(PETSC_NULL_CHARACTER,perr)<br class="gmail_msg"></div></div><div class="gmail_msg">[...]</div><div class="gmail_msg"><br class="gmail_msg"></div><div class="gmail_msg">subroutine MyMatMult(TheShellMatrix,T,AT,ierr)<br class="gmail_msg"></div><div class="gmail_msg">[...]</div><div class="gmail_msg">Vec T, AT<br class="gmail_msg"></div><div class="gmail_msg">Mat TheShellMatrix</div><div class="gmail_msg">PetscReal fval(1), u0(1)</div><div class="gmail_msg">[...]</div><div class="gmail_msg">call PetscInitialize(PETSC_NULL_CHARACTER,ierr)<br class="gmail_msg"></div><div class="gmail_msg"><div class="gmail_msg">ind(1) = 10</div><div class="gmail_msg">call VecCreate(PETSC_COMM_WORLD,AT,ierr)</div><div class="gmail_msg">call VecSetSizes(AT,PETSC_DECIDE,ind,ierr)</div><div class="gmail_msg">call VecSetFromOptions(AT,ierr) </div></div></div></blockquote><div class="gmail_msg"><br class="gmail_msg"></div></div></div></div><div class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg"><div class="gmail_msg">Its not your job to create AT. We are passing it in, so just use it.</div></div></div></div><div class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg"><div class="gmail_msg"> </div><blockquote class="gmail_quote gmail_msg" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="gmail_msg"><div class="gmail_msg"><div class="gmail_msg">do i =0,9</div><div class="gmail_msg"><span class="m_6294238755145690588gmail-m_-4755789698870337980gmail-Apple-tab-span gmail_msg" style="white-space:pre-wrap"> </span>ind(1) = i</div><div class="gmail_msg"><span class="m_6294238755145690588gmail-m_-4755789698870337980gmail-Apple-tab-span gmail_msg" style="white-space:pre-wrap"> </span>call VecGetValues(T,1,ind,u0(1),ierr)</div><div class="gmail_msg"><span class="m_6294238755145690588gmail-m_-4755789698870337980gmail-Apple-tab-span gmail_msg" style="white-space:pre-wrap"> </span>fval(1) = u0(1)</div><div class="gmail_msg"><span class="m_6294238755145690588gmail-m_-4755789698870337980gmail-Apple-tab-span gmail_msg" style="white-space:pre-wrap"> </span>call VecSetValues(AT,1,ind,fval(1),INSERT_VALUES,ierr)</div></div></div></blockquote><div class="gmail_msg"><br class="gmail_msg"></div></div></div></div><div class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg"><div class="gmail_msg">You can do it this way, but its easier to use</div><div class="gmail_msg"><br class="gmail_msg"></div><div class="gmail_msg"> <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecGetArray.html" class="gmail_msg" target="_blank">http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecGetArray.html</a></div><div class="gmail_msg"><br class="gmail_msg"></div><div class="gmail_msg">outside the loop for both vectors.</div><div class="gmail_msg"><br class="gmail_msg"></div><div class="gmail_msg"> Matt</div></div></div></div><div class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg"><div class="gmail_msg"> </div><blockquote class="gmail_quote gmail_msg" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="gmail_msg"><div class="gmail_msg"><div class="gmail_msg">enddo</div><div class="gmail_msg">call VecAssemblyBegin(AT,ierr)</div><div class="gmail_msg">call VecAssemblyEnd(AT,ierr)</div><div class="gmail_msg">return</div><div class="gmail_msg">end subroutine MyMatMult</div></div><div class="gmail_msg"><br class="gmail_msg"></div><div class="gmail_msg">The output of this code is something completely invented but in some way related to the actual solution:</div><div class="gmail_msg"><div class="gmail_msg">5.0964719143762542E-002</div><div class="gmail_msg">0.10192943828752508 </div><div class="gmail_msg">0.15289415743128765 </div><div class="gmail_msg">0.20385887657505017 </div><div class="gmail_msg">0.25482359571881275 </div><div class="gmail_msg">0.30578831486257529 </div><div class="gmail_msg">0.35675303400633784 </div><div class="gmail_msg">0.40771775315010034 </div><div class="gmail_msg">0.45868247229386289 </div><div class="gmail_msg">0.50964719143762549 </div></div><div class="gmail_msg"><br class="gmail_msg"></div><div class="gmail_msg">Instead, if I use MatMult in MyMatMult I get the right solution. Here's the code</div><div class="gmail_msg"><br class="gmail_msg"></div><div class="gmail_msg"><div class="gmail_msg">subroutine MyMatMult(TheShellMatrix,T,AT,ierr)</div><div class="gmail_msg">[...]</div><div class="gmail_msg">Vec T, AT</div><div class="gmail_msg">Mat TheShellMatrix, IDEN</div><div class="gmail_msg">PetscReal fval(1)</div><div class="gmail_msg">[...]</div><div class="gmail_msg">call PetscInitialize(PETSC_NULL_CHARACTER,ierr)</div><div class="gmail_msg">ind(1) = 10</div><div class="gmail_msg">call MatCreate(PETSC_COMM_WORLD,IDEN,ierr)</div><div class="gmail_msg">call MatSetSizes(IDEN,PETSC_DECIDE,PETSC_DECIDE,ind,ind,ierr)</div><div class="gmail_msg">call MatSetUp(IDEN,ierr)</div><div class="gmail_msg">do i =0,9</div><div class="gmail_msg"><span class="m_6294238755145690588gmail-m_-4755789698870337980gmail-Apple-tab-span gmail_msg" style="white-space:pre-wrap"> </span>ind(1) = i</div><div class="gmail_msg"><span class="m_6294238755145690588gmail-m_-4755789698870337980gmail-Apple-tab-span gmail_msg" style="white-space:pre-wrap"> </span>fval(1) = 1</div><div class="gmail_msg"><span class="m_6294238755145690588gmail-m_-4755789698870337980gmail-Apple-tab-span gmail_msg" style="white-space:pre-wrap"> </span>call MatSetValues(IDEN,1,ind,1,ind,fval(1),INSERT_VALUES,ierr)</div><div class="gmail_msg">enddo</div><div class="gmail_msg">call MatAssemblyBegin(IDEN,MAT_FINAL_ASSEMBLY,ierr)</div><div class="gmail_msg">call MatAssemblyEnd(IDEN,MAT_FINAL_ASSEMBLY,ierr)</div><div class="gmail_msg">call MatMult(IDEN,T,AT,ierr)</div><div class="gmail_msg">return</div><div class="gmail_msg">end subroutine MyMatMult</div></div><div class="gmail_msg"><br class="gmail_msg"></div><div class="gmail_msg">Thanks in advance for any answer!</div><span class="m_6294238755145690588gmail-HOEnZb gmail_msg"><font color="#888888" class="gmail_msg"><div class="gmail_msg">Francesco Migliorini</div><div class="gmail_msg"><br class="gmail_msg"></div></font></span></div>
</blockquote></div></div></div><div class="gmail_msg"><div class="gmail_extra gmail_msg"><br class="gmail_msg"><br clear="all" class="gmail_msg"><div class="gmail_msg"><br class="gmail_msg"></div>-- <br class="gmail_msg"><div class="m_6294238755145690588gmail_signature gmail_msg">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br class="gmail_msg">-- Norbert Wiener</div>
</div></div></blockquote></div></div>