[petsc-users] About MatGetOwnershipRange()

[Joon Hee Choi]

Hi Jed,

Thank you for your reply. I tried to follow your method, but I think my code still has something wrong because "localsize" is the same as global row size("I"). Could you let me know what I am missing? My new code is as follows:

PetscErrorCode SetUpMatrix(Mat *X, vector< std::tr1::tuple< PetscInt, PetscInt, PetscInt > >, PetscInt I, PetscInt J)
{
   PetscInt i, j;
   PetscScalar val;

   MatCreate(PETSC_COMM_WORLD, X);
   MatSetType(*X, MATMPIAIJ);
   MatSetSizes(*X, PETSC_DECIDE, PETSC_DECIDE, I, J);

   PetscInt begin, end, localsize=PETSC_DECIDE;
   PetscSplitOwnership(PETSC_COMM_WORLD, &localsize, &I);
   MPI_Scan(&localsize, &end, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD);
   begin = end - localsize;
   PetscPrintf(PETSC_COMM_WORLD, "Local Size: %d, begin: %d, end: %d \n", localsize, begin, end);

   ...
}

Also, tuples are the total pairs of (row dimension, column dimension, element), values of sparse matrix, which are read from a file. The tuples are sorted and distributed. Thus, I tried to set up d_nnz and o_nnz. I implemented d_nnz[i-begin]++ only if begin<="i"&&"j"<end, otherwise o_nnz[i-begin]++ because "i" means the row number and "j" means the column number. Please let me know if there is something wrong in the following code:

   ...
   for (int x=0; x<tups.size(); x++) {
      i = std::tr1::get<0>(tups[x]);
      if (i>=end) break;
      if (i>=begin) {
         j = std::tr1::get<1>(tups[x]);
         if (j>=begin && j<end)  d_nnz[i-begin]++;
         else                    o_nnz[i-begin]++;
      }
   }
   MatMPIAIJSetPreallocation(*X, 0, d_nnz, 0, o_nnz);
   for (int x=0; x<tups.size(); x++) {
      i = std::tr1::get<0>(tups[x]);
      if (i>=end) break;
      if (i>=begin) {
         j = std::tr1::get<1>(tups[x]);
         val = std::tr1::get<2>(tups[x]);
         MatSetValues(*X, 1, &i, 1, &j, &val, INSERT_VALUES);
      }
   }
   MatAssemblyBegin(*X, MAT_FINAL_ASSEMBLY);
   MatAssemblyEnd(*X, MAT_FINAL_ASSEMBLY);
}


----- Original Message -----
From: "Jed Brown" <jedbrown at mcs.anl.gov>
To: "Joon Hee Choi" <choi240 at purdue.edu>, petsc-users at mcs.anl.gov
Sent: Thursday, September 12, 2013 4:04:17 AM
Subject: Re: [petsc-users] About MatGetOwnershipRange()

Joon Hee Choi <choi240 at purdue.edu> writes:

>> Hello,
>>
>> I am trying to set up a large MATMPIAIJ matrix from tuples. The tuples
>> is pairs of (row, column, value). I will use 12~20 machines for this
>> implementation. However, I think I need to know rows on each process
>> before preallocation, because for large matrices the exact
>> preallocation is important for good performance. Thus, I tried to use
>> MatGetOwnershipRange(), but I found out that I cannot use it before
>> preallocation. Could you let me know how to change my code? I am
>> attaching the code.
>>
>> Thank you,
>> Joon
>>
>>
>> PetscErrorCode SetUpMatrix(Mat *X, vector< std::tr1::tuple< PetscInt, PetscInt, PetscInt > >, PetscInt I, PetscInt J)
>> {
>>    PetscInt i, j;
>>    PetscScalar val;
>>
>>    MatCreate(PETSC_COMM_WORLD, X);
>>    MatSetType(*X, MATMPIAIJ);
>>    MatSetSizes(*X, PETSC_DECIDE, PETSC_DECIDE, I, J);
>>
>>    PetscInt begin, end;
>>    MatGetOwnershipRange(*X, &begin, &end);

>Use PetscSplitOwnership to determine sizes rows/columns, then MPI_Scan
>(or MPI_Exscan) to find your starting offset.  Alternatively, use the
>following sequence (from PetscLayoutCreate).
>
>    Notes: Typical calling sequence
>       PetscLayoutCreate(MPI_Comm,PetscLayout *);
>       PetscLayoutSetBlockSize(PetscLayout,1);
>       PetscLayoutSetSize(PetscLayout,n) or PetscLayoutSetLocalSize(PetscLayout,N);
>       PetscLayoutSetUp(PetscLayout);
>       Optionally use any of the following:
>          PetscLayoutGetSize(PetscLayout,PetscInt *); or PetscLayoutGetLocalSize(PetscLayout,PetscInt *;)
>          PetscLayoutGetRange(PetscLayout,PetscInt *rstart,PetscInt *rend); or PetscLayoutGetRanges(PetscLayout,const PetscInt *range[])
>       PetscLayoutDestroy(PetscLayout*);

>>    PetscInt localsize = end - begin;
>>
>>    PetscInt d_nnz[localsize];
>>    PetscInt o_nnz[localsize];

>This allocation will overflow your stack.  You probably want to put it
>on the heap.

>>    for (int x=0; x<tups.size(); x++) {
>>       i = std::tr1::get<0>(tups[x]);
>>       if (i>=end) break;

>Are your tuples sorted to begin with?  Are they specified redundantly or
>distributed?  Are any of the tuples nonlocal?  If it is distributed
>without matching the row distribution, then you need to communicate
>d_nnz and o_nnz.

>>       if (i>=begin) {
>>          j = std::tr1::get<1>(tups[x]);
>>          if (j>=begin && j<end)  d_nnz[i-begin]++;
>>          else                    o_nnz[i-begin]++;
>>       }
>>    }
>>    MatMPIAIJSetPreallocation(*X, 0, d_nnz, 0, o_nnz);
>>
>>    for (int x=0; x<tups.size(); x++) {
>>       i = std::tr1::get<0>(tups[x]);
>>       if (i>=end) break;
>>       if (i>=begin) {
>>          j = std::tr1::get<1>(tups[x]);
>>          val = std::tr1::get<2>(tups[x]);
>>          MatSetValues(*X, 1, &i, 1, &j, &val, INSERT_VALUES);
>>       }
>>    }
>>    MatAssemblyBegin(*X, MAT_FINAL_ASSEMBLY);
>>    MatAssemblyEnd(*X, MAT_FINAL_ASSEMBLY);
>> }