<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Wed, Jun 11, 2014 at 6:47 AM, Jonas Mairhofer <span dir="ltr"><<a href="mailto:mairhofer@itt.uni-stuttgart.de" target="_blank">mairhofer@itt.uni-stuttgart.de</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Dear PETSc-Team,<br>
<br>
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.<br>
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.<br>
<br>
This is the output I get (first vectors X and Y, then matrices A and B):<br>
<br>
Vector Object:Vec_0x84000000_0 1 MPI processes<br>
type: mpi<br>
Process [0]<br>
1<br>
1<br>
1<br>
Vector Object:Vec_0x84000000_1 1 MPI processes<br>
type: mpi<br>
Process [0]<br>
2<br>
2<br>
2<br>
Matrix Object: 1 MPI processes<br>
type: seqdense<br>
2.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00<br>
2.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00<br>
2.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00<br>
Matrix Object: 1 MPI processes<br>
type: seqdense<br>
2.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00<br>
2.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00<br>
2.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00<br>
<br>
I also attached the Fortran code:<br>
<br>
<br>
<br>
<br>
program main<br>
implicit none<br>
<br>
#include <finclude/petscsys.h><br>
#include <finclude/petscvec.h><br>
#include <finclude/petscdmda.h><br>
#include <finclude/petscis.h><br>
#include <finclude/petscmat.h><br>
#include <finclude/petscksp.h><br>
#include <finclude/petscpc.h><br>
#include <finclude/petscsnes.h><br>
#include <finclude/petscvec.h90><br>
#include <finclude/petscdmda.h90><br>
<br>
PetscInt i,j,m,n<br>
PetscErrorCode ierr<br>
<br>
PetscScalar one<br>
Mat A<br>
Mat B<br>
<br>
PetscInt row(3)<br>
DM da<br>
Vec X<br>
Vec Y<br>
PetscScalar,pointer :: X_pointer(:)<br>
PetscScalar,pointer :: Y_pointer(:)<br>
PetscInt xs,xm<br>
<br>
!Initialize<br>
call PetscInitialize(PETSC_NULL_<u></u>CHARACTER,ierr)<br>
<br>
m = 3 !col<br>
n = 3 !row<br>
<br>
<br>
call DMDACreate1d(PETSC_COMM_WORLD, & !MPI communicator<br>
& DMDA_BOUNDARY_NONE, & !Boundary type at boundary of physical domain<br>
& n, & !global dimension of array (if negative number, then value can be changed by user via command line!)<br>
& 1, & !number of degrees of freedom per grid point (number of unknowns at each grid point)<br>
& 0, & !number of ghost points accessible to local vectors<br>
& PETSC_NULL_INTEGER, & !could be an array to specify the number of grid points per processor<br>
& da, & !the resulting distributed array object<br>
& ierr)<br>
<br>
call DMDAGetCorners(da, & !the distributed array<br>
& xs, & !corner index in x direction<br>
& PETSC_NULL_INTEGER, & !corner index in y direction<br>
& PETSC_NULL_INTEGER, & !corner index in z direction<br>
& xm, & !width of locally owned part in x direction<br>
& PETSC_NULL_INTEGER, & !width of locally owned part in y direction<br>
& PETSC_NULL_INTEGER, & !width of locally owned part in z direction<br>
& ierr) !error check<br>
<br>
<br>
<br>
one = 1.0<br>
<br>
<br>
!set indices of matrix rows to be set<br>
Do i = xs,xs+xm-1<br>
row(i) = i<br>
END DO<br>
<br>
<br>
!Set up vector X and Matrix A<br>
call DMCreateGlobalVector(da,X,<u></u>ierr)<br>
call VecSet(X,one,ierr) !set all entries of global array to 1<br>
call VecGetArrayF90(X,X_pointer,<u></u>ierr)<br>
<br>
call MatCreateDense(PETSC_COMM_<u></u>WORLD,PETSC_DECIDE,PETSC_<u></u>DECIDE,n,m,PETSC_NULL_<u></u>CHARACTER,A,ierr) <br></blockquote><div><br></div><div>The bug is here (and in the other create call). You should not be using PETSC_NULL_CHARACTER, but PETSC_NULL_SCALAR.</div>
<div><br></div><div>This is a good reason for reconsidering the use of Fortran. It just cannot warn you about damaging errors like this.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
call MatSetValues(A,xm,row(xs:xs+<u></u>xm-1),1,0,X_pointer(1:xm),<u></u>INSERT_VALUES,ierr)<br>
call MatAssemblyBegin(A,MAT_FINAL_<u></u>ASSEMBLY,ierr)<br>
call MatAssemblyEnd(A,MAT_FINAL_<u></u>ASSEMBLY,ierr)<br>
<br>
<br>
!Set up vector Y and matrix B<br>
call DMCreateGlobalVector(da,Y,<u></u>ierr)<br>
call VecSet(Y,2.d0*one,ierr) !set all entries of global array to 2<br>
call VecGetArrayF90(Y,Y_pointer,<u></u>ierr)<br>
<br>
call MatCreateDense(PETSC_COMM_<u></u>WORLD,PETSC_DECIDE,PETSC_<u></u>DECIDE,n,m,PETSC_NULL_<u></u>CHARACTER,B,ierr) <br>
call MatSetValues(B,xm,row(xs:xs+<u></u>xm-1),1,0,Y_pointer(1:xm),<u></u>INSERT_VALUES,ierr)<br>
call MatAssemblyBegin(B,MAT_FINAL_<u></u>ASSEMBLY,ierr)<br>
call MatAssemblyEnd(B,MAT_FINAL_<u></u>ASSEMBLY,ierr)<br>
<br>
<br>
!print out matrices and vectors to check entries<br>
call VecView(X,PETSC_VIEWER_STDOUT_<u></u>WORLD,ierr)<br>
call VecView(Y,PETSC_VIEWER_STDOUT_<u></u>WORLD,ierr)<br>
call MatView(A,PETSC_VIEWER_STDOUT_<u></u>WORLD,ierr)<br>
call MatView(B,PETSC_VIEWER_STDOUT_<u></u>WORLD,ierr)<br>
<br>
!free objects<br>
call VecRestoreArrayF90(X,X_<u></u>pointer,ierr)<br>
call VecDestroy(X,ierr)<br>
call MatDestroy(A,ierr)<br>
call VecRestoreArrayF90(Y,Y_<u></u>pointer,ierr)<br>
call VecDestroy(Y,ierr)<br>
call MatDestroy(B,ierr)<br>
<br>
<br>
call PetscFinalize(ierr)<br>
end<br>
<br>
<br>
<br>
Am I doing something wrong in my code?<br>
Thank you for your help!<span class="HOEnZb"><font color="#888888"><br>
Jonas<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
</font></span></blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div>