[petsc-users] Fortran 77 VecGetArray problems

[Matthew Knepley]

On Thu, Nov 20, 2014 at 9:19 PM, Valerio Grazioso <graziosov at libero.it>
wrote:

> Dear Petsc Users,
> I’m trying to write a simple parallel linear solver (A*x=q) in fortran 77
> with Petsc.
> I managed to Setup a MPIAIJ matrix A  and a VecMPI qvec
> then the KSP and the solution.
> Everything seems fine (I’ve checked both A and q with MatView and VecView)
> but I’m having problems getting the solution in a “normal” local double
> precision fortran vector.
>
> To be more than sure I’ve tried to repeat the necessary steps in order to
> achieve this result on qvec :
>
> After creating and calling a VecScatter in order to bring a sequential
> version on all the processors (as described in
>
> http://www.mcs.anl.gov/petsc/petsc-as/documentation/faq.html#mpi-vec-to-mpi-vec
> )
> I’ve tried to use VecGetValue as described in the User manual in the
> Fortran users chapter (pag. 147-148)
> It seems to fail in both the local xx_v vector (that has wrong values) and
> in calculating the i_x index offset that has a very large negative value
> that, when it is used, generates an obvious “Aut of bound” error.
> I’ve tried to use VecGetValue also on qvec directly : same result!
>
> There must be something wrong that I’m doing.
> Can anybody point me in the right direction?
>

The code looks correct to me. As a first step, can you try running an
example, to make sure its
not a compiler problem. This example


http://www.mcs.anl.gov/petsc/petsc-current/src/vec/vec/examples/tutorials/ex4f.F.html

should be enough to check.

  Thanks,

     Matt


> Thanks
>
> Regards
> Valerio Grazioso
>
>
> Here it is the relevant code:
>
>
>
> #define xx_a(ib)  xx_v(i_x + (ib))
>
>  …………………..
>
>       Vec         xSEQ
>       PetscScalar xx_v(1)
>       PetscOffset i_x
>
>       ……………………..
>      ………………..
>
>
>       write(6,'(a,i2,a)'),"Processor  ",mpi_rank,": petsc_solve --->
> VecCreate"
>       call VecCreate(comm,qvec,ierr); CHKERRQ(ierr)
>
>       write(6,'(a,i2,a)'),"Processor  ",mpi_rank,": petsc_solve --->
> VecSetSizes"
>       call VecSetSizes(qvec,PETSC_DECIDE,NCELL,ierr); CHKERRQ(ierr)
>
>       write(6,'(a,i2,a)'),"Processor  ",mpi_rank,": petsc_solve --->
> VecSetFromOptions"
>       call VecSetFromOptions(qvec,ierr); CHKERRQ(ierr)
>
>       write(6,'(a,i2,a)'),"Processor  ",mpi_rank,": petsc_solve --->
> VecDuplicate"
>       call VecDuplicate(qvec,xvec,ierr); CHKERRQ(ierr)
>
>       write(6,'(a,i2,a)'),"Processor  ",mpi_rank,": petsc_solve --->
> VecGetOwnershipRange"
>       call VecGetOwnershipRange(qvec,Istart,Iend,ierr); CHKERRQ(ierr)
>
>       write(6,'(a,i2,a,2i8)'),"Processor  ",mpi_rank,": petsc_solve --->
> Istart,Iend (Vet)=", Istart,Iend
>
>       do Iproc = 0,(mpi_size-1)
>         if (Iproc == mpi_rank) then
>           do I= (Istart+1), (Iend)
> C            write(6,'(a,i2,a,2i8)'),"Processor  ",mpi_rank,": petsc_solve
> ---> MatSetValues",I,I
>             call VecSetValues(qvec,1,I-1,SU(I),INSERT_VALUES,ierr);
> CHKERRQ(ierr)
>             call VecSetValues(xvec,1,I-1,0.E0,INSERT_VALUES,ierr);
> CHKERRQ(ierr)
>           enddo
>         endif
>       enddo
>
>       write(6,'(a,i2,a)'),"Processor  ",mpi_rank,": petsc_solve --->
> VecAssemblyBegin"
>       call VecAssemblyBegin(qvec,ierr); CHKERRQ(ierr)
>       call VecAssemblyEnd(qvec,ierr); CHKERRQ(ierr)
>       call VecAssemblyBegin(xvec,ierr); CHKERRQ(ierr)
>       call VecAssemblyEnd(xvec,ierr); CHKERRQ(ierr)
>       write(6,'(a,i2,a)'),"Processor  ",mpi_rank,": petsc_solve --->
> VecAssemblyEnd"
>
> C ---------------------------------  DEBUG
> -------------------------------------
>       if (1 == 1) then
>
>          call VecScatterCreateToAll(qvec,ctx,xSEQ,ierr); CHKERRQ(ierr)
>          call
> VecScatterBegin(ctx,qvec,xSEQ,INSERT_VALUES,SCATTER_FORWARD,ierr)
>          CHKERRQ(ierr)
>          call
> VecScatterEnd(ctx,qvec,xSEQ,INSERT_VALUES,SCATTER_FORWARD,ierr)
>          CHKERRQ(ierr)
>          call VecScatterDestroy(ctx,ierr); CHKERRQ(ierr)
>
> C         call PetscViewerCreate(comm, view, ierr); CHKERRQ(ierr)
> C         call PetscViewerSetType(view, PETSCVIEWERASCII, ierr);
> CHKERRQ(ierr)
> C         call PetscViewerSetFormat(view, PETSC_VIEWER_ASCII_INDEX, ierr)
> C         CHKERRQ(ierr)
> C         call PetscViewerFileSetName(view, 'qvec.m', ierr)
> C         call VecView(xSEQ, view, ierr); CHKERRQ(ierr)
> C         call PetscViewerDestroy(view, ierr); CHKERRQ(ierr)
>
>          call VecGetArray(xSEQ,xx_v,i_x,ierr); CHKERRQ(ierr)
>
>             write(6, *) (xx_a(I), I=1,NCELL)
>             write(6,'(a,i20)')"i_x = ", i_x
>
>
>          call VecRestoreArray(xSEQ,xx_v,i_x,ierr); CHKERRQ(ierr)
>
>          call PetscFinalize(ierr); CHKERRQ(ierr)
>          STOP
>       endif
> C ---------------------------------- END DEBUG
> ------------------------------
>
>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20141121/f0df4e0c/attachment.html>