[petsc-users] Why MatMult(A, x, Ax) != A*x?
Jed Brown
jed at 59A2.org
Sun Nov 28 14:18:59 CST 2010
On Sun, Nov 28, 2010 at 21:11, Rongliang Chen <rongliang.chan at gmail.com>wrote:
> The format of the matrix A is AIJ and is obtained by function
> SNESDefaultComputeJacobian.
> I compute A*x in Matlab with the following A and x and compare it with Ax.
>
> .......................
> ierr = MatMult(*A, x, Ax);CHKERRQ(ierr);
>
> sprintf(filename,"x.m");
> ierr =
> PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&view->viewer);CHKERRQ(ierr);
> ierr = PetscViewerSetFormat(view->viewer,
> PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);
> ierr = VecView(x, view->viewer);CHKERRQ(ierr);
>
> sprintf(filename,"Ax.m");
> ierr =
> PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&view->viewer);CHKERRQ(ierr);
> ierr = PetscViewerSetFormat(view->viewer,
> PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);
> ierr = VecView(Ax, view->viewer);CHKERRQ(ierr);
>
Are the vectors and matrices obtained from a DA (DACreateGlobalVector,
DAGetMatrix)? How are you comparing the matrix that PETSc sees with the
matrix that Matlab sees? I suggest writing the matrix and both vectors to a
PETSc binary file
PetscViewer viewer = PETSC_VIEWER_BINARY(PETSC_COMM_WORLD);
MatView(A,viewer);
VecView(x,viewer);
MatMult(A,x,y);
VecView(y,viewer);
and read them into Matlab with
[A,x,y] = PetscBinaryRead('binaryoutput')
norm(A*x - y) % This should be small
Jed
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20101128/d97eefaa/attachment.htm>
More information about the petsc-users
mailing list