On 01/13/2015 11:32 AM, Jørgen Kvalsvik wrote: > Hi, > > 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.
Oh, and a quick addendum: I'm not sure if this is the right project for it, but a bolder choice would be to say "screw it" to iterators and completely redesign all matrix operations in the spirit of matrix algebra.
signature.asc
Description: OpenPGP digital signature
_______________________________________________ Opm mailing list [email protected] http://www.opm-project.org/mailman/listinfo/opm
