Matrixfree modification

[Matthew Knepley]

On Tue, Oct 28, 2008 at 4:35 PM, Johann Schlamp <schlamp at in.tum.de> wrote:
> Barry Smith wrote:
>>    If you have the "matrix-based" version that you can run on the same
>> problem then
>> look at the residuals computed in the 0th and 1st step of the Krylov
>> method
>> to see how they are different in the two runs.
> The residuals in the 0th and 1st step of the linear solver are
> 0.00363413
> 0.00189276
> for the "matrix-based" version, and
> 0.00363413
> 1.27858e-17

No this looks wrong. Shouldn't these be identical? It looks like you
are wiping out the input vector instead.

   Matt

> for the matrix-free version. That's definitely smaller than epsilon, so
> it converges. By the way, the "matrix-based" version doesn't converge
> either, as I was not using a preconditioner for getting comparable results.
>
> Simply thought: the residual is in the magnitude of machine accuracy, so
> I would have concluded that the calculation of the residual (y-A*x)
> results in zero with respect to some rounding errors. Unfortunately, I
> don't completely understand the PETSc code for calculating the residual
> and therfore cannot verify it for my new matrix structure.
>
>>    Perhaps your matrix-free is corrupting memory? Run with -malloc_debug
>> and put a CHKMEMQ; at the end of your matrix free multiply. Or better
>> run through
>> valgrind,. www.valgrind.org
> Interesting thought! I will check this tomorrow.
>
>
> Johann
>
>
>> On Oct 28, 2008, at 3:39 PM, Johann Schlamp wrote:
>>
>>> Thanks for your reply, Barry!
>>>
>>> Barry Smith wrote:
>>>>   Your general approach seems fine. I would put a break point in
>>>> your MatMult routine for a tiny problem and verify that the input
>>>> vector for u and
>>>> x are what you expect (and then at the end of the function make sure
>>>> the
>>>> output y is what you expect).
>>> I have already done this, everything is as expected.
>>>
>>>>   Here is my guess: The matrix vector product is J(u)*x; when your
>>>> calculation
>>>> is done you need to use the correct u value. This is the vector that is
>>>> passed into your "empty method as user-provided Jacobian method".
>>>> In your "empty method as user-provided Jacobian method" you should make
>>>> a copy of this vector so that you have it for each of the matrix
>>>> vector products.
>>>> At each Newton step your "empty method as user-provided Jacobian
>>>> method"
>>>> will be called and you will copy the new value of u over.
>>> It took me some time to understand what you mean. But after that I got
>>> excited trying it out. :)
>>> It definitely had some effect on the nonlinear solution, but the linear
>>> solver still finishes after one iteration with the way too small
>>> residual norm.
>>>
>>> Anyway, thanks for the anticipated bugfix!
>>>
>>> Do you have any further suggestions?
>>>
>>>
>>> Best regards,
>>> Johann
>>>
>>>
>>>> On Oct 28, 2008, at 10:08 AM, Johann Schlamp <schlamp at in.tum.de> wrote:
>>>>
>>>>> Hello folks,
>>>>>
>>>>> I have to implement some advanced matrixfree calculation.
>>>>>
>>>>> Usually, we use our own analytical method for calculating the Jacobian
>>>>> matrix and provide it via SNESSetJacobian(). For complex problems, the
>>>>> global matrix takes too much memory. So here's the idea:
>>>>>
>>>>> First, I set some empty method as user-provided Jacobian method.
>>>>> After a
>>>>> corresponding call, SNES hopefully thinks the Jacobian matrix got
>>>>> calculated right. Then it will use it in some multiplication like
>>>>> y=A*x.
>>>>> For that I created the matrix A with MatCreateShell() and set up my
>>>>> own
>>>>> MatMult method. This new method iterates over my grid, calculates
>>>>> local
>>>>> Jacobians, multiplies them with the corresponding part of the vector x
>>>>> and writes them into y. After this full iteration, it should look
>>>>> like I
>>>>> had the full global Jacobian and multiplied it with the full vector x.
>>>>> In y, the result will be the same.
>>>>> The matrix A and the empty Jacobian method are set through
>>>>> SNESSetJacobian().
>>>>>
>>>>> I have implemented some unittests for generally proofing the idea,
>>>>> seems
>>>>> to work.
>>>>>
>>>>> But if I run a complete simulation, the KSP converges after the first
>>>>> iteration with a residual norm of about 1.0e-18, which is
>>>>> definitely not
>>>>> right.
>>>>>
>>>>>
>>>>> Now my question: does this procedure make sense at all, and if so - is
>>>>> it possible that just the calculation of the residual norm goes wrong
>>>>> due to my new matrix structure? I searched the PETSc code, but I
>>>>> wasn't
>>>>> able to find a proof or solution for that.
>>>>>
>>>>> Any help would be appreciated.
>>>>>
>>>>>
>>>>> Best regards,
>>>>> Johann Schlamp
>>>>>
>>>>>
>>>>
>>>
>>
>
>



-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which
their experiments lead.
-- Norbert Wiener