armadillo drawbacks

[jbrezmorf]

Armadillo provides very fast algebra for dens matrices and objects with very simple syntax. How ever there are two drawbacks for its usage in Flow123d:

1. Its fixed size objects (namely vectors) are very big which makes them inconvenient for storing e.g. coordinates in 3D space. Possible solutions:
   • fork or patch Armadillo to have true fixed objects. This could be possibly done by turning arma::Mat<> into more complex template with some parameters with default value.
   • make our own library
2. Armadillo matrices are stored in column first ordering. This is inconvenient for matrix assembly since local matrices should be in row first ordering. This also should be faster. But this speed argument has to be checked (e.g. PETSC allows both assembly types)
3. Fixed size objects are not literals, so could not be used as in const expr variables and functions.

[jbrezmorf]

Proposed concept for the first drawback:

1. • add our own test that measure: a. performance b. memory for fixed Mat type
2. • add one more template parameter to arma::Mat<>, Shape.
   • Define Shape class template. Generic and Shape<0,0> specialization.
   • Define Mat as Mat<eT, Shape<0,0>> Test.
   • Define Mat::fixed<N,M> as derived from Mat<eT, Shape<N,M>> Test.
3. • add Shape member to Mat<>
   • move protected mem_local member into Shape. Test.
   • use fixed space (see Mat<>::fixed<>::mem_local_extra) in generic Shape<> template, use preallocation, in Shape<0,0> specialization.
4. • declare n_rows, n_cols, n_elem, vec_state, mem_state as const references
   • declare mem using a typedef from State
   • set them in all constructors to refer to equivalent Shape members.
   • Test.
5. • Define Mat<eT>::fixed<N,M> as Mat<eT, Shape<N,M>> (or derived)
   • Introduce specialization Shape<0,0> with const static members for n_rows, n_cols, n_elem, ... (reuse from Mat::fixed)

   for both fixed dimensions nonzero we try to set all types into auxiliary classes which contains no data.

6. introduce Mat<...>::nrows(), ncols() methods and use them through the library instead of nrows, ncols attributes

7. derive final classes exposed to users, the fixed variant defines nrows, ncols attributes as static consts, while the free size variant defines them as current implementation and set them in constructors (these have to be defined)

Question: how can expression templates work with Mat_fixed derived from the Mat implementation?

[jbrezmorf]

Point 1. resolved by introduction of Armor interface classes. Point 2. can be resolved by own implementation of matrices in the LocalSystem. Point 3. ... not clear why this is an issue. Armadillo fixed Mat are not literals probably because of its derivation from general Mat class. So possibly Armor Mat can resolve this issue as well.

[jbrezmorf]

Resolution:

• Continue current approach : storing data in ArmorArray and provide temporary construction of Aramdilo objects with these data
• For critical parts of code implement vectorized operations in ArmorArray: FieldModel, FEValues