<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Thu, Nov 20, 2014 at 9:19 PM, Valerio Grazioso <span dir="ltr"><<a href="mailto:graziosov@libero.it" target="_blank">graziosov@libero.it</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">Dear Petsc Users,<br>
I’m trying to write a simple parallel linear solver (A*x=q) in fortran 77 with Petsc.<br>
I managed to Setup a MPIAIJ matrix A and a VecMPI qvec<br>
then the KSP and the solution.<br>
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.<br>
<br>
To be more than sure I’ve tried to repeat the necessary steps in order to achieve this result on qvec :<br>
<br>
After creating and calling a VecScatter in order to bring a sequential version on all the processors (as described in<br>
<a href="http://www.mcs.anl.gov/petsc/petsc-as/documentation/faq.html#mpi-vec-to-mpi-vec" target="_blank">http://www.mcs.anl.gov/petsc/petsc-as/documentation/faq.html#mpi-vec-to-mpi-vec</a>)<br>
I’ve tried to use VecGetValue as described in the User manual in the Fortran users chapter (pag. 147-148)<br>
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<br>
that, when it is used, generates an obvious “Aut of bound” error.<br>
I’ve tried to use VecGetValue also on qvec directly : same result!<br>
<br>
There must be something wrong that I’m doing.<br>
Can anybody point me in the right direction?<br></blockquote><div><br></div><div>The code looks correct to me. As a first step, can you try running an example, to make sure its</div><div>not a compiler problem. This example</div><div><br></div><div> <a href="http://www.mcs.anl.gov/petsc/petsc-current/src/vec/vec/examples/tutorials/ex4f.F.html">http://www.mcs.anl.gov/petsc/petsc-current/src/vec/vec/examples/tutorials/ex4f.F.html</a></div><div><br></div><div>should be enough to check.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
Thanks<br>
<br>
Regards<br>
Valerio Grazioso<br>
<br>
<br>
Here it is the relevant code:<br>
<br>
<br>
<br>
#define xx_a(ib) xx_v(i_x + (ib))<br>
<br>
…………………..<br>
<br>
Vec xSEQ<br>
PetscScalar xx_v(1)<br>
PetscOffset i_x<br>
<br>
……………………..<br>
………………..<br>
<br>
<br>
write(6,'(a,i2,a)'),"Processor ",mpi_rank,": petsc_solve ---> VecCreate"<br>
call VecCreate(comm,qvec,ierr); CHKERRQ(ierr)<br>
<br>
write(6,'(a,i2,a)'),"Processor ",mpi_rank,": petsc_solve ---> VecSetSizes"<br>
call VecSetSizes(qvec,PETSC_DECIDE,NCELL,ierr); CHKERRQ(ierr)<br>
<br>
write(6,'(a,i2,a)'),"Processor ",mpi_rank,": petsc_solve ---> VecSetFromOptions"<br>
call VecSetFromOptions(qvec,ierr); CHKERRQ(ierr)<br>
<br>
write(6,'(a,i2,a)'),"Processor ",mpi_rank,": petsc_solve ---> VecDuplicate"<br>
call VecDuplicate(qvec,xvec,ierr); CHKERRQ(ierr)<br>
<br>
write(6,'(a,i2,a)'),"Processor ",mpi_rank,": petsc_solve ---> VecGetOwnershipRange"<br>
call VecGetOwnershipRange(qvec,Istart,Iend,ierr); CHKERRQ(ierr)<br>
<br>
write(6,'(a,i2,a,2i8)'),"Processor ",mpi_rank,": petsc_solve ---> Istart,Iend (Vet)=", Istart,Iend<br>
<br>
do Iproc = 0,(mpi_size-1)<br>
if (Iproc == mpi_rank) then<br>
do I= (Istart+1), (Iend)<br>
C write(6,'(a,i2,a,2i8)'),"Processor ",mpi_rank,": petsc_solve ---> MatSetValues",I,I<br>
call VecSetValues(qvec,1,I-1,SU(I),INSERT_VALUES,ierr); CHKERRQ(ierr)<br>
call VecSetValues(xvec,1,I-1,0.E0,INSERT_VALUES,ierr); CHKERRQ(ierr)<br>
enddo<br>
endif<br>
enddo<br>
<br>
write(6,'(a,i2,a)'),"Processor ",mpi_rank,": petsc_solve ---> VecAssemblyBegin"<br>
call VecAssemblyBegin(qvec,ierr); CHKERRQ(ierr)<br>
call VecAssemblyEnd(qvec,ierr); CHKERRQ(ierr)<br>
call VecAssemblyBegin(xvec,ierr); CHKERRQ(ierr)<br>
call VecAssemblyEnd(xvec,ierr); CHKERRQ(ierr)<br>
write(6,'(a,i2,a)'),"Processor ",mpi_rank,": petsc_solve ---> VecAssemblyEnd"<br>
<br>
C --------------------------------- DEBUG -------------------------------------<br>
if (1 == 1) then<br>
<br>
call VecScatterCreateToAll(qvec,ctx,xSEQ,ierr); CHKERRQ(ierr)<br>
call VecScatterBegin(ctx,qvec,xSEQ,INSERT_VALUES,SCATTER_FORWARD,ierr)<br>
CHKERRQ(ierr)<br>
call VecScatterEnd(ctx,qvec,xSEQ,INSERT_VALUES,SCATTER_FORWARD,ierr)<br>
CHKERRQ(ierr)<br>
call VecScatterDestroy(ctx,ierr); CHKERRQ(ierr)<br>
<br>
C call PetscViewerCreate(comm, view, ierr); CHKERRQ(ierr)<br>
C call PetscViewerSetType(view, PETSCVIEWERASCII, ierr); CHKERRQ(ierr)<br>
C call PetscViewerSetFormat(view, PETSC_VIEWER_ASCII_INDEX, ierr)<br>
C CHKERRQ(ierr)<br>
C call PetscViewerFileSetName(view, 'qvec.m', ierr)<br>
C call VecView(xSEQ, view, ierr); CHKERRQ(ierr)<br>
C call PetscViewerDestroy(view, ierr); CHKERRQ(ierr)<br>
<br>
call VecGetArray(xSEQ,xx_v,i_x,ierr); CHKERRQ(ierr)<br>
<br>
write(6, *) (xx_a(I), I=1,NCELL)<br>
write(6,'(a,i20)')"i_x = ", i_x<br>
<br>
<br>
call VecRestoreArray(xSEQ,xx_v,i_x,ierr); CHKERRQ(ierr)<br>
<br>
call PetscFinalize(ierr); CHKERRQ(ierr)<br>
STOP<br>
endif<br>
C ---------------------------------- END DEBUG ------------------------------<br>
<br>
<br>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</div></div>