[petsc-users] error in calling VecGetArrayf90()

Ethan Coon ecoon at lanl.gov
Wed Jan 5 13:52:20 CST 2011


See the below program and output (done with dev, but I don't think this
has changed).  In F90, the local arrays, gotten by VecGetArrayF90, are
all indexed:

xx_v(1:(Iend-Istart)).  

This should really be moot -- It's way easier to follow the example code
for ex50f90.F and pass their arrays into a subroutine which casts the
shape in the "PETSc way", indexed from (Istart:Iend-1).  

Ethan


---


#define PETSC_USE_FORTRAN_MODULES 1
#include "finclude/petscsysdef.h"
#include "finclude/petscvecdef.h"

program test
  use petsc
  implicit none

  Vec x
  PetscInt Istart, Iend
  PetscScalar, pointer:: xx_v(:)
  PetscInt i
  PetscErrorCode ierr
  PetscInt twenty
  PetscInt rank
  
  twenty = 20

  call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
  call mpi_comm_rank(PETSC_COMM_WORLD,rank,ierr)
  call VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, twenty, x, ierr)

  call VecGetArrayF90(x, xx_v, ierr)
  xx_v = rank+1
  call VecRestoreArrayF90(x, xx_v, ierr)

  call VecView(x, PETSC_VIEWER_STDOUT_WORLD, ierr)

  call VecGetOwnershipRange(x, Istart, Iend)
  write (*,*) 'range on rank', rank, ':', Istart, Iend

  call VecGetArrayF90(x, xx_v, ierr)
  CHKMEMQ
  do i = 1,(Iend-Istart)
     write (*,*) i, xx_v(i), rank
  enddo
  CHKMEMQ
  call VecRestoreArrayF90(x, xx_v, ierr)

  call VecDestroy(x, ierr)
  call PetscFinalize(ierr)
end program test

-----

$> ${PETSC_DIR}/${PETSC_ARCH}/bin/mpiexec -n 2 ./test
Vector Object:
  type: mpi
Process [0]
1
1
1
1
1
1
1
1
1
1
Process [1]
2
2
2
2
2
2
2
2
2
2
 range on rank           1 :          10          20
           1   2.0000000000000000                1
           2   2.0000000000000000                1
           3   2.0000000000000000                1
           4   2.0000000000000000                1
           5   2.0000000000000000                1
           6   2.0000000000000000                1
           7   2.0000000000000000                1
           8   2.0000000000000000                1
           9   2.0000000000000000                1
          10   2.0000000000000000                1
 range on rank           0 :           0          10
           1  1.00000000000000000                0
           2  1.00000000000000000                0
           3  1.00000000000000000                0
           4  1.00000000000000000                0
           5  1.00000000000000000                0
           6  1.00000000000000000                0
           7  1.00000000000000000                0
           8  1.00000000000000000                0
           9  1.00000000000000000                0
          10  1.00000000000000000                0

---


