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 ------------------------------ >> >> >
