Hi, I am having trouble with using the function MatMult in PETSc. The result of MatMult(A, x, Ax) (*using multiprocessors*) is not equal to the result of A*x computed in Matlab. 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); ..................... What kind of bug can cause this problem? Thanks. Best, Rongliang
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
On Sun, Nov 28, 2010 at 21:18, Jed Brown <[email protected]> wrote:
PetscViewer viewer = PETSC_VIEWER_BINARY(PETSC_COMM_WORLD);
This line should read: PetscViewer viewer = PETSC_VIEWER_BINARY_(PETSC_COMM_WORLD); You can use PetscViewerBinaryOpen for more control, and to specify the filename. Jed
Hi, I know that I should call PetscViewerDestroy() when I have done with each viewer. But if I called this I just can view the last Matrix and vector. What I want to do is that I let it do three steps of SNES, view the jacobian of each step. I found that A*x is correct at the first step and it's wrong for the second and third step. I also tried using the binary output and input, and it has the same results. But If I destroy the matrix A after the first SNES, recreate it and then do the second and third SNES, it will be alright. Thanks. Best, Rongliang
------------------------------
Message: 3 Date: Tue, 30 Nov 2010 09:29:17 +0100 From: Jed Brown <[email protected]> Subject: Re: [petsc-users] petsc-users Digest, Vol 23, Issue 32 To: PETSc users list <[email protected]> Message-ID: <[email protected]> Content-Type: text/plain; charset="utf-8"
On Tue, Nov 30, 2010 at 03:26, Rongliang Chen <[email protected]
wrote:
[A,x,y] = PetscBinaryRead('binaryoutput'), can this command load the binary files into matlab matrix and vector?
Yes, this reads the binary files. You will need to add $PETSC_DIR/bin/matlab to your path (see "help addpath").
You should also call PetscViewerDestroy() when you are done with each viewer. In the source you quoted, this is never called.
Jed
On Tue, Nov 30, 2010 at 21:42, Rongliang Chen <[email protected]>wrote:
I know that I should call PetscViewerDestroy() when I have done with each viewer. But if I called this I just can view the last Matrix and vector. What I want to do is that I let it do three steps of SNES, view the jacobian of each step.
Maybe view it to different files. The code you posted for viewing is certainly not correct, especially if called in a loop.
I found that A*x is correct at the first step and it's wrong for the second and third step. I also tried using the binary output and input, and it has the same results.
The problem is most likely that you are picking these pieces out of the files differently than they are put in. MatMult being incorrect is by far the least likely scenario. I don't know why you want to do this, but you need to find a way to manage the files so that you can be certain that you are reading (from Matlab) exactly what you think you are reading. One way to do this is to use a binary file that gets all three items together. Jed
participants (2)
-
Jed Brown -
Rongliang Chen