[petsc-users] extract Mat data

[Barry Smith]

On Mar 31, 2011, at 12:38 PM, gouarin wrote:

> On 31/03/2011 19:06, Matthew Knepley wrote:
>> On Thu, Mar 31, 2011 at 10:10 AM, gouarin <loic.gouarin at math.u-psud.fr> wrote:
>> Thanks for the reply.
>> 
>> 
>> On 31/03/2011 16:48, Hong Zhang wrote:
>> Loic,
>> 
>> I want to do the same for MATMPI. This is what I do:
>> 
>> MatGetLocalMat(A,MAT_INITIAL_MATRIX,&Aloc);
>> MatGetRowIJ(Aloc, 1, PETSC_FALSE, PETSC_FALSE,&(Acsr.n),&(Acsr.ia),
>> &(Acsr.ja),&done);
>> MatGetArray(Aloc,&(Acsr.a));
>> 
>> Now, I need to know on what process the non local points are.
>> Each process get its own local submatrix. See
>> http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Mat/MatCreateMPIAIJ.html
>> on data structure of MPIAIJ format. Applying your statements above to
>> the example listed on the page (a 8x8 matrix), Proc0 will get
>> Aloc=A(0:2,:), ..., Porc2 gets Aloc=A(6:7,:).
>> Is there a petsc function that gives this partition ?
>> MatGetOwnershipRange()
>> I use this function to know the local points index. In the same way, with this function I know the non local index.
>> 
>> But this function doesn't give the process. I give you an example. I have this matrix:
>> 
>> 1 0 | 0 1
>> 1 2 | 1 0
>> -----------
>> 0 0 | 1 0
>> 0 2 | 0 1
>> 
>> I want this kind of format
>> 
>> for proc 0
>> a = {1,1,1,2,1}
>> ia = {0,2,5}
>> ja = {0,3,0,1,2}
>> rank = {0,0,1,1}
>> 
>> for proc 1
>> a = {1,2,1}
>> ia = {0,1,3}
>> ja = {0,2,1}
>> rank = {1,1,0}
>> 
>> rank give me the process for each non zero column. Local rows should have numbers 1,...N and offdiagonal entries with j>N.
>> 
>> The ownership range is the row partition. However, MPIAIJ matrices are not stored as you indicate above.
>> 
>> Different storage formats make the process you describe fragile. WHy do you need the CSR form of the matrix? If you really
>> do, use MatGetSubMatrix() to get a  SEQAIJ matrix on each process and pull out the values.
>> 
> I want to use this format because I try to create a new multigrid PC with this method:
> 
> http://homepages.ulb.ac.be/~ynotay/AGMG/userguide/index.html
> 
> I don't understand why it's not possible to have the process partition. When you create a MPIAIJ matrix and you use a KSP to solve the system, you know what values to exchange. No? Where is this information ?

   The point is not that the information is not available, it is just that it is not trivially available since we don't use the data in that format.

   Here is what you have to do. Include "src/mat/impls/aij/mpi/mpiaij.h" directly in your interface files and access the data directly, do not try to find some PETSc function that provides exactly what you want directly, that cannot exist because we cannot know what format each other package wants. In particular inside Mat_MPIAJ are the two submatrices called A and B which are Mat_SeqAIJ matrices. The data in
#if defined (PETSC_USE_CTABLE)
  PetscTable    colmap;
#else
  PetscInt      *colmap;                /* local col number of off-diag col */
#endif
  PetscInt      *garray;                /* global index of all off-processor columns */
provides the information about the "off process" columns for the matrices stored in B. 

   But you will have to understand the data structures somewhat in order to get the data in the form you want.

  Barry





> 
> I can also use MatGetOwnershipRange() and MPI_Allgather to have this information.
> 
> Loic
> 
> 
> 
> 
>>    Matt
>>  
>> I suppose that I have no information about the structure of MPIAIJ.
>> 
>> Loic
>> 
>> 
>> 
>> 
>> -- 
>> 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
> 
> 
> -- 
> Loic Gouarin
> Laboratoire de Mathématiques 
> Université Paris-Sud  
> Bâtiment 425 
> 91405 Orsay Cedex
> France 
> Tel: (+33) 1 69 15 60 14
> Fax: (+33) 1 69 15 67 18
>