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

[Barry Smith]

On Nov 30, 2010, at 7:13 PM, Rongliang Chen 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.
> 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.

  You are going about this all wrong. Forget about checking that A*x = A*x. That is not the bug and not the way to find the bug. 

  Run with -snes_monitor -ksp_monitor_true_residual -snes_converged_reason -ksp_converged_reason  -snes_ls_monitor 

likely either the Jacobian is wrong or the linear solver is not working. Look at the output and if confused send the output to petsc-maint at mcs.anl.gov

  Barry

> 
> 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>
> 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
> 
>