Dear fellow PETSc users,

I've been trying to subtract one matrix from another, i.e. Y = Y - X, using 
`call MatAXPY(Y, minus_one, X, SAME_NONZERO_PATTERN, ierr)`. In my example,
some entries come out correct (e.g. at [0,0]), some come out wrong 
(e.g. at [29,29]). Here is the complete code, in Fortran:

program main

    implicit none

#include "finclude/petscsys.h"
#include "finclude/petscvec.h"
#include "finclude/petscmat.h"
#include "finclude/petscviewer.h"

    Mat :: X, Y
    PetscReal :: minus_one = -1.d0
    PetscErrorCode :: ierr
    PetscViewer :: X_viewer, Y_viewer

    call PetscInitialize(PETSC_NULL_CHARACTER,ierr); CHKERRQ(ierr)

    call MatCreate(PETSC_COMM_WORLD,X,ierr); CHKERRQ(ierr)
    call MatCreate(PETSC_COMM_WORLD,Y,ierr); CHKERRQ(ierr)
    
    call PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mat1.dat", FILE_MODE_READ, 
X_viewer, ierr); CHKERRQ(ierr)
    call PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mat2.dat", FILE_MODE_READ, 
Y_viewer, ierr); CHKERRQ(ierr)

    call MatLoad(X, X_viewer, ierr); CHKERRQ(ierr)
    call MatLoad(Y, Y_viewer, ierr); CHKERRQ(ierr)
    call PetscViewerDestroy(X_viewer, ierr); CHKERRQ(ierr)
    call PetscViewerDestroy(Y_viewer, ierr); CHKERRQ(ierr)

    call MatView(X, PETSC_VIEWER_STDOUT_SELF, ierr); CHKERRQ(ierr)
    call MatView(Y, PETSC_VIEWER_STDOUT_SELF, ierr); CHKERRQ(ierr)
    call MatAXPY(Y, minus_one, X, SAME_NONZERO_PATTERN, ierr); CHKERRQ(ierr)
    call MatView(Y, PETSC_VIEWER_STDOUT_SELF, ierr); CHKERRQ(ierr)

    call MatDestroy(X, ierr); CHKERRQ(ierr)
    call MatDestroy(Y, ierr); CHKERRQ(ierr)

    call PetscFinalize(ierr); CHKERRQ(ierr)
end program

Attachment: mat1.dat
Description: Netscape Proxy Auto Config

Attachment: mat2.dat
Description: Netscape Proxy Auto Config

Reply via email to