Hi Matt, thanks for the reply. I've solved the problem thanks to Barry Smith’s suggestion… Cheers Valerio
On 21/nov/2014, at 12:18, Matthew Knepley <[email protected]> wrote: > On Thu, 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. > 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
