26 #ifndef TILEDARRAY_CONVERSIONS_EIGEN_H__INCLUDED
27 #define TILEDARRAY_CONVERSIONS_EIGEN_H__INCLUDED
41 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
43 typedef Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
45 typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic,
48 typedef Eigen::Matrix<std::complex<float>, Eigen::Dynamic, Eigen::Dynamic,
51 typedef Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
53 typedef Eigen::Matrix<long, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
57 typedef Eigen::Matrix<std::complex<double>, 1, Eigen::Dynamic>
EigenVectorXcd;
58 typedef Eigen::Matrix<std::complex<float>, 1, Eigen::Dynamic>
EigenVectorXcf;
72 template <
typename T, Eigen::StorageOptions Storage = Eigen::RowMajor,
73 std::enable_if_t<detail::is_contiguous_tensor_v<T>>* =
nullptr>
74 inline Eigen::Map<
const Eigen::Matrix<
typename T::value_type, Eigen::Dynamic,
75 Eigen::Dynamic, Storage>,
77 eigen_map(
const T& tensor,
const std::size_t m,
const std::size_t n) {
80 return Eigen::Map<
const Eigen::Matrix<
typename T::value_type, Eigen::Dynamic,
81 Eigen::Dynamic, Storage>,
82 Eigen::AutoAlign>(tensor.data(), m, n);
95 template <
typename T, Eigen::StorageOptions Storage = Eigen::RowMajor,
96 std::enable_if_t<detail::is_contiguous_tensor_v<T>>* =
nullptr>
97 inline Eigen::Map<
Eigen::Matrix<
typename T::value_type, Eigen::Dynamic,
98 Eigen::Dynamic, Storage>,
100 eigen_map(T& tensor,
const std::size_t m,
const std::size_t n) {
103 return Eigen::Map<
Eigen::Matrix<
typename T::value_type, Eigen::Dynamic,
104 Eigen::Dynamic, Storage>,
105 Eigen::AutoAlign>(tensor.data(), m, n);
115 template <
typename T,
116 std::enable_if_t<detail::is_contiguous_tensor_v<T>>* =
nullptr>
118 const Eigen::Matrix<typename T::value_type, Eigen::Dynamic, 1>,
124 const Eigen::Matrix<typename T::value_type, Eigen::Dynamic, 1>,
125 Eigen::AutoAlign>(tensor.data(), n);
135 template <
typename T,
136 std::enable_if_t<detail::is_contiguous_tensor_v<T>>* =
nullptr>
137 inline Eigen::Map<Eigen::Matrix<typename T::value_type, Eigen::Dynamic, 1>,
142 return Eigen::Map<Eigen::Matrix<typename T::value_type, Eigen::Dynamic, 1>,
143 Eigen::AutoAlign>(tensor.data(), n);
156 template <
typename T, Eigen::StorageOptions Storage = Eigen::RowMajor,
157 std::enable_if_t<detail::is_contiguous_tensor_v<T>>* =
nullptr>
158 inline Eigen::Map<
const Eigen::Matrix<
typename T::value_type, Eigen::Dynamic,
159 Eigen::Dynamic, Storage>,
162 TA_ASSERT((tensor.range().rank() == 2u) || (tensor.range().rank() == 1u));
163 const auto* MADNESS_RESTRICT
const tensor_extent =
164 tensor.range().extent_data();
165 return eigen_map<T, Storage>(
166 tensor, tensor_extent[0],
167 (tensor.range().rank() == 2u ? tensor_extent[1] : 1ul));
179 template <
typename T, Eigen::StorageOptions Storage = Eigen::RowMajor,
180 std::enable_if_t<detail::is_contiguous_tensor_v<T>>* =
nullptr>
181 inline Eigen::Map<
Eigen::Matrix<
typename T::value_type, Eigen::Dynamic,
182 Eigen::Dynamic, Storage>,
185 TA_ASSERT((tensor.range().rank() == 2u) || (tensor.range().rank() == 1u));
186 const auto* MADNESS_RESTRICT
const tensor_extent =
187 tensor.range().extent_data();
188 return eigen_map<T, Storage>(
189 tensor, tensor_extent[0],
190 (tensor.range().rank() == 2u ? tensor_extent[1] : 1ul));
205 template <
typename T,
typename Derived,
206 std::enable_if_t<detail::is_contiguous_tensor_v<T>>* =
nullptr>
209 [[maybe_unused]]
typedef typename T::index1_type size_type;
210 TA_ASSERT((tensor.range().rank() == 2u) || (tensor.range().rank() == 1u));
213 const auto* MADNESS_RESTRICT
const tensor_lower =
214 tensor.range().lobound_data();
215 const auto* MADNESS_RESTRICT
const tensor_upper =
216 tensor.range().upbound_data();
217 const auto* MADNESS_RESTRICT
const tensor_extent =
218 tensor.range().extent_data();
220 if (tensor.range().rank() == 2u) {
222 const std::size_t tensor_lower_0 = tensor_lower[0];
223 const std::size_t tensor_lower_1 = tensor_lower[1];
224 [[maybe_unused]]
const std::size_t tensor_upper_0 = tensor_upper[0];
225 [[maybe_unused]]
const std::size_t tensor_upper_1 = tensor_upper[1];
226 const std::size_t tensor_extent_0 = tensor_extent[0];
227 const std::size_t tensor_extent_1 = tensor_extent[1];
229 TA_ASSERT(tensor_upper_0 <= std::size_t(matrix.rows()));
230 TA_ASSERT(tensor_upper_1 <= std::size_t(matrix.cols()));
233 eigen_map(tensor, tensor_extent_0, tensor_extent_1) = matrix.block(
234 tensor_lower_0, tensor_lower_1, tensor_extent_0, tensor_extent_1);
237 const std::size_t tensor_lower_0 = tensor_lower[0];
238 [[maybe_unused]]
const std::size_t tensor_upper_0 = tensor_upper[0];
239 const std::size_t tensor_extent_0 = tensor_extent[0];
242 TA_ASSERT((matrix.rows() == 1) || (matrix.cols() == 1));
244 if (matrix.rows() == 1) {
245 TA_ASSERT(tensor_upper_0 <= std::size_t(matrix.cols()));
249 matrix.block(0, tensor_lower_0, 1, tensor_extent_0);
251 TA_ASSERT(tensor_upper_0 <= std::size_t(matrix.rows()));
255 matrix.block(tensor_lower_0, 0, tensor_extent_0, 1);
272 template <
typename T,
typename Derived,
273 std::enable_if_t<detail::is_contiguous_tensor_v<T>>* =
nullptr>
275 Eigen::MatrixBase<Derived>& matrix) {
276 [[maybe_unused]]
typedef typename T::index1_type size_type;
277 TA_ASSERT((tensor.range().rank() == 2u) || (tensor.range().rank() == 1u));
280 const auto* MADNESS_RESTRICT
const tensor_lower =
281 tensor.range().lobound_data();
282 const auto* MADNESS_RESTRICT
const tensor_upper =
283 tensor.range().upbound_data();
284 const auto* MADNESS_RESTRICT
const tensor_extent =
285 tensor.range().extent_data();
287 if (tensor.range().rank() == 2) {
289 const std::size_t tensor_lower_0 = tensor_lower[0];
290 const std::size_t tensor_lower_1 = tensor_lower[1];
291 [[maybe_unused]]
const std::size_t tensor_upper_0 = tensor_upper[0];
292 [[maybe_unused]]
const std::size_t tensor_upper_1 = tensor_upper[1];
293 const std::size_t tensor_extent_0 = tensor_extent[0];
294 const std::size_t tensor_extent_1 = tensor_extent[1];
296 TA_ASSERT(tensor_upper_0 <= std::size_t(matrix.rows()));
297 TA_ASSERT(tensor_upper_1 <= std::size_t(matrix.cols()));
300 matrix.block(tensor_lower_0, tensor_lower_1, tensor_extent_0,
302 eigen_map(tensor, tensor_extent_0, tensor_extent_1);
305 const std::size_t tensor_lower_0 = tensor_lower[0];
306 [[maybe_unused]]
const std::size_t tensor_upper_0 = tensor_upper[0];
307 const std::size_t tensor_extent_0 = tensor_extent[0];
309 TA_ASSERT((matrix.rows() == 1) || (matrix.cols() == 1));
311 if (matrix.rows() == 1) {
312 TA_ASSERT(tensor_upper_0 <= std::size_t(matrix.cols()));
315 matrix.block(0, tensor_lower_0, 1, tensor_extent_0) =
318 TA_ASSERT(tensor_upper_0 <= std::size_t(matrix.rows()));
321 matrix.block(tensor_lower_0, 0, tensor_extent_0, 1) =
337 template <
typename A,
typename Derived>
340 const typename A::ordinal_type i,
341 madness::AtomicInt* counter) {
342 typename A::value_type tensor(array->trange().make_tile_range(i));
344 array->set(i, tensor);
355 template <
typename Derived,
typename T>
357 Eigen::MatrixBase<Derived>* matrix,
358 madness::AtomicInt* counter) {
408 template <
typename A,
typename Derived>
410 const Eigen::MatrixBase<Derived>& matrix,
411 bool replicated =
false,
412 std::shared_ptr<typename A::pmap_interface> pmap = {}) {
413 typedef typename A::index1_type size_type;
415 const auto rank = trange.tiles_range().rank();
418 "TiledArray::eigen_to_array(): The number of dimensions in "
419 "trange must match that of the Eigen matrix.");
423 trange.elements_range().extent(0) == size_type(matrix.rows()) &&
424 "TiledArray::eigen_to_array(): The number of rows in trange is not "
425 "equal to the number of rows in the Eigen matrix.");
427 trange.elements_range().extent(1) == size_type(matrix.cols()) &&
428 "TiledArray::eigen_to_array(): The number of columns in trange is not "
429 "equal to the number of columns in the Eigen matrix.");
432 trange.elements_range().extent(0) == size_type(matrix.size()) &&
433 "TiledArray::eigen_to_array(): The size of trange must be equal to the "
438 if (replicated && (world.size() > 1))
439 pmap = std::static_pointer_cast<typename A::pmap_interface>(
440 std::make_shared<detail::ReplicatedPmap>(
441 world, trange.tiles_range().volume()));
442 A array = (pmap ? A(world, trange, pmap) : A(world, trange));
445 madness::AtomicInt counter;
448 for (std::size_t i = 0; i < array.size(); ++i) {
449 if (array.is_local(i)) {
450 world.taskq.add(&detail::counted_eigen_submatrix_to_tensor<A, Derived>,
451 &matrix, &array, i, &counter);
457 array.world().await([&counter, n]() {
return counter == n; });
492 template <
typename Tile,
typename Policy,
493 unsigned int EigenStorageOrder = Eigen::ColMajor>
498 Eigen::Dynamic, EigenStorageOrder>
501 const auto rank = array.
trange().tiles_range().rank();
506 "TiledArray::array_to_eigen(): The array dimensions must be "
511 if (!array.
pmap()->is_replicated())
513 "TiledArray::array_to_eigen(): non-replicated Array cannot "
514 "be assigned to an Eigen::Matrix when the number of World "
515 "ranks is greater than 1.");
518 const auto* MADNESS_RESTRICT
const array_extent =
519 array.
trange().elements_range().extent_data();
522 EigenMatrix::Zero(array_extent[0], (
rank == 2 ? array_extent[1] : 1));
525 madness::AtomicInt counter;
528 for (std::size_t i = 0; i < array.
size(); ++i) {
530 array.
world().taskq.add(
533 array.
find(i), &matrix, &counter);
540 array.
world().await([&counter, n]() {
return counter == n; });
592 template <
typename A>
594 World& world,
const typename A::trange_type& trange,
595 const typename A::value_type::value_type* buffer,
const std::size_t m,
596 const std::size_t n,
const bool replicated =
false,
597 std::shared_ptr<typename A::pmap_interface> pmap = {}) {
598 TA_ASSERT(trange.elements_range().extent(0) == m &&
599 "TiledArray::eigen_to_array(): The number of rows in trange "
600 "is not equal to m.");
601 TA_ASSERT(trange.elements_range().extent(1) == n &&
602 "TiledArray::eigen_to_array(): The number of columns in "
603 "trange is not equal to n.");
605 typedef Eigen::Matrix<
typename A::value_type::value_type, Eigen::Dynamic,
606 Eigen::Dynamic, Eigen::RowMajor>
610 Eigen::Map<const matrix_type, Eigen::AutoAlign>(buffer, m, n), replicated,
661 template <
typename A>
663 World& world,
const typename A::trange_type& trange,
664 const typename A::value_type::value_type* buffer,
const std::size_t m,
665 const std::size_t n,
const bool replicated =
false,
666 std::shared_ptr<typename A::pmap_interface> pmap = {}) {
667 TA_ASSERT(trange.elements_range().extent(0) == m &&
668 "TiledArray::eigen_to_array(): The number of rows in trange "
669 "is not equal to m.");
670 TA_ASSERT(trange.elements_range().extent(1) == n &&
671 "TiledArray::eigen_to_array(): The number of columns in "
672 "trange is not equal to n.");
674 typedef Eigen::Matrix<
typename A::value_type::value_type, Eigen::Dynamic,
675 Eigen::Dynamic, Eigen::ColMajor>
679 Eigen::Map<const matrix_type, Eigen::AutoAlign>(buffer, m, n), replicated,
685 #endif // TILEDARRAY_CONVERSIONS_EIGEN_H__INCLUDED