On Wed, 2011-01-05 at 13:21 -0600, Barry Smith wrote:
> I don't think this is correct. You are suppose to use the local indexing for each process. With the strange index starting at 1 instead of 0.
> 
> 
>    Barry
> 
> 
> On Jan 5, 2011, at 1:03 PM, Ethan Coon wrote:
> 
> > On all processors, the array you get is indexed:
> > 
> > xx_v(1:(Iend-Istart))
> > 
> > while the Istart and Iend are global indices into the global Vec.
> > 
> > It's only by luck that the values on proc 3 are correct (this code
> > should probably seg fault as it is accessing memory outside of the
> > bounds of xx_v).
> > 
> > To access the array, you'll want:
> > 
> > do i=1, (Iend-Istart)
> >   write(6,*)'check xx_v',i,xx_v(i),myid
> > enddo
> > 
> > or (better yet) pass it into a subroutine to get the array indexed
> > correctly, like demonstrated in the example.
> > 
> > Ethan
> > 
> > On Wed, 2011-01-05 at 11:23 -0600, Peter Wang wrote:
> >> Thanks, Satish,
> >> 
> >>     The index of the array is modified to i+1:
> >> !===================
> >>         do i=Istart,Iend-1
> >>              write(6,*)'check xx_v',i+1,xx_v(i+1),myid
> >>         enddo  
> >> !===================
> >> 
> >>    However, only the elements on root process (process 0) and the
> >> last process (process 3) are corrent, is there any ohter logical
> >> error?
> >> 
> >> check xx_v           1   3999.9999999999982                0
> >> check xx_v           2   3999.9999999999982                0
> >> check xx_v           3   3999.9999999999982                0
> >> check xx_v           4   3999.9999999999982                0
> >> check xx_v           5   3999.9999999999982                0
> >> check xx_v           6   3000.0000000000005                0
> >> check xx_v           7   3000.0000000000005                0
> >> check xx_v           8  2.61360726650019422E-321           1 
> >> check xx_v           9  7.90505033345994471E-323           1
> >> check xx_v          10  1.69759663277221785E-312           1
> >> check xx_v          11  6.16840020108069212E-317           1
> >> check xx_v          12  6.16840316547456717E-317           1
> >> check xx_v          13  6.16832658529946177E-317           1
> >> check xx_v          14  1.99665037664579820E-314           2
> >> check xx_v          15  6.19784009071943448E-317           2
> >> check xx_v          16  6.20249221284067566E-317           2
> >> check xx_v          17  6.20218737433719161E-317           2
> >> check xx_v          18  6.18236051996958238E-317           2
> >> check xx_v          19  6.16840316547456717E-317           2
> >> check xx_v          20  6.18107199676522841E-317           3
> >> check xx_v          21   0.0000000000000000                3
> >> check xx_v          22   0.0000000000000000                3
> >> check xx_v          23   0.0000000000000000                3
> >> check xx_v          24   0.0000000000000000                3
> >> check xx_v          25   0.0000000000000000                3
> >> 
> >> 
> >> 
> >>> Date: Tue, 4 Jan 2011 22:49:50 -0600
> >>> From: balay at mcs.anl.gov
> >>> To: petsc-users at mcs.anl.gov
> >>> Subject: Re: [petsc-users] error in calling VecGetArrayf90()
> >>> 
> >>> The global index starts at Istart - but the array index starts at 1
> >> [for fortran array]
> >>> 
> >>> Satish
> >>> 
> >>> On Tue, 4 Jan 2011, Peter Wang wrote:
> >>> 
> >>>> 
> >>>> In last question, the pointer xx_v is local data. However, if
> >> write them to the monitor or assign them to another array, the value
> >> is incorrect.
> >>>> 
> >>>> The protion of the code to display them on the monitor is like as
> >> following: 
> >>>> call MatGetOwnershipRange(A,Istart,Iend,ierr)
> >>>> call VecGetArrayF90(x,xx_v,ierr) ! Vector x is matched with Matrix
> >> A in the same communicator
> >>>> 
> >>>> write(*,*)xx_v,myid ! write the poiner array together
> >>>> 
> >>>> do i=Istart,Iend-1
> >>>> write(6,*)'check xx_v',i,xx_v(i),myid !write the element of the
> >> array one by one with local range (Istart to Iend-1)
> >>>> enddo 
> >>>> 
> >>>> 
> >>>> =========The result is as following: ( the values of the elements
> >> from 7 to 20 are not correct !!)
> >>>> 
> >>>> 3999.9999999999982 3999.9999999999982 3999.9999999999982
> >> 3999.9999999999982 3999.9999999999982 3000.0000000000005
> >> 3000.0000000000005 0
> >>>> 
> >>>> 3000.0000000000009 3000.0000000000009 3000.0000000000009
> >> 2000.0000000000011 2000.0000000000011 2000.0000000000000 1
> >>>> 
> >>>> 2000.0000000000009 2000.0000000000009 1000.0000000000003
> >> 1000.0000000000003 1000.0000000000003 999.99999999999989 2
> >>>> 
> >>>> 1000.0000000000000 0.0000000000000000 0.0000000000000000
> >> 0.0000000000000000 0.0000000000000000 0.0000000000000000 3
> >>>> 
> >>>> 
> >>>> check xx_v 0 0.0000000000000000 0.0000000000000000 0
> >>>> check xx_v 1 3999.9999999999982 3999.9999999999982 0
> >>>> check xx_v 2 3999.9999999999982 3999.9999999999982 0
> >>>> check xx_v 3 3999.9999999999982 3999.9999999999982 0
> >>>> check xx_v 4 3999.9999999999982 3999.9999999999982 0
> >>>> check xx_v 5 3999.9999999999982 3999.9999999999982 0
> >>>> check xx_v 6 3000.0000000000005 3000.0000000000005 0
> >>>> check xx_v 7 1.99665037664579820E-314 1.99665037664579820E-314 1
> >>>> check xx_v 8 2.61360726650019422E-321 2.61360726650019422E-321 1
> >>>> check xx_v 9 7.90505033345994471E-323 7.90505033345994471E-323 1
> >>>> check xx_v 10 1.69759663277221785E-312 1.69759663277221785E-312 1
> >>>> check xx_v 11 6.16846344148335980E-317 6.16846344148335980E-317 1
> >>>> check xx_v 12 6.16846640587723485E-317 6.16846640587723485E-317 1
> >>>> check xx_v 13 6.16838982570212945E-317 6.16838982570212945E-317 2
> >>>> check xx_v 14 1.99665037664579820E-314 1.99665037664579820E-314 2
> >>>> check xx_v 15 6.19790333112210216E-317 6.19790333112210216E-317 2
> >>>> check xx_v 16 6.20255545324334334E-317 6.20255545324334334E-317 2
> >>>> check xx_v 17 6.20225061473985929E-317 6.20225061473985929E-317 2
> >>>> check xx_v 18 6.18242376037225006E-317 6.18242376037225006E-317 2
> >>>> check xx_v 19 6.16846640587723485E-317 6.16846640587723485E-317 3
> >>>> check xx_v 20 6.18113523716789609E-317 6.18113523716789609E-317 3
> >>>> check xx_v 21 0.0000000000000000 0.0000000000000000 3
> >>>> check xx_v 22 0.0000000000000000 0.0000000000000000 3
> >>>> check xx_v 23 0.0000000000000000 0.0000000000000000 3
> >>>> check xx_v 24 0.0000000000000000 0.0000000000000000 3
> >>>> 
> >>>> ======The vector x is :
> >>>> Process [0]
> >>>> 4000
> >>>> 4000
> >>>> 4000
> >>>> 4000
> >>>> 4000
> >>>> 3000
> >>>> 3000
> >>>> Process [1]
> >>>> 3000
> >>>> 3000
> >>>> 3000
> >>>> 2000
> >>>> 2000
> >>>> 2000
> >>>> Process [2]
> >>>> 2000
> >>>> 2000
> >>>> 1000
> >>>> 1000
> >>>> 1000
> >>>> 1000
> >>>> Process [3]
> >>>> 1000
> >>>> 0
> >>>> 0
> >>>> 0
> >>>> 0
> >>>> 0
> >>>> 
> >>>> 
> >>>> 
> >>>> 
> >>>>> Date: Tue, 4 Jan 2011 17:50:11 -0600
> >>>>> From: balay at mcs.anl.gov
> >>>>> To: petsc-users at mcs.anl.gov
> >>>>> Subject: Re: [petsc-users] error in calling VecGetArrayf90()
> >>>>> 
> >>>>> Did you included "finclude/petscvec.h90" in your code - as the
> >> example did?
> >>>>> 
> >>>>> satish
> >>>>> 
> >>>>> On Tue, 4 Jan 2011, Peter Wang wrote:
> >>>>> 
> >>>>>> 
> >>>>>> I am trying to obtain the value of each element of a solution
> >> Vector by KSPsolve(). 
> >>>>>> 
> >>>>>> The variables are defined according the example of ex4f90.F in
> >> \petsc-3.1-p5\src\snes\examples\tutorials\ as following,
> >>>>>> 
> >>>>>> PetscScalar, pointer :: xx_v(:)
> >>>>>> 
> >>>>>> ...
> >>>>>> call KSPSolve(ksp,b,x,ierr)
> >>>>>> call VecView(x,PETSC_VIEWER_STDOUT_WORLD,ierr)
> >>>>>> 
> >>>>>> call VecGetArrayF90(x,xx_v,ierr)
> >>>>>> call VecRestoreArrayF90(x,xx_v,ierr)
> >>>>>> 
> >>>>>> ...
> >>>>>> 
> >>>>>> But, the error keeps coming out when call
> >> VecGetArrayF90(x,xx_v,ierr) and call VecRestoreArrayF90(x,xx_v,ierr)
> >> are not commented off.
> >>>>>> 
> >>>>>> 
> >>>>>> The error information shows:
> >>>>>> Caught signal number 11 SEGV: Segmentation Violation, probably
> >> memory access out of range
> >>>>>> 
> >>>>>> [0]PETSC ERROR: --------------------- Stack Frames
> >> ------------------------------------
> >>>>>> [0]PETSC ERROR: Note: The EXACT line numbers in the stack are
> >> not available,
> >>>>>> [0]PETSC ERROR: INSTEAD the line number of the start of the
> >> function
> >>>>>> [0]PETSC ERROR: is given.
> >>>>>> [0]PETSC ERROR: [0] F90Array1dCreate line 52
> >> src/sys/f90-src/f90_cwrap.c
> >>>>>> [0]PETSC ERROR: --------------------- Error Message
> >> ------------------------------------
> >>>>>> 
> >>>>>> I checked the code according the example, but cannot see any
> >> difference to that. Just don't know why the pointer array xx_v doesn't
> >> work here? Thanks.
> >>>>>> 
> >>>>>> 
> >>>>> 
> >>>> 
> >>> 
> > 
> > -- 
> > -------------------------------------
> > Ethan Coon
> > Post-Doctoral Researcher
> > Mathematical Modeling and Analysis
> > Los Alamos National Laboratory
> > 505-665-8289
> > 
> > http://www.ldeo.columbia.edu/~ecoon/
> > -------------------------------------
> > 
> 

-- 
-------------------------------------
Ethan Coon
Post-Doctoral Researcher
Mathematical Modeling and Analysis
Los Alamos National Laboratory
505-665-8289

http://www.ldeo.columbia.edu/~ecoon/
-------------------------------------



More information about the petsc-users mailing list