[petsc-users] Problem with assembling 2 matrices

[Jonas Mairhofer]

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