# Debugging Shell matrix

I am trying to solve complex valued Poisson equation $$(C + \nabla. D \nabla )u = f \text{ ;where C, D, u and f are complex numbers.}$$

I am breaking this eqn into real valued problem, which is of the form

$$\underbrace{\begin{pmatrix} C_r + \nabla. D_r \nabla & -C_i - \nabla. D_i \nabla \\ C_i + \nabla. D_i \nabla & C_r + \nabla. D_r \nabla \end{pmatrix} }_{A} \begin{pmatrix} u_r \\ u_i \end{pmatrix} = \begin{pmatrix} f_r \\ f_i \end{pmatrix}$$

where 'r' and 'i' denote the real and imaginary parts. I wrote a Preconditioner (PC) which has the same form as the above equation by defining MAT A and Vec. I verified the PC was working fine by solving the above equation for an example. I have also defined a shell matrix operator which computes the action of A operator on the vector and registered the shell matrix and PC with ksp. When I try to solve equations using shell matrix and PC I am getting wrong solution. I am using Dirichlet boundary conditions(bcs) for $u_r$ and $u_i$ and homogeneous bcs for residual equation passed to PC. As a sanity check I did -ksp_type preonly with Dirichlet bcs for PC and I do get a correct solution. The ksp for both PC and Shell matrix is gmres.

Any suggestions on what could be going wrong?

• We resolved the issue offline. The shell operator had implemented boundary conditions incorrectly. You could turn this question into a debugging question for shell matrices (the part about equivalent real formulation is distracting), though most of the techniques are covered by scicomp.stackexchange.com/q/513/119. – Jed Brown Jul 15 '12 at 15:41

Run with -ksp_monitor_true_residual. Does the residual go to zero? If so, then the definition of the matrix is incorrect (or the right hand side or your meaning of "correct"). If it does not converge, try a small well-conditioned case with -pc_type none to rule out bugs in the preconditioner. If the problem is in the preconditioner, try -ksp_type fgmres to see if the preconditioner is nonlinear (FGMRES converging when GMRES does not usually means the preconditioner is nonlinear).
1. Check that matrix-multiply is correct by setting the "correct solution" in $x$, then computing $b \gets A x$, then solving $A y = b$. If $x = y$, the matrix multiply is likely correct.
2. Compare the true residual to the preconditioned residual in -ksp_monitor_true_residual. Are they similar? Try the same with -ksp_norm_type unpreconditioned.
• -ksp_monitor_true_residual is giving me residual close enough to zero with the PC. With -pc_type none it shows that residual is going to zero but it might take tens of thousands of iterations to converge. With both gmres' and fgmres the solver converges ; but the solution is not correct. With an another example I do get correct and converged solution with ksp+PC combination. – amneet Jul 15 '12 at 2:06
• Also if I build my PC as a sparse matrix A defined in the above equation, doesn't that make my PC linear only? It should not be nonlinear right? Do I need to specify nullspace for this problem? – amneet Jul 15 '12 at 2:12
• I added a couple more things to check. If that doesn't work, send code snippets and the full output to petsc-maint@mcs.anl.gov. – Jed Brown Jul 15 '12 at 2:21
• 1) I have checked my linear operator. Ax does indeed give me correct b. I took known x and b`. I am also getting second order convergence on $r = b - Ax$ (I am using second order finite difference scheme for discretization.) – amneet Jul 15 '12 at 3:19