<html>
<head>
<meta content="text/html; charset=UTF-8" http-equiv="Content-Type">
</head>
<body bgcolor="#FFFFFF" text="#000000">
Thank you! <br>
And of course you are right, the more i work with PETSc the more I
think about changing from Fortran to C...<br>
<br>
<br>
<br>
<br>
<div class="moz-cite-prefix">Am 11.06.2014 14:43, schrieb Matthew
Knepley:<br>
</div>
<blockquote
cite="mid:CAMYG4GnG1A+YMeZBRBWgDZJ_gN5-k4yhp+0EkcvcnL=N3b8suw@mail.gmail.com"
type="cite">
<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
moz-do-not-send="true"
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_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,ierr)<br>
call VecSet(X,one,ierr) !set all entries of
global array to 1<br>
call VecGetArrayF90(X,X_pointer,ierr)<br>
<br>
call MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,m,PETSC_NULL_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+xm-1),1,0,X_pointer(1:xm),INSERT_VALUES,ierr)<br>
call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)<br>
call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)<br>
<br>
<br>
!Set up vector Y and matrix B<br>
call DMCreateGlobalVector(da,Y,ierr)<br>
call VecSet(Y,2.d0*one,ierr) !set all entries
of global array to 2<br>
call VecGetArrayF90(Y,Y_pointer,ierr)<br>
<br>
call MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,m,PETSC_NULL_CHARACTER,B,ierr)
<br>
call MatSetValues(B,xm,row(xs:xs+xm-1),1,0,Y_pointer(1:xm),INSERT_VALUES,ierr)<br>
call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr)<br>
call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr)<br>
<br>
<br>
!print out matrices and vectors to check entries<br>
call VecView(X,PETSC_VIEWER_STDOUT_WORLD,ierr)<br>
call VecView(Y,PETSC_VIEWER_STDOUT_WORLD,ierr)<br>
call MatView(A,PETSC_VIEWER_STDOUT_WORLD,ierr)<br>
call MatView(B,PETSC_VIEWER_STDOUT_WORLD,ierr)<br>
<br>
!free objects<br>
call VecRestoreArrayF90(X,X_pointer,ierr)<br>
call VecDestroy(X,ierr)<br>
call MatDestroy(A,ierr)<br>
call VecRestoreArrayF90(Y,Y_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>
</blockquote>
<br>
</body>
</html>