Flow123d  master-92884d111
armor.hh
Go to the documentation of this file.
1 #ifndef ARMOR_HH
2 #define ARMOR_HH
3 
4 /**
5  * The Mat class template is used for small vectors and matricies.
6  * It is just a wrapper class for an actual storage, for calculations the armadillo library
7  * is used, in combination with inlining and expression templates all auxiliary constructors
8  * and data copy are usually optimized out leaving just optimal code for the desired calculation.
9  *
10  * When deducing template parameters for the function templates the compiler consider only
11  * exact type match with possible const conversion and conversion to the parent. In particular
12  * implicit conversions are not considered. This makes transition between Armor and armadillo
13  * a bit complicated. In order to make everythink as smooth as possible we use several
14  * tricks:
15  * 1) defining a friend function inside of the class template creates an explicit function
16  * instance. As the instance is already created the implicit conversion is considered.
17  * See: https://stackoverflow.com/questions/9787593/implicit-type-conversion-with-template
18  * At least with GCC it seems, that this approach works for operators. Argument dependent lookup (ADL)
19  * finds the operator or function if at least one of its parameters is Mat template. However
20  * for operators the implicit conversion of the other argument is applied, while for the function it is not.
21  *
22  *
23  * 2) In order to consturct Mat<Type, nr, 1> and Mat<Type, 1, 1> also from arma::Col and from
24  * the Type variable, we use specialization that is derived from the generic case with nr>1, nc>1.
25  *
26  * 3) In order to prevent ambiguous method resolution, we can have just single
27  * TODO: try to transpose storage format used in Armor (prefered row first) and Armadillo (column first).
28  */
29 
30 //#define ARMA_DONT_USE_WRAPPER
31 //#define ARMA_NO_DEBUG
32 #include <armadillo>
33 #include <array>
34 #include "system/asserts.hh"
35 #include "system/logger.hh"
36 
37 
38 //
39 //template<typename T1, typename T2>
40 //inline bool is_close(
41 // const arma::BaseCube<typename T1::elem_type,T1>& a,
42 // const arma::BaseCube<typename T1::elem_type,T2>& b,
43 // const typename T1::pod_type abs_tol=1e-10, const typename T1::pod_type rel_tol=1e-6)
44 //{
45 // double a_diff = arma::norm(a.get_ref() - b.get_ref(), 1);
46 // if (a_diff < abs_tol) return true;
47 // if (a_diff < rel_tol * arma::norm(a.get_ref(), 1)) return true;
48 // return false;
49 //}
50 //
51 //template<typename T1, typename T2>
52 //inline bool is_close(
53 // const arma::Base<typename T1::elem_type,T1>& a,
54 // const arma::Base<typename T1::elem_type,T2>& b,
55 // const typename T1::pod_type abs_tol=1e-10, const typename T1::pod_type rel_tol=1e-6)
56 //{
57 // double a_diff = arma::norm(a.get_ref() - b.get_ref(), 1);
58 // if (a_diff < abs_tol) return true;
59 // if (a_diff < rel_tol * arma::norm(a.get_ref(), 1)) return true;
60 // return false;
61 //}
62 
63 
64 namespace Armor {
65 
66 //template <class type, uint nr, uint nc>
67 //struct _MatSimpleType {
68 // typedef typename arma::Mat<type>::template fixed<nr, nc> AType;
69 // inline static AType convert(const type *begin) {return AType(begin);}
70 //};
71 //
72 //template <class type, uint nr>
73 //struct _MatSimpleType<type, nr, 1> {
74 // typedef typename arma::Col<type>::template fixed<nr> AType;
75 // inline static AType convert(const type *begin) {return AType(begin);}
76 //};
77 //
78 //template <class type>
79 //struct _MatSimpleType<type, 1, 1> {
80 // typedef type AType;
81 // inline static AType convert(const type *begin) {return *begin;}
82 //};
83 
84 
85 //
86 //
87 // template <class Type, uint nRows, uint nCols>
88 // class _Mat
89 // {
90 // protected:
91 // Type * const data;
92 //
93 // public:
94 // static const uint n_rows = nRows;
95 // static const uint n_cols = nCols;
96 //
97 // typedef typename arma::Mat<Type>::template fixed<nRows, nCols> Arma;
98 // // typedef typename arma::_Mat<Type>::template fixed<nRows, nCols> BaseAr_Matype;
99 // // typedef typename __MatSimpleType<Type, nRows, nCols>::AType Ar_Matype;
100 // // typedef typename arma::Col<Type>::template fixed<nRows> VecType;
101 //
102 // _Mat(Type * const mem)
103 // : data(mem)
104 // {}
105 //
106 // inline _Mat(const Armor::_Mat<Type, nRows, nCols> & other)
107 // : data(other.data)
108 // {}
109 //
110 // inline _Mat(const Arma &arma_mat)
111 // : data(const_cast<Type *>(arma_mat.memptr()))
112 // // TODO: Here possible problem with creating non-constant Mat interface to
113 // // a temporary Arma object.
114 // {}
115 //
116 //
117 // inline operator Arma() const {return arma();}
118 //
119 // // inline operator Arma() {return arma();}
120 //
121 // // inline const Arma &arma() const {
122 // // return Arma(begin());
123 // // }
124 //
125 // inline Arma arma() const {
126 // // See proper way how to deal with l-value wrapper and r-value wrapper distinction:
127 // // https://stackoverflow.com/questions/20928044/best-way-to-write-a-value-wrapper-class
128 // // Currently we provide only r-value convertion to Arma types.
129 // // Use Index access and assignement to modify values.
130 // return Arma(begin());
131 // }
132 //
133 //
134 // inline const Type * begin() const {
135 // return data;
136 // }
137 //
138 // inline Type * begin() {
139 // return data;
140 // }
141 //
142 // inline const Type * end() const {
143 // return begin() + nCols * nRows;
144 // }
145 //
146 // inline Type * end() {
147 // return begin() + nCols * nRows;
148 // }
149 //
150 // inline uint size() const {
151 // return nRows * nCols;
152 // }
153 //
154 // inline Type * memptr() {
155 // return begin();
156 // }
157 //
158 // inline const Type & operator[](uint index) const {
159 // return data[index];
160 // }
161 //
162 // inline Type & operator[](uint index) {
163 // return data[index];
164 // }
165 //
166 // inline const Type & operator()(uint row, uint col) const {
167 // return data[row+col*nRows];
168 // }
169 //
170 // inline Type & operator()(uint row, uint col) {
171 // return data[row+col*nRows];
172 // }
173 //
174 // inline const Type & at(uint row, uint col) const {
175 // return data[row+col*nRows];
176 // }
177 //
178 // inline Type & at(uint row, uint col) {
179 // return data[row+col*nRows];
180 // }
181 //
182 //
183 // inline const _Mat<Type, nRows, nCols> & operator=(const _Mat<Type, nRows, nCols> & other) {
184 // _data_copy(other.data);
185 // return *this;
186 // }
187 //
188 // // inline const _Mat<Type, nRows, nCols> & operator=(const Arma & other) {
189 // // _data_copy(other.memptr());
190 // // return *this;
191 // // }
192 //
193 // private:
194 //
195 // inline void _data_copy(Type *other_data) {
196 // // DebugOut() << "data: " << data << "\n";
197 // // DebugOut() << "other: " << other_data << "\n";
198 // //DebugOut() << "dc size : " << n_rows * n_cols << "\n";
199 // for (uint i = 0; i < nRows * nCols; ++i) {
200 // data[i] = other_data[i];
201 // }
202 // }
203 //
204 // public:
205 // // inline const _Mat<Type, nRows, nCols> & operator=(const VecType & other) {
206 // // for (uint i = 0; i < nRows * nCols; ++i) {
207 // // data[i] = other[i];
208 // // }
209 // // return *this;
210 // // }
211 //
212 //
213 // // inline const _Mat<Type, nRows, nCols> & operator=(std::initializer_list<std::initializer_list<Type>> list) {
214 // // *this = (Arma(list));
215 // //
216 // //// const auto * listIt = list.begin();
217 // //// const Type * it;
218 // //// for (uint i = 0; i < nRows; ++i) {
219 // //// it = (listIt + i)->begin();
220 // //// for (uint j = 0; j < nCols; ++j) {
221 // //// data[i+j*nRows] = *(it + j);
222 // //// }
223 // //// }
224 // //// return *this;
225 // // }
226 // //
227 // // inline const _Mat<Type, nRows, nCols> & operator=(std::initializer_list<Type> list) {
228 // // *this = (Arma(list));
229 // //
230 // //// const Type * it = list.begin();
231 // //// for (uint i = 0; i < nRows * nCols; ++i) {
232 // //// data[i] = *(it + i);
233 // //// }
234 // //// return *this;
235 // // }
236 //
237 //
238 //
239 // // inline const _Mat<Type, nRows, nCols> & init(std::initializer_list<Type> list) {
240 // // *this = (Arma(list));
241 // // }
242 //
243 //
244 // // inline bool operator==(const Arma & other) {
245 // // const Type * first1 = begin();
246 // // const Type * last1 = end();
247 // // const Type * first2 = other.begin();
248 // // for (; first1 != last1; ++first1, ++first2) {
249 // // if (*first1 != *first2) {
250 // // return false;
251 // // }
252 // // }
253 // // return true;
254 // // }
255 // //
256 // // inline bool operator==(const _Mat<Type, nRows, nCols> & other) {
257 // // const Type * first1 = begin();
258 // // const Type * last1 = end();
259 // // const Type * first2 = other.begin();
260 // // for (; first1 != last1; ++first1, ++first2) {
261 // // if (*first1 != *first2) {
262 // // return false;
263 // // }
264 // // }
265 // // return true;
266 // // }
267 //
268 // //
269 // // inline friend bool operator==(const _Mat& a, const _Mat& b) {
270 // // return a.arma() == b.arma();
271 // // }
272 //
273 //
274 //
275 // inline friend Arma operator+=(_Mat& a, const _Mat& b) {
276 // return (a = (a.arma() + b.arma()));
277 // }
278 //
279 // inline friend Arma operator-(const _Mat& a, const _Mat& b) {
280 // return a.arma() - b.arma();
281 // }
282 //
283 // inline friend Arma operator-=(_Mat& a, const _Mat& b) {
284 // return (a = (a.arma() - b.arma()));
285 // }
286 //
287 //
288 // // // Matrix multiplication, matrix-vector multiplication.
289 // // inline friend Arma operator*=(_Mat& a, const _Mat& b) {
290 // // return (a = (a.arma() * b.arma()));
291 // // }
292 //
293 // // Elementwise multiplication.
294 // inline friend Arma operator%(const _Mat& a, const _Mat& b) {
295 // return a.arma() % b.arma();
296 // }
297 //
298 // // Elementwise multiplication.
299 // inline friend Arma operator%=(_Mat& a, const _Mat& b) {
300 // return (a = (a.arma() % b.arma()));
301 // }
302 //
303 // // Elementwise division.
304 // inline friend Arma operator/(const _Mat& a, const _Mat& b) {
305 // return a.arma() / b.arma();
306 // }
307 //
308 // // Elementwise division.
309 // inline friend Arma operator/=(_Mat& a, const _Mat& b) {
310 // return (a = (a.arma() / b.arma()));
311 // }
312 //
313 //
314 //
315 // inline friend Type dot(const _Mat& a, const _Mat& b) {
316 // return arma::dot(a.arma(), b.arma());
317 // }
318 //
319 // inline friend bool approx_equal(const _Mat &a, const _Mat &b,
320 // const char* method,
321 // double a_tol = 1e-10, double r_tol = 1e-6)
322 // {
323 // return arma::approx_equal(a.arma(), b.arma(), method, a_tol, r_tol);
324 // }
325 //
326 //
327 // //template <class Type, uint nRows, uint nCols>
328 // inline friend bool is_close(const _Mat &a, const _Mat &b, double a_tol = 1e-10, double r_tol = 1e-6) {
329 // is_close(a.arma(), b.arma(), a_tol, r_tol);
330 // }
331 //
332 //
333 //
334 //
335 // };
336 //
337 //
338 //
339 //
340 
341 // // Matrix multiplication, matrix-vector multiplication.
342 // template<class Type, uint nrA, uint ncA, uint nrB, uint ncB>
343 // inline auto operator*(
344 // const _Mat<Type, nrA, ncA>& a,
345 // const _Mat<Type, nrB, ncB>& b)
346 // -> decltype(a.arma() * b.arma())
347 // {
348 // return a.arma() * b.arma();
349 // }
350 //
351 // template<class Type, uint nrA, uint ncA, uint nrB, uint ncB>
352 // inline auto operator*(
353 // const typename _Mat<Type, nrA, ncA>::Arma& a,
354 // const _Mat<Type, nrB, ncB>& b)
355 // -> decltype(a * b.arma())
356 // {
357 // return a * b.arma();
358 // }
359 //
360 // template<class Type, uint nrA, uint ncA, uint nrB, uint ncB>
361 // inline auto operator*(
362 // const _Mat<Type, nrA, ncA>& a,
363 // const typename _Mat<Type, nrB, ncB>::Arma& b)
364 // -> decltype(a.arma() * b)
365 // {
366 // return a.arma() * b;
367 // }
368 //
369 // /**
370 // * Check for close vectors or matrices.
371 // * |a - b| < a_tol or |a-b| < r_tol * |a|
372 // *
373 // * Note: Similar function arma::approx_equal.
374 // */
375 //
376 //
377 // /** dot ********************************/
378 //
379 //
380 // /** operator + ********************************/
381 // template <class Type, uint nRows, uint nCols>
382 // inline auto operator+(const _Mat<Type, nRows, nCols> & a, const _Mat<Type, nRows, nCols> & b)
383 // ->decltype(a.arma() + b.arma())
384 // {
385 // return a.arma() + b.arma();
386 // }
387 //
388 // template <class Type, uint nRows, uint nCols, class TB>
389 // inline auto operator+(const _Mat<Type, nRows, nCols> & a, const TB & b)
390 // ->decltype(a.arma() + b)
391 // {
392 // return a.arma() + b;
393 // }
394 //
395 // template <class Type, uint nRows, uint nCols, class TB>
396 // inline auto operator+(const TB& a, const _Mat<Type, nRows, nCols> & b)
397 // ->decltype(a + b.arma())
398 // {
399 // return a + b.arma();
400 // }
401 //
402 //
403 //
404 //
405 // /** operator - ********************************/
406 // //
407 // //template <class Type, uint nRows, uint nCols>
408 // //inline typename Mat<Type, nRows, nCols>::ArmaType operator-(const Mat<Type, nRows, nCols> & a, const Mat<Type, nRows, nCols> & b) {
409 // // return a.arma() - b.arma();
410 // //}
411 // //
412 // //template <class Type, uint nRows, uint nCols>
413 // //inline typename Mat<Type, nRows, nCols>::ArmaType operator-(const Mat<Type, nRows, nCols> & a, const typename Mat<Type, nRows, nCols>::ArmaType & b) {
414 // // return a.arma() - b;
415 // //}
416 // //
417 // //template <class Type, uint nRows, uint nCols>
418 // //inline typename Mat<Type, nRows, nCols>::ArmaType operator-(const typename Mat<Type, nRows, nCols>::ArmaType & a, const Mat<Type, nRows, nCols> & b) {
419 // // return a - b.arma();
420 // //}
421 // //
422 // //template <class Type, uint nRows, uint nCols>
423 // //inline typename Mat<Type, nRows, nCols>::ArmaType operator-=(typename Mat<Type, nRows, nCols>::ArmaType & a, const Mat<Type, nRows, nCols> & b) {
424 // // return (a = (a - b.arma()));
425 // //}
426 //
427 // /** operator * ********************************/
428 // //template <class Type, uint nRows, uint nCols>
429 // //inline typename Mat<Type, nRows, nCols>::ArmaType operator*(const Mat<Type, nRows, nCols> & a, const Mat<Type, nRows, nCols> & b) {
430 // // return a.arma() * b.arma();
431 // //}
432 // //
433 // //template <class Type, uint nRows, uint nCols>
434 // //inline typename Mat<Type, nRows, nCols>::ArmaType operator*(const Mat<Type, nRows, nCols> & a, const typename Mat<Type, nRows, nCols>::ArmaType & b) {
435 // // return a.arma() * b;
436 // //}
437 // //
438 // //template <class Type, uint nRows, uint nCols>
439 // //inline typename Mat<Type, nRows, nCols>::ArmaType operator*(const typename Mat<Type, nRows, nCols>::ArmaType & a, const Mat<Type, nRows, nCols> & b) {
440 // // return a * b.arma();
441 // //}
442 //
443 //
444 // /*
445 //
446 // template <class Type, uint nRows, uint nCols>
447 // inline typename Mat<Type, nRows, nCols>::ArmaType operator%(const Mat<Type, nRows, nCols> & a, const Mat<Type, nRows, nCols> & b) {
448 // return a.arma() % b.arma();
449 // }
450 //
451 // template <class Type, uint nRows, uint nCols>
452 // inline typename Mat<Type, nRows, nCols>::ArmaType operator%(const Mat<Type, nRows, nCols> & a, const typename Mat<Type, nRows, nCols>::ArmaType & b) {
453 // return a.arma() % b;
454 // }
455 //
456 // template <class Type, uint nRows, uint nCols>
457 // inline typename Mat<Type, nRows, nCols>::ArmaType operator%(const typename Mat<Type, nRows, nCols>::ArmaType & a, const Mat<Type, nRows, nCols> & b) {
458 // return a % b.arma();
459 // }
460 // */
461 //
462 //
463 //
464 // /** scalar * and / **/
465 // //template <class Type, uint nRows, uint nCols>
466 // //inline typename Mat<Type, nRows, nCols>::ArmaType operator*(Type number, const Mat<Type, nRows, nCols> & a) {
467 // // return number * a.arma();
468 // //}
469 // //
470 // //template <class Type, uint nRows, uint nCols>
471 // //inline typename Mat<Type, nRows, nCols>::ArmaType operator/(const Mat<Type, nRows, nCols> & a, Type number) {
472 // // return a.arma() / number;
473 // //}
474 //
475 //
476 //
477 // template <uint N>
478 // using vec = Mat<double, N, 1>;
479 //
480 // template <uint N, uint M>
481 // using mat = Mat<double, N, M>;
482 //
483 //
484 // template<class Type, int nr, int nc>
485 // Mat<Type, nr, nc> mat_type(typename arma::Mat<Type>::template fixed<nr, nc> &x) {return Mat<Type, nr, nc>(x.memptr());}
486 //
487 // template<class Type>
488 // Mat<Type, 1, 1> mat_type(double &x) {return Mat<Type, 1, 1>(&x);}
489 //
490 // template<class eT, int nr>
491 // Mat<eT, nr, 1> mat_type(typename arma::Col<eT>::template fixed<nr> &x) {return Mat<eT, nr, 1>(x.memptr());}
492 //
493 //
494 //
495 // //template<Type, int nr>
496 // //class Mat<Type, nr, 1> : public Mat<Type, nr, 1>{
497 // //
498 // //}
499 //
500 
501 template <class Type, uint nr, uint nc>
502 using ArmaMat = typename arma::Mat<Type>::template fixed<nr, nc>;
503 
504 template <class Type, uint nr>
505 using ArmaVec = typename arma::Col<Type>::template fixed<nr>;
506 
507 //// assignment wrapper
508 //
509 //template <class Type, uint nRows, uint nCols>
510 //class Mat
511 //{
512 // Type *ptr;
513 //public:
514 // typedef typename ArmaMat<Type, nRows, nCols> Arma;
515 //
516 // // inherit constructors
517 // //using ArmaMat<Type, nRows, nCols>::ArmaMat;
518 //
519 // Mat(Type *ptr)
520 // : ptr(ptr)
521 // {}
522 //
523 // const Arma &operator=(const Arma &other)
524 // {
525 // for (uint i = 0; i < nRows * nCols; ++i) {
526 // *(ptr + i) = *(other.mem + i);
527 // }
528 // return *this;
529 // }
530 //
531 //};
532 //
533 //
534 //template <class Type, uint nRows>
535 //class Mat<Type, nRows, 1> : public ArmaMat<Type, nRows, 1>
536 //{
537 //public:
538 // typedef typename arma::Col<Type>::template fixed<nRows> ArmaVec;
539 // // inherit constructors
540 // using ArmaMat<Type, nRows, nCols>::ArmaMat;
541 //
542 // // Add Col constructor
543 // Mat<Type, nRows, 1>(const ArmaVec &other)
544 // : ArmaMat<Type, nRows, 1>(other.memptr())
545 // {}
546 //
547 //// const Mat<Type, nRows, 1> &operator=(const ArmaVec &other)
548 //// {
549 //// this->_data_copy(other.memptr());
550 //// return *this;
551 //// }
552 //};
553 //
554 //
555 //template <class Type>
556 //class Mat<Type, 1, 1> : public _Mat<Type, 1, 1>
557 //{
558 //public:
559 // typedef typename _Mat<Type, 1, 1>::Arma Arma;
560 // typedef typename arma::Col<Type>::template fixed<1> ArmaVec;
561 // typedef Type Scalar;
562 //
563 // // inherit constructors
564 // using _Mat<Type, 1, 1>::_Mat;
565 //
566 // // Add Col and Scalar constructor
567 // Mat<Type, 1, 1>(const ArmaVec &other)
568 // : _Mat<Type, 1, 1>(other.memptr())
569 // {}
570 //
571 // Mat<Type, 1, 1>(const Scalar &other)
572 // : _Mat<Type, 1, 1>(&other)
573 // {}
574 //
575 //// const Mat<Type, 1, 1> &operator=(const ArmaVec &other)
576 //// {
577 //// this->_data_copy(other.memptr());
578 //// return *this;
579 //// }
580 ////
581 //// const Mat<Type, 1, 1> &operator=(const Scalar &other)
582 //// {
583 //// this->_data_copy(&other);
584 //// return *this;
585 //// }
586 //
587 //};
588 //
589 
590 
591 /**
592  * Array of Armor::Mat with given shape. Provides contiguous storage for the data and access to the array elements.
593  * The shape of the matrices is specified at run time, so the class Array is independent of additional template parameters.
594  * However, to access the array elements, one must use the templated method get().
595  */
596 template<class Type>
597 class Array {
598 public:
599  class ArrayMatSet {
600  Type * ptr_;
601  uint n_rows_, n_cols_;
602  uint reserved_; ///< Reserved size of Armor::Array
603  public:
604  inline ArrayMatSet(Type *ptr, uint n_rows, uint n_cols, uint reserved)
605  : ptr_(ptr), n_rows_(n_rows), n_cols_(n_cols), reserved_(reserved) {}
606 
607 // template<class T>
608 // ArrayMatSet &operator=(const typename arma::Base<Type, T>& arma_x)
609 // {
610 // const T &derived = arma_x.get_ref();
611 // ASSERT_EQ(n_rows_, T::n_rows);
612 // ASSERT_EQ(n_cols_, T::n_cols);
613 // copy<T::n_rows, T::n_cols>(derived.memptr());
614 // }
615 
616  template<long long unsigned int nr, long long unsigned int nc>
617  ArrayMatSet &operator=(const ArmaMat<Type, nr, nc>& arma_x)
618  {
619  ASSERT_EQ(n_rows_, nr);
620  ASSERT_EQ(n_cols_, nc);
621  copy<nr, nc>(arma_x.memptr());
622  return *this;
623  }
624 
625  template<long long unsigned int nr>
626  ArrayMatSet &operator=(const ArmaVec<Type, nr>& arma_x)
627  {
628  ASSERT_EQ(n_rows_, nr);
629  ASSERT_EQ(n_cols_, 1);
630  copy<nr, 1>(arma_x.memptr());
631  return *this;
632  }
633 
634 
635  ArrayMatSet &operator=(const Type& arma_x)
636  {
637  ASSERT_EQ(n_rows_, 1);
638  ASSERT_EQ(n_cols_, 1);
639  copy<1, 1>(&arma_x);
640  return *this;
641  }
642 
643 
644  template <uint nr, uint nc>
645  void copy(const Type *other_ptr) {
646  for (uint i = 0; i < nr * nc; ++i) {
647  *(ptr_ + i * reserved_) = *(other_ptr + i);
648  }
649  }
650  };
651 
652 public:
653  /**
654  * Construct array of Armor matrices.
655  * @param nv Number of matrices in the array.
656  * @param nr Number of rows in each matrix.
657  * @param nc Number of columns in each matrix.
658  */
659  Array(uint nr, uint nc = 1, uint size = 0)
660  : data_(new Type[nr * nc * size]),
661  n_rows_(nr),
662  n_cols_(nc),
663  size_(size),
664  reserved_(size)
665  {
666  }
667 
668  Array(const Array &other)
669  : Array(other.n_rows_, other.n_cols_, other.size_)
670  {
671  for(uint i = 0; i < n_rows_ * n_cols_ * size(); i++) {
672  data_[i] = other.data_[i];
673  }
674  }
675 
676  ~Array() {
677  delete [] data_;
678  data_ = nullptr;
679  }
680 
681  Array &operator=(const Array &other)
682  {
683  ASSERT( (n_rows_ == other.n_rows_) && (n_cols_ == other.n_cols_) );
684  reinit(other.size());
685  resize(other.size());
686 
687  for(uint i = 0; i < n_rows_ * n_cols_ * size(); i++) {
688  data_[i] = other.data_[i];
689  }
690  return *this;
691  }
692 
693  /**
694  * Drop all data and allocate new space of given size.
695  * Change number of elements in the array, while keeping the shape of arrays.
696  * @param size New size of array.
697  */
698  void reinit(uint size) {
699  delete [] data_;
700  data_ = nullptr;
701  reserved_ = size;
702  size_ = 0;
703  data_ = new Type[n_rows_ * n_cols_ * reserved_];
704  }
705 
706 
707  /**
708  * Resize active part of the allocated space.
709  */
710  void resize(uint size) {
711  ASSERT_LE(size, reserved_);
712  size_ = size;
713  }
714 
715  inline uint n_rows() const
716  {
717  return n_rows_;
718  }
719 
720  inline uint n_cols() const
721  {
722  return n_cols_;
723  }
724 
725  /**
726  * Get size of active space.
727  */
728  inline unsigned int size() const {
729  return size_;
730  }
731 
732  /**
733  * Increase active space by 1 and store given Mat value to the end of the active space.
734  */
735  template<unsigned long long int nr, unsigned long long int nc = 1>
736  inline void append(const ArmaMat<Type,nr,nc> &item)
737  {
738  ASSERT_LT(size_, reserved_);
739  size_ += 1;
740  set(size_ - 1) = item;
741 
742  }
743 
744  template<unsigned long long int nr>
745  inline void append(const ArmaVec<Type,nr> &item)
746  {
747  append<nr, 1>(ArmaMat<Type,nr,1>(item));
748  }
749 
750  /**
751  * Return matrix at given position in array. The returned object is a Armor::Mat
752  * pointing to the respective data_ block in the Array's storage.
753  * One can assign to the Armor::Mat which performs postponed evaluation and storing the result to the array.
754  *
755  * @param i Index of the matrix.
756  * TODO: Should be renamed to item(), but we have compilation problem in Field::loc_point_value
757  */
758 // template<uint nr, uint nc = 1>
759 // inline Mat<Type,nr,nc> get(uint mat_index) const
760 // {
761 // ASSERT( (nr == n_rows_) && (nc == n_cols_) );
762 // return Mat<Type,nr,nc>( data_ + mat_index * n_rows_ * n_cols_ );
763 // }
764 //
765 // template<uint nr, uint nc = 1>
766 // inline Mat<Type,nr,nc> get(uint mat_index)
767 // {
768 // ASSERT( (nr == n_rows_) && (nc == n_cols_) );
769 // return Mat<Type,nr,nc>( data_ + mat_index * n_rows_ * n_cols_ );
770 // }
771 
772 
773 // template<uint nr, uint nc = 1>
774 // inline Mat<Type,nr,nc> get(uint mat_index) const
775 // {
776 // ASSERT( (nr == n_rows_) && (nc == n_cols_) );
777 // return Mat<Type,nr,nc>( data_ + mat_index * n_rows_ * n_cols_ );
778 // }
779 
780  template<uint nr, uint nc = 1>
781  inline ArmaMat<Type,nr,nc> mat(uint mat_index) const
782  {
783  ASSERT( (nr == n_rows_) && (nc == n_cols_) );
784 // ArmaMat<Type,nr,nc> m;
785 // double ** ptr = const_cast<double **>(&(m.mem));
786 // *ptr = data_ + mat_index * n_rows_ * n_cols_;
787 // return m;
788  ArmaMat<Type,nr,nc> mat;
789  for (uint row=0; row<n_rows_; ++row)
790  for (uint col=0; col<n_cols_; ++col)
791  mat(row,col) = data_[mat_index + (col*n_rows_+row) * reserved_];
792  return mat;
793  }
794 
795 // template<long long unsigned int nr, long long unsigned int nc = 1>
796 // inline void set(uint mat_index, const ArmaMat<Type,nr,nc> &mat)
797 // {
798 // ASSERT( (nr == n_rows_) && (nc == n_cols_) );
799 //
800 // for (uint i = 0; i < nr * nc; ++i) {
801 // *(data_ + i) = *(mat.mem + i);
802 // }
803 //
804 //
805 //// ArmaMat<Type, nr, nc> m;
806 //// double ** ptr = const_cast<double **>(&(m.mem));
807 //// *ptr = data_ + mat_index * n_rows_ * n_cols_;
808 //// m = mat;
809 // }
810 
811 
812 
813 // template<uint nr, uint nc = 1>
814 // inline Vec<Type,nr,nc> get_vec(uint mat_index) const
815 // {
816 // ASSERT( (nr == n_rows_) && (nc == n_cols_) );
817 // return Vec<Type,nr,nc>( data_ + mat_index * n_rows_ * n_cols_ );
818 // }
819 
820  template<uint nr>
821  inline ArmaVec<Type, nr> vec(uint mat_index) const
822  {
823  ASSERT( (nr == n_rows_) && (1 == n_cols_) )(n_rows_)(n_cols_);
824  ASSERT_LT(mat_index, size());
825  ArmaVec<Type, nr> vec;
826  for (uint i=0; i<n_rows_; ++i)
827  vec(i) = data_[mat_index + i * reserved_];
828  return vec;
829  }
830 
831  inline Type scalar(uint mat_index) const
832  {
833  ASSERT( (1 == n_rows_) && (1 == n_cols_) )(n_rows_)(n_cols_);
834  ASSERT_LT(mat_index, size());
835  return data_[mat_index];
836  }
837 
838  inline ArrayMatSet set(uint index) {
839  ASSERT_LT(index, size());
840  return ArrayMatSet(data_ + index, n_rows_, n_cols_, reserved_);
841  }
842 
843 
844  /**
845  * Return armadillo matrix at given position in array.
846  * @param i Index of matrix.
847  */
848  inline arma::mat arma_mat(uint i) const
849  {
850  ASSERT_LT(i, size());
851  arma::mat mat(n_rows_, n_cols_);
852  for (uint row=0; row<n_rows_; ++row)
853  for (uint col=0; col<n_cols_; ++col)
854  mat(row,col) = data_[i + (col*n_rows_+row) * reserved_];
855  return mat;
856  }
857 
858  /**
859  * Return armadillo vector at given position in array.
860  * Warning! Method can be used only if nCols == 1.
861  * @param i Index of matrix.
862  */
863  inline arma::vec arma_vec(uint i) const
864  {
865  ASSERT_LT(i, size());
866  ASSERT_EQ(n_cols_, 1);
867  arma::vec vec(n_rows_);
868  for (uint row=0; row<n_rows_; ++row)
869  vec(i) = data_[i + row * reserved_];
870  return vec;
871  }
872 
873  Type * data_;
874 
875 private:
876  inline uint space_() { return n_rows_ * n_cols_ * reserved_; }
877  uint n_rows_;
878  uint n_cols_;
879  uint size_;
880  uint reserved_;
881 };
882 
883 
884 template <uint N>
885 using vec = ArmaVec<double, N>;
886 
887 template <uint N, uint M>
888 using mat = ArmaMat<double, N, M>;
889 
890 using array = Array<double>;
891 }
892 
893 #endif
Armor::Array::append
void append(const ArmaVec< Type, nr > &item)
Definition: armor.hh:745
Armor::Array::ArrayMatSet::operator=
ArrayMatSet & operator=(const ArmaMat< Type, nr, nc > &arma_x)
Definition: armor.hh:617
Armor::vec
ArmaVec< double, N > vec
Definition: armor.hh:885
Armor::Array::ArrayMatSet
Definition: armor.hh:599
Armor::Array::arma_vec
arma::vec arma_vec(uint i) const
Definition: armor.hh:863
Armor::Array::ArrayMatSet::ptr_
Type * ptr_
Definition: armor.hh:600
ASSERT
#define ASSERT(expr)
Definition: asserts.hh:350
Armor::Array::vec
ArmaVec< Type, nr > vec(uint mat_index) const
Definition: armor.hh:821
asserts.hh
Definitions of ASSERTS.
Armor::Array::ArrayMatSet::reserved_
uint reserved_
Reserved size of Armor::Array.
Definition: armor.hh:602
Armor::Array::ArrayMatSet::n_rows_
uint n_rows_
Definition: armor.hh:601
Armor::Array::ArrayMatSet::operator=
ArrayMatSet & operator=(const ArmaVec< Type, nr > &arma_x)
Definition: armor.hh:626
Armor::Array::Array
Array(const Array &other)
Definition: armor.hh:668
Armor::Array::ArrayMatSet::n_cols_
uint n_cols_
Definition: armor.hh:601
Armor::Array::resize
void resize(uint size)
Definition: armor.hh:710
uint
unsigned int uint
Definition: mh_dofhandler.hh:101
Armor::Array::n_cols_
uint n_cols_
Definition: armor.hh:878
Armor::Array::reserved_
uint reserved_
Definition: armor.hh:880
Armor::Array::mat
ArmaMat< Type, nr, nc > mat(uint mat_index) const
Definition: armor.hh:781
ASSERT_LT
#define ASSERT_LT(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
Definition: asserts.hh:300
Armor::Array::n_rows
uint n_rows() const
Definition: armor.hh:715
Armor::Array::Array
Array(uint nr, uint nc=1, uint size=0)
Definition: armor.hh:659
Armor
Definition: armor.hh:64
Armor::mat
ArmaMat< double, N, M > mat
Definition: armor.hh:888
Armor::Array::operator=
Array & operator=(const Array &other)
Definition: armor.hh:681
Armor::Array::arma_mat
arma::mat arma_mat(uint i) const
Definition: armor.hh:848
Armor::Array::ArrayMatSet::operator=
ArrayMatSet & operator=(const Type &arma_x)
Definition: armor.hh:635
ASSERT_EQ
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:332
logger.hh
Armor::Array::n_cols
uint n_cols() const
Definition: armor.hh:720
Armor::Array::set
ArrayMatSet set(uint index)
Definition: armor.hh:838
Armor::Array::size
unsigned int size() const
Definition: armor.hh:728
Armor::Array::size_
uint size_
Definition: armor.hh:879
Armor::Array::reinit
void reinit(uint size)
Definition: armor.hh:698
Armor::Array::append
void append(const ArmaMat< Type, nr, nc > &item)
Definition: armor.hh:736
ASSERT_LE
#define ASSERT_LE(a, b)
Definition of comparative assert macro (Less or Equal) only for debug mode.
Definition: asserts.hh:308
Armor::Array
Definition: armor.hh:597
Armor::Array::data_
Type * data_
Definition: armor.hh:873
Armor::Array::n_rows_
uint n_rows_
Definition: armor.hh:877
Armor::Array::scalar
Type scalar(uint mat_index) const
Definition: armor.hh:831
Armor::Array::space_
uint space_()
Definition: armor.hh:876
Armor::ArmaVec
typename arma::Col< Type >::template fixed< nr > ArmaVec
Definition: armor.hh:505
Armor::Array::~Array
~Array()
Definition: armor.hh:676
Armor::ArmaMat
typename arma::Mat< Type >::template fixed< nr, nc > ArmaMat
Definition: armor.hh:502
Armor::Array::ArrayMatSet::copy
void copy(const Type *other_ptr)
Definition: armor.hh:645
Armor::Array::ArrayMatSet::ArrayMatSet
ArrayMatSet(Type *ptr, uint n_rows, uint n_cols, uint reserved)
Definition: armor.hh:604
Generated by   doxygen 1.8.17