Hello Jørgen, I think the deal.II example is a good example how to incorporate PETSc bindings into a finite element code. I don't think it is a good idea to support an operator [][] for matrices for the reasons that you mentioned below. You usually get away with the kind of SetMatrixEntry(i,j,value) interface that PETSc provides. You can also take a look into the DUNE module dune-fem (http://dune.mathematik.uni-freiburg.de/) where we have implemented parallel PETSc bindings was well as bindings to the BCRSMatrix from dune-istl. In dune-fem we support local matrix views (element wise) and then only allow to set or add values to matrix entries. So far this covers all application scenarios without making it impossible to interface to such different packages as PETSc or dune-istl or other matrix implementations.
Best, Robert -- I'm currently playing with implementing (proper) support for petsc in OPM. For this we need better containers and interfaces than petsc's own. The thing is, petsc doesn't support direct memory access. While this is good object oriented practice, it creates some headaches when implementing references to the scalar in an arbitrary cell. The question is: do we want at()/operator[][] for matrices at all? Or is it sufficient to say "hey pal, if you want to access the value in a specific cell you need to go through iterators"?. While it is possible to implement reference-esque features towards petsc, it will be subject to data races. I looked up deal.II, a FEM library that can use petsc, and they do not support random access for matrices. http://www.dealii.org/8.2.0/doxygen/deal.II/classPETScWrappers_1_1BlockSparseMatrix.html Regarding data races: Iterators are needed, that's a no-brainer. I think of doing this the deal.II way, i.e. reading per-row into a buffer which in turn writes to the matrix in the case of operator*. This again can either lead to a data race if there is no global look-up mechanism for uniquely spawned instances to share buffers (which again provides races :)), or we can simply declare it undefined behaviour to have two live iterators to the same row. Since petsc uses much MPI this isn't too unreasonable. Of course, we can extend the same "this-is-undefined" statement to also concern scalar references, but I personally think that this much undefined behaviour of unrelated means of access is just bad design, and I'd rather omit it. So: do we need or want per-cell references to matrices? Finally, some code examples. Keep in mind that petsc ONLY supports GetRow(). Consider this case of iterators. Perfectly defined and as you expect. matrix m; for( auto itr = m.begin(); itr != m.end(); ++itr ) *itr += 3; Now, another well-formed example. Notice that other_itr was created through itr. They share underlying buffers. for( auto itr = m.begin(); itr != m.end(); ++itr ) { auto other_itr = itr; while( other_itr++ != m.end() ) { *other_itr += 3; } } An undefined example: auto itr1 = m.begin(); *itr1 += 3; auto itr2 = m.begin() + 2; *itr2 += 4; ++itr1; *itr1 += 2; // POSSIBLY UNDEFINED. *itr1 += 2 may or may not be undefined (by us). It's ok in the sense that the underlying matrix isn't written in this cell by itr2, which is ok. However, assuming writes aren't deferred itr2's buffers up until line 4 are identical to itr1's. Should itr1 increment once more it will not have received writes from itr2 and reading would definitely be undefined. This can be partially solved with subscription and (global) cache invalidation, but personally I think this adds more complexity than gain, so I'm with non-enforced read-on-write being undefined. Unfortunately petsc introduces the issues of concurrent read and write even in the single-CPU single-threaded case, and provides no elegant means of preventing it. _______________________________________________ Opm mailing list [email protected] http://www.opm-project.org/mailman/listinfo/opm [attachment "signature.asc" removed by Robert Klöfkorn/IRIS] "This e-mail with attached documents is only for the addressee indicated above. The e-mail and documents may contain confidential information. If you are not the correct addressee or are responsible for transfer of this e-mail with attached documents to the correct addressee, any copying or further transfer of information contained herein is expressly forbidden. If you have received this e-mail by mistake, you are kindly requested to immediately inform me thereof by e-mail and delete the e-mail and the received documents." http://www.iris.no _______________________________________________ Opm mailing list [email protected] http://www.opm-project.org/mailman/listinfo/opm
