<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class="">Hi,<div class="">I am having trouble with the implementation of SOR preconditioner with matshell. I am using Fortran 90.</div><div class=""><br class=""></div><div class="">(1) First, I implement a function called MatSOR_shell which has the same signature as MatSOR function, i.e.</div><div class="">MatSOR_wrapper(A, b, omega, flag, shift, its, lits, x, ierr)</div><div class="">This function returns a PETSc vector x which is the next iteration of the input vector x using the SOR algorithm.</div><div class=""><br class=""></div><div class="">(2) Next, I set matrix operation for my shell matrix, i.e.</div><div class="">CALL MatShellSetOperation(my_matshell, MATOP_SOR, MatSOR_wrapper, ierr)</div><div class=""><br class=""></div><div class="">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 class=""><br class=""></div><div class="">SUBROUTINE MatSOR_wrapper(A, b, omega, flag, shift, its, lits, x, ierr)</div><div class=""> Create, set, and assemble PETSc matrix <b class="">aa</b></div><div class=""> Create, set, and assemble PETSc rhs vector <b class="">bb</b></div><div class=""> CALL MatSOR(<b class="">aa</b>, <b class="">bb</b>, omega, flag, shift, its, lits, x, ierr)</div><div class="">END SUBROUTINE MatSOR_wrapper</div><div class=""><br class=""></div><div class="">With the above wrapper function, I suppose that the MatSOR function computes and return the vector x based on the input matrix <b class="">aa</b> and rhs vector <b class="">bb </b>using SOR algorithm.</div><div class=""><br class=""></div><div class="">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 class="">-pc_type sor -pc_sor_local_forward -ksp_monitor</div><div class=""><br class=""></div><div class="">The solving converts after exactly 2 iterations, with the residual norm as</div><div class=""><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class=""><span style="font-variant-ligatures: no-common-ligatures" class=""> 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;" class=""><span style="font-variant-ligatures: no-common-ligatures" class=""> </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;" class=""><span style="font-variant-ligatures: no-common-ligatures" class="">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;" class=""><br class=""></div><div style="margin: 0px; font-stretch: normal; line-height: normal;" class=""><span style="font-size: 11px; font-family: Menlo;" class=""><span style="font-family: Helvetica; font-size: 12px;" class="">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;" class=""><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class=""><span style="font-variant-ligatures: no-common-ligatures" class=""> </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;" class=""><span style="font-variant-ligatures: no-common-ligatures" class=""> 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;" class=""><span style="font-variant-ligatures: no-common-ligatures" class=""> 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;" class=""><span style="font-variant-ligatures: no-common-ligatures" class=""> 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;" class=""><span style="font-variant-ligatures: no-common-ligatures" class=""> 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;" class=""><span style="font-variant-ligatures: no-common-ligatures" class=""> 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;" class=""><span style="font-variant-ligatures: no-common-ligatures" class="">Linear solve converged due to CONVERGED_RTOL iterations 5</span></div><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class=""><span style="font-variant-ligatures: no-common-ligatures" class=""><br class=""></span></div><div style="margin: 0px; font-stretch: normal; line-height: normal;" class=""><span style="font-variant-ligatures: no-common-ligatures;" class="">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;" class=""><br class=""></div><div style="margin: 0px; font-stretch: normal; line-height: normal;" class="">My questions are:</div><div style="margin: 0px; font-stretch: normal; line-height: normal;" class="">(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;" class=""><font class="">(2) For parallel, as shown in PETSc manual at <a href="https://petsc.org/release/docs/manualpages/PC/PCSOR.html" class="">https://petsc.org/release/docs/manualpages/PC/PCSOR.html</a> as </font>"<span style="background-color: rgb(255, 255, 255);" class="">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;" class=""><span style="background-color: rgb(255, 255, 255);" class=""><br class=""></span></div><div style="margin: 0px; font-stretch: normal; line-height: normal;" class=""><span style="background-color: rgb(255, 255, 255);" class="">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;" class=""><span style="background-color: rgb(255, 255, 255);" class=""><br class=""></span></div><div style="margin: 0px; font-stretch: normal; line-height: normal;" class=""><span style="background-color: rgb(255, 255, 255);" class="">Han Tran</span></div><div style="margin: 0px; font-stretch: normal; line-height: normal;" class=""><span style="font-size: 11px;" class=""><span style="background-color: rgb(255, 255, 255);" class=""><br class=""></span></span></div></div><div class=""><br class=""></div></body></html>