Great! BTW: since most compilers (including gfortran) now support F90 - its best to use VecGetArrayF90 - instead of VecGetArray() [eventhough the rest of the code is F77]
Satish On Sat, 22 Nov 2014, Valerio Grazioso wrote: > Ehi Barry, > thanks a lot… that was the problem: I had the -fbounds-check flag of gfortran > in my Makefile…removing it everything is ok ! > (Of course …it was in the manual… !! :) ) > > Valerio > > > > > On 21/nov/2014, at 14:18, Barry Smith <[email protected]> wrote: > > > > >> On Nov 20, 2014, at 9:19 PM, Valerio Grazioso <[email protected]> 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. > > > > It is completely possible that i_x is a very large negative value, in > > fact it usually is. > > > > From the users manual: Note: If using VecGetArray(), MatDenseGetArray(), > > or ISGetIndices(), from Fortran, the user must not compile the Fortran code > > with options to check for “array entries out of bounds” (e.g., on the IBM > > RS/6000 this is done with the -C compiler option, so never use the -C > > option with this). > > > > My guess is that your Fortran compiler is checking for array out of bounds > > and (of course) finding them since xx_v() is declared as xx_v(1). You need > > to turn of this array bounds checking when compiling. Or you can use > > VecGetArrayF90(). > > > > Barry > > > > > > > >> 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? > >> 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 > >> ------------------------------ > >> > >> > > > >
