Andrey's revamp of generic algebra

CMakeLists.txt

@@ -59,9 +59,9 @@ endif(TILEDARRAY_PRERELEASE_ID)
 
 # Set install paths ============================================================
 
-set(TILEDARRAY_INSTALL_BINDIR "bin" 
+set(TILEDARRAY_INSTALL_BINDIR "bin"
     CACHE PATH "TiledArray binary install directory")
-set(TILEDARRAY_INSTALL_INCLUDEDIR "include" 
+set(TILEDARRAY_INSTALL_INCLUDEDIR "include"
     CACHE PATH "TiledArray INCLUDE install directory")
 set(TILEDARRAY_INSTALL_LIBDIR "lib"
     CACHE PATH "TiledArray LIB install directory")
@@ -162,16 +162,15 @@ else ()
 endif()
 redefaultable_option(CMAKE_POSITION_INDEPENDENT_CODE "Default value for POSITION_INDEPENDENT_CODE of targets" ${default_CMAKE_POSITION_INDEPENDENT_CODE})
 
+set(BLA_STATIC FALSE CACHE BOOL "Whether to use static linkage for BLAS, LAPACK, and related libraries")
 if(BUILD_SHARED_LIBS)
-  set(BLA_STATIC FALSE CACHE BOOL "Whether to use static linkage for BLAS, LAPACK, and related libraries")
   set(CMAKE_MACOSX_RPATH TRUE)
 else()
-  set(BLA_STATIC TRUE CACHE BOOL "Whether to use static linkage for BLAS, LAPACK, and related libraries")
   set(CMAKE_MACOSX_RPATH FALSE)
 endif()
 
 # miscellaneous cmake platform-neutral and platform-specific configuration =============================
-set(CMAKE_FIND_NO_INSTALL_PREFIX TRUE)  # do not search in CMAKE_INSTALL_PREFIX 
+set(CMAKE_FIND_NO_INSTALL_PREFIX TRUE)  # do not search in CMAKE_INSTALL_PREFIX
 set(CMAKE_SKIP_RPATH FALSE)
 set(CMAKE_SKIP_BUILD_RPATH  FALSE)
 set(CMAKE_SKIP_INSTALL_RPATH FALSE)
@@ -394,19 +393,19 @@ export(EXPORT tiledarray
 configure_package_config_file(cmake/tiledarray-config.cmake.in
     "${PROJECT_BINARY_DIR}/tiledarray-config.cmake"
   INSTALL_DESTINATION "${TILEDARRAY_INSTALL_CMAKEDIR}"
-  PATH_VARS CMAKE_INSTALL_PREFIX TILEDARRAY_INSTALL_BINDIR 
+  PATH_VARS CMAKE_INSTALL_PREFIX TILEDARRAY_INSTALL_BINDIR
             TILEDARRAY_INSTALL_INCLUDEDIR TILEDARRAY_INSTALL_LIBDIR
             TILEDARRAY_INSTALL_DOCDIR TILEDARRAY_INSTALL_CMAKEDIR)
 
 # Install config, version, and target files
 install(EXPORT tiledarray
     FILE "tiledarray-targets.cmake"
-    DESTINATION "${TILEDARRAY_INSTALL_CMAKEDIR}" 
+    DESTINATION "${TILEDARRAY_INSTALL_CMAKEDIR}"
     COMPONENT tiledarray)
 install(FILES
     "${PROJECT_BINARY_DIR}/tiledarray-config.cmake"
     "${PROJECT_BINARY_DIR}/tiledarray-config-version.cmake"
-    DESTINATION "${TILEDARRAY_INSTALL_CMAKEDIR}" 
+    DESTINATION "${TILEDARRAY_INSTALL_CMAKEDIR}"
     COMPONENT tiledarray)
 
 

src/CMakeLists.txt

@@ -60,15 +60,17 @@ TiledArray/zero_tensor.h
 TiledArray/algebra/conjgrad.h
 TiledArray/algebra/diis.h
 TiledArray/algebra/utils.h
-TiledArray/algebra/chol.h
+TiledArray/algebra/cholesky.h
 TiledArray/algebra/heig.h
 TiledArray/algebra/lu.h
 TiledArray/algebra/svd.h
 TiledArray/algebra/types.h
-TiledArray/algebra/lapack/chol.h
+TiledArray/algebra/lapack/cholesky.h
 TiledArray/algebra/lapack/heig.h
 TiledArray/algebra/lapack/util.h
-TiledArray/algebra/scalapack/chol.h
+TiledArray/algebra/lapack/lu.h
+TiledArray/algebra/lapack/svd.h
+TiledArray/algebra/scalapack/cholesky.h
 TiledArray/algebra/scalapack/heig.h
 TiledArray/algebra/scalapack/lu.h
 TiledArray/algebra/scalapack/svd.h
@@ -207,7 +209,7 @@ TiledArray/array_impl.cpp
 TiledArray/dist_array.cpp
 TiledArray/util/backtrace.cpp
 TiledArray/util/bug.cpp
-TiledArray/algebra/lapack/lapack.cc
+TiledArray/algebra/lapack/lapack.cpp
 )
 
 # the list of libraries on which TiledArray depends on, will be cached later

src/TiledArray/algebra/chol.h → src/TiledArray/algebra/cholesky.h

@@ -26,9 +26,9 @@
 
 #include <TiledArray/config.h>
 #if TILEDARRAY_HAS_SCALAPACK
-#include <TiledArray/algebra/scalapack/chol.h>
+#include <TiledArray/algebra/scalapack/cholesky.h>
 #endif
-#include <TiledArray/algebra/lapack/chol.h>
+#include <TiledArray/algebra/lapack/cholesky.h>
 
 namespace TiledArray {
 
@@ -37,19 +37,17 @@ auto cholesky(const Array& A, TiledRange l_trange = TiledRange()) {
 #if TILEDARRAY_HAS_SCALAPACK
   if (A.world().size() > 1 && A.range().volume() > 10000000)
     return scalapack::cholesky<Array>(A, l_trange);
-  else
 #endif
-    return lapack::cholesky<Array>(A, l_trange);
+  return lapack::cholesky<Array>(A, l_trange);
 }
 
 template <typename Array, bool RetL = false>
 auto cholesky_linv(const Array& A, TiledRange l_trange = TiledRange()) {
 #if TILEDARRAY_HAS_SCALAPACK
   if (A.world().size() > 1 && A.range().volume() > 10000000)
     return scalapack::cholesky_linv<Array, RetL>(A, l_trange);
-  else
 #endif
-    return lapack::cholesky_linv<Array, RetL>(A, l_trange);
+  return lapack::cholesky_linv<Array, RetL>(A, l_trange);
 }
 
 template <typename Array>
@@ -58,9 +56,8 @@ auto cholesky_solve(const Array& A, const Array& B,
 #if TILEDARRAY_HAS_SCALAPACK
   if (A.world().size() > 1 && A.range().volume() > 10000000)
     return scalapack::cholesky_solve<Array>(A, B, x_trange);
-  else
 #endif
-    return lapack::cholesky_solve<Array>(A, B, x_trange);
+  return lapack::cholesky_solve<Array>(A, B, x_trange);
 }
 
 template <typename Array>
@@ -71,9 +68,8 @@ auto cholesky_lsolve(TransposeFlag transpose, const Array& A, const Array& B,
   if (A.world().size() > 1 && A.range().volume() > 10000000)
     return scalapack::cholesky_lsolve<Array>(transpose, A, B, l_trange,
                                              x_trange);
-  else
 #endif
-    return lapack::cholesky_lsolve<Array>(transpose, A, B, l_trange, x_trange);
+  return lapack::cholesky_lsolve<Array>(transpose, A, B, l_trange, x_trange);
 }
 
 }  // namespace TiledArray

src/TiledArray/algebra/heig.h

@@ -27,16 +27,30 @@
 #include <TiledArray/config.h>
 #if TILEDARRAY_HAS_SCALAPACK
 #include <TiledArray/algebra/scalapack/heig.h>
-#else
-// eigen
 #endif
+#include <TiledArray/algebra/lapack/heig.h>
 
 namespace TiledArray {
 
+template <typename Array>
+auto heig(const Array& A, TiledRange evec_trange = TiledRange()) {
 #if TILEDARRAY_HAS_SCALAPACK
-using scalapack::heig;
+  if (A.world().size() > 1 && A.range().volume() > 10000000) {
+    return scalapack::heig(A, evec_trange);
+  }
 #endif
+  return lapack::heig(A, evec_trange);
+}
+
+template <typename ArrayA, typename ArrayB, typename EVecType = ArrayA>
+auto heig(const ArrayA& A, const ArrayB& B, TiledRange evec_trange = TiledRange()) {
+#if TILEDARRAY_HAS_SCALAPACK
+  if (A.world().size() > 1 && A.range().volume() > 10000000) {
+    return scalapack::heig(A, B, evec_trange);
+  }
 #endif
+  return lapack::heig(A, B, evec_trange);
+}
 
 }  // namespace TiledArray
 

src/TiledArray/algebra/lapack/chol.h → src/TiledArray/algebra/lapack/cholesky.h

@@ -24,9 +24,9 @@
 #ifndef TILEDARRAY_ALGEBRA_LAPACK_CHOL_H__INCLUDED
 #define TILEDARRAY_ALGEBRA_LAPACK_CHOL_H__INCLUDED
 
-#include <TiledArray/config.h>
 #include <TiledArray/algebra/lapack/lapack.h>
 #include <TiledArray/algebra/lapack/util.h>
+#include <TiledArray/config.h>
 #include <TiledArray/conversions/eigen.h>
 
 namespace TiledArray {
@@ -45,7 +45,7 @@ auto make_L_eig(const DistArray<Tile, Policy>& A) {
   World& world = A.world();
   auto A_eig = detail::to_eigen(A);
   if (world.rank() == 0) {
-    chol_eig(A_eig);
+    lapack::cholesky(A_eig);
   }
   world.gop.broadcast_serializable(A_eig, 0);
   return A_eig;

src/TiledArray/algebra/lapack/heig.h

@@ -24,15 +24,16 @@
 #ifndef TILEDARRAY_ALGEBRA_LAPACK_HEIG_H__INCLUDED
 #define TILEDARRAY_ALGEBRA_LAPACK_HEIG_H__INCLUDED
 
-#include <TiledArray/algebra/lapack/util.h>
 #include <TiledArray/config.h>
+
+#include <TiledArray/algebra/lapack/lapack.h>
+#include <TiledArray/algebra/lapack/util.h>
 #include <TiledArray/conversions/eigen.h>
 
-namespace TiledArray {
-namespace lapack {
+namespace TiledArray::lapack {
 
 /**
- *  @brief Solve the standard eigenvalue problem with ScaLAPACK
+ *  @brief Solve the standard eigenvalue problem with LAPACK
  *
  *  A(i,k) X(k,j) = X(i,j) E(j)
  *
@@ -51,80 +52,24 @@ namespace lapack {
  */
 template <typename Array>
 auto heig(const Array& A, TiledRange evec_trange = TiledRange()) {
-  using scalar_type = typename Array::scalar_type;
-  using numeric_type = typename Array::numeric_type;
-  constexpr const bool is_real = std::is_same_v<scalar_type, numeric_type>;
-  static_assert(std::is_same_v<numeric_type, typename Array::element_type>,
-                "TA::lapack::{cholesky*} are only usable with a DistArray of "
-                "scalar types");
-
+  using numeric_type = typename lapack::array_traits<Array>::numeric_type;
   World& world = A.world();
   auto A_eig = detail::to_eigen(A);
-  std::vector<scalar_type> evals;
-//   if (world.rank() == 0) {
-//     char jobz = 'V';
-//     char uplo = 'L';
-//     integer n = A_eig.rows();
-//     numeric_type* a = A_eig.data();
-//     integer lda = n;
-//     integer info = 0;
-//     evals.resize(n);
-//     integer lwork = -1;
-//     std::vector<numeric_type> work(1);
-//     // run once to query, then to compute
-//     while (lwork != static_cast<integer>(work.size())) {
-//       if (lwork > 0) {
-//         work.resize(lwork);
-//       }
-//       if constexpr (is_real) {
-// #if defined(MADNESS_LINALG_USE_LAPACKE)
-//         MADNESS_DISPATCH_LAPACK_FN(syev, &jobz, &uplo, &n, a, &lda,
-//                                    evals.data(), work.data(), &lwork, &info);
-// #else
-//         MADNESS_DISPATCH_LAPACK_FN(syev, &jobz, &uplo, &n, a, &lda,
-//                                    evals.data(), work.data(), &lwork, &info,
-//                                    sizeof(char), sizeof(char));
-// #endif
-//       } else {
-//         std::vector<scalar_type> rwork;
-//         if (lwork == static_cast<integer>(work.size())) rwork.resize(3 * n - 2);
-// #if defined(MADNESS_LINALG_USE_LAPACKE)
-//         MADNESS_DISPATCH_LAPACK_FN(heev, &jobz, &uplo, &n, a, &lda,
-//                                    evals.data(), work.data(), &lwork,
-//                                    &rwork.data(), &info);
-// #else
-//         MADNESS_DISPATCH_LAPACK_FN(
-//             heev, &jobz, &uplo, &n, a, &lda, evals.data(), work.data(), &lwork,
-//             &rwork.data(), &info, sizeof(char), sizeof(char));
-// #endif
-//       }
-//       if (lwork == -1) {
-//         if constexpr (is_real) {
-//           lwork = static_cast<integer>(work[0]);
-//         } else {
-//           lwork = static_cast<integer>(work[0].real());
-//         }
-//         TA_ASSERT(lwork > 1);
-//       }
-//     };
-
-//     if (info != 0) {
-//       if (is_real)
-//         TA_EXCEPTION("LAPACK::syev failed");
-//       else
-//         TA_EXCEPTION("LAPACK::heev failed");
-//     }
-//   }
-
+  std::vector<numeric_type> evals;
+  if (world.rank() == 0) {
+    lapack::heig(A_eig, evals);
+  }
   world.gop.broadcast_serializable(A_eig, 0);
   world.gop.broadcast_serializable(evals, 0);
   if (evec_trange.rank() == 0) evec_trange = A.trange();
-  return std::tuple(evals,
-                    eigen_to_array<Array>(A.world(), evec_trange, A_eig));
+  return std::tuple(
+    evals,
+    eigen_to_array<Array>(world, evec_trange, A_eig)
+  );
 }
 
 /**
- *  @brief Solve the generalized eigenvalue problem with ScaLAPACK
+ *  @brief Solve the generalized eigenvalue problem with LAPACK
  *
  *  A(i,k) X(k,j) = B(i,k) X(k,j) E(j)
  *
@@ -142,26 +87,30 @@ auto heig(const Array& A, TiledRange evec_trange = TiledRange()) {
  *  @param[in] B           Positive-definite matrix
  *  @param[in] evec_trange TiledRange for resulting eigenvectors. If left empty,
  *                         will default to array.trange()
- *  @param[in] NB          ScaLAPACK block size. Defaults to 128
  *
  *  @returns A tuple containing the eigenvalues and eigenvectors of input array
  *  as std::vector and in TA format, respectively.
  */
 template <typename ArrayA, typename ArrayB, typename EVecType = ArrayA>
-auto heig(const ArrayA& A, const ArrayB& B,
-          TiledRange evec_trange = TiledRange()) {
-  using scalar_type = typename ArrayA::scalar_type;
-  using numeric_type = typename ArrayA::numeric_type;
-  constexpr const bool is_real = std::is_same_v<scalar_type, numeric_type>;
-  static_assert(std::is_same_v<numeric_type, typename ArrayA::element_type>,
-                "TA::lapack::{cholesky*} are only usable with a DistArray of "
-                "scalar types");
-
-  abort();
-  return std::tuple(std::vector<scalar_type>{}, EVecType{});
+auto heig(const ArrayA& A, const ArrayB& B, TiledRange evec_trange = TiledRange()) {
+  using numeric_type = typename lapack::array_traits<ArrayA>::numeric_type;
+  (void)lapack::array_traits<ArrayB>{};
+  World& world = A.world();
+  auto A_eig = detail::to_eigen(A);
+  auto B_eig = detail::to_eigen(B);
+  std::vector<numeric_type> evals;
+  if (world.rank() == 0) {
+    lapack::heig(A_eig, B_eig, evals);
+  }
+  world.gop.broadcast_serializable(A_eig, 0);
+  world.gop.broadcast_serializable(evals, 0);
+  if (evec_trange.rank() == 0) evec_trange = A.trange();
+  return std::tuple(
+    evals,
+    eigen_to_array<ArrayA>(A.world(), evec_trange, A_eig)
+  );
 }
 
-}  // namespace lapack
-}  // namespace TiledArray
+}  // namespace TiledArray::lapack
 
-#endif  // TILEDARRAY_ALGEBRA_SCALAPACK_HEIG_H__INCLUDED
+#endif  // TILEDARRAY_ALGEBRA_LAPACK_HEIG_H__INCLUDED