On Sun, Nov 28, 2010 at 21:11, Rongliang Chen
<rongliang.chan@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