heig.h
Go to the documentation of this file.
1 /*
2  * This file is a part of TiledArray.
3  * Copyright (C) 2020 Virginia Tech
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * David Williams-Young
19  * Computational Research Division, Lawrence Berkeley National Laboratory
20  *
21  * heig.h
22  * Created: 13 May, 2020
23  * Edited: 8 June, 2020
24  *
25  */
26 #ifndef TILEDARRAY_MATH_LINALG_SCALAPACK_HEIG_H__INCLUDED
27 #define TILEDARRAY_MATH_LINALG_SCALAPACK_HEIG_H__INCLUDED
28 
29 #include <TiledArray/config.h>
30 #if TILEDARRAY_HAS_SCALAPACK
31 
33 
34 #include <scalapackpp/eigenvalue_problem/gevp.hpp>
35 #include <scalapackpp/eigenvalue_problem/sevp.hpp>
36 
38 
58 template <typename Array>
59 auto heig(const Array& A, TiledRange evec_trange = TiledRange(),
60  size_t NB = default_block_size()) {
61  using value_type = typename Array::element_type;
62  using real_type = scalapackpp::detail::real_t<value_type>;
63 
64  auto& world = A.world();
65  auto world_comm = world.mpi.comm().Get_mpi_comm();
66  // auto world_comm = MPI_COMM_WORLD;
67  blacspp::Grid grid = blacspp::Grid::square_grid(world_comm);
68 
69  world.gop.fence(); // stage ScaLAPACK execution
70  auto matrix = scalapack::array_to_block_cyclic(A, grid, NB, NB);
71  world.gop.fence(); // stage ScaLAPACK execution
72 
73  auto [M, N] = matrix.dims();
74  if (M != N) TA_EXCEPTION("Matrix must be square for EVP");
75 
76  auto [Mloc, Nloc] = matrix.dist().get_local_dims(N, N);
77  auto desc = matrix.dist().descinit_noerror(N, N, Mloc);
78 
79  std::vector<real_type> evals(N);
80  scalapack::BlockCyclicMatrix<value_type> evecs(world, grid, N, N, NB, NB);
81 
82  auto info = scalapackpp::hereig(
83  scalapackpp::VectorFlag::Vectors, blacspp::Triangle::Lower, N,
84  matrix.local_mat().data(), 1, 1, desc, evals.data(),
85  evecs.local_mat().data(), 1, 1, desc);
86  if (info) TA_EXCEPTION("EVP Failed");
87 
88  if (evec_trange.rank() == 0) evec_trange = A.trange();
89 
90  world.gop.fence();
91  auto evecs_ta = scalapack::block_cyclic_to_array<Array>(evecs, evec_trange);
92  world.gop.fence();
93 
94  return std::tuple(evals, evecs_ta);
95 }
96 
121 template <typename ArrayA, typename ArrayB, typename EVecType = ArrayA>
122 auto heig(const ArrayA& A, const ArrayB& B,
123  TiledRange evec_trange = TiledRange(),
124  size_t NB = default_block_size()) {
125  using value_type = typename ArrayA::element_type;
126  static_assert(std::is_same_v<typename ArrayB::element_type, value_type>);
127  using real_type = scalapackpp::detail::real_t<value_type>;
128 
129  auto& world = A.world();
130  auto world_comm = world.mpi.comm().Get_mpi_comm();
131  // auto world_comm = MPI_COMM_WORLD;
132  blacspp::Grid grid = blacspp::Grid::square_grid(world_comm);
133 
134  world.gop.fence(); // stage ScaLAPACK execution
135  auto A_sca = scalapack::array_to_block_cyclic(A, grid, NB, NB);
136  auto B_sca = scalapack::array_to_block_cyclic(B, grid, NB, NB);
137  world.gop.fence(); // stage ScaLAPACK execution
138 
139  auto [M, N] = A_sca.dims();
140  if (M != N) TA_EXCEPTION("Matrix must be square for EVP");
141 
142  auto [B_M, B_N] = B_sca.dims();
143  if (B_M != M or B_N != N)
144  TA_EXCEPTION("A and B must have the same dimensions");
145 
146  auto [Mloc, Nloc] = A_sca.dist().get_local_dims(N, N);
147  auto desc = A_sca.dist().descinit_noerror(N, N, Mloc);
148 
149  std::vector<real_type> evals(N);
150  scalapack::BlockCyclicMatrix<value_type> evecs(world, grid, N, N, NB, NB);
151 
152  auto info = scalapackpp::hereig_gen(
153  scalapackpp::VectorFlag::Vectors, blacspp::Triangle::Lower, N,
154  A_sca.local_mat().data(), 1, 1, desc, B_sca.local_mat().data(), 1, 1,
155  desc, evals.data(), evecs.local_mat().data(), 1, 1, desc);
156  if (info) TA_EXCEPTION("EVP Failed");
157 
158  if (evec_trange.rank() == 0) evec_trange = A.trange();
159 
160  world.gop.fence();
161  auto evecs_ta =
162  scalapack::block_cyclic_to_array<EVecType>(evecs, evec_trange);
163  world.gop.fence();
164 
165  return std::tuple(evals, evecs_ta);
166 }
167 
168 } // namespace TiledArray::math::linalg::scalapack
169 
170 #endif // TILEDARRAY_HAS_SCALAPACK
171 #endif // TILEDARRAY_MATH_LINALG_SCALAPACK_HEIG_H__INCLUDED
auto heig(const Array &A, TiledRange evec_trange=TiledRange(), size_t NB=default_block_size())
Solve the standard eigenvalue problem with ScaLAPACK.
Definition: heig.h:59
std::size_t default_block_size()
Definition: util.h:88
#define TA_EXCEPTION(m)
Definition: error.h:83
Range data of a tiled array.
Definition: tiled_range.h:32
const trange_type & trange() const
Tiled range accessor.
Definition: dist_array.h:917
Forward declarations.
Definition: dist_array.h:57
value_type::value_type element_type
Definition: dist_array.h:102
World & world() const
World accessor.
Definition: dist_array.h:1007
BlockCyclicMatrix< typename std::remove_cv_t< Array >::element_type > array_to_block_cyclic(const Array &array, const blacspp::Grid &grid, size_t MB, size_t NB)
Convert a dense DistArray to block-cyclic storage format.
Definition: block_cyclic.h:308