1 #ifndef _MARRAY_DPD_MARRAY_BASE_HPP_ 2 #define _MARRAY_DPD_MARRAY_BASE_HPP_ 3 4 #include "marray_view.hpp" 5 #include "dpd_range.hpp" 6 7 namespace MArray 8 { 9 10 template <typename Type, unsigned NDim, typename Derived, bool Owner> 11 class dpd_marray_base; 12 13 template <typename Type, unsigned NDim> 14 class dpd_marray_view; 15 16 template <typename Type, unsigned NDim, typename Allocator=std::allocator<Type>> 17 class dpd_marray; 18 19 template <typename Type, typename Derived, bool Owner> 20 class dpd_varray_base; 21 22 template <typename Type> 23 class dpd_varray_view; 24 25 template <typename Type, typename Allocator=std::allocator<Type>> 26 class dpd_varray; 27 28 namespace detail 29 { 30 31 template <typename Derived> 32 struct dpd_base 33 { default_depthMArray::detail::dpd_base34 static dim_vector default_depth(dpd_layout layout, len_type ndim) 35 { 36 dim_vector depth(ndim); 37 38 if (layout == BALANCED_ROW_MAJOR || 39 layout == BALANCED_COLUMN_MAJOR) 40 { 41 unsigned dl = sizeof(unsigned)*8 - __builtin_clz(ndim) - 1; 42 unsigned du = dl + (ndim&(ndim-1) ? 1 : 0); 43 unsigned o = ndim - (1<<dl); 44 unsigned k = 2*((o+1)/2); 45 unsigned l = 2*(o/2); 46 unsigned h = (ndim+1)/2; 47 48 for (unsigned i = 0;i < k;i++) depth[i] = du; 49 for (unsigned i = k;i < h;i++) depth[i] = dl; 50 for (unsigned i = h;i < h+l;i++) depth[i] = du; 51 for (unsigned i = h+l;i < ndim;i++) depth[i] = dl; 52 } 53 else if (layout == PREFIX_ROW_MAJOR || 54 layout == PREFIX_COLUMN_MAJOR) 55 { 56 depth[0] = ndim-1; 57 for (unsigned i = 1;i < ndim;i++) depth[i] = ndim-i; 58 } 59 else //if (layout == BLOCKED_ROW_MAJOR || 60 // layout == BLOCKED_COLUMN_MAJOR) 61 { 62 for (unsigned i = 0;i < ndim-1;i++) depth[i] = i+1; 63 depth[ndim-1] = ndim-1; 64 } 65 66 return depth; 67 } 68 set_treeMArray::detail::dpd_base69 void set_tree() 70 { 71 auto& leaf = derived().leaf_; 72 auto& parent = derived().parent_; 73 auto nirrep = derived().nirrep_; 74 unsigned ndim = leaf.size(); 75 dim_vector depth(derived().depth_.begin(), derived().depth_.end()); 76 dim_vector node(ndim); 77 len_vector leaf_idx = range(ndim); 78 79 MARRAY_ASSERT(depth.size() == ndim); 80 for (unsigned i = 0;i < depth.size();i++) 81 MARRAY_ASSERT(depth[i] < ndim); 82 83 unsigned pos = 0; 84 for (unsigned d = ndim;d --> 0;) 85 { 86 for (unsigned i = 0;i < depth.size();i++) 87 { 88 /* 89 * If we encounter a pair of nodes with depth greater than the 90 * current level, combine them into an interior node and 91 * assign the parent pointers. 92 */ 93 if (depth[i] == d+1) 94 { 95 MARRAY_ASSERT(i < depth.size()-1 && depth[i+1] == d+1); 96 parent[node[i]] = parent[node[i+1]] = pos; 97 depth.erase(depth.begin()+i+1); 98 depth[i]--; 99 node.erase(node.begin()+i+1); 100 node[i] = pos++; 101 leaf_idx.erase(leaf_idx.begin()+i+1); 102 leaf_idx[i] = -1; 103 } 104 /* 105 * For nodes on the current depth level, assign a position 106 * in the linearized tree. 107 */ 108 else if (depth[i] == d) 109 { 110 node[i] = pos++; 111 if (leaf_idx[i] != -1) leaf[leaf_idx[i]] = node[i]; 112 } 113 } 114 } 115 116 MARRAY_ASSERT(pos == 2*ndim-1); 117 MARRAY_ASSERT(depth.size() == 1); 118 MARRAY_ASSERT(depth[0] == 0); 119 } 120 set_sizeMArray::detail::dpd_base121 void set_size() 122 { 123 auto& perm = derived().perm_; 124 auto& size = derived().size_; 125 auto& len = derived().len_; 126 const auto& leaf = derived().leaf_; 127 const auto& parent = derived().parent_; 128 auto nirrep = derived().nirrep_; 129 auto layout = derived().layout_; 130 unsigned ndim = perm.size(); 131 132 if (layout == COLUMN_MAJOR) 133 { 134 // Column major 135 for (unsigned i = 0;i < ndim;i++) 136 { 137 std::copy_n(&len[i][0], nirrep, &size[leaf[i]][0]); 138 perm[i] = i; 139 } 140 } 141 else 142 { 143 // Row major: reverse the dimensions and treat as 144 // permuted column major 145 for (unsigned i = 0;i < ndim;i++) 146 { 147 std::copy_n(&len[i][0], nirrep, &size[leaf[ndim-1-i]][0]); 148 perm[i] = ndim-1-i; 149 } 150 151 for (unsigned i = 0;i < ndim/2;i++) 152 for (unsigned j = 0;j < nirrep;j++) 153 std::swap(len[i][j], len[ndim-1-i][j]); 154 } 155 156 for (unsigned i = 0;i < ndim-1;i++) 157 { 158 unsigned next = parent[2*i]; 159 160 for (unsigned irr1 = 0;irr1 < nirrep;irr1++) 161 { 162 size[next][irr1] = 0; 163 for (unsigned irr2 = 0;irr2 < nirrep;irr2++) 164 { 165 size[next][irr1] += size[2*i][irr1^irr2]*size[2*i+1][irr2]; 166 } 167 } 168 } 169 } 170 171 template <typename U, typename V, typename W, typename X> get_blockMArray::detail::dpd_base172 void get_block(const U& irreps, V& len, W& data, X& stride) const 173 { 174 auto& perm = derived().perm_; 175 auto& size = derived().size_; 176 auto& extent = derived().len_; 177 auto& off = derived().off_; 178 auto& leap = derived().stride_; 179 auto& leaf = derived().leaf_; 180 auto& parent = derived().parent_; 181 unsigned ndim = perm.size(); 182 183 short_vector<unsigned, 2*MARRAY_OPT_NDIM-1> dpd_irrep(2*ndim-1); 184 short_vector<stride_type, 2*MARRAY_OPT_NDIM-1> dpd_stride(2*ndim-1); 185 dpd_stride[2*ndim-2] = 1; 186 187 auto it = irreps.begin(); 188 for (unsigned i = 0;i < ndim;i++) 189 { 190 dpd_irrep[leaf[perm[i]]] = *it; 191 ++it; 192 } 193 194 for (unsigned i = 0;i < ndim-1;i++) 195 { 196 dpd_irrep[parent[2*i]] = dpd_irrep[2*i]^dpd_irrep[2*i+1]; 197 } 198 199 for (unsigned i = ndim-1;i --> 0;) 200 { 201 unsigned irrep = dpd_irrep[parent[2*i]]; 202 203 dpd_stride[2*i] = dpd_stride[parent[2*i]]; 204 dpd_stride[2*i+1] = dpd_stride[2*i]*size[2*i][dpd_irrep[2*i]]; 205 206 stride_type offset = 0; 207 for (unsigned irr1 = 0;irr1 < dpd_irrep[2*i+1];irr1++) 208 { 209 offset += size[2*i][irr1^irrep]*size[2*i+1][irr1]; 210 } 211 212 data += offset*dpd_stride[2*i]; 213 } 214 215 it = irreps.begin(); 216 auto l = len.begin(); 217 auto s = stride.begin(); 218 for (unsigned i = 0;i < ndim;i++) 219 { 220 auto stride = dpd_stride[leaf[perm[i]]]*leap[perm[i]][dpd_irrep[leaf[perm[i]]]]; 221 *s = stride; 222 *l = extent[perm[i]][*it]; 223 data += stride*off[perm[i]][*it]; 224 ++it; 225 ++l; 226 ++s; 227 } 228 } 229 derivedMArray::detail::dpd_base230 const Derived& derived() const { return static_cast<const Derived&>(*this); } 231 derivedMArray::detail::dpd_base232 Derived& derived() { return static_cast<Derived&>(*this); } 233 }; 234 235 } 236 237 template <typename Type, unsigned NDim, typename Derived, bool Owner> 238 class dpd_marray_base : protected detail::dpd_base<dpd_marray_base<Type, NDim, Derived, Owner>> 239 { 240 static_assert(NDim > 0, "NDim must be positive"); 241 242 template <typename> friend struct detail::dpd_base; 243 template <typename, unsigned, typename, bool> friend class dpd_marray_base; 244 template <typename, unsigned> friend class dpd_marray_view; 245 template <typename, unsigned, typename> friend class dpd_marray; 246 247 public: 248 typedef Type value_type; 249 typedef Type* pointer; 250 typedef const Type* const_pointer; 251 typedef Type& reference; 252 typedef const Type& const_reference; 253 254 typedef typename std::conditional<Owner,const Type,Type>::type ctype; 255 typedef ctype& cref; 256 typedef ctype* cptr; 257 258 protected: 259 std::array<std::array<stride_type,8>, 2*NDim-1> size_ = {}; 260 std::array<std::array<len_type,8>, NDim> len_ = {}; 261 std::array<std::array<len_type,8>, NDim> off_ = {}; 262 std::array<std::array<stride_type,8>, NDim> stride_ = {}; 263 std::array<unsigned, NDim> leaf_ = {}; 264 std::array<unsigned, 2*NDim-1> parent_ = {}; 265 std::array<unsigned, NDim> perm_ = {}; 266 std::array<unsigned, NDim> depth_ = {}; 267 pointer data_ = nullptr; 268 unsigned irrep_ = 0; 269 unsigned nirrep_ = 0; 270 layout layout_ = DEFAULT; 271 272 /*********************************************************************** 273 * 274 * Reset 275 * 276 **********************************************************************/ 277 reset()278 void reset() 279 { 280 size_ = {}; 281 len_ = {}; 282 off_ = {}; 283 stride_ = {}; 284 leaf_ = {}; 285 parent_ = {}; 286 perm_ = {}; 287 depth_ = {}; 288 data_ = nullptr; 289 irrep_ = 0; 290 nirrep_ = 0; 291 layout_ = DEFAULT; 292 } 293 294 template <typename U, bool O, typename D, 295 typename=detail::enable_if_convertible_t< 296 typename dpd_marray_base<U, NDim, D, O>::cptr,pointer>> reset(const dpd_marray_base<U,NDim,D,O> & other)297 void reset(const dpd_marray_base<U, NDim, D, O>& other) 298 { 299 reset(const_cast<dpd_marray_base<U, NDim, D, O>&>(other)); 300 } 301 302 template <typename U, bool O, typename D, 303 typename=detail::enable_if_convertible_t< 304 typename dpd_marray_base<U, NDim, D, O>::pointer,pointer>> reset(dpd_marray_base<U,NDim,D,O> & other)305 void reset(dpd_marray_base<U, NDim, D, O>& other) 306 { 307 size_ = other.size_; 308 len_ = other.len_; 309 off_ = other.off_; 310 stride_ = other.stride_; 311 leaf_ = other.leaf_; 312 parent_ = other.parent_; 313 perm_ = other.perm_; 314 depth_ = other.depth_; 315 data_ = other.data_; 316 irrep_ = other.irrep_; 317 nirrep_ = other.nirrep_; 318 layout_ = other.layout_; 319 } 320 reset(unsigned irrep,unsigned nirrep,const detail::array_2d<len_type> & len,pointer ptr,dpd_layout layout=DEFAULT)321 void reset(unsigned irrep, unsigned nirrep, 322 const detail::array_2d<len_type>& len, pointer ptr, 323 dpd_layout layout = DEFAULT) 324 { 325 reset(irrep, nirrep, len, ptr, 326 this->default_depth(layout, NDim), layout.base()); 327 } 328 reset(unsigned irrep,unsigned nirrep,const detail::array_2d<len_type> & len,pointer ptr,const detail::array_1d<unsigned> & depth,layout layout=DEFAULT)329 void reset(unsigned irrep, unsigned nirrep, 330 const detail::array_2d<len_type>& len, pointer ptr, 331 const detail::array_1d<unsigned>& depth, layout layout = DEFAULT) 332 { 333 MARRAY_ASSERT(nirrep == 1 || nirrep == 2 || 334 nirrep == 4 || nirrep == 8); 335 336 MARRAY_ASSERT(len.length(0) == NDim); 337 MARRAY_ASSERT(len.length(1) >= nirrep); 338 MARRAY_ASSERT(depth.size() == NDim); 339 340 data_ = ptr; 341 irrep_ = irrep; 342 nirrep_ = nirrep; 343 layout_ = layout; 344 depth.slurp(depth_); 345 len.slurp(len_); 346 std::fill_n(&stride_[0][0], NDim*8, 1); 347 348 this->set_tree(); 349 this->set_size(); 350 } 351 352 /*********************************************************************** 353 * 354 * Private helper functions 355 * 356 **********************************************************************/ 357 358 template <typename View, typename Func, unsigned... I> for_each_block(Func && f,detail::integer_sequence<unsigned,I...>) const359 void for_each_block(Func&& f, detail::integer_sequence<unsigned, I...>) const 360 { 361 typedef typename View::pointer Ptr; 362 363 std::array<unsigned, NDim-1> nirrep; 364 nirrep.fill(nirrep_); 365 366 const_pointer cptr; 367 std::array<unsigned, NDim> irreps; 368 std::array<len_type, NDim> len; 369 std::array<stride_type, NDim> stride; 370 371 miterator<NDim-1, 0> it(nirrep); 372 while (it.next()) 373 { 374 irreps[0] = irrep_; 375 for (unsigned i = 1;i < NDim;i++) 376 { 377 irreps[0] ^= irreps[i] = it.position()[i-1]; 378 } 379 380 bool empty = false; 381 for (unsigned i = 0;i < NDim;i++) 382 { 383 if (length(i, irreps[i]) == 0) empty = true; 384 } 385 if (empty) continue; 386 387 cptr = data(); 388 this->get_block(irreps, len, cptr, stride); 389 390 detail::call(std::forward<Func>(f), 391 View(len, const_cast<Ptr>(cptr), stride), 392 irreps[I]...); 393 } 394 } 395 396 template <typename Tp, typename Func, unsigned... I> for_each_element(Func && f,detail::integer_sequence<unsigned,I...>) const397 void for_each_element(Func&& f, detail::integer_sequence<unsigned, I...>) const 398 { 399 typedef Tp* Ptr; 400 401 std::array<unsigned, NDim-1> nirrep; 402 nirrep.fill(nirrep_); 403 404 const_pointer cptr; 405 std::array<unsigned, NDim> irreps; 406 std::array<len_type, NDim> len; 407 std::array<stride_type, NDim> stride; 408 409 miterator<NDim-1, 0> it1(nirrep); 410 while (it1.next()) 411 { 412 irreps[0] = irrep_; 413 for (unsigned i = 1;i < NDim;i++) 414 { 415 irreps[0] ^= irreps[i] = it1.position()[i-1]; 416 } 417 418 bool empty = false; 419 for (unsigned i = 0;i < NDim;i++) 420 { 421 if (length(i, irreps[i]) == 0) empty = true; 422 } 423 if (empty) continue; 424 425 cptr = data(); 426 this->get_block(irreps, len, cptr, stride); 427 428 miterator<NDim, 1> it2(len, stride); 429 Ptr ptr = const_cast<Ptr>(cptr); 430 while (it2.next(ptr)) detail::call(std::forward<Func>(f), *ptr, 431 irreps[I]..., it2.position()[I]...); 432 } 433 } 434 swap(dpd_marray_base & other)435 void swap(dpd_marray_base& other) 436 { 437 using std::swap; 438 swap(size_, other.size_); 439 swap(len_, other.len_); 440 swap(off_, other.off_); 441 swap(stride_, other.stride_); 442 swap(leaf_, other.leaf_); 443 swap(parent_, other.parent_); 444 swap(perm_, other.perm_); 445 swap(depth_, other.depth_); 446 swap(data_, other.data_); 447 swap(irrep_, other.irrep_); 448 swap(nirrep_, other.nirrep_); 449 swap(layout_, other.layout_); 450 } 451 452 public: 453 454 /*********************************************************************** 455 * 456 * Static helper functions 457 * 458 **********************************************************************/ 459 size(unsigned irrep,const detail::array_2d<len_type> & len_)460 static stride_type size(unsigned irrep, const detail::array_2d<len_type>& len_) 461 { 462 if (len_.length(0) == 0) return 1; 463 464 //TODO: add alignment option 465 466 matrix<len_type> len; 467 len_.slurp(len); 468 469 unsigned nirrep = len.length(1); 470 unsigned ndim = len.length(0); 471 472 MARRAY_ASSERT(nirrep == 1 || nirrep == 2 || 473 nirrep == 4 || nirrep == 8); 474 475 std::array<stride_type, 8> size; 476 477 for (unsigned irr = 0;irr < nirrep;irr++) 478 size[irr] = len[0][irr]; 479 480 for (unsigned i = 1;i < ndim;i++) 481 { 482 std::array<stride_type, 8> new_size = {}; 483 484 for (unsigned irr1 = 0;irr1 < nirrep;irr1++) 485 { 486 for (unsigned irr2 = 0;irr2 < nirrep;irr2++) 487 { 488 new_size[irr1] += size[irr1^irr2]*len[i][irr2]; 489 } 490 } 491 492 size = new_size; 493 } 494 495 return size[irrep]; 496 } 497 498 /*********************************************************************** 499 * 500 * Operators 501 * 502 **********************************************************************/ 503 operator =(const dpd_marray_base & other)504 Derived& operator=(const dpd_marray_base& other) 505 { 506 return operator=<>(other); 507 } 508 509 template <typename U, typename D, bool O, 510 typename=detail::enable_if_t<std::is_assignable<reference,U>::value>> operator =(const dpd_marray_base<U,NDim,D,O> & other)511 Derived& operator=(const dpd_marray_base<U, NDim, D, O>& other) 512 { 513 MARRAY_ASSERT(nirrep_ == other.nirrep_); 514 MARRAY_ASSERT(irrep_ == other.irrep_); 515 516 for (unsigned i = 0;i < NDim;i++) 517 { 518 MARRAY_ASSERT(lengths(i) == other.lengths(i)); 519 } 520 521 unsigned mask = nirrep_-1; 522 unsigned shift = (nirrep_>1) + (nirrep_>2) + (nirrep_>4); 523 524 unsigned nblocks = 1u << (shift*(NDim-1)); 525 std::array<unsigned, NDim> irreps; 526 for (unsigned block = 0;block < nblocks;block++) 527 { 528 unsigned b = block; 529 irreps[0] = irrep_; 530 for (unsigned i = 1;i < NDim;i++) 531 { 532 irreps[0] ^= irreps[i] = b & mask; 533 b >>= shift; 534 } 535 536 (*this)(irreps) = other(irreps); 537 } 538 539 return static_cast<Derived&>(*this); 540 } 541 operator =(const Type & value)542 Derived& operator=(const Type& value) 543 { 544 unsigned mask = nirrep_-1; 545 unsigned shift = (nirrep_>1) + (nirrep_>2) + (nirrep_>4); 546 547 unsigned nblocks = 1u << (shift*(NDim-1)); 548 std::array<unsigned, NDim> irreps; 549 for (unsigned block = 0;block < nblocks;block++) 550 { 551 unsigned b = block; 552 irreps[0] = irrep_; 553 for (unsigned i = 1;i < NDim;i++) 554 { 555 irreps[0] ^= irreps[i] = b & mask; 556 b >>= shift; 557 } 558 559 (*this)(irreps) = value; 560 } 561 562 return static_cast<Derived&>(*this); 563 } 564 565 /*********************************************************************** 566 * 567 * Views 568 * 569 **********************************************************************/ 570 cview() const571 dpd_marray_view<const Type, NDim> cview() const 572 { 573 return const_cast<dpd_marray_base&>(*this).view(); 574 } 575 view() const576 dpd_marray_view<ctype, NDim> view() const 577 { 578 return const_cast<dpd_marray_base&>(*this).view(); 579 } 580 view()581 dpd_marray_view<Type, NDim> view() 582 { 583 return *this; 584 } 585 cview(const dpd_marray_base & x)586 friend dpd_marray_view<const Type, NDim> cview(const dpd_marray_base& x) 587 { 588 return x.view(); 589 } 590 view(const dpd_marray_base & x)591 friend dpd_marray_view<ctype, NDim> view(const dpd_marray_base& x) 592 { 593 return x.view(); 594 } 595 view(dpd_marray_base & x)596 friend dpd_marray_view<Type, NDim> view(dpd_marray_base& x) 597 { 598 return x.view(); 599 } 600 601 /*********************************************************************** 602 * 603 * Permutation 604 * 605 **********************************************************************/ 606 permuted(const detail::array_1d<unsigned> & perm) const607 dpd_marray_view<ctype,NDim> permuted(const detail::array_1d<unsigned>& perm) const 608 { 609 return const_cast<dpd_marray_base&>(*this).permuted(perm); 610 } 611 permuted(const detail::array_1d<unsigned> & perm)612 dpd_marray_view<Type,NDim> permuted(const detail::array_1d<unsigned>& perm) 613 { 614 dpd_marray_view<Type,NDim> r(*this); 615 r.permute(perm); 616 return r; 617 } 618 619 template <unsigned N=NDim, typename=detail::enable_if_t<N==2>> transposed() const620 dpd_marray_view<ctype, NDim> transposed() const 621 { 622 return const_cast<dpd_marray_base&>(*this).transposed(); 623 } 624 625 template <unsigned N=NDim, typename=detail::enable_if_t<N==2>> transposed()626 dpd_marray_view<Type, NDim> transposed() 627 { 628 return permuted({1, 0}); 629 } 630 631 template <unsigned N=NDim, typename=detail::enable_if_t<N==2>> T() const632 dpd_marray_view<ctype, NDim> T() const 633 { 634 return const_cast<dpd_marray_base&>(*this).T(); 635 } 636 637 template <unsigned N=NDim, typename=detail::enable_if_t<N==2>> T()638 dpd_marray_view<Type, NDim> T() 639 { 640 return transposed(); 641 } 642 643 /*********************************************************************** 644 * 645 * Indexing 646 * 647 **********************************************************************/ 648 649 template <typename... Irreps> 650 detail::enable_if_t<detail::are_assignable<unsigned&, Irreps...>::value && 651 sizeof...(Irreps) == NDim, marray_view<ctype, NDim>> operator ()(const Irreps &...irreps) const652 operator()(const Irreps&... irreps) const 653 { 654 return const_cast<dpd_marray_base&>(*this)(irreps...); 655 } 656 657 template <typename... Irreps> 658 detail::enable_if_t<detail::are_assignable<unsigned&, Irreps...>::value && 659 sizeof...(Irreps) == NDim, marray_view<Type, NDim>> operator ()(const Irreps &...irreps)660 operator()(const Irreps&... irreps) 661 { 662 return const_cast<dpd_marray_base&>(*this)({(unsigned)irreps...}); 663 } 664 operator ()(const detail::array_1d<unsigned> & irreps) const665 marray_view<ctype, NDim> operator()(const detail::array_1d<unsigned>& irreps) const 666 { 667 return const_cast<dpd_marray_base&>(*this)(irreps); 668 } 669 operator ()(const detail::array_1d<unsigned> & irreps_)670 marray_view<Type, NDim> operator()(const detail::array_1d<unsigned>& irreps_) 671 { 672 MARRAY_ASSERT(irreps_.size() == NDim); 673 674 std::array<unsigned, NDim> irreps; 675 std::array<len_type, NDim> len; 676 std::array<stride_type, NDim> stride; 677 678 irreps_.slurp(irreps); 679 680 unsigned irrep = 0; 681 for (auto& i: irreps) irrep ^= i; 682 MARRAY_ASSERT(irrep == irrep_); 683 684 pointer ptr = data(); 685 this->get_block(irreps, len, ptr, stride); 686 687 return marray_view<Type, NDim>(len, ptr, stride); 688 } 689 690 /*********************************************************************** 691 * 692 * Slicing 693 * 694 **********************************************************************/ 695 696 template <typename... Slices> 697 detail::enable_if_t<detail::are_dpd_indices_or_slices<Slices...>::value && 698 (detail::sliced_dimension<Slices...>::value > 0) && 699 (detail::sliced_dimension<Slices...>::value < NDim), 700 dpd_marray_view<ctype, detail::sliced_dimension<Slices...>::value>> operator ()(const Slices &...slices) const701 operator()(const Slices&... slices) const 702 { 703 return const_cast<dpd_marray_base&>(*this)(slices...); 704 } 705 706 template <typename... Slices> 707 detail::enable_if_t<detail::are_dpd_indices_or_slices<Slices...>::value && 708 (detail::sliced_dimension<Slices...>::value > 0) && 709 (detail::sliced_dimension<Slices...>::value < NDim), 710 dpd_marray_view<Type, detail::sliced_dimension<Slices...>::value>> operator ()(const Slices &...slices)711 operator()(const Slices&... slices) 712 { 713 constexpr unsigned NDim2 = detail::sliced_dimension<Slices...>::value; 714 715 abort(); 716 //TODO 717 } 718 719 template <typename... Slices> 720 detail::enable_if_t<detail::are_assignable<dpd_range_t&, Slices...>::value && 721 sizeof...(Slices) == NDim, dpd_marray_view<ctype, NDim>> operator ()(const Slices &...slices) const722 operator()(const Slices&... slices) const 723 { 724 return const_cast<dpd_marray_base&>(*this)(slices...); 725 } 726 727 template <typename... Slices> 728 detail::enable_if_t<detail::are_assignable<dpd_range_t&, Slices...>::value && 729 sizeof...(Slices) == NDim, dpd_marray_view<Type, NDim>> operator ()(const Slices &...slices)730 operator()(const Slices&... slices) 731 { 732 return const_cast<dpd_marray_base&>(*this)({slices...}); 733 } 734 operator ()(const detail::array_1d<dpd_range_t> & slices) const735 dpd_marray_view<ctype, NDim> operator()(const detail::array_1d<dpd_range_t>& slices) const 736 { 737 return const_cast<dpd_marray_base&>(*this)(slices); 738 } 739 operator ()(const detail::array_1d<dpd_range_t> & slices_)740 dpd_marray_view<Type, NDim> operator()(const detail::array_1d<dpd_range_t>& slices_) 741 { 742 MARRAY_ASSERT(slices_.size() == NDim); 743 744 std::array<dpd_range_t, NDim> slices; 745 slices_.slurp(slices); 746 747 dpd_marray_view<Type, NDim> v = view(); 748 749 for (unsigned i = 0;i < NDim;i++) 750 { 751 for (unsigned j = 0;j < nirrep_;j++) 752 { 753 v.off_[i][j] = slices[i].from(j); 754 v.len_[i][j] = slices[i].size(j); 755 } 756 } 757 758 return v; 759 } 760 761 /*********************************************************************** 762 * 763 * Iteration 764 * 765 **********************************************************************/ 766 767 template <typename Func> for_each_block(Func && f) const768 void for_each_block(Func&& f) const 769 { 770 for_each_block<marray_view<ctype, NDim>>(std::forward<Func>(f), detail::static_range<unsigned, NDim>{}); 771 } 772 773 template <typename Func> for_each_block(Func && f)774 void for_each_block(Func&& f) 775 { 776 for_each_block<marray_view<Type, NDim>>(std::forward<Func>(f), detail::static_range<unsigned, NDim>{}); 777 } 778 779 template <typename Func> for_each_element(Func && f) const780 void for_each_element(Func&& f) const 781 { 782 for_each_element<ctype>(std::forward<Func>(f), detail::static_range<unsigned, NDim>{}); 783 } 784 785 template <typename Func> for_each_element(Func && f)786 void for_each_element(Func&& f) 787 { 788 for_each_element<Type>(std::forward<Func>(f), detail::static_range<unsigned, NDim>{}); 789 } 790 791 /*********************************************************************** 792 * 793 * Basic getters 794 * 795 **********************************************************************/ 796 cdata() const797 const_pointer cdata() const 798 { 799 return const_cast<dpd_marray_base&>(*this).data(); 800 } 801 data() const802 cptr data() const 803 { 804 return const_cast<dpd_marray_base&>(*this).data(); 805 } 806 data()807 pointer data() 808 { 809 return data_; 810 } 811 812 template <unsigned Dim> length(unsigned irrep) const813 len_type length(unsigned irrep) const 814 { 815 static_assert(Dim < NDim, "Dim out of range"); 816 return length(Dim, irrep); 817 } 818 length(unsigned dim,unsigned irrep) const819 len_type length(unsigned dim, unsigned irrep) const 820 { 821 MARRAY_ASSERT(irrep < nirrep_); 822 return lengths(dim)[irrep]; 823 } 824 825 template <unsigned Dim> lengths() const826 const std::array<len_type,8>& lengths() const 827 { 828 static_assert(Dim < NDim, "Dim out of range"); 829 return lengths(Dim); 830 } 831 lengths(unsigned dim) const832 const std::array<len_type,8>& lengths(unsigned dim) const 833 { 834 MARRAY_ASSERT(dim < NDim); 835 return len_[perm_[dim]]; 836 } 837 lengths() const838 std::array<std::array<len_type,8>, NDim> lengths() const 839 { 840 std::array<std::array<len_type,8>, NDim> len = {}; 841 for (unsigned i = 0;i < NDim;i++) len[i] = lengths(i); 842 return len; 843 } 844 irrep() const845 unsigned irrep() const 846 { 847 return irrep_; 848 } 849 num_irreps() const850 unsigned num_irreps() const 851 { 852 return nirrep_; 853 } 854 permutation() const855 const std::array<unsigned, NDim>& permutation() const 856 { 857 return perm_; 858 } 859 dimension()860 static constexpr unsigned dimension() 861 { 862 return NDim; 863 } 864 }; 865 866 } 867 868 #endif 869