<div dir="ltr"><div dir="ltr">On Wed, Nov 10, 2021 at 11:16 PM Han Tran <<a href="mailto:hantran@cs.utah.edu">hantran@cs.utah.edu</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div style="overflow-wrap: break-word;">Hi,<div>I am having trouble with the implementation of SOR preconditioner with matshell. I am using Fortran 90.</div><div><br></div><div>(1) First, I implement a function called MatSOR_shell which has the same signature as MatSOR function, i.e.</div><div>MatSOR_wrapper(A, b, omega, flag, shift, its, lits, x, ierr)</div><div>This function returns a PETSc vector x which is the next iteration of the input vector x using the SOR algorithm.</div><div><br></div><div>(2) Next, I set matrix operation for my shell matrix, i.e.</div><div>CALL MatShellSetOperation(my_matshell, MATOP_SOR, MatSOR_wrapper, ierr)</div><div><br></div><div>To test if it works in sequential (before implementing my own way to compute the vector x using SOR algorithm), I want to use PETSc function MatSOR inside my MatSOR_wrapper as follows</div><div><br></div><div>SUBROUTINE MatSOR_wrapper(A, b, omega, flag, shift, its, lits, x, ierr)</div><div>   Create, set, and assemble PETSc matrix <b>aa</b></div><div>   Create, set, and assemble PETSc rhs vector <b>bb</b></div><div>   CALL MatSOR(<b>aa</b>, <b>bb</b>, omega, flag, shift, its, lits, x, ierr)</div><div>END SUBROUTINE MatSOR_wrapper</div><div><br></div><div>With the above wrapper function, I suppose that the MatSOR function computes and return the vector x based on the input matrix <b>aa</b> and rhs vector <b>bb </b>using SOR algorithm.</div><div><br></div><div>I run my code to solve a linear system Ax = b using my matshell SOR preconditioner. I used GMRES for KSPSolve, and set preconditioner at runtime as:</div><div>-pc_type sor -pc_sor_local_forward -ksp_monitor</div><div><br></div><div>The solving converts after exactly 2 iterations, with the residual norm as</div><div><div style="margin:0px;font-stretch:normal;font-size:11px;line-height:normal;font-family:Menlo"><span style="font-variant-ligatures:no-common-ligatures">  0 KSP Residual norm 1.147215489446e+06 </span></div></div><div style="margin:0px;font-stretch:normal;font-size:11px;line-height:normal;font-family:Menlo"><span style="font-variant-ligatures:no-common-ligatures">  </span>1 KSP Residual norm 3.798348914902e-11 </div><div style="margin:0px;font-stretch:normal;font-size:11px;line-height:normal;font-family:Menlo"><span style="font-variant-ligatures:no-common-ligatures">Linear solve converged due to CONVERGED_RTOL iterations 1</span></div><div style="margin:0px;font-stretch:normal;font-size:11px;line-height:normal;font-family:Menlo"><br></div><div style="margin:0px;font-stretch:normal;line-height:normal"><span style="font-size:11px;font-family:Menlo"><span style="font-family:Helvetica;font-size:12px">However, the converged solution is wrong. For comparison, I solved the same above system without using matshell (i.e. using standard PETSc assembled matrix) and also using -pc_type sor </span>-pc_sor_local_forward, with the following results</span></div><div style="margin:0px;font-stretch:normal;line-height:normal"><div style="margin:0px;font-stretch:normal;font-size:11px;line-height:normal;font-family:Menlo"><span style="font-variant-ligatures:no-common-ligatures">  </span>0 KSP Residual norm 1.147215489446e+06 </div><div style="margin:0px;font-stretch:normal;font-size:11px;line-height:normal;font-family:Menlo"><span style="font-variant-ligatures:no-common-ligatures">  1 KSP Residual norm 7.296837190230e+05 </span></div><div style="margin:0px;font-stretch:normal;font-size:11px;line-height:normal;font-family:Menlo"><span style="font-variant-ligatures:no-common-ligatures">  2 KSP Residual norm 4.881295368957e+05 </span></div><div style="margin:0px;font-stretch:normal;font-size:11px;line-height:normal;font-family:Menlo"><span style="font-variant-ligatures:no-common-ligatures">  3 KSP Residual norm 3.176644114439e+05 </span></div><div style="margin:0px;font-stretch:normal;font-size:11px;line-height:normal;font-family:Menlo"><span style="font-variant-ligatures:no-common-ligatures">  4 KSP Residual norm 1.941135089932e+05 </span></div><div style="margin:0px;font-stretch:normal;font-size:11px;line-height:normal;font-family:Menlo"><span style="font-variant-ligatures:no-common-ligatures">  5 KSP Residual norm 2.870541410834e-10 </span></div><div style="margin:0px;font-stretch:normal;font-size:11px;line-height:normal;font-family:Menlo"><span style="font-variant-ligatures:no-common-ligatures">Linear solve converged due to CONVERGED_RTOL iterations 5</span></div></div></div></blockquote><div><br></div><div>I would run a problem small enough that you can print all the vectors to the screen. You know that the residual diverges at iterate 1, so I would just print vectors</div><div>until I found the difference. For example, print b and x for each iterate of your SOR and the original SOR.</div><div><br></div><div>  Thanks,</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 style="overflow-wrap: break-word;"><div style="margin:0px;font-stretch:normal;line-height:normal"><div style="margin:0px;font-stretch:normal;line-height:normal"><span style="font-variant-ligatures:no-common-ligatures">It can be seen that with my matshell SOR, the solver computes correctly the residual norm of the first iteration (i.e. both methods give the same value </span>1.147215489446e+06) but not the second iteration. I do not understand why the solver quickly gives a very small residual norm and stops the solving because the relative decrease of residual norm is smaller then RTOL. I am quite sure that the matrix aa and rhs vector bb provided inside MatSOR_wrapper are correct because I did compare with the assembled matrix used in second method (i.e. the assembled-matrix method). I also verified that the MatSOR_wrapper is called 2 times during the solving (because 2 iterations are done) and it returns a correct vector x with an input vector x=0 using forward SOR algorithm.</div><div style="margin:0px;font-stretch:normal;line-height:normal"><br></div><div style="margin:0px;font-stretch:normal;line-height:normal">My questions are:</div><div style="margin:0px;font-stretch:normal;line-height:normal">(1) Did I miss providing something else inside my function MatSOR_wrapper, or any other options for the solver?</div><div style="margin:0px;font-stretch:normal;line-height:normal"><font>(2) For parallel, as shown in PETSc manual at <a href="https://petsc.org/release/docs/manualpages/PC/PCSOR.html" target="_blank">https://petsc.org/release/docs/manualpages/PC/PCSOR.html</a> as </font>"<span style="background-color:rgb(255,255,255)">Not a true parallel SOR, in parallel this implementation corresponds to block Jacobi with SOR on each block”. If I use the same wrapper as shown above, my matrix aa should be the diagonal block (sequential matrix, diagonal submatrix owned by my proc), and my vector bb should be the local part of the global rhs vector. Is it correct?</span></div><div style="margin:0px;font-stretch:normal;line-height:normal"><span style="background-color:rgb(255,255,255)"><br></span></div><div style="margin:0px;font-stretch:normal;line-height:normal"><span style="background-color:rgb(255,255,255)">Thank you very much for your help. I am very sorry for my lengthy email. I just want to explain the problem as much as I can. I look forward to your answers.</span></div><div style="margin:0px;font-stretch:normal;line-height:normal"><span style="background-color:rgb(255,255,255)"><br></span></div><div style="margin:0px;font-stretch:normal;line-height:normal"><span style="background-color:rgb(255,255,255)">Han Tran</span></div><div style="margin:0px;font-stretch:normal;line-height:normal"><span style="font-size:11px"><span style="background-color:rgb(255,255,255)"><br></span></span></div></div><div><br></div></div></blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>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><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>