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

[Jed Brown]

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>