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

[Rongliang Chen]

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