[petsc-users] About MatGetOwnershipRange()

[Joon Hee Choi]

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);
   PetscInt localsize = end - begin;

   PetscInt d_nnz[localsize];
   PetscInt o_nnz[localsize];

   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);
}