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)
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
