On Wed, Jun 11, 2014 at 6:47 AM, Jonas Mairhofer
<[email protected]
<mailto:[email protected]>> wrote:
Dear PETSc-Team,
I am trying to set a global array X as column 0 of a dense matrix
A and another global array Y as the first column of a second dense
matrix B.
Everything works perfectly fine in serial and parallel as long as
I do this only for one of the matrices. However, when I try to do
it for both matrices, the entries of matrix A seem to be
overwritten by the values of matrix B.
This is the output I get (first vectors X and Y, then matrices A
and B):
Vector Object:Vec_0x84000000_0 1 MPI processes
type: mpi
Process [0]
1
1
1
Vector Object:Vec_0x84000000_1 1 MPI processes
type: mpi
Process [0]
2
2
2
Matrix Object: 1 MPI processes
type: seqdense
2.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
2.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
2.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
Matrix Object: 1 MPI processes
type: seqdense
2.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
2.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
2.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
I also attached the Fortran code:
program main
implicit none
#include <finclude/petscsys.h>
#include <finclude/petscvec.h>
#include <finclude/petscdmda.h>
#include <finclude/petscis.h>
#include <finclude/petscmat.h>
#include <finclude/petscksp.h>
#include <finclude/petscpc.h>
#include <finclude/petscsnes.h>
#include <finclude/petscvec.h90>
#include <finclude/petscdmda.h90>
PetscInt i,j,m,n
PetscErrorCode ierr
PetscScalar one
Mat A
Mat B
PetscInt row(3)
DM da
Vec X
Vec Y
PetscScalar,pointer :: X_pointer(:)
PetscScalar,pointer :: Y_pointer(:)
PetscInt xs,xm
!Initialize
call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
m = 3 !col
n = 3 !row
call DMDACreate1d(PETSC_COMM_WORLD, & !MPI communicator
& DMDA_BOUNDARY_NONE, & !Boundary type at
boundary of physical domain
& n, & !global dimension of array (if negative
number, then value can be changed by user via command line!)
& 1, & !number of degrees of freedom per grid
point (number of unknowns at each grid point)
& 0, & !number of ghost points accessible to
local vectors
& PETSC_NULL_INTEGER, & !could be an array to
specify the number of grid points per processor
& da, & !the resulting distributed array object
& ierr)
call DMDAGetCorners(da, & !the distributed array
& xs, & !corner index in x direction
& PETSC_NULL_INTEGER, & !corner index in y direction
& PETSC_NULL_INTEGER, & !corner index in z direction
& xm, & !width of locally owned part in x
direction
& PETSC_NULL_INTEGER, & !width of locally owned part
in y direction
& PETSC_NULL_INTEGER, & !width of locally owned part
in z direction
& ierr) !error check
one = 1.0
!set indices of matrix rows to be set
Do i = xs,xs+xm-1
row(i) = i
END DO
!Set up vector X and Matrix A
call DMCreateGlobalVector(da,X,ierr)
call VecSet(X,one,ierr) !set all entries of global
array to 1
call VecGetArrayF90(X,X_pointer,ierr)
call
MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,m,PETSC_NULL_CHARACTER,A,ierr)
The bug is here (and in the other create call). You should not be
using PETSC_NULL_CHARACTER, but PETSC_NULL_SCALAR.
This is a good reason for reconsidering the use of Fortran. It just
cannot warn you about damaging errors like this.
Thanks,
Matt
call
MatSetValues(A,xm,row(xs:xs+xm-1),1,0,X_pointer(1:xm),INSERT_VALUES,ierr)
call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
!Set up vector Y and matrix B
call DMCreateGlobalVector(da,Y,ierr)
call VecSet(Y,2.d0*one,ierr) !set all entries of global
array to 2
call VecGetArrayF90(Y,Y_pointer,ierr)
call
MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,m,PETSC_NULL_CHARACTER,B,ierr)
call
MatSetValues(B,xm,row(xs:xs+xm-1),1,0,Y_pointer(1:xm),INSERT_VALUES,ierr)
call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr)
call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr)
!print out matrices and vectors to check entries
call VecView(X,PETSC_VIEWER_STDOUT_WORLD,ierr)
call VecView(Y,PETSC_VIEWER_STDOUT_WORLD,ierr)
call MatView(A,PETSC_VIEWER_STDOUT_WORLD,ierr)
call MatView(B,PETSC_VIEWER_STDOUT_WORLD,ierr)
!free objects
call VecRestoreArrayF90(X,X_pointer,ierr)
call VecDestroy(X,ierr)
call MatDestroy(A,ierr)
call VecRestoreArrayF90(Y,Y_pointer,ierr)
call VecDestroy(Y,ierr)
call MatDestroy(B,ierr)
call PetscFinalize(ierr)
end
Am I doing something wrong in my code?
Thank you for your help!
Jonas
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which
their experiments lead.
-- Norbert Wiener
