CSR Sparse Matrix Query

[Randall Mackie]

Hi Tim,

The way I understand it, Distributed Arrays are used for structured problems
and they automatically take care of ghost values and so forth, and for me
were much more efficient and easier to use (once I got them figured out, but
then I'm not an expert). There is a good section on Distributed Arrays in the
users manual.

So, now in my more recent code, I replaced what I outlined in the first email with calls
like:

       call DACreate3d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,
      .     l+1,m+1,nzn,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,i3,i1,
      .     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,
      .     da,ierr)

       call DACreateGlobalVector(da,b,ierr)
       call VecDuplicate(b,xsol,ierr)

       ngrow=3*(l+1)*(m+1)*(n+1)

       call MatCreateMPIAIJ(PETSC_COMM_WORLD,mloc,mloc,ngrow,ngrow,
      .     i15,PETSC_NULL_INTEGER,i7,PETSC_NULL_INTEGER,A,ierr)

       call DAGetInfo(da,PETSC_NULL_INTEGER,mx,my,mz,
      .               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,
      .               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,
      .               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,
      .               PETSC_NULL_INTEGER,ierr)

       call DAGetCorners(da,xs,ys,zs,xm,ym,zm,ierr)
       call DAGetGhostCorners(da,gxs,gys,gzs,gxm,gym,gzm,ierr)

[and then to set the matrix elements, calls like]:

       do kk=zs,zs+zm-1
         do jj=ys,ys+ym-1
           do ii=xs,xs+xm-1

             row=ii-gxs + (jj-gys)*gxm + (kk-gzs)*gxm*gym

             i=ii+1
             j=jj+1
             k=kk+1

!--Begin-Hx--====================================================================

             grow=ltog(idltog + 3*row + 1)

             IF (jj == 0 .or. jj == my-1 .or. kk == 0 .or.
      .          ii == mx-1) THEN

               v(1) = 1.d0
               col(1) = grow
               call MatSetValues(A,i1,grow,i1,col,v,INSERT_VALUES,
      .             ierr)
             ELSE

               fact=-rijcy*dx*dy/z(k-1)         ! Hxijc
               ic=ic+1
               v(ic)=fact
               col(ic)=ltog(idltog + 3*(row-gxm*gym) +1)


[etc]


Regards,

Randy


Tim Stitt wrote:
> Randall,
> 
> Many thanks for that...very helpful indeed.
> 
> Just out of curiosity...what do you mean by "Distributed Arrays". This 
> is not the same as PETSc distributed arrays which I thought are implicit 
> in the Matrix and Vector generation routines (depending on whether 
> multiple nodes are available or not)?
> 
> Regards,
> 
> Tim.
> 
> Randall Mackie wrote:
>> Hi Tim,
>>
>> There's no need to generate a dense matrix if you're dealing with sparse
>> matrices. There are many approaches, but if you just want to get your
>> existing matrix into PETSc, you can use something like the following 
>> (written
>> in Fortran, but almost the same in C):
>>
>> Assume the global vector length is np:
>>
>>       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
>>
>>
>>       ALLOCATE (d_nnz(mloc), o_nnz(mloc), STAT=istat)
>>
>> Note: the d_nnz and o_nnz are used to determine the pre-allocation of 
>> the global matrix
>> according to the PETSc manual (basically just looping through and 
>> figuring out which
>> elements are in d_nnz and which are in o_nnz).
>>
>> Then you can create the matrix:
>>
>>
>>       call MatCreateMPIAIJ(PETSC_COMM_WORLD,mloc,mloc,np,np,
>>      .     PETSC_NULL_INTEGER, d_nnz, PETSC_NULL_INTEGER,
>>      .     o_nnz,P,ierr)
>>
>>
>> And then I had a subroutine to fill it like this (the critical part is 
>> the check
>> whether that part of the matrix is owned by the current processor):
>>
>>       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
>>
>>               if(k.ne.2) then
>>                 fact=-rijcy*dx*dy/z(k-1)       ! Hxijc
>>                 ic=ic+1
>>                 v(ic)=fact
>>                 col(ic)=ijkhx(i,j,k-1)-1
>>               end if
>>
>> [snip code]
>>
>>               call MatSetValues(A,i1,row,ic,col,v,INSERT_VALUES,
>>      .             ierr)
>>
>>
>> [and finally at the end you assemble the matrix]
>>
>>       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
>>       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
>>
>>
>> Since these early days, I've switched over to using Distributed Arrays,
>> which in my opinion are much better and much easier to use.
>>
>> Hope this helps,
>>
>> Randy
>>
>>
>> Dr. Timothy Stitt wrote:
>>> Hi all,
>>>
>>> I have just begun to use PETSc and I have already ran into a problem. 
>>> Apologies if this has been dealt with before.
>>>
>>> My code currently generates its own non-distributed Hamiltonian 
>>> sparse matrix in CSR format (due to the large size of the dense 
>>> representation). I want to exploit PETSc's (actually SLEPc's) 
>>> eigensolver routines and hence I need to convert my CSR 
>>> representation into PETSc CSR representation for use by the eigensolver.
>>>  From reading the documentation I am under the impression that you 
>>> initially supply the dense matrix to the PETSc matrix routines for 
>>> conversion to PETSc CSR representation. In my case I do not generate 
>>> the dense matrix. I would be grateful if someone could guide me on to 
>>> how to create a PETSC sparse matrix from my own in-code generated 
>>> sparse matrix.
>>>
>>> Thanks in advance for any assistance that can be provided.
>>>
>>> Regards,
>>>
>>> Tim.
>>>
>>
> 
> 

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