Hi Jed,<br><br>Thank you for your reply.<br>The Matrix A and Vector x are created by MatCreateMPIAIJ and VecCreateMPI and<br>matrix A is obtained by SNESDefaultComputeJacobian. <br>I compared the results using the following routine:<br>
<br> ..........................................................<br> sprintf(filename,"A.m");<br> ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&view->viewer);CHKERRQ(ierr);<br>
ierr = PetscViewerSetFormat(view->viewer, PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr); <br> ierr = MatView(A, view->viewer);CHKERRQ(ierr); <br><br> Vec y;<br> ierr = VecDuplicate(x, &y);CHKERRQ(ierr);<br>
ierr = MatMult(A, x, y);CHKERRQ(ierr);<br><br> sprintf(filename,"x.m");<br> ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&view->viewer);CHKERRQ(ierr);<br> ierr = PetscViewerSetFormat(view->viewer, PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr); <br>
ierr = VecView(x, view->viewer);CHKERRQ(ierr);<br><br> sprintf(filename,"y.m");<br> ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&view->viewer);CHKERRQ(ierr);<br> ierr = PetscViewerSetFormat(view->viewer, PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr); <br>
ierr = VecView(y, view->viewer);CHKERRQ(ierr); <br> ...........................................................<br><br>In the Matlab, I first load the files A.m, x.m and y.m into A, x and y, and then do "error=norm(A*x - y)".<br>
The 'error' is very big(about 200). <br><br>Following your suggestion, I can write the matrix and vectors into binary files, but I do not how to load <br>these binary files into matlab matrix and vector. <br>I do not understand the command: [A,x,y] = PetscBinaryRead('binaryoutput'), can this command load <br>
the binary files into matlab matrix and vector? Thanks.<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>
Message: 2<br>
Date: Sun, 28 Nov 2010 21:18:59 +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">petsc-users@mcs.anl.gov</a>><br>
Message-ID:<br>
<AANLkTikGhD+2WcBmaZOR8xfT=<a href="mailto:D2Vh5xeK4HYHY_bG2vM@mail.gmail.com">D2Vh5xeK4HYHY_bG2vM@mail.gmail.com</a>><br>
Content-Type: text/plain; charset="utf-8"<br>
<br>
On Sun, Nov 28, 2010 at 21:11, Rongliang Chen <<a href="mailto:rongliang.chan@gmail.com">rongliang.chan@gmail.com</a>>wrote:<br>
<br>
> The format of the matrix A is AIJ and is obtained by function<br>
> SNESDefaultComputeJacobian.<br>
> I compute A*x in Matlab with the following A and x and compare it with Ax.<br>
><br>
> .......................<br>
> ierr = MatMult(*A, x, Ax);CHKERRQ(ierr);<br>
><br>
> sprintf(filename,"x.m");<br>
> ierr =<br>
> PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&view->viewer);CHKERRQ(ierr);<br>
> ierr = PetscViewerSetFormat(view->viewer,<br>
> PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);<br>
> ierr = VecView(x, view->viewer);CHKERRQ(ierr);<br>
><br>
> sprintf(filename,"Ax.m");<br>
> ierr =<br>
> PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&view->viewer);CHKERRQ(ierr);<br>
> ierr = PetscViewerSetFormat(view->viewer,<br>
> PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);<br>
> ierr = VecView(Ax, view->viewer);CHKERRQ(ierr);<br>
><br>
<br>
Are the vectors and matrices obtained from a DA (DACreateGlobalVector,<br>
DAGetMatrix)? How are you comparing the matrix that PETSc sees with the<br>
matrix that Matlab sees? I suggest writing the matrix and both vectors to a<br>
PETSc binary file<br>
<br>
PetscViewer viewer = PETSC_VIEWER_BINARY(PETSC_COMM_WORLD);<br>
MatView(A,viewer);<br>
VecView(x,viewer);<br>
MatMult(A,x,y);<br>
VecView(y,viewer);<br>
<br>
and read them into Matlab with<br>
<br>
[A,x,y] = PetscBinaryRead('binaryoutput')<br>
norm(A*x - y) % This should be small<br>
<br>
Jed<br>
-------------- next part --------------<br>
An HTML attachment was scrubbed...<br>
URL: <<a href="http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20101128/d97eefaa/attachment-0001.htm" target="_blank">http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20101128/d97eefaa/attachment-0001.htm</a>><br>
<br></blockquote></div>