[petsc-users] petsc-users Digest, Vol 23, Issue 32

[Rongliang Chen]

Hi Jed,

Thank you for your reply.
The Matrix A and Vector x are created by MatCreateMPIAIJ and VecCreateMPI
and
matrix A is obtained by SNESDefaultComputeJacobian.
I compared the results using the following routine:

        ..........................................................
        sprintf(filename,"A.m");
        ierr =
PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&view->viewer);CHKERRQ(ierr);
        ierr = PetscViewerSetFormat(view->viewer,
PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);
        ierr = MatView(A, view->viewer);CHKERRQ(ierr);

        Vec y;
        ierr = VecDuplicate(x, &y);CHKERRQ(ierr);
        ierr = MatMult(A, x, y);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,"y.m");
        ierr =
PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&view->viewer);CHKERRQ(ierr);
        ierr = PetscViewerSetFormat(view->viewer,
PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);
        ierr = VecView(y, view->viewer);CHKERRQ(ierr);
        ...........................................................

In the Matlab, I first load the files A.m, x.m and y.m into A, x and y, and
then do "error=norm(A*x - y)".
The 'error' is very big(about 200).

Following your suggestion, I can write the matrix and vectors into binary
files, but I do not how to load
these binary files into matlab matrix and vector.
I do not understand the command:   [A,x,y] =
PetscBinaryRead('binaryoutput'), can this command load
the binary files into matlab matrix and vector? Thanks.

Best,

Rongliang


> Message: 2
> Date: Sun, 28 Nov 2010 21:18:59 +0100
> From: Jed Brown <jed at 59A2.org>
> Subject: Re: [petsc-users] Why MatMult(A, x, Ax) != A*x?
> To: PETSc users list <petsc-users at mcs.anl.gov>
> Message-ID:
>        <AANLkTikGhD+2WcBmaZOR8xfT=D2Vh5xeK4HYHY_bG2vM at mail.gmail.com>
> Content-Type: text/plain; charset="utf-8"
>
> 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-0001.htm
> >
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20101129/1f0c733b/attachment-0001.htm>