37 #ifndef BGEOT_TENSOR_H__
38 #define BGEOT_TENSOR_H__
53 class multi_index :
public std::vector<size_type> {
56 void incrementation(
const multi_index &m) {
57 iterator it = begin(), ite = end();
58 const_iterator itm = m.begin();
61 while (*it >= *itm && it != (ite-1)) { *it = 0; ++it; ++itm; ++(*it); }
65 void reset() { std::fill(begin(), end(), 0); }
67 inline bool finished(
const multi_index &m) {
71 return ((*
this)[size()-1] >= m[size()-1]);
74 multi_index(
size_t n) : std::vector<
size_type>(n)
75 { std::fill(begin(), end(),
size_type(0)); }
78 { (*this)[0] = i; (*this)[1] = j; }
81 { (*this)[0] = i; (*this)[1] = j; (*this)[2] = k; }
84 { (*this)[0] = i; (*this)[1] = j; (*this)[2] = k; (*this)[3] = l; }
88 bool is_equal(
const multi_index &m)
const {
89 if (this->size() != m.size())
return false;
91 if (m[i] != (*
this)[i])
return false;
97 for (
size_type k = 0; k < this->size(); ++k) s *= (*
this)[k];
102 return std::vector<size_type>::capacity()*
sizeof(
size_type) +
108 const multi_index& mi) {
109 multi_index::const_iterator it = mi.begin(), ite = mi.end();
112 for ( ; it != ite; ++it)
113 {
if (!f) o <<
", "; o << *it; f =
false; }
118 template<
class T>
class tensor :
public std::vector<T> {
127 typedef typename std::vector<T>::iterator iterator;
128 typedef typename std::vector<T>::const_iterator const_iterator;
130 template<
class CONT>
inline const T& operator ()(
const CONT &c)
const {
131 typename CONT::const_iterator it = c.begin();
132 multi_index::const_iterator q = coeff_.begin(), e = coeff_.end();
134 multi_index::const_iterator qv = sizes_.begin();
137 for ( ; q != e; ++q, ++it) {
139 GMM_ASSERT2(*it < *qv++,
"Index out of range.");
141 return *(this->begin() + d);
146 GMM_ASSERT2(order() == 4,
"Bad tensor order.");
147 size_type d = coeff_[0]*i + coeff_[1]*j + coeff_[2]*k + coeff_[3]*l;
148 GMM_ASSERT2(d < size(),
"Index out of range.");
149 return *(this->begin() + d);
153 GMM_ASSERT2(order() == 3,
"Bad tensor order.");
154 size_type d = coeff_[0]*i + coeff_[1]*j + coeff_[2]*k;
155 GMM_ASSERT2(d < size(),
"Index out of range.");
156 return *(this->begin() + d);
160 GMM_ASSERT2(order() == 2,
"Bad tensor order");
162 GMM_ASSERT2(d < size(),
"Index out of range.");
163 return *(this->begin() + d);
168 GMM_ASSERT2(order() == 4,
"Bad tensor order.");
169 size_type d = coeff_[0]*i + coeff_[1]*j + coeff_[2]*k + coeff_[3]*l;
170 GMM_ASSERT2(d < size(),
"Index out of range.");
171 return *(this->begin() + d);
176 GMM_ASSERT2(order() == 3,
"Bad tensor order.");
177 size_type d = coeff_[0]*i + coeff_[1]*j + coeff_[2]*k;
178 GMM_ASSERT2(d < size(),
"Index out of range.");
179 return *(this->begin() + d);
183 GMM_ASSERT2(order() == 2,
"Bad tensor order.");
185 GMM_ASSERT2(d < size(),
"Index out of range.");
186 return *(this->begin() + d);
189 template<
class CONT>
inline T& operator ()(
const CONT &c) {
190 typename CONT::const_iterator it = c.begin();
191 multi_index::iterator q = coeff_.begin(), e = coeff_.end();
193 for ( ; q != e; ++q, ++it) d += (*q) * (*it);
194 GMM_ASSERT2(d < size(),
"Index out of range.");
195 return *(this->begin() + d);
198 inline size_type size()
const {
return std::vector<T>::size(); }
200 inline const multi_index &sizes()
const {
return sizes_; }
201 inline size_type order()
const {
return sizes_.size(); }
203 void init(
const multi_index &c) {
206 sizes_ = c; coeff_.resize(c.size());
207 auto p = coeff_.begin(), pe = coeff_.end();
208 for ( ; p != pe; ++p, ++it) { *p = d; d *= *it; }
212 inline void init() { sizes_.resize(0); coeff_.resize(0); this->
resize(1); }
215 sizes_.resize(1); sizes_[0] = i; coeff_.resize(1); coeff_[0] = 1;
220 sizes_.resize(2); sizes_[0] = i; sizes_[1] = j;
221 coeff_.resize(2); coeff_[0] = 1; coeff_[1] = i;
226 sizes_.resize(3); sizes_[0] = i; sizes_[1] = j; sizes_[2] = k;
227 coeff_.resize(3); coeff_[0] = 1; coeff_[1] = i; coeff_[2] = i*j;
233 sizes_[0] = i; sizes_[1] = j; sizes_[2] = k; sizes_[3] = k;
235 coeff_[0] = 1; coeff_[1] = i; coeff_[2] = i*j; coeff_[3] = i*j*k;
239 inline void adjust_sizes(
const multi_index &mi) { init(mi); }
240 inline void adjust_sizes() { init(); }
241 inline void adjust_sizes(
size_type i) { init(i); }
246 { init(i, j, k, l); }
249 const multi_index &mi = t.sizes_;
size_type d = mi.size();
250 sizes_.resize(d); coeff_.resize(d);
252 std::copy(mi.begin(), mi.end(), sizes_.begin());
253 std::copy(t.coeff_.begin(), t.coeff_.end(), coeff_.begin());
264 inline void remove_unit_dim() {
267 for (; i < sizes_.size(); ++i)
268 if (sizes_[i] != 1) { sizes_[j]=sizes_[i]; coeff_[j]=coeff_[i]; ++j; }
278 void mat_reduction(
const tensor &t,
const gmm::dense_matrix<T> &m,
int ni);
279 void mat_transp_reduction(
const tensor &t,
const gmm::dense_matrix<T> &m,
282 void mat_mult(
const gmm::dense_matrix<T> &m, gmm::dense_matrix<T> &mm);
285 void product(
const tensor &t2, tensor &tt);
287 void dot_product(
const tensor &t2, tensor &tt);
288 void dot_product(
const gmm::dense_matrix<T> &m, tensor &tt);
290 void double_dot_product(
const tensor &t2, tensor &tt);
291 void double_dot_product(
const gmm::dense_matrix<T> &m, tensor &tt);
294 return sizeof(T) * this->size()
295 +
sizeof(*this) + sizes_.memsize() + coeff_.memsize();
298 std::vector<T> &as_vector() {
return *
this; }
299 const std::vector<T> &as_vector()
const {
return *
this; }
302 tensor<T>& operator +=(
const tensor<T>& w)
303 { gmm::add(w.as_vector(), this->as_vector());
return *
this; }
305 tensor<T>& operator -=(
const tensor<T>& w) {
306 gmm::add(gmm::scaled(w.as_vector(), T(-1)), this->as_vector());
310 tensor<T>& operator *=(
const scalar_type w)
311 { gmm::scale(this->as_vector(), w);
return *
this; }
313 tensor<T>& operator /=(
const scalar_type w)
314 { gmm::scale(this->as_vector(), scalar_type(1)/w);
return *
this; }
316 tensor &operator =(
const tensor &t) {
317 if (this->size() != t.size()) this->
resize(t.size());
318 std::copy(t.begin(), t.end(), this->begin());
319 if (sizes_.size() != t.sizes_.size()) sizes_.resize(t.sizes_.size());
320 std::copy(t.sizes_.begin(), t.sizes_.end(), sizes_.begin());
321 if (coeff_.size() != t.coeff_.size()) coeff_.resize(t.coeff_.size());
322 std::copy(t.coeff_.begin(), t.coeff_.end(), coeff_.begin());
326 tensor(
const tensor &t)
327 : std::vector<T>(t), sizes_(t.sizes_), coeff_(t.coeff_) { }
328 tensor(
const multi_index &c) { init(c); }
333 { init(i, j, k, l); }
337 template<
class T>
void tensor<T>::mat_transp_reduction
338 (
const tensor &t,
const gmm::dense_matrix<T> &m,
int ni) {
341 THREAD_SAFE_STATIC std::vector<T> tmp;
342 THREAD_SAFE_STATIC multi_index mi;
345 size_type dimt = mi[ni], dim = m.nrows();
347 GMM_ASSERT2(dimt,
"Inconsistent dimension.");
348 GMM_ASSERT2(dimt == m.ncols(),
"Dimensions mismatch.");
349 GMM_ASSERT2(&t !=
this,
"Does not work when t and *this are the same.");
352 if (tmp.size() < dimt) tmp.resize(dimt);
355 const_iterator pft = t.begin();
356 iterator pf = this->begin();
357 size_type dd = coeff_[ni]*( sizes()[ni]-1)-1, co = coeff_[ni];
358 size_type ddt = t.coeff_[ni]*(t.sizes()[ni]-1)-1, cot = t.coeff_[ni];
359 std::fill(mi.begin(), mi.end(), 0);
360 for (;!mi.finished(sizes()); mi.incrementation(sizes()), ++pf, ++pft) {
364 pf += dd; pft += ddt;
366 const_iterator pl = pft; iterator pt = tmp.begin();
368 for(
size_type k = 1; k < dimt; ++k, ++pt) { pl += cot; *pt = *pl;}
373 *pff = T(0); pt = tmp.begin(); pl = m.begin() + k;
374 *pff += (*pl) * (*pt); ++pt;
375 for (
size_type l = 1; l < dimt; ++l, ++pt) {
377 *pff += (*pl) * (*pt);
384 template<
class T>
void tensor<T>::mat_mult(
const gmm::dense_matrix<T> &m,
385 gmm::dense_matrix<T> &mm) {
386 GMM_ASSERT2(order() == 4,
387 "This operation is for order four tensors only.");
388 GMM_ASSERT2(sizes_[2] == gmm::mat_nrows(m) &&
389 sizes_[3] == gmm::mat_ncols(m),
"Dimensions mismatch.");
390 mm.base_resize(sizes_[0], sizes_[1]);
393 const_iterator pt = this->begin();
394 const_iterator pm = m.begin();
395 for (
size_type l = 0; l < sizes_[3]; ++l)
396 for (
size_type k = 0; k < sizes_[2]; ++k) {
397 iterator pmm = mm.begin();
398 for (
size_type j = 0; j < sizes_[1]; ++j)
399 for (
size_type i = 0; i < sizes_[0]; ++i)
400 *pmm++ += *pt++ * (*pm);
405 template<
class T>
void tensor<T>::mat_reduction
406 (
const tensor &t,
const gmm::dense_matrix<T> &m,
int ni) {
408 THREAD_SAFE_STATIC std::vector<T> tmp;
409 THREAD_SAFE_STATIC multi_index mi;
412 size_type dimt = mi[ni], dim = m.ncols();
413 GMM_ASSERT2(dimt,
"Inconsistent dimension.");
414 GMM_ASSERT2(dimt == m.nrows(),
"Dimensions mismatch.");
415 GMM_ASSERT2(&t !=
this,
"Does not work when t and *this are the same.");
418 if (tmp.size() < dimt) tmp.resize(dimt);
420 const_iterator pft = t.begin();
421 iterator pf = this->begin();
422 size_type dd = coeff_[ni]*( sizes()[ni]-1)-1, co = coeff_[ni];
423 size_type ddt = t.coeff_[ni]*(t.sizes()[ni]-1)-1, cot = t.coeff_[ni];
424 std::fill(mi.begin(), mi.end(), 0);
425 for (;!mi.finished(sizes()); mi.incrementation(sizes()), ++pf, ++pft) {
429 pf += dd; pft += ddt;
432 const_iterator pl = pft; iterator pt = tmp.begin();
434 for(
size_type k = 1; k < dimt; ++k, ++pt) { pl += cot; *pt = *pl; }
436 iterator pff = pf; pl = m.begin();
439 *pff = T(0); pt = tmp.begin();
440 for (
size_type l = 0; l < dimt; ++l, ++pt, ++pl)
441 *pff += (*pl) * (*pt);
448 template<
class T>
void tensor<T>::product(
const tensor<T> &t2,
450 size_type res_order = order() + t2.order();
451 multi_index res_size(res_order);
452 for (
size_type i = 0 ; i < this->order(); ++i) res_size[i] = this->size(i);
453 for (
size_type i = 0 ; i < t2.order(); ++i) res_size[order() + i] = t2.size(i);
454 tt.adjust_sizes(res_size);
459 const_iterator pt2 = t2.begin();
460 iterator ptt = tt.begin();
461 for (
size_type j = 0; j < size2; ++j, ++pt2) {
462 const_iterator pt1 = this->begin();
463 for (
size_type i = 0; i < size1; ++i, ++pt1, ++ptt)
464 *ptt += *pt1 * (*pt2);
469 template<
class T>
void tensor<T>::dot_product(
const tensor<T> &t2,
471 GMM_ASSERT2(size(order()-1) == t2.size(0),
472 "Dimensions mismatch between last dimension of first tensor "
473 "and first dimension of second tensor.");
474 size_type res_order = order() + t2.order() - 2;
475 multi_index res_size(res_order);
476 for (
size_type i = 0 ; i < this->order() - 1; ++i) res_size[i] = this->size(i);
477 for (
size_type i = 0 ; i < t2.order() - 1; ++i) res_size[order() - 1 + i] = t2.size(i);
478 tt.adjust_sizes(res_size);
484 const_iterator pt2 = t2.begin();
485 iterator ptt = tt.begin();
487 const_iterator pt1 = this->begin();
489 for (
size_type q = 0; q < size0; ++q, ++pt2) {
491 for (
size_type i = 0; i < size1; ++i, ++pt1, ++ptt)
492 *ptt += *pt1 * (*pt2);
497 template<
class T>
void tensor<T>::dot_product(
const gmm::dense_matrix<T> &m,
499 GMM_ASSERT2(size(order()-1) == gmm::mat_nrows(m),
500 "Dimensions mismatch between last dimensions of tensor "
501 "and rows of the matrix.");
502 tensor<T> t2(multi_index(gmm::mat_nrows(m),gmm::mat_ncols(m)));
503 gmm::copy(m.as_vector(), t2.as_vector());
508 template<
class T>
void tensor<T>::double_dot_product(
const tensor<T> &t2,
510 GMM_ASSERT2(order() >= 2 && t2.order() >= 2,
511 "Tensors of wrong size. Tensors of order two or higher are required.");
512 GMM_ASSERT2(size(order()-2) == t2.size(0) && size(order()-1) == t2.size(1),
513 "Dimensions mismatch between last two dimensions of first tensor "
514 "and first two dimensions of second tensor.");
515 size_type res_order = order() + t2.order() - 4;
516 multi_index res_size(res_order);
517 for (
size_type i = 0 ; i < this->order() - 2; ++i) res_size[i] = this->size(i);
518 for (
size_type i = 0 ; i < t2.order() - 2; ++i) res_size[order() - 2 + i] = t2.size(i);
519 tt.adjust_sizes(res_size);
525 const_iterator pt2 = t2.begin();
526 iterator ptt = tt.begin();
528 const_iterator pt1 = this->begin();
530 for (
size_type q = 0; q < size0; ++q, ++pt2) {
532 for (
size_type i = 0; i < size1; ++i, ++pt1, ++ptt)
533 *ptt += *pt1 * (*pt2);
538 template<
class T>
void tensor<T>::double_dot_product(
const gmm::dense_matrix<T> &m,
540 GMM_ASSERT2(order() >= 2,
541 "Tensor of wrong size. Tensor of order two or higher is required.");
542 GMM_ASSERT2(size(order()-2) == gmm::mat_nrows(m) &&
543 size(order()-1) == gmm::mat_ncols(m),
544 "Dimensions mismatch between last two dimensions of tensor "
545 "and dimensions of the matrix.");
546 tensor<T> t2(multi_index(gmm::mat_nrows(m),gmm::mat_ncols(m)));
547 gmm::copy(m.as_vector(), t2.as_vector());
548 double_dot_product(t2, tt);
552 template<
class T> std::ostream &
operator <<
553 (std::ostream &o,
const tensor<T>& t) {
554 o <<
"sizes " << t.sizes() <<
" " << vref(t.as_vector());
558 typedef tensor<scalar_type> base_tensor;
559 typedef tensor<complex_type> base_complex_tensor;