Matrix input file format

[Matthew Knepley]

The Istart and Iend are the rows present on this particular processor,
which can be obtain from MatGetOwnershipRange().

  Matt

On 6/15/06, Evrim Dizemen <gudik at ae.metu.edu.tr> wrote:
>
> Dear Randall,
>
> I guess i began to understand the consepts but there is still a missing
> point that i do not know when and how we define Istart and Iend. I'll be
> glade if you can send me the pre_op3d routine so i can see the algorithm
> which is a black box for me now.
>
> Thanks a lot
>
> EVRIM
>
>
> Randall Mackie wrote:
> > Hi Evrim,
> >
> > It's quite easy to modify your Fortran code to do what you want. I
> > thought
> > I had written it all out before, but I'll try again. There are many ways
> > to do this, but I'll start with the easiest, at least if you're going to
> > just modify your current sequential code.
> >
> > Let's say that your matrix has np global rows. Then
> >
> >       call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,np,b,ierr)
> >       call VecDuplicate(b,xsol,ierr)
> >       call VecGetLocalSize(b,mloc,ierr)
> >       call VecGetOwnershipRange(b,Istart,Iend,ierr)
> >
> >       do i=Istart+1,Iend
> >         loc(i)=i-1
> >       end do
> >
> > These statements create parallel vectors for the solution (xsol) and
> > the right hand side (b). The vector loc(i)is used to set values in the
> > vectors later.
> >
> > Then
> >
> > ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> -
> > !         Create the linear solver context
> > ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> -
> >
> >       call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
> >
> > ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> -
> > !         Create the scatter context for getting results back to
> > !         each node.
> > ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> -
> >
> >       call VecScatterCreateToAll(xsol,xToLocalAll,xseq,ierr)
> >
> > ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> -
> > !     Create the matrix that defines the linear system, Ax = b,
> > !     for the EM problem.
> > ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> -
> >
> >       call pre_op3d(l,m,nzn,Istart,Iend,ijkhx,ijkhy,ijkhz,d_nnz,o_nnz)
> >
> >       call MatCreateMPIAIJ(PETSC_COMM_WORLD,mloc,mloc,np,np,
> >      .     PETSC_NULL_INTEGER, d_nnz, PETSC_NULL_INTEGER,
> >      .     o_nnz,A,ierr)
> >
> >       call set_op3d(A,l,m,nzn,period,resist,x,y,z,Istart,Iend,
> >      .           ijkhx,ijkhy,ijkhz)
> >
> >
> > The subroutine pre_op3d is an important routine and it figures out
> > how to pre-allocate space for your parallel matrix. This will be
> > the difference between near-instanteous assembly and 2.5 hours that
> > you experienced. Basically, it just computes the global column numbers,
> > and figures out if they are between Istart and Iend. I can send you
> > my subroutine if you'd like.
> >
> > The subroutine set_op3d.F actually assembles the parallel matrix and
> > goes like this:
> >
> >
> >       jj=0
> >
> >       do i=1,l
> >         do k=2,n
> >           do j=2,m
> >
> >             jj=jj+1
> >             row = jj-1
> >
> >             IF (jj >= Istart+1 .and. jj <= Iend) THEN
> >
> >            compute elements...
> >
> >               call MatSetValues(A,i1,row,ic,col,v,INSERT_VALUES,
> >      .             ierr)
> >
> >             END IF
> >
> >           end do
> >         end do
> >       end do
> >
> > At the end,
> >
> >       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
> >       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
> >
> >
> > Again, because you compute the pre-allocation, this is near-instanteous,
> > even for large models (larger than you're using).
> >
> >
> > Once you do that, you're golden:
> >
> >         call KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN,ierr)
> >
> > etc.
> >
> > Randy M.
> > San Francisco
> >
> >
> > Evrim Dizemen wrote:
> >> Hi all,
> >>
> >> Again thanks for your comments. I guess i can not define the problem
> >> correctly. I have a sequantial fortran code giving me the global
> >> matrix of 200000x200000. The code writes the matrix in a binary file
> >> at little endian mode (i can write only the nonzero terms or the
> >> entire matrix). I tried to change the binary mode to big endian and
> >> read the global matrix by a c program as in the example
> >> /src/mat/example/tests/ex31.c. However the program reads binary file
> >> wrong and gives the following error message but the true value of
> >> no-nonzero in the binary file is 6 (for the test case 3x3 matrix) :
> >>
> >> reading matrix in binary from matrix.dat ...
> >>
> ------------------------------------------------------------------------
> >> Petsc Release Version 2.3.1, Patch 13, Wed May 10 11:08:35 CDT 2006
> >> BK revision: balay at asterix.mcs.anl.gov|ChangeSet|20060510160640|13832
> >> See docs/changes/index.html for recent updates.
> >> See docs/faq.html for hints about trouble shooting.
> >> See docs/index.html for manual pages.
> >>
> ------------------------------------------------------------------------
> >> ./ex31 on a linux named akbaba.ae.metu.edu.tr by evrim Thu Jun 15
> >> 09:26:36 2006
> >> Libraries linked from /home/evrim/petsc-2.3.1-p13/lib/linux
> >> Configure run at Tue May 30 10:26:48 2006
> >> Configure options --with-scalar-type=complex --with-shared=0
> >>
> ------------------------------------------------------------------------
> >> [0]PETSC ERROR: MatLoad_SeqAIJ() line 3055 in
> >> src/mat/impls/aij/seq/aij.c
> >> [0]PETSC ERROR: Read from file failed!
> >> [0]PETSC ERROR: Inconsistant matrix data in file. no-nonzeros =
> >> 100663296, sum-row-lengths = 234300
> >> !
> >> [0]PETSC ERROR: MatLoad() line 149 in src/mat/utils/matio.c
> >> [0]PETSC ERROR: main() line 37 in src/mat/examples/tests/ex31.c
> >>
> >> I want to send the global matrix at once to Petsc by a written input
> >> file (as i'm working on now) or by sending the matrix array from my
> >> fortran code and then partition it and solve iteratively. After the
> >> solution i also want to get the solution vector back to my fortran
> >> code. As i told in the previous mails i tried to send the matrix
> >> array to Petsc and used MatSetValues as reading one value at time in
> >> a do loop but it took about 2,5 hours to read the global matrix.
> >> Additionally i tried to read a row at  a time but can not figure out
> >> a algorithm for this. Hence i do not prefer to create the matrix
> >> again in Petsc by MatSetValues.
> >>
> >> Aside i figured out that the binary files written in fortran and c
> >> are completely different from each other (fortran adds the size of
> >> the characters to the beginning and end of each character) so i wrote
> >> a c interface code to get the matrix array from the fortran code and
> >> write it to a binary file in c format. By this code i avoided from
> >> the additional information in the binary file but i still have the
> >> endianness problem.
> >>
> >> I know that i asked so much but since i'm a rookie in parallel
> >> programming, c language and library using i really need your comments
> >> on my problem. Sorry for this long mail and thanks a lot for your
> >> kind effort on guiding me.
> >>
> >> Thanks
> >>
> >> EVRIM
> >>
> >
>
>


-- 
"Failure has a thousand explanations. Success doesn't need one" -- Sir Alec
Guiness
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20060615/ea53e332/attachment.htm>