[petsc-users] Why MatMult(A, x, Ax) != A*x?

[Matthew Knepley]

On Tue, Nov 30, 2010 at 7:13 PM, Rongliang Chen <rongliang.chan at gmail.com>wrote:

> Hi,
>
> I have tried the following method:
>
>         Vec Ax;
>         ierr = VecDuplicate(x, &Ax);CHKERRQ(ierr);
>         ierr = MatMult(*A, x, Ax);CHKERRQ(ierr);
>         PetscViewer viewer = PETSC_VIEWER_BINARY_(PETSC_COMM_WORLD);
>         ierr = MatView(*A,viewer);CHKERRQ(ierr);
>         ierr = VecView(x,viewer);CHKERRQ(ierr);
>         ierr = VecView(Ax,viewer);CHKERRQ(ierr);
> //      ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);   (If I use this
> line, I just can obtain A,x,y and I can not get the A1, x1, y1 (see
> below.)).
>
> After do two SNES, I use the following commands in matlab:
>
> >> [A,x,y,A1,x1,y1] = PetscBinaryRead('binaryoutput');
> >> norm(A*x-y)
>
> ans =
>
>    1.2000e-13
>
> >> norm(A1*x1-y1)
>
> ans =
>
>    2.6612e+03
>
> The reason that I want to do this is that I doubt that there is a bug in my
> code
> and I want to following these information to find out the bug.
>

When you create the Jacobian, do you remember to zero out the matrix? This
is a common
bug when creating the matrix after the initial iterate.

   Matt


> There is a very strange thing in my code:(I use GMRES with ASM
> preconditioner to solve the linear system.)
>
> Solving Shape Optimization problem
>   0 SNES norm 1.1359606379e+02, 0 KSP its last norm 0.0000000000e+00.
>   1 SNES norm 4.1335342946e+01, 1 KSP its last norm 2.3855348597e-03.
>   2 SNES norm *4.1085343652e+01, 13 KSP* its last norm 3.4528582271e-02.
> *SNES diverged DIVERGED_LS_FAILURE.*
>
> I think there is some thing wrong with the second SNES, because its KSP
> should
> be 1 or 2. If I destroy the Jacobian matrix here and recreate it, then
> continue the
> SNES, it will be alright:
>
>   0 SNES norm *4.1085343652e+01,* 0 KSP its last norm 0.0000000000e+00.
>   1 SNES norm 6.9950501590e+00, 1 KSP its last norm 2.8278185399e-03.
>   2 SNES norm 1.3445253365e-01,  1 KSP its last norm 2.3243944231e-04.
>   3 SNES norm 5.0914645029e-04,  1 KSP its last norm 5.0954271956e-06.
> .................
> So I output the matrix and vectors into matlab and do the matrix and vector
>
> multiplication and I found that the result is different from the MatMult in
> Petsc.
>
> Best,
>
> Rongliang
>
> ------------------------------
>>
>> Message: 2
>> Date: Tue, 30 Nov 2010 21:52:45 +0100
>> From: Jed Brown <jed at 59A2.org>
>> Subject: Re: [petsc-users] Why MatMult(A, x, Ax) != A*x?
>> To: PETSc users list <petsc-users at mcs.anl.gov>
>> Message-ID:
>>        <AANLkTiks1Go4J56qYDuFabb+QQBtx0DmfaythUkzgJik at mail.gmail.com<AANLkTiks1Go4J56qYDuFabb%2BQQBtx0DmfaythUkzgJik at mail.gmail.com>
>> >
>> Content-Type: text/plain; charset="utf-8"
>>
>> On Tue, Nov 30, 2010 at 21:42, Rongliang Chen <rongliang.chan at gmail.com
>> >wrote:
>>
>> > I know that I should call PetscViewerDestroy() when I have done with
>> each
>> > viewer. But if I called this I just can view the last Matrix and vector.
>> > What I want
>> > to do is that I let it do three steps of SNES, view the jacobian of each
>> > step.
>> >
>>
>> Maybe view it to different files.  The code you posted for viewing is
>> certainly not correct, especially if called in a loop.
>>
>>
>> > I found that A*x is correct at the first step and it's wrong for the
>> second
>> > and
>> > third step.
>> > I also tried using the binary output and input, and it has the same
>> > results.
>> >
>>
>> The problem is most likely that you are picking these pieces out of the
>> files differently than they are put in.  MatMult being incorrect is by far
>> the least likely scenario.  I don't know why you want to do this, but you
>> need to find a way to manage the files so that you can be certain that you
>> are reading (from Matlab) exactly what you think you are reading.  One way
>> to do this is to use a binary file that gets all three items together.
>>
>> Jed
>>
>>
>>


-- 
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20101130/d96960c7/attachment.htm>