Matrix input file format

[Randall Mackie]

Evrim,

I was in your spot several years ago, wanting to convert a serial
code to solve complex systems to parallel, and PETSc seemed like a
good solution. It was a steep learning curve for me, but well worth
it. Barry told you how to read a file with PETSc, but maybe you want
to consider to generate your matrix and assemble it in parallel?
This has many advantages, one being that no single node
has to store the entire matrix in memory, or that no single node
has to do all the calculations to generate the matrix.

There are different ways to do this in PETSc, but a simple way that
doesn't take too much effort would be as follows:

Say your system is np x np in size, first create a parallel vector
for this system for the right hand side (b) and the solution xsol:

       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

Istart and Iend determine the rows in the global system that are owned by
a particular node. The variable loc(i) just converts to 0-based indexing.
I use it for setting the vector values:

           call VecSetValues(xsol,mloc,loc(Istart+1),xvec(Istart+1),
      .         INSERT_VALUES,ierr)
           call VecAssemblyBegin(xsol,ierr)
           call VecAssemblyEnd(xsol,ierr)


Then, you create the parallel matrix (but you need to know something
about it's structure in order to do an efficient matrix assembly,
and for that you'll need to read the manual):

       call MatCreateMPIAIJ(PETSC_COMM_WORLD,mloc,mloc,np,np,
      .     PETSC_NULL_INTEGER, d_nnz, PETSC_NULL_INTEGER,
      .     o_nnz,A,ierr)


Then, in your program that computes the elements of the coefficient matrix,
each node only computes those rows that it owns (from Istart and Iend), something
like:

	do jj=1,np

             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

After that, you assemble the matrix:

       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

You can do something similar to set up the right hand side (b), and
then it's a rather simple matter to have PETSc solve the system.


 From this, you should be able to fill in the details by reading the
manual and looking at the examples.

Good luck,

Randy Mackie


Evrim Dizemen wrote:
> Dear all,
> 
> I'm using Petsc 2.3.1 for my master thesis. I want to solve a linear 
> problem Ax=b in parallel programming using MPICH with complex variables. 
> I can get the A matrix from my fortran code but i can not find in which 
> format i must write it for reading by Petsc. I'm really not familiar 
> with c programming  so  I'll be happy for a solution in fortran.
> 
> THANKS,
> 
> Evrim
> 

-- 
Randall Mackie
GSY-USA, Inc.
PMB# 643
2261 Market St.,
San Francisco, CA 94114-1600
Tel (415) 469-8649
Fax (415) 469-5044

California Registered Geophysicist
License No. GP 1034