Iterating over all the entries in a PETScWrapper::MPI::SparseMatrix in parallel

[Feimi Yu]

Hi,

I'm using PETScWrapper to parallelize my code. In my preconditioner for the GMRES solver, there is one step that requires a matrix copied from the system matrix, and set all the elements to be the absolute value. It was fine in serial because I can iterator over all the entries simply using the iterator, but it does not work in parallel since each process only owns local entries. When I run the code using iteration the error message will tell me "argument out of range; Not local rows". It also tells me "Object is in wrong state; Not for unassembled matrix" despite the fact that I called compress(VectorOperation:add) right before the iteration.

Is there any way to bypass the iteration to get an absolute value matrix or to iterate locally and not trigger the "unassembled matrix" error?

Thanks,

Feimi

[Feimi Yu]

An update:

I tried to use a iteration below to iterate over local entries: (The reason I use local_range() for only (0, 0) block and iterator for the entire block matrix is because I only need the block(0, 0), and sparse matrix class does not have a non-const iterator, I have to call the local range in the (0, 0) block and initialize an iterator from the blockmatrix)

for (auto r = Abs_A_matrix->block(0, 0).local_range().first;

r < Abs_A_matrix->block(0, 0).local_range().second; ++r)

{

for (auto itr = Abs_A_matrix->begin(r);

itr != Abs_A_matrix->end(r);

++itr)

{

itr->set_value(std::abs(itr->value()));

}

}

However, it doesn't work. The error message says:

dealii/install/include/deal.II/lac/matrix_iterator.h:45:41:   required from ‘class dealii::MatrixIterator<dealii::BlockMatrixIterators::Accessor<dealii::BlockMatrixBase<dealii::PETScWrappers::MPI::SparseMatrix>, false> >’

dealii/install/include/deal.II/lac/block_matrix_base.h:178:51: error: no type named ‘iterator’ in ‘dealii::BlockMatrixBase<dealii::PETScWrappers::MPI::SparseMatrix>::BlockType {aka class dealii::PETScWrappers::MPI::SparseMatrix}’
     typename BlockMatrixType::BlockType::iterator base_iterator;
                                                   ^~~~~~~~~~~~~

Thanks for any help!
Feimi

[Feimi Yu]

Second update (sorry for so many updates)

I changed my strategy to use set(r, c, v) function to set the values so that I can use the const iterators. also called compress after every add:

for (auto r = Abs_A_matrix->block(0, 0).local_range().first;

r < Abs_A_matrix->block(0, 0).local_range().second; ++r)

{

for (auto itr = Abs_A_matrix->block(0, 0).begin(r);

itr != Abs_A_matrix->block(0, 0).end(r);

++itr)

{

Abs_A_matrix->block(0, 0).set(itr->row(), itr->column(), std::abs(itr->value()));

Abs_A_matrix->block(0, 0).compress(VectorOperation::add);

}

}

This time it says:
*** Error in ' ' :free(): invalid size: 0x0000556443d63e50 ***

This seems like a out of bounds access. This happens even when I run with 1 process.

Thanks,

Feimi

On Monday, March 5, 2018 at 12:21:46 PM UTC-5, Feimi Yu wrote:

[Wolfgang Bangerth]

On 03/05/2018 02:54 PM, Feimi Yu wrote:
>
> I changed my strategy to use set(r, c, v) function to set the values so that I
> can use the const iterators. also called compress after every add:

> for(autor =Abs_A_matrix->block(0, 0).local_range().first;

> r <Abs_A_matrix->block(0, 0).local_range().second; ++r)
> {
> for(autoitr =Abs_A_matrix->block(0, 0).begin(r);

> itr !=Abs_A_matrix->block(0, 0).end(r);

> ++itr)
> {
> Abs_A_matrix->block(0, 0).set(itr->row(), itr->column(), std::abs(itr->value()));
> Abs_A_matrix->block(0, 0).compress(VectorOperation::add);
> }
> }
>
> This time it says:
> *** Error in ' ' :free(): invalid size: 0x0000556443d63e50 ***
>
> This seems like a out of bounds access. This happens even when I run with 1
> process.

This is the right approach -- you do need to restrict work to the locally
owned range of rows, but if you do that on every processor, then you cover all
rows.

But you can't call 'compress()' after every set operation. That's because
compress() is a *collective* operation, i.e., it has to be called on all
processors at the same time. Since different processors may have different
numbers of entries, what you are doing here leads to different numbers of
calls on different processors.

In addition, it is *very* inefficient to call the function so many times. Just
move the call to after the loop. (And change the argument to
VectorOperation::set -- you are *setting* elements, after all, not adding to
them.)

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Feimi Yu]

Hi Wolfgang,

This time I used VectorOperation::insert and it didn't happen the memory error that I posted before. 

However, if I put the compress function after the loop, it only sets one entry on each rank then throws the exception

 "Object is in wrong state, not for unassembled matrix." 

Putting the compress() in the loop does not throw exceptions, instead, it iterates forever. 

