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

28 Nov 2010


      On Sun, Nov 28, 2010 at 21:11, Rongliang Chen <[email protected]>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

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

Jed Brown