cleanup complex API

INSTALL.md

@@ -40,7 +40,7 @@ Both methods are supported. However, for most users we _strongly_ recommend to b
   - Boost.Container: header-only
   - Boost.Test: header-only or (optionally) as a compiled library, *only used for unit testing*
   - Boost.Range: header-only, *only used for unit testing*
-- [BTAS](http://github.com/ValeevGroup/BTAS), tag 6fcb6451bc7ca46a00534a30c51dc5c230c39ac3 . If usable BTAS installation is not found, TiledArray will download and compile
+- [BTAS](http://github.com/ValeevGroup/BTAS), tag 561fe1bff7f3374814111a15e28c7a141ab9b67a . If usable BTAS installation is not found, TiledArray will download and compile
   BTAS from source. *This is the recommended way to compile BTAS for all users*.
 - [MADNESS](https://github.com/m-a-d-n-e-s-s/madness), tag 91fff76deba20c751d0646c54f2f1c1e07bd6156 .
   Only the MADworld runtime and BLAS/LAPACK C API component of MADNESS is used by TiledArray.

external/versions.cmake

@@ -24,8 +24,8 @@ set(TA_TRACKED_MADNESS_PREVIOUS_TAG 0b44ef319643cb9721fbe17d294987c146e6460e)
 set(TA_TRACKED_MADNESS_VERSION 0.10.1)
 set(TA_TRACKED_MADNESS_PREVIOUS_VERSION 0.10.1)
 
-set(TA_TRACKED_BTAS_TAG 6fcb6451bc7ca46a00534a30c51dc5c230c39ac3)
-set(TA_TRACKED_BTAS_PREVIOUS_TAG 474ddc095cbea12a1d28aca5435703dd9f69b166)
+set(TA_TRACKED_BTAS_TAG 561fe1bff7f3374814111a15e28c7a141ab9b67a)
+set(TA_TRACKED_BTAS_PREVIOUS_TAG 6fcb6451bc7ca46a00534a30c51dc5c230c39ac3)
 
 set(TA_TRACKED_LIBRETT_TAG 68abe31a9ec6fd2fd9ffbcd874daa80457f947da)
 set(TA_TRACKED_LIBRETT_PREVIOUS_TAG 7e27ac766a9038df6aa05613784a54a036c4b796)

src/TiledArray/dist_array.h

@@ -136,6 +136,27 @@ class DistArray : public madness::archive::ParallelSerializableObject {
       std::is_same_v<std::decay_t<Value>, Future<value_type>> ||
       std::is_same_v<std::decay_t<Value>, value_type>;
 
+  /// compute type of DistArray with different Policy and/or Tile
+  template <typename TileU, typename PolicyU = Policy>
+  using rebind_t = DistArray<TileU, PolicyU>;
+
+ private:
+  template <typename Numeric, typename = void>
+  struct rebind_numeric;
+  template <typename Numeric>
+  struct rebind_numeric<
+      Numeric, std::enable_if_t<detail::has_rebind_numeric_v<Tile, Numeric>>> {
+    using type =
+        DistArray<typename Tile::template rebind_numeric_t<Numeric>, Policy>;
+  };
+
+ public:
+  /// compute type of DistArray with Tile's rebound numeric type
+  /// @note this is SFINAE-disabled if `Tile::rebind_numeric_t<Numeric>` is not
+  /// defined
+  template <typename Numeric>
+  using rebind_numeric_t = typename rebind_numeric<Numeric>::type;
+
  private:
   pimpl_type pimpl_;  ///< managed ptr to Array implementation
   bool defer_deleter_to_next_fence_ =
@@ -1838,6 +1859,22 @@ DistArray<T, P> replicated(const DistArray<T, P>& a) {
   return result;
 }
 
+namespace detail {
+
+template <typename Tile, typename Policy>
+struct real_t_impl<DistArray<Tile, Policy>> {
+  using type = typename DistArray<Tile, Policy>::template rebind_numeric_t<
+      typename Tile::scalar_type>;
+};
+
+template <typename Tile, typename Policy>
+struct complex_t_impl<DistArray<Tile, Policy>> {
+  using type = typename DistArray<Tile, Policy>::template rebind_numeric_t<
+      std::complex<typename Tile::scalar_type>>;
+};
+
+}  // namespace detail
+
 }  // namespace TiledArray
 
 // serialization

src/TiledArray/external/btas.h

@@ -841,6 +841,20 @@ struct ordinal_traits<btas::RangeNd<_Order, _Index, _Ordinal>> {
                                          : OrdinalType::ColMajor;
 };
 
+template <typename T, typename Range, typename Storage>
+struct real_t_impl<btas::Tensor<T, Range, Storage>> {
+  using type =
+      typename btas::Tensor<T, Range, Storage>::template rebind_numeric_t<
+          typename btas::Tensor<T, Range, Storage>::scalar_type>;
+};
+
+template <typename T, typename Range, typename Storage>
+struct complex_t_impl<btas::Tensor<T, Range, Storage>> {
+  using type =
+      typename btas::Tensor<T, Range, Storage>::template rebind_numeric_t<
+          std::complex<typename btas::Tensor<T, Range, Storage>::scalar_type>>;
+};
+
 }  // namespace detail
 }  // namespace TiledArray
 

src/TiledArray/math/linalg/non-distributed/heig.h

@@ -52,10 +52,10 @@ namespace TiledArray::math::linalg::non_distributed {
  */
 template <typename Array>
 auto heig(const Array& A, TiledRange evec_trange = TiledRange()) {
-  using numeric_type = typename detail::array_traits<Array>::numeric_type;
+  using scalar_type = typename detail::array_traits<Array>::scalar_type;
   World& world = A.world();
   auto A_eig = detail::make_matrix(A);
-  std::vector<numeric_type> evals;
+  std::vector<scalar_type> evals;
   if (world.rank() == 0) {
     linalg::rank_local::heig(A_eig, evals);
   }
@@ -93,12 +93,12 @@ auto heig(const Array& A, TiledRange evec_trange = TiledRange()) {
 template <typename ArrayA, typename ArrayB, typename EVecType = ArrayA>
 auto heig(const ArrayA& A, const ArrayB& B,
           TiledRange evec_trange = TiledRange()) {
-  using numeric_type = typename detail::array_traits<ArrayA>::numeric_type;
+  using scalar_type = typename detail::array_traits<ArrayA>::scalar_type;
   (void)detail::array_traits<ArrayB>{};
   World& world = A.world();
   auto A_eig = detail::make_matrix(A);
   auto B_eig = detail::make_matrix(B);
-  std::vector<numeric_type> evals;
+  std::vector<scalar_type> evals;
   if (world.rank() == 0) {
     linalg::rank_local::heig(A_eig, B_eig, evals);
   }

src/TiledArray/math/linalg/non-distributed/svd.h

@@ -27,9 +27,9 @@
 
 #include <TiledArray/config.h>
 
-#include <TiledArray/math/linalg/util.h>
-#include <TiledArray/math/linalg/rank-local.h>
 #include <TiledArray/conversions/eigen.h>
+#include <TiledArray/math/linalg/rank-local.h>
+#include <TiledArray/math/linalg/util.h>
 
 namespace TiledArray::math::linalg::non_distributed {
 
@@ -52,13 +52,14 @@ namespace TiledArray::math::linalg::non_distributed {
  *  @param[in] vt_trange   TiledRange for resulting right singular vectors
  * (transposed).
  *
- *  @returns A tuple containing the eigenvalues and eigenvectors of input array
- *  as std::vector and in TA format, respectively.
+ *  @returns A tuple containing the singular values and singular vectors of
+ * input array as std::vector and in TA format, respectively.
  */
-template<SVD::Vectors Vectors, typename Array>
-auto svd(const Array& A, TiledRange u_trange = TiledRange(), TiledRange vt_trange = TiledRange()) {
-
+template <SVD::Vectors Vectors, typename Array>
+auto svd(const Array& A, TiledRange u_trange = TiledRange(),
+         TiledRange vt_trange = TiledRange()) {
   using T = typename Array::numeric_type;
+  using TS = typename Array::scalar_type;
   using Matrix = linalg::rank_local::Matrix<T>;
 
   World& world = A.world();
@@ -68,7 +69,7 @@ auto svd(const Array& A, TiledRange u_trange = TiledRange(), TiledRange vt_trang
   constexpr bool need_u = (Vectors == SVD::LeftVectors) or svd_all_vectors;
   constexpr bool need_vt = (Vectors == SVD::RightVectors) or svd_all_vectors;
 
-  std::vector<T> S;
+  std::vector<TS> S;
   std::unique_ptr<Matrix> U, VT;
 
   if constexpr (need_u) U = std::make_unique<Matrix>();
@@ -82,7 +83,7 @@ auto svd(const Array& A, TiledRange u_trange = TiledRange(), TiledRange vt_trang
   if (U) world.gop.broadcast_serializable(*U, 0);
   if (VT) world.gop.broadcast_serializable(*VT, 0);
 
-  auto make_array = [&world](auto && ... args) {
+  auto make_array = [&world](auto&&... args) {
     return eigen_to_array<Array>(world, args...);
   };
 
@@ -97,7 +98,6 @@ auto svd(const Array& A, TiledRange u_trange = TiledRange(), TiledRange vt_trang
   }
 
   if constexpr (!need_u && !need_vt) return S;
-
 }
 
 }  // namespace TiledArray::math::linalg::non_distributed

src/TiledArray/math/linalg/rank-local.cpp

@@ -113,19 +113,23 @@ void cholesky_lsolve(Op transpose, Matrix<T>& A, Matrix<T>& X) {
 }
 
 template <typename T>
-void heig(Matrix<T>& A, std::vector<T>& W) {
+void heig(Matrix<T>& A, std::vector<TiledArray::detail::real_t<T>>& W) {
   auto jobz = lapack::Job::Vec;
   auto uplo = lapack::Uplo::Lower;
   integer n = A.rows();
   T* a = A.data();
   integer lda = A.rows();
   W.resize(n);
-  T* w = W.data();
-  TA_LAPACK(syev, jobz, uplo, n, a, lda, w);
+  auto* w = W.data();
+  if constexpr (TiledArray::detail::is_complex_v<T>)
+    TA_LAPACK(heev, jobz, uplo, n, a, lda, w);
+  else
+    TA_LAPACK(syev, jobz, uplo, n, a, lda, w);
 }
 
 template <typename T>
-void heig(Matrix<T>& A, Matrix<T>& B, std::vector<T>& W) {
+void heig(Matrix<T>& A, Matrix<T>& B,
+          std::vector<TiledArray::detail::real_t<T>>& W) {
   integer itype = 1;
   auto jobz = lapack::Job::Vec;
   auto uplo = lapack::Uplo::Lower;
@@ -135,53 +139,60 @@ void heig(Matrix<T>& A, Matrix<T>& B, std::vector<T>& W) {
   T* b = B.data();
   integer ldb = B.rows();
   W.resize(n);
-  T* w = W.data();
-  TA_LAPACK(sygv, itype, jobz, uplo, n, a, lda, b, ldb, w);
+  auto* w = W.data();
+  if constexpr (TiledArray::detail::is_complex_v<T>)
+    TA_LAPACK(hegv, itype, jobz, uplo, n, a, lda, b, ldb, w);
+  else
+    TA_LAPACK(sygv, itype, jobz, uplo, n, a, lda, b, ldb, w);
 }
 
 template <typename T>
-void svd(Job jobu, Job jobvt, Matrix<T>& A, std::vector<T>& S, Matrix<T>* U, Matrix<T>* VT) {
+void svd(Job jobu, Job jobvt, Matrix<T>& A,
+         std::vector<TiledArray::detail::real_t<T>>& S, Matrix<T>* U,
+         Matrix<T>* VT) {
   integer m = A.rows();
   integer n = A.cols();
   integer k = std::min(m, n);
   T* a = A.data();
   integer lda = A.rows();
 
   S.resize(k);
-  T* s = S.data();
+  auto* s = S.data();
 
-  T* u  = nullptr;
+  T* u = nullptr;
   T* vt = nullptr;
   integer ldu = 1, ldvt = 1;
-  if( (jobu == Job::SomeVec or jobu == Job::AllVec) and (not U) ) 
-    TA_LAPACK_ERROR("Requested out-of-place right singular vectors with null U input");
-  if( (jobvt == Job::SomeVec or jobvt == Job::AllVec) and (not VT) ) 
-    TA_LAPACK_ERROR("Requested out-of-place left singular vectors with null VT input");
+  if ((jobu == Job::SomeVec or jobu == Job::AllVec) and (not U))
+    TA_LAPACK_ERROR(
+        "Requested out-of-place right singular vectors with null U input");
+  if ((jobvt == Job::SomeVec or jobvt == Job::AllVec) and (not VT))
+    TA_LAPACK_ERROR(
+        "Requested out-of-place left singular vectors with null VT input");
 
-  if( jobu == Job::SomeVec ) {
+  if (jobu == Job::SomeVec) {
     U->resize(m, k);
     u = U->data();
     ldu = m;
   }
 
-  if( jobu == Job::AllVec ) {
+  if (jobu == Job::AllVec) {
     U->resize(m, m);
     u = U->data();
     ldu = m;
   }
 
-  if( jobvt == Job::SomeVec ) {
+  if (jobvt == Job::SomeVec) {
     VT->resize(k, n);
     vt = VT->data();
     ldvt = k;
   }
 
-  if( jobvt == Job::AllVec ) {
+  if (jobvt == Job::AllVec) {
     VT->resize(n, n);
     vt = VT->data();
     ldvt = n;
   }
-    
+
   TA_LAPACK(gesvd, jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt);
 }
 
@@ -208,47 +219,51 @@ void lu_inv(Matrix<T>& A) {
 }
 
 template <bool QOnly, typename T>
-void householder_qr( Matrix<T> &V, Matrix<T> &R ) {
+void householder_qr(Matrix<T>& V, Matrix<T>& R) {
   integer m = V.rows();
   integer n = V.cols();
-  integer k = std::min(m,n);
-  integer ldv = V.rows(); // Col Major
+  integer k = std::min(m, n);
+  integer ldv = V.rows();  // Col Major
   T* v = V.data();
   std::vector<T> tau(k);
-  lapack::geqrf( m, n, v, ldv, tau.data() );
+  lapack::geqrf(m, n, v, ldv, tau.data());
 
   // Extract R
-  if constexpr ( not QOnly ) {
+  if constexpr (not QOnly) {
     // Resize R just in case
-    R.resize(k,n);
+    R.resize(k, n);
     R.fill(0.);
     // Extract Upper triangle into R
     integer ldr = R.rows();
     T* r = R.data();
-    lapack::lacpy( lapack::MatrixType::Upper, k, n, v, ldv, r, ldr );
+    lapack::lacpy(lapack::MatrixType::Upper, k, n, v, ldv, r, ldr);
   }
 
   // Explicitly form Q
   // TODO: This is wrong for complex, but it doesn't look like R/C is caught
   //       anywhere else either...
-  lapack::orgqr( m, n, k, v, ldv, tau.data() );
-
+  if constexpr (TiledArray::detail::is_complex_v<T>)
+    lapack::ungqr(m, n, k, v, ldv, tau.data());
+  else
+    lapack::orgqr(m, n, k, v, ldv, tau.data());
 }
 
-#define TA_LAPACK_EXPLICIT(MATRIX, VECTOR)                       \
-  template void cholesky(MATRIX&);                               \
-  template void cholesky_linv(MATRIX&);                          \
-  template void cholesky_solve(MATRIX&, MATRIX&);                \
-  template void cholesky_lsolve(Op, MATRIX&, MATRIX&);           \
-  template void heig(MATRIX&, VECTOR&);                          \
-  template void heig(MATRIX&, MATRIX&, VECTOR&);                 \
-  template void svd(Job,Job,MATRIX&, VECTOR&, MATRIX*, MATRIX*); \
-  template void lu_solve(MATRIX&, MATRIX&);                      \
-  template void lu_inv(MATRIX&);                                 \
-  template void householder_qr<true>(MATRIX&,MATRIX&);           \
-  template void householder_qr<false>(MATRIX&,MATRIX&);
+#define TA_LAPACK_EXPLICIT(MATRIX, VECTOR)                         \
+  template void cholesky(MATRIX&);                                 \
+  template void cholesky_linv(MATRIX&);                            \
+  template void cholesky_solve(MATRIX&, MATRIX&);                  \
+  template void cholesky_lsolve(Op, MATRIX&, MATRIX&);             \
+  template void heig(MATRIX&, VECTOR&);                            \
+  template void heig(MATRIX&, MATRIX&, VECTOR&);                   \
+  template void svd(Job, Job, MATRIX&, VECTOR&, MATRIX*, MATRIX*); \
+  template void lu_solve(MATRIX&, MATRIX&);                        \
+  template void lu_inv(MATRIX&);                                   \
+  template void householder_qr<true>(MATRIX&, MATRIX&);            \
+  template void householder_qr<false>(MATRIX&, MATRIX&);
 
 TA_LAPACK_EXPLICIT(Matrix<double>, std::vector<double>);
 TA_LAPACK_EXPLICIT(Matrix<float>, std::vector<float>);
+TA_LAPACK_EXPLICIT(Matrix<std::complex<double>>, std::vector<double>);
+TA_LAPACK_EXPLICIT(Matrix<std::complex<float>>, std::vector<float>);
 
 }  // namespace TiledArray::math::linalg::rank_local

src/TiledArray/math/linalg/rank-local.h

@@ -42,17 +42,20 @@ template <typename T>
 void cholesky_lsolve(Op transpose, Matrix<T> &A, Matrix<T> &X);
 
 template <typename T>
-void heig(Matrix<T> &A, std::vector<T> &W);
+void heig(Matrix<T> &A, std::vector<TiledArray::detail::real_t<T>> &W);
 
 template <typename T>
-void heig(Matrix<T> &A, Matrix<T> &B, std::vector<T> &W);
+void heig(Matrix<T> &A, Matrix<T> &B,
+          std::vector<TiledArray::detail::real_t<T>> &W);
 
 template <typename T>
-void svd(Job jobu, Job jobvt, Matrix<T> &A, std::vector<T> &S, Matrix<T> *U,
+void svd(Job jobu, Job jobvt, Matrix<T> &A,
+         std::vector<TiledArray::detail::real_t<T>> &S, Matrix<T> *U,
          Matrix<T> *VT);
 
 template <typename T>
-void svd(Matrix<T> &A, std::vector<T> &S, Matrix<T> *U, Matrix<T> *VT) {
+void svd(Matrix<T> &A, std::vector<TiledArray::detail::real_t<T>> &S,
+         Matrix<T> *U, Matrix<T> *VT) {
   svd(U ? Job::SomeVec : Job::NoVec, VT ? Job::SomeVec : Job::NoVec, A, S, U,
       VT);
 }