So just as you said it is very inefficient and not an option.

[Wolfgang Bangerth]

> This time I used VectorOperation::insert and it didn't happen the memory error
> that I posted before.

OK, so that then clearly helps :-)

> However, if I put the compress function after the loop, it only sets one entry
> on each rank then throws the exception
> "Object is in wrong state, not for unassembled matrix."
> Putting the compress() in the loop does not throw exceptions, instead, it
> iterates forever.
> So just as you said it is very inefficient and not an option.

I see. Yes, this is an effect of the PETSc memory model then. I'm not sure how
to solve this other than to not iterate over the entries of the matrix after
assembly :-(

[Feimi Yu]

Hi Wolfgang,

Fortunately, I managed to solve that problem. I found that every single operation on the iterator needs to call row_length(), which
requires the assembled status of the matrix, and apparently set() operation will break this status. My solution is to iterate and cache
the rows, columns, values information in a big vector and use another loop to set them all in a time without the use of iterator.

Now I have another issue, which is I think is not so much reasonable in the way of implementation: (but I cannot think of a better way,
:(

The problem is that I still encounter the "out of range" problem even when I do iterator over the local rows. I debugged my code and
checked the source code, and found where the problem is:
The end iterator of each row is pointing to the first entry of the next row (or true end if the row is the end row).
However when we have multiple ranks, when one wants to get the end iterator of the last LOCAL row, the end
iterator will point to the next row, which is not belong to local.
Let me take an example to be clear:
I have 1223 rows, where rank 0 has the information of row 0 to 633 and rank 1 has 634 to 1222.
When I iterate over the entries on rank 1, everything works as it is.
When I iterate over the entries on rank 0, especially when it reaches row 633, I will call matrix.end(633)
However this end is essentially matrix.begin(634)! So PETSc tells me I am out of range, because rank 0
does not have information of row 634.

I have attached the source code here for your reference:
in petsc__matrix__base.h file
1449  inline

1450  MatrixBase::const_iterator

1451  MatrixBase::end(const size_type r) const

1452  {

1453  Assert (r < m(), ExcIndexRange(r, 0, m()));

1454 

1455  // place the iterator on the first entry

1456  // past this line, or at the end of the

1457  // matrix

1458  for (size_type i=r+1; i<m(); ++i)

1459  if (row_length(i) > 0)

1460  return const_iterator(this, i, 0);

1461 

1462  // if there is no such line, then take the

1463  // end iterator of the matrix

1464  return end();

1465  }

Thanks,
Feimi

[Wolfgang Bangerth]

On 03/08/2018 02:55 PM, Feimi Yu wrote:
>
> The problem is that I still encounter the "out of range" problem even when I
> do iterator over the local rows. I debugged my code and
> checked the source code, and found where the problem is:
> The end iterator of each row is pointing to the first entry of the next row
> (or true end if the row is the end row).
> However when we have multiple ranks, when one wants to get the end iterator of
> the last LOCAL row, the end
> iterator will point to the next row, which is not belong to local.
> Let me take an example to be clear:
> I have 1223 rows, where rank 0 has the information of row 0 to 633 and rank 1
> has 634 to 1222.
> When I iterate over the entries on rank 1, everything works as it is.
> When I iterate over the entries on rank 0, especially when it reaches row 633,
> I will call matrix.end(633)
> However this end is essentially matrix.begin(634)! So PETSc tells me I am out
> of range, because rank 0
> does not have information of row 634.

Ah, yes, I see how this can happen.

Do you think you could come up with a small testcase that shows this? This
would make it much easier to debug the problem. In essence, all you'd have to
do is create a matrix on two processors, put a few values into it (say, one
value per row) and create the iterators. The whole thing doesn't have to do
anything useful -- no need to create a mesh, assemble a linear system using
FEValues, etc --, it only has to demonstrate the exception you get.

[Feimi Yu]

I'm sorry for the late reply. Here is my small testcase:

I just add one line to call the end iterator of the last local row in the built-in test case reinit_preconditioner_01.cc
under tests/petsc folder:
auto itr = mat.end(mat.local_range().second);

It produces the same error as I mentioned before.

I added the file for your reference.

Thanks!
Feimi

[Feimi Yu]

Please ignore my last post. I made a mistake there.
Attached is the revised version to better illustrate the problem.

Thanks,
Feimi

[Wolfgang Bangerth]

On 03/18/2018 04:41 PM, Feimi Yu wrote:
> Please ignore my last post. I made a mistake there.
> Attached is the revised version to better illustrate the problem.

Great, much appreciated -- a most excellent testcase! I can reproduce
the problem and have something that may work. Will finish this tomorrow!

[Wolfgang Bangerth]

On 03/18/2018 04:41 PM, Feimi Yu wrote:

> Please ignore my last post. I made a mistake there.
> Attached is the revised version to better illustrate the problem.

Patch is now here:
https://github.com/dealii/dealii/pull/6087

[Feimi Yu]

Got it. Thank you so much!

Thanks,
Feimi