On Tue, Nov 30, 2010 at 7:13 PM, Rongliang Chen <span dir="ltr"><<a href="mailto:rongliang.chan@gmail.com">rongliang.chan@gmail.com</a>></span> wrote:<br><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
Hi,<br><br>I have tried the following method: <br><br> Vec Ax;<br> ierr = VecDuplicate(x, &Ax);CHKERRQ(ierr);<br> ierr = MatMult(*A, x, Ax);CHKERRQ(ierr); <br> PetscViewer viewer = PETSC_VIEWER_BINARY_(PETSC_COMM_WORLD);<br>
ierr = MatView(*A,viewer);CHKERRQ(ierr);<br> ierr = VecView(x,viewer);CHKERRQ(ierr);<br> ierr = VecView(Ax,viewer);CHKERRQ(ierr);<br>// 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.)).<br>
<br>After do two SNES, I use the following commands in matlab:<br><br>>> [A,x,y,A1,x1,y1] = PetscBinaryRead('binaryoutput');<br>>> norm(A*x-y)<br><br>ans =<br><br> 1.2000e-13<br><br>>> norm(A1*x1-y1)<br>
<br>ans =<br><br> 2.6612e+03<br> <br>The reason that I want to do this is that I doubt that there is a bug in my code <br>and I want to following these information to find out the bug.<br></blockquote><div><br></div><div>
When you create the Jacobian, do you remember to zero out the matrix? This is a common</div><div>bug when creating the matrix after the initial iterate.</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
There is a very strange thing in my code:(I use GMRES with ASM preconditioner to solve the linear system.)<br>
<br>Solving Shape Optimization problem<br> 0 SNES norm 1.1359606379e+02, 0 KSP its last norm 0.0000000000e+00.<br> 1 SNES norm 4.1335342946e+01, 1 KSP its last norm 2.3855348597e-03.<br> 2 SNES norm <b>4.1085343652e+01, 13 KSP</b> its last norm 3.4528582271e-02.<br>
<i>SNES diverged DIVERGED_LS_FAILURE.</i><br><br>I think there is some thing wrong with the second SNES, because its KSP should <br>be 1 or 2. If I destroy the Jacobian matrix here and recreate it, then continue the <br>
SNES, it will be alright:<br>
<br> 0 SNES norm <b>4.1085343652e+01,</b> 0 KSP its last norm 0.0000000000e+00.<br> 1 SNES norm 6.9950501590e+00, 1 KSP its last norm 2.8278185399e-03.<br> 2 SNES norm 1.3445253365e-01, 1 KSP its last norm 2.3243944231e-04.<br>
3 SNES norm 5.0914645029e-04, 1 KSP its last norm 5.0954271956e-06.<br>.................<br>So I output the matrix and vectors into matlab and do the matrix and vector <br>multiplication and I found that the result is different from the MatMult in Petsc.<br>
<br>Best,<br><br>Rongliang<br><br><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0pt 0pt 0pt 0.8ex;border-left:1px solid rgb(204, 204, 204);padding-left:1ex">
------------------------------<br>
<br>
Message: 2<br>
Date: Tue, 30 Nov 2010 21:52:45 +0100<br>
From: Jed Brown <jed@59A2.org><br>
Subject: Re: [petsc-users] Why MatMult(A, x, Ax) != A*x?<br>
To: PETSc users list <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>><br>
Message-ID:<br>
<<a href="mailto:AANLkTiks1Go4J56qYDuFabb%2BQQBtx0DmfaythUkzgJik@mail.gmail.com" target="_blank">AANLkTiks1Go4J56qYDuFabb+QQBtx0DmfaythUkzgJik@mail.gmail.com</a>><br>
Content-Type: text/plain; charset="utf-8"<br>
<br>
On Tue, Nov 30, 2010 at 21:42, Rongliang Chen <<a href="mailto:rongliang.chan@gmail.com" target="_blank">rongliang.chan@gmail.com</a>>wrote:<br>
<br>
> I know that I should call PetscViewerDestroy() when I have done with each<br>
> viewer. But if I called this I just can view the last Matrix and vector.<br>
> What I want<br>
> to do is that I let it do three steps of SNES, view the jacobian of each<br>
> step.<br>
><br>
<br>
Maybe view it to different files. The code you posted for viewing is<br>
certainly not correct, especially if called in a loop.<br>
<br>
<br>
> I found that A*x is correct at the first step and it's wrong for the second<br>
> and<br>
> third step.<br>
> I also tried using the binary output and input, and it has the same<br>
> results.<br>
><br>
<br>
The problem is most likely that you are picking these pieces out of the<br>
files differently than they are put in. MatMult being incorrect is by far<br>
the least likely scenario. I don't know why you want to do this, but you<br>
need to find a way to manage the files so that you can be certain that you<br>
are reading (from Matlab) exactly what you think you are reading. One way<br>
to do this is to use a binary file that gets all three items together.<br>
<br>
Jed<br>
<br>
<br></blockquote></div>
</blockquote></div><br><br clear="all"><br>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener<br>