[petsc-users] Why MatMult(A, x, Ax) != A*x?
Rongliang Chen
rongliang.chan at gmail.com
Tue Nov 30 19:13:35 CST 2010
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>
More information about the petsc-users
mailing list