From e7409bcb382bb77040ded5984fe1230a9f4e4dbb Mon Sep 17 00:00:00 2001 From: asadchev Date: Sun, 29 Nov 2020 17:49:01 -0500 Subject: [PATCH 1/7] Replace f77 lapack calls with lapacke --- src/TiledArray/algebra/lapack/lapack.cpp | 171 ++++++++++++----------- src/TiledArray/algebra/lapack/lapack.h | 78 ----------- 2 files changed, 86 insertions(+), 163 deletions(-) diff --git a/src/TiledArray/algebra/lapack/lapack.cpp b/src/TiledArray/algebra/lapack/lapack.cpp index c29f0d350f..014f08f209 100644 --- a/src/TiledArray/algebra/lapack/lapack.cpp +++ b/src/TiledArray/algebra/lapack/lapack.cpp @@ -24,45 +24,60 @@ #include #include -#include + +#include + +#define TA_LAPACKE_THROW(F) throw std::runtime_error("lapacke::" #F " failed") + +#define TA_LAPACKE_CALL(F, ARGS...) \ + (LAPACKE_##F(LAPACK_COL_MAJOR, ARGS) == 0 || (TA_LAPACKE_THROW(F), false)) + +/// TA_LAPACKE(fn,args) can be called only from template context, with `T` +/// defining the element type +#define TA_LAPACKE(name, args...) { \ + using numeric_type = T; \ + if constexpr (std::is_same_v) \ + TA_LAPACKE_CALL(d##name, args); \ + else if constexpr (std::is_same_v) \ + TA_LAPACKE_CALL(s##name, args); \ + else if constexpr (std::is_same_v>) \ + TA_LAPACKE_CALL(z##name, args); \ + else if constexpr (std::is_same_v>) \ + TA_LAPACKE_CALL(c##name, args); \ + else std::abort(); \ + } namespace TiledArray::lapack { template void cholesky(Matrix& A) { char uplo = 'L'; - integer n = A.rows(); + lapack_int n = A.rows(); auto* a = A.data(); - integer lda = n; - integer info = 0; - TA_LAPACK_POTRF(&uplo, &n, a, &lda, &info); - if (info != 0) TA_EXCEPTION("lapack::cholesky failed"); + lapack_int lda = n; + TA_LAPACKE(potrf, uplo, n, a, lda); } template void cholesky_linv(Matrix& A) { char uplo = 'L'; char diag = 'N'; - integer n = A.rows(); + lapack_int n = A.rows(); auto* l = A.data(); - integer lda = n; - integer info = 0; - TA_LAPACK_TRTRI(&uplo, &diag, &n, l, &lda, &info); - if (info != 0) TA_EXCEPTION("lapack::cholesky_linv failed"); + lapack_int lda = n; + TA_LAPACKE(trtri, uplo, diag, n, l, lda); } template void cholesky_solve(Matrix& A, Matrix& X) { char uplo = 'L'; - integer n = A.rows(); - integer nrhs = X.cols(); + lapack_int n = A.rows(); + lapack_int nrhs = X.cols(); auto* a = A.data(); auto* b = X.data(); - integer lda = n; - integer ldb = n; - integer info = 0; - TA_LAPACK_POSV(&uplo, &n, &nrhs, a, &lda, b, &ldb, &info); - if (info != 0) TA_EXCEPTION("lapack::cholesky_solve failed"); + lapack_int lda = n; + lapack_int ldb = n; + TA_LAPACKE(posv, uplo, n, nrhs, a, lda, b, ldb); } template @@ -72,73 +87,65 @@ void cholesky_lsolve(TransposeFlag transpose, Matrix& A, Matrix& X) { ? 'T' : (transpose == TransposeFlag::NoTranspose ? 'N' : 'C'); char diag = 'N'; - integer n = A.rows(); - integer nrhs = X.cols(); + lapack_int n = A.rows(); + lapack_int nrhs = X.cols(); auto* a = A.data(); auto* b = X.data(); - integer lda = n; - integer ldb = n; - integer info = 0; - TA_LAPACK_TRTRS(&uplo, &trans, &diag, &n, &nrhs, a, &lda, b, &ldb, &info); - if (info != 0) TA_EXCEPTION("lapack::cholesky_lsolve failed"); + lapack_int lda = n; + lapack_int ldb = n; + TA_LAPACKE(trtrs, uplo, trans, diag, n, nrhs, a, lda, b, ldb); } template void heig(Matrix& A, std::vector& W) { char jobz = 'V'; char uplo = 'L'; - integer n = A.rows(); + lapack_int n = A.rows(); T* a = A.data(); - integer lda = A.rows(); + lapack_int lda = A.rows(); W.resize(n); T* w = W.data(); - integer lwork = -1; - integer info; - T lwork_dummy; - TA_LAPACK_SYEV(&jobz, &uplo, &n, a, &lda, w, &lwork_dummy, &lwork, &info); - lwork = integer(lwork_dummy); - std::vector work(lwork); - TA_LAPACK_SYEV(&jobz, &uplo, &n, a, &lda, w, work.data(), &lwork, &info); - if (info != 0) TA_EXCEPTION("lapack::heig failed"); + lapack_int lwork = -1; + std::vector work(1); + TA_LAPACKE(syev_work, jobz, uplo, n, a, lda, w, work.data(), lwork); + lwork = lapack_int(work[0]); + work.resize(lwork); + TA_LAPACKE(syev_work, jobz, uplo, n, a, lda, w, work.data(), lwork); } template void heig(Matrix& A, Matrix& B, std::vector& W) { - integer itype = 1; + lapack_int itype = 1; char jobz = 'V'; char uplo = 'L'; - integer n = A.rows(); + lapack_int n = A.rows(); T* a = A.data(); - integer lda = A.rows(); + lapack_int lda = A.rows(); T* b = B.data(); - integer ldb = B.rows(); + lapack_int ldb = B.rows(); W.resize(n); T* w = W.data(); - integer lwork = -1; - integer info; - T lwork_dummy; - TA_LAPACK_SYGV(&itype, &jobz, &uplo, &n, a, &lda, b, &ldb, w, &lwork_dummy, - &lwork, &info); - lwork = integer(lwork_dummy); - std::vector work(lwork); - TA_LAPACK_SYGV(&itype, &jobz, &uplo, &n, a, &lda, b, &ldb, w, work.data(), - &lwork, &info); - if (info != 0) TA_EXCEPTION("lapack::heig failed"); + std::vector work(1); + lapack_int lwork = -1; + TA_LAPACKE(sygv_work, itype, jobz, uplo, n, a, lda, b, ldb, w, work.data(), lwork); + lwork = lapack_int(work[0]); + work.resize(lwork); + TA_LAPACKE(sygv_work, itype, jobz, uplo, n, a, lda, b, ldb, w, work.data(), lwork); } template void svd(Matrix& A, std::vector& S, Matrix* U, Matrix* VT) { - integer m = A.rows(); - integer n = A.cols(); + lapack_int m = A.rows(); + lapack_int n = A.cols(); T* a = A.data(); - integer lda = A.rows(); + lapack_int lda = A.rows(); S.resize(std::min(m, n)); T* s = S.data(); char jobu = 'N'; T* u = nullptr; - integer ldu = m; + lapack_int ldu = m; if (U) { jobu = 'A'; U->resize(m, n); @@ -148,7 +155,7 @@ void svd(Matrix& A, std::vector& S, Matrix* U, Matrix* VT) { char jobvt = 'N'; T* vt = nullptr; - integer ldvt = n; + lapack_int ldvt = n; if (VT) { jobvt = 'A'; VT->resize(n, m); @@ -156,49 +163,43 @@ void svd(Matrix& A, std::vector& S, Matrix* U, Matrix* VT) { ldvt = VT->rows(); } - integer lwork = -1; - integer info; - T lwork_dummy; - - TA_LAPACK_GESVD(&jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, - &lwork_dummy, &lwork, &info); - lwork = integer(lwork_dummy); - std::vector work(lwork); - TA_LAPACK_GESVD(&jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, - work.data(), &lwork, &info); - if (info != 0) TA_EXCEPTION("lapack::svd failed"); + std::vector work(1); + lapack_int lwork = -1; + + TA_LAPACKE(gesvd_work, jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, + work.data(), lwork); + lwork = lapack_int(work[0]); + work.resize(lwork); + TA_LAPACKE(gesvd_work, jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, + work.data(), lwork); + } template void lu_solve(Matrix& A, Matrix& B) { - integer n = A.rows(); - integer nrhs = B.cols(); + lapack_int n = A.rows(); + lapack_int nrhs = B.cols(); T* a = A.data(); - integer lda = A.rows(); + lapack_int lda = A.rows(); T* b = B.data(); - integer ldb = B.rows(); - std::vector ipiv(n); - integer info; - TA_LAPACK_GESV(&n, &nrhs, a, &lda, ipiv.data(), b, &ldb, &info); - if (info != 0) TA_EXCEPTION("lapack::lu_solve failed"); + lapack_int ldb = B.rows(); + std::vector ipiv(n); + TA_LAPACKE(gesv, n, nrhs, a, lda, ipiv.data(), b, ldb); } template void lu_inv(Matrix& A) { - integer n = A.rows(); + lapack_int n = A.rows(); T* a = A.data(); - integer lda = A.rows(); - integer lwork = -1; + lapack_int lda = A.rows(); + std::vector ipiv(n); + TA_LAPACKE(getrf, n, n, a, lda, ipiv.data()); std::vector work(1); - std::vector ipiv(n); - integer info; - TA_LAPACK_GETRF(&n, &n, a, &lda, ipiv.data(), &info); - if (info != 0) TA_EXCEPTION("lapack::lu_inv failed"); - TA_LAPACK_GETRI(&n, a, &lda, ipiv.data(), work.data(), &lwork, &info); - lwork = (integer)work[0]; + lapack_int lwork = -1; + TA_LAPACKE(getri_work, n, a, lda, ipiv.data(), work.data(), lwork); + lwork = (lapack_int)work[0]; work.resize(lwork); - TA_LAPACK_GETRI(&n, a, &lda, ipiv.data(), work.data(), &lwork, &info); - if (info != 0) TA_EXCEPTION("lapack::lu_inv failed"); + TA_LAPACKE(getri_work, n, a, lda, ipiv.data(), work.data(), lwork); } #define TA_LAPACK_EXPLICIT(MATRIX, VECTOR) \ @@ -213,6 +214,6 @@ void lu_inv(Matrix& A) { template void lu_inv(MATRIX&); TA_LAPACK_EXPLICIT(lapack::Matrix, std::vector); -// TA_LAPACK_EXPLICIT(lapack::Matrix, std::vector); +TA_LAPACK_EXPLICIT(lapack::Matrix, std::vector); } // namespace TiledArray::lapack diff --git a/src/TiledArray/algebra/lapack/lapack.h b/src/TiledArray/algebra/lapack/lapack.h index 738543370f..f3d2b47f0e 100644 --- a/src/TiledArray/algebra/lapack/lapack.h +++ b/src/TiledArray/algebra/lapack/lapack.h @@ -26,86 +26,8 @@ #include #include -#include #include -/// TA_LAPACK_CALL(fn,args) can be called only from template context, with `T` -/// defining the element type -#define TA_LAPACK_CALL(name, args...) \ - using numeric_type = T; \ - if constexpr (std::is_same_v) \ - d##name##_(args); \ - else if constexpr (std::is_same_v) \ - s##name##_(args); \ - else if constexpr (std::is_same_v>) \ - z##name##_(args); \ - else if constexpr (std::is_same_v>) \ - c##name##_(args); \ - else \ - std::abort(); - -#define TA_LAPACK_GESV(...) TA_LAPACK_CALL(gesv, __VA_ARGS__) -#define TA_LAPACK_GETRF(...) TA_LAPACK_CALL(getrf, __VA_ARGS__) -#define TA_LAPACK_GETRI(...) TA_LAPACK_CALL(getri, __VA_ARGS__) - -#ifdef MADNESS_LINALG_USE_LAPACKE - -#define TA_LAPACK_POTRF(...) TA_LAPACK_CALL(potrf, __VA_ARGS__) -#define TA_LAPACK_POTRS(...) TA_LAPACK_CALL(potrs, __VA_ARGS__) -#define TA_LAPACK_POSV(...) TA_LAPACK_CALL(posv, __VA_ARGS__) -#define TA_LAPACK_GESVD(...) TA_LAPACK_CALL(gesvd, __VA_ARGS__) -#define TA_LAPACK_TRTRI(...) TA_LAPACK_CALL(trtri, __VA_ARGS__) -#define TA_LAPACK_TRTRS(...) TA_LAPACK_CALL(trtrs, __VA_ARGS__) -#define TA_LAPACK_SYEV(...) TA_LAPACK_CALL(syev, __VA_ARGS__) -#define TA_LAPACK_SYEVR(...) TA_LAPACK_CALL(syevr, __VA_ARGS__) -#define TA_LAPACK_SYGV(...) TA_LAPACK_CALL(sygv, __VA_ARGS__) - -#else - -#ifdef FORTRAN_LINKAGE_LCU -#define dtrtrs dtrtrs_ -#define dposv dposv_ -#define dpotrs dpotrs_ -#define dsyevr dsyevr_ -#endif - -extern "C" { // these arent in madness/clapack_fortran.h -void dtrtrs(const char* uplo, const char* trans, const char* diag, - const integer* n, const integer* nrhs, const real8* a, - const integer* lda, const real8* b, const integer* ldb, - integer* info, char_len, char_len, char_len); -void dposv(const char* uplo, const integer* n, const integer* nrhs, - const real8* a, const integer* lda, const real8* b, - const integer* ldb, integer* info, char_len); -void dpotrs(char* uplo, integer* n, integer* nrhs, const real8* a, integer* lda, - real8* b, integer* ldb, integer* info, char_len); -void dsyevr(char* jobz, char* range, char* uplo, integer* n, real8* a, - integer* lda, real8* vl, real8* vu, integer* il, integer* iu, - real8* abstol, integer* m, real8* w, real8* z, integer* ldz, - integer* isuppz, real8* work, integer* lwork, integer* iwork, - integer* liwork, integer* info, char_len, char_len, char_len); -} - -#define TA_LAPACK_POTRF(...) TA_LAPACK_CALL(potrf, __VA_ARGS__, sizeof(char)) -#define TA_LAPACK_POTRS(...) TA_LAPACK_CALL(potrs, __VA_ARGS__, sizeof(char)) -#define TA_LAPACK_POSV(...) TA_LAPACK_CALL(posv, __VA_ARGS__, sizeof(char)) -#define TA_LAPACK_GESVD(...) \ - TA_LAPACK_CALL(gesvd, __VA_ARGS__, sizeof(char), sizeof(char)) -#define TA_LAPACK_TRTRI(...) \ - TA_LAPACK_CALL( \ - trtri, \ - __VA_ARGS__) // see clapack_fortran.h: does not pass sizeof(char) -#define TA_LAPACK_TRTRS(...) \ - TA_LAPACK_CALL(trtrs, __VA_ARGS__, sizeof(char), sizeof(char), sizeof(char)) -#define TA_LAPACK_SYEV(...) \ - TA_LAPACK_CALL(syev, __VA_ARGS__, sizeof(char), sizeof(char)) -#define TA_LAPACK_SYEVR(...) \ - TA_LAPACK_CALL(syevr, __VA_ARGS__, sizeof(char), sizeof(char), sizeof(char)) -#define TA_LAPACK_SYGV(...) \ - TA_LAPACK_CALL(sygv, __VA_ARGS__, sizeof(char), sizeof(char)) - -#endif // MADNESS_LINALG_USE_LAPACKE - namespace TiledArray::lapack { template From b0c8303200a0ea071b6ca5f687f5c48afe1c58a0 Mon Sep 17 00:00:00 2001 From: asadchev Date: Wed, 2 Dec 2020 11:44:51 -0500 Subject: [PATCH 2/7] Refactor rank-local Matrix, blas, and linalg. - Reorganize rank-local matrix, blas, and linalg - Replace madness::cblas::Transpose by TA::blas::TransposeFlag - Move non-distributed DistArray linalg to TiledArray/algebra/non-distributed --- cmake/toolchains/travis.cmake | 6 - examples/vector_tests/vector.cpp | 2 +- external/eigen.cmake | 18 +-- src/CMakeLists.txt | 38 +++--- src/TiledArray/algebra/cholesky.h | 10 +- src/TiledArray/algebra/diis.h | 1 - src/TiledArray/algebra/heig.h | 6 +- src/TiledArray/algebra/lapack/lapack.cc | 0 src/TiledArray/algebra/lapack/lapack.h | 75 ------------ src/TiledArray/algebra/lu.h | 2 +- .../{lapack => non-distributed}/cholesky.h | 25 ++-- .../{lapack => non-distributed}/heig.h | 22 ++-- .../algebra/{lapack => non-distributed}/lu.h | 22 ++-- .../algebra/{lapack => non-distributed}/svd.h | 23 ++-- .../{lapack => non-distributed}/util.h | 23 ++-- .../{lapack/lapack.cpp => rank-local.cpp} | 11 +- src/TiledArray/algebra/rank-local.h | 50 ++++++++ src/TiledArray/algebra/svd.h | 4 +- src/TiledArray/conversions/eigen.h | 2 +- src/TiledArray/cuda/btas_cublas.h | 2 +- src/TiledArray/{math => cuda}/cublas.h | 0 src/TiledArray/expressions/cont_engine.h | 16 +-- src/TiledArray/expressions/permopt.h | 6 +- src/TiledArray/external/btas.h | 16 ++- src/TiledArray/external/eigen.h | 14 +-- src/TiledArray/initialize.h | 2 +- src/TiledArray/math/blas.h | 53 +++++---- src/TiledArray/math/eigen.h | 112 ------------------ src/TiledArray/math/gemm_helper.h | 23 ++-- src/TiledArray/math/parallel_gemm.h | 6 +- src/TiledArray/proc_grid.h | 1 - src/TiledArray/tensor/kernels.h | 1 - src/TiledArray/tensor/tensor.h | 42 ++++--- tests/linear_algebra.cpp | 82 ++++++------- tests/math_blas.cpp | 12 +- 35 files changed, 292 insertions(+), 436 deletions(-) delete mode 100644 src/TiledArray/algebra/lapack/lapack.cc delete mode 100644 src/TiledArray/algebra/lapack/lapack.h rename src/TiledArray/algebra/{lapack => non-distributed}/cholesky.h (92%) rename src/TiledArray/algebra/{lapack => non-distributed}/heig.h (84%) rename src/TiledArray/algebra/{lapack => non-distributed}/lu.h (76%) rename src/TiledArray/algebra/{lapack => non-distributed}/svd.h (82%) rename src/TiledArray/algebra/{lapack => non-distributed}/util.h (87%) rename src/TiledArray/algebra/{lapack/lapack.cpp => rank-local.cpp} (95%) create mode 100644 src/TiledArray/algebra/rank-local.h rename src/TiledArray/{math => cuda}/cublas.h (100%) delete mode 100644 src/TiledArray/math/eigen.h diff --git a/cmake/toolchains/travis.cmake b/cmake/toolchains/travis.cmake index cf4abd498e..3a3d9f5497 100644 --- a/cmake/toolchains/travis.cmake +++ b/cmake/toolchains/travis.cmake @@ -17,12 +17,6 @@ set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "-O2 -g -Wall" CACHE STRING "Inital C++ relea set(BLAS_LINKER_FLAGS "-L/usr/lib/libblas" "-lblas" "-L/usr/lib/lapack" "-llapack" "-L/usr/lib" "-llapacke" CACHE STRING "BLAS linker flags") set(LAPACK_LIBRARIES ${BLAS_LINKER_FLAGS} CACHE STRING "LAPACK linker flags") set(LAPACK_INCLUDE_DIRS "/usr/include" CACHE STRING "LAPACK include directories") -set( - LAPACK_COMPILE_DEFINITIONS - MADNESS_LINALG_USE_LAPACKE - TILEDARRAY_EIGEN_USE_LAPACKE - CACHE STRING "LAPACK preprocessor definitions" - ) set(INTEGER4 TRUE CACHE BOOL "Set Fortran integer size to 4 bytes") set(BUILD_SHARED_LIBS OFF CACHE BOOL "Build shared libraries") diff --git a/examples/vector_tests/vector.cpp b/examples/vector_tests/vector.cpp index 712af309e0..aa10372819 100644 --- a/examples/vector_tests/vector.cpp +++ b/examples/vector_tests/vector.cpp @@ -391,7 +391,7 @@ int main(int argc, char** argv) { start = madness::wall_time(); for (std::size_t r = 0ul; r < repeat; ++r) { - madness::cblas::scal(n, 3.0, c, 1); + TiledArray::blas::scale(n, 3.0, c); } stop = madness::wall_time(); diff --git a/external/eigen.cmake b/external/eigen.cmake index f62afa7462..f7499711c1 100644 --- a/external/eigen.cmake +++ b/external/eigen.cmake @@ -12,19 +12,6 @@ else(ENABLE_CUDA) set(_tiledarray_required_eigen_version ${TA_TRACKED_EIGEN_VERSION}) endif(ENABLE_CUDA) -set(_tiledarray_eigen_use_lapacke FALSE) -if ("${LAPACK_COMPILE_DEFINITIONS}" MATCHES "TILEDARRAY_EIGEN_USE_LAPACKE") - set(_tiledarray_eigen_use_lapacke TRUE) - if (_tiledarray_required_eigen_version VERSION_LESS 3.3.7) - message( - WARNING - "Eigen3 version => 3.3.7 is required if TILEDARRAY_EIGEN_USE_LAPACKE is set. " - "Prior Eigen3 with LAPACKE enabled may give incorrect eigenvalue results" - ) - set(_tiledarray_required_eigen_version 3.3.7) - endif() -endif () - # Check for existing Eigen # prefer CMake-configured-and-installed instance # re:NO_CMAKE_PACKAGE_REGISTRY: eigen3 registers its *build* tree with the user package registry ... @@ -168,11 +155,8 @@ if (TARGET TiledArray_Eigen) target_compile_definitions(TiledArray_Eigen INTERFACE EIGEN_USE_MKL_ALL) else(MADNESS_HAS_MKL) # Eigen's prototypes for non-MKL (i.e. F77) BLAS interface libraries do not match those in MADNESS (and are not const correct) - # thus can't use non-MKL BLAS, only LAPACK + # thus can't use non-MKL BLAS # target_compile_definitions(TiledArray_Eigen INTERFACE EIGEN_USE_BLAS) - if (_tiledarray_eigen_use_lapacke) - target_compile_definitions(TiledArray_Eigen INTERFACE EIGEN_USE_LAPACKE EIGEN_USE_LAPACKE_STRICT) - endif () endif(MADNESS_HAS_MKL) install(TARGETS TiledArray_Eigen EXPORT tiledarray COMPONENT tiledarray) endif(TARGET TiledArray_Eigen) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 7893ab40ea..21bc79749b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -65,11 +65,6 @@ TiledArray/algebra/heig.h TiledArray/algebra/lu.h TiledArray/algebra/svd.h TiledArray/algebra/types.h -TiledArray/algebra/lapack/cholesky.h -TiledArray/algebra/lapack/heig.h -TiledArray/algebra/lapack/util.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 @@ -120,7 +115,6 @@ TiledArray/expressions/unary_expr.h TiledArray/expressions/index_list.h TiledArray/external/btas.h TiledArray/math/blas.h -TiledArray/math/eigen.h TiledArray/math/gemm_helper.h TiledArray/math/outer.h TiledArray/math/parallel_gemm.h @@ -180,6 +174,10 @@ TiledArray/util/random.h TiledArray/util/singleton.h TiledArray/util/time.h TiledArray/util/vector.h + +TiledArray/algebra/rank-local.h +TiledArray/algebra/rank-local.cpp + ) if(CUDA_FOUND) @@ -187,19 +185,19 @@ if(CUDA_FOUND) list(APPEND TILEDARRAY_HEADER_FILES TiledArray/external/cuda.h TiledArray/external/cutt.h - TiledArray/math/cublas.h - TiledArray/cuda/btas_cublas.h - TiledArray/cuda/btas_um_tensor.h - TiledArray/cuda/cpu_cuda_vector.h - TiledArray/cuda/cuda_task_fn.h - TiledArray/cuda/kernel/mult_kernel.h - TiledArray/cuda/kernel/mult_kernel_impl.h - TiledArray/cuda/kernel/reduce_kernel.h - TiledArray/cuda/kernel/reduce_kernel_impl.h - TiledArray/cuda/platform.h - TiledArray/cuda/thrust.h - TiledArray/cuda/um_allocator.h - TiledArray/cuda/um_storage.h) + TiledArray/cuda/cublas.h + TiledArray/cuda/btas_cublas.h + TiledArray/cuda/btas_um_tensor.h + TiledArray/cuda/cpu_cuda_vector.h + TiledArray/cuda/cuda_task_fn.h + TiledArray/cuda/kernel/mult_kernel.h + TiledArray/cuda/kernel/mult_kernel_impl.h + TiledArray/cuda/kernel/reduce_kernel.h + TiledArray/cuda/kernel/reduce_kernel_impl.h + TiledArray/cuda/platform.h + TiledArray/cuda/thrust.h + TiledArray/cuda/um_allocator.h + TiledArray/cuda/um_storage.h) endif(CUDA_FOUND) @@ -211,7 +209,7 @@ TiledArray/array_impl.cpp TiledArray/dist_array.cpp TiledArray/util/backtrace.cpp TiledArray/util/bug.cpp -TiledArray/algebra/lapack/lapack.cpp +TiledArray/algebra/rank-local.cpp ) # the list of libraries on which TiledArray depends on, will be cached later diff --git a/src/TiledArray/algebra/cholesky.h b/src/TiledArray/algebra/cholesky.h index 5a06ea3821..7f3a4e716d 100644 --- a/src/TiledArray/algebra/cholesky.h +++ b/src/TiledArray/algebra/cholesky.h @@ -28,7 +28,7 @@ #if TILEDARRAY_HAS_SCALAPACK #include #endif -#include +#include namespace TiledArray { @@ -38,7 +38,7 @@ auto cholesky(const Array& A, TiledRange l_trange = TiledRange()) { if (A.world().size() > 1 && A.range().volume() > 10000000) return scalapack::cholesky(A, l_trange); #endif - return lapack::cholesky(A, l_trange); + return non_distributed::cholesky(A, l_trange); } template @@ -47,7 +47,7 @@ auto cholesky_linv(const Array& A, TiledRange l_trange = TiledRange()) { if (A.world().size() > 1 && A.range().volume() > 10000000) return scalapack::cholesky_linv(A, l_trange); #endif - return lapack::cholesky_linv(A, l_trange); + return non_distributed::cholesky_linv(A, l_trange); } template @@ -57,7 +57,7 @@ auto cholesky_solve(const Array& A, const Array& B, if (A.world().size() > 1 && A.range().volume() > 10000000) return scalapack::cholesky_solve(A, B, x_trange); #endif - return lapack::cholesky_solve(A, B, x_trange); + return non_distributed::cholesky_solve(A, B, x_trange); } template @@ -69,7 +69,7 @@ auto cholesky_lsolve(TransposeFlag transpose, const Array& A, const Array& B, return scalapack::cholesky_lsolve(transpose, A, B, l_trange, x_trange); #endif - return lapack::cholesky_lsolve(transpose, A, B, l_trange, x_trange); + return non_distributed::cholesky_lsolve(transpose, A, B, l_trange, x_trange); } } // namespace TiledArray diff --git a/src/TiledArray/algebra/diis.h b/src/TiledArray/algebra/diis.h index 537b39b44e..5b50ede9f3 100644 --- a/src/TiledArray/algebra/diis.h +++ b/src/TiledArray/algebra/diis.h @@ -27,7 +27,6 @@ #define TILEDARRAY_ALGEBRA_DIIS_H__INCLUDED #include -#include #include #include #include "../dist_array.h" diff --git a/src/TiledArray/algebra/heig.h b/src/TiledArray/algebra/heig.h index 68e39c9de9..23d28329bc 100644 --- a/src/TiledArray/algebra/heig.h +++ b/src/TiledArray/algebra/heig.h @@ -28,7 +28,7 @@ #if TILEDARRAY_HAS_SCALAPACK #include #endif -#include +#include namespace TiledArray { @@ -39,7 +39,7 @@ auto heig(const Array& A, TiledRange evec_trange = TiledRange()) { return scalapack::heig(A, evec_trange); } #endif - return lapack::heig(A, evec_trange); + return non_distributed::heig(A, evec_trange); } template @@ -49,7 +49,7 @@ auto heig(const ArrayA& A, const ArrayB& B, TiledRange evec_trange = TiledRange( return scalapack::heig(A, B, evec_trange); } #endif - return lapack::heig(A, B, evec_trange); + return non_distributed::heig(A, B, evec_trange); } } // namespace TiledArray diff --git a/src/TiledArray/algebra/lapack/lapack.cc b/src/TiledArray/algebra/lapack/lapack.cc deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/src/TiledArray/algebra/lapack/lapack.h b/src/TiledArray/algebra/lapack/lapack.h deleted file mode 100644 index f3d2b47f0e..0000000000 --- a/src/TiledArray/algebra/lapack/lapack.h +++ /dev/null @@ -1,75 +0,0 @@ -/* - * This file is a part of TiledArray. - * Copyright (C) 2020 Virginia Tech - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - * - * Eduard Valeyev - * - * lapack.h - * Created: 16 October, 2020 - * - */ -#ifndef TILEDARRAY_ALGEBRA_LAPACK_LAPACK_H__INCLUDED -#define TILEDARRAY_ALGEBRA_LAPACK_LAPACK_H__INCLUDED - -#include -#include -#include - -namespace TiledArray::lapack { - -template -struct array_traits { - using scalar_type = typename A::scalar_type; - using numeric_type = typename A::numeric_type; - static const bool complex = !std::is_same_v; - static_assert(std::is_same_v, - "TA::lapack is only usable with a DistArray of scalar types"); -}; - -template -using Matrix = - Eigen::Matrix; - -template -void cholesky(Matrix &A); - -template -void cholesky_linv(Matrix &A); - -template -void cholesky_solve(Matrix &A, Matrix &X); - -template -void cholesky_lsolve(TransposeFlag transpose, Matrix &A, Matrix &X); - -template -void heig(Matrix &A, std::vector &W); - -template -void heig(Matrix &A, Matrix &B, std::vector &W); - -template -void svd(Matrix &A, std::vector &S, Matrix *U, Matrix *VT); - -template -void lu_solve(Matrix &A, Matrix &B); - -template -void lu_inv(Matrix &A); - -} // namespace TiledArray::lapack - -#endif // TILEDARRAY_ALGEBRA_LAPACK_LAPACK_H__INCLUDED diff --git a/src/TiledArray/algebra/lu.h b/src/TiledArray/algebra/lu.h index 330f49c8f7..bedb2df59a 100644 --- a/src/TiledArray/algebra/lu.h +++ b/src/TiledArray/algebra/lu.h @@ -28,7 +28,7 @@ #if TILEDARRAY_HAS_SCALAPACK #include #endif -#include +#include namespace TiledArray { diff --git a/src/TiledArray/algebra/lapack/cholesky.h b/src/TiledArray/algebra/non-distributed/cholesky.h similarity index 92% rename from src/TiledArray/algebra/lapack/cholesky.h rename to src/TiledArray/algebra/non-distributed/cholesky.h index 05af51fc81..037685430a 100644 --- a/src/TiledArray/algebra/lapack/cholesky.h +++ b/src/TiledArray/algebra/non-distributed/cholesky.h @@ -21,16 +21,16 @@ * Created: 16 October, 2020 * */ -#ifndef TILEDARRAY_ALGEBRA_LAPACK_CHOL_H__INCLUDED -#define TILEDARRAY_ALGEBRA_LAPACK_CHOL_H__INCLUDED +#ifndef TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_CHOL_H__INCLUDED +#define TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_CHOL_H__INCLUDED -#include -#include #include + +#include #include +#include -namespace TiledArray { -namespace lapack { +namespace TiledArray::non_distributed { namespace detail { @@ -45,7 +45,7 @@ auto make_L_eig(const DistArray& A) { World& world = A.world(); auto A_eig = detail::to_eigen(A); if (world.rank() == 0) { - lapack::cholesky(A_eig); + algebra::rank_local::cholesky(A_eig); } world.gop.broadcast_serializable(A_eig, 0); return A_eig; @@ -102,7 +102,7 @@ template >> auto cholesky(const ContiguousTensor& A) { auto A_eig = detail::to_eigen(A); - lapack::cholesky(A_eig); + algebra::rank_local::cholesky(A_eig); detail::zero_out_upper_triangle(A_eig); return detail::from_eigen(A_eig, A.range()); } @@ -141,7 +141,7 @@ auto cholesky_linv(const Array& A, TiledRange l_trange = TiledRange()) { if (world.rank() == 0) { if (RetL) L_inv_eig = L_eig; auto& L_inv_eig_ref = RetL ? L_inv_eig : L_eig; - cholesky_linv(L_inv_eig_ref); + algebra::rank_local::cholesky_linv(L_inv_eig_ref); detail::zero_out_upper_triangle(L_inv_eig_ref); } world.gop.broadcast_serializable(RetL ? L_inv_eig : L_eig, 0); @@ -167,7 +167,7 @@ auto cholesky_solve(const Array& A, const Array& B, auto X_eig = detail::to_eigen(B); World& world = A.world(); if (world.rank() == 0) { - cholesky_solve(A_eig, X_eig); + algebra::rank_local::cholesky_solve(A_eig, X_eig); } world.gop.broadcast_serializable(X_eig, 0); if (x_trange.rank() == 0) x_trange = B.trange(); @@ -190,7 +190,7 @@ auto cholesky_lsolve(TransposeFlag transpose, const Array& A, const Array& B, auto X_eig = detail::to_eigen(B); if (world.rank() == 0) { - cholesky_lsolve(transpose, L_eig, X_eig); + algebra::rank_local::cholesky_lsolve(transpose, L_eig, X_eig); } world.gop.broadcast_serializable(X_eig, 0); if (l_trange.rank() == 0) l_trange = A.trange(); @@ -199,7 +199,6 @@ auto cholesky_lsolve(TransposeFlag transpose, const Array& A, const Array& B, eigen_to_array(world, x_trange, X_eig)); } -} // namespace lapack } // namespace TiledArray -#endif // TILEDARRAY_ALGEBRA_LAPACK_CHOL_H__INCLUDED +#endif // TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_CHOL_H__INCLUDED diff --git a/src/TiledArray/algebra/lapack/heig.h b/src/TiledArray/algebra/non-distributed/heig.h similarity index 84% rename from src/TiledArray/algebra/lapack/heig.h rename to src/TiledArray/algebra/non-distributed/heig.h index 5fe4270c8f..484f0bdb20 100644 --- a/src/TiledArray/algebra/lapack/heig.h +++ b/src/TiledArray/algebra/non-distributed/heig.h @@ -21,16 +21,16 @@ * Created: 19 October, 2020 * */ -#ifndef TILEDARRAY_ALGEBRA_LAPACK_HEIG_H__INCLUDED -#define TILEDARRAY_ALGEBRA_LAPACK_HEIG_H__INCLUDED +#ifndef TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_HEIG_H__INCLUDED +#define TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_HEIG_H__INCLUDED #include -#include -#include +#include +#include #include -namespace TiledArray::lapack { +namespace TiledArray::non_distributed { /** * @brief Solve the standard eigenvalue problem with LAPACK @@ -52,12 +52,12 @@ namespace TiledArray::lapack { */ template auto heig(const Array& A, TiledRange evec_trange = TiledRange()) { - using numeric_type = typename lapack::array_traits::numeric_type; + using numeric_type = typename detail::array_traits::numeric_type; World& world = A.world(); auto A_eig = detail::to_eigen(A); std::vector evals; if (world.rank() == 0) { - lapack::heig(A_eig, evals); + algebra::rank_local::heig(A_eig, evals); } world.gop.broadcast_serializable(A_eig, 0); world.gop.broadcast_serializable(evals, 0); @@ -93,14 +93,14 @@ auto heig(const Array& A, TiledRange evec_trange = TiledRange()) { */ template auto heig(const ArrayA& A, const ArrayB& B, TiledRange evec_trange = TiledRange()) { - using numeric_type = typename lapack::array_traits::numeric_type; - (void)lapack::array_traits{}; + using numeric_type = typename detail::array_traits::numeric_type; + (void)detail::array_traits{}; World& world = A.world(); auto A_eig = detail::to_eigen(A); auto B_eig = detail::to_eigen(B); std::vector evals; if (world.rank() == 0) { - lapack::heig(A_eig, B_eig, evals); + algebra::rank_local::heig(A_eig, B_eig, evals); } world.gop.broadcast_serializable(A_eig, 0); world.gop.broadcast_serializable(evals, 0); @@ -113,4 +113,4 @@ auto heig(const ArrayA& A, const ArrayB& B, TiledRange evec_trange = TiledRange( } // namespace TiledArray::lapack -#endif // TILEDARRAY_ALGEBRA_LAPACK_HEIG_H__INCLUDED +#endif // TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_HEIG_H__INCLUDED diff --git a/src/TiledArray/algebra/lapack/lu.h b/src/TiledArray/algebra/non-distributed/lu.h similarity index 76% rename from src/TiledArray/algebra/lapack/lu.h rename to src/TiledArray/algebra/non-distributed/lu.h index 389af28942..76beba06ac 100644 --- a/src/TiledArray/algebra/lapack/lu.h +++ b/src/TiledArray/algebra/non-distributed/lu.h @@ -22,27 +22,29 @@ * Created: 19 June, 2020 * */ -#ifndef TILEDARRAY_ALGEBRA_LAPACK_LU_H__INCLUDED -#define TILEDARRAY_ALGEBRA_LAPACK_LU_H__INCLUDED +#ifndef TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_LU_H__INCLUDED +#define TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_LU_H__INCLUDED #include -#include +#include +#include +#include -namespace TiledArray::lapack { +namespace TiledArray::non_distributed { /** * @brief Solve a linear system via LU factorization */ template auto lu_solve(const ArrayA& A, const ArrayB& B, TiledRange x_trange = TiledRange()) { - (void)lapack::array_traits{}; - (void)lapack::array_traits{}; + (void)detail::array_traits{}; + (void)detail::array_traits{}; auto& world = A.world(); auto A_eig = detail::to_eigen(A); auto B_eig = detail::to_eigen(B); if (world.rank() == 0) { - lapack::lu_solve(A_eig, B_eig); + algebra::rank_local::lu_solve(A_eig, B_eig); } world.gop.broadcast_serializable(B_eig, 0); if (x_trange.rank() == 0) x_trange = B.trange(); @@ -54,11 +56,11 @@ auto lu_solve(const ArrayA& A, const ArrayB& B, TiledRange x_trange = TiledRange */ template auto lu_inv(const Array& A, TiledRange ainv_trange = TiledRange()) { - (void)lapack::array_traits{}; + (void)detail::array_traits{}; auto& world = A.world(); auto A_eig = detail::to_eigen(A); if (world.rank() == 0) { - lapack::lu_inv(A_eig); + algebra::rank_local::lu_inv(A_eig); } world.gop.broadcast_serializable(A_eig, 0); if (ainv_trange.rank() == 0) ainv_trange = A.trange(); @@ -67,4 +69,4 @@ auto lu_inv(const Array& A, TiledRange ainv_trange = TiledRange()) { } // namespace TiledArray::lapack -#endif // TILEDARRAY_ALGEBRA_LAPACK_LU_H__INCLUDED +#endif // TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_LU_H__INCLUDED diff --git a/src/TiledArray/algebra/lapack/svd.h b/src/TiledArray/algebra/non-distributed/svd.h similarity index 82% rename from src/TiledArray/algebra/lapack/svd.h rename to src/TiledArray/algebra/non-distributed/svd.h index 8ea0afd3da..01f3ef8eee 100644 --- a/src/TiledArray/algebra/lapack/svd.h +++ b/src/TiledArray/algebra/non-distributed/svd.h @@ -22,14 +22,16 @@ * Created: 12 June, 2020 * */ -#ifndef TILEDARRAY_ALGEBRA_LAPACK_SVD_H__INCLUDED -#define TILEDARRAY_ALGEBRA_LAPACK_SVD_H__INCLUDED +#ifndef TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_SVD_H__INCLUDED +#define TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_SVD_H__INCLUDED #include -#include +#include +#include +#include -namespace TiledArray::lapack { +namespace TiledArray::non_distributed { /** * @brief Compute the singular value decomposition (SVD) via ScaLAPACK @@ -58,6 +60,7 @@ template ; World& world = A.world(); auto A_eig = detail::to_eigen(A); @@ -67,13 +70,13 @@ auto svd(const Array& A, TiledRange u_trange = TiledRange(), TiledRange vt_trang constexpr bool need_vt = std::is_same_v or svd_all_vectors; std::vector S; - std::unique_ptr< Matrix > U, VT; + std::unique_ptr U, VT; - if constexpr (need_u) U = std::make_unique< Matrix >(); - if constexpr (need_vt) VT = std::make_unique< Matrix >(); + if constexpr (need_u) U = std::make_unique(); + if constexpr (need_vt) VT = std::make_unique(); if (world.rank() == 0) { - lapack::svd(A_eig, S, U.get(), VT.get()); + algebra::rank_local::svd(A_eig, S, U.get(), VT.get()); } world.gop.broadcast_serializable(S, 0); @@ -98,6 +101,6 @@ auto svd(const Array& A, TiledRange u_trange = TiledRange(), TiledRange vt_trang } -} // namespace scalapack::TiledArray +} // namespace TiledArray::non_distributed -#endif // TILEDARRAY_ALGEBRA_LAPACK_SVD_H__INCLUDED +#endif // TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_SVD_H__INCLUDED diff --git a/src/TiledArray/algebra/lapack/util.h b/src/TiledArray/algebra/non-distributed/util.h similarity index 87% rename from src/TiledArray/algebra/lapack/util.h rename to src/TiledArray/algebra/non-distributed/util.h index daf822ff8e..ec72ce963d 100644 --- a/src/TiledArray/algebra/lapack/util.h +++ b/src/TiledArray/algebra/non-distributed/util.h @@ -21,15 +21,23 @@ * Created: 19 October, 2020 * */ -#ifndef TILEDARRAY_ALGEBRA_LAPACK_UTIL_H__INCLUDED -#define TILEDARRAY_ALGEBRA_LAPACK_UTIL_H__INCLUDED +#ifndef TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_UTIL_H__INCLUDED +#define TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_UTIL_H__INCLUDED #include #include -namespace TiledArray { -namespace lapack { -namespace detail { +namespace TiledArray::non_distributed::detail { + +template +struct array_traits { + using scalar_type = typename A::scalar_type; + using numeric_type = typename A::numeric_type; + static const bool complex = !std::is_same_v; + static_assert(std::is_same_v, + "TA::linalg is only usable with a DistArray of scalar types"); +}; + template auto to_eigen(const DistArray& A) { @@ -101,9 +109,6 @@ void zero_out_upper_triangle(Eigen::MatrixBase& A) { A.template triangularView().setZero(); } -} // namespace detail - -} // namespace lapack } // namespace TiledArray -#endif // TILEDARRAY_ALGEBRA_LAPACK_UTIL_H__INCLUDED +#endif // TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_UTIL_H__INCLUDED diff --git a/src/TiledArray/algebra/lapack/lapack.cpp b/src/TiledArray/algebra/rank-local.cpp similarity index 95% rename from src/TiledArray/algebra/lapack/lapack.cpp rename to src/TiledArray/algebra/rank-local.cpp index 014f08f209..20fd9dcaa0 100644 --- a/src/TiledArray/algebra/lapack/lapack.cpp +++ b/src/TiledArray/algebra/rank-local.cpp @@ -22,8 +22,7 @@ * */ -#include -#include +#include #include @@ -47,7 +46,7 @@ else std::abort(); \ } -namespace TiledArray::lapack { +namespace TiledArray::algebra::rank_local { template void cholesky(Matrix& A) { @@ -213,7 +212,7 @@ void lu_inv(Matrix& A) { template void lu_solve(MATRIX&, MATRIX&); \ template void lu_inv(MATRIX&); -TA_LAPACK_EXPLICIT(lapack::Matrix, std::vector); -TA_LAPACK_EXPLICIT(lapack::Matrix, std::vector); +TA_LAPACK_EXPLICIT(Matrix, std::vector); +TA_LAPACK_EXPLICIT(Matrix, std::vector); -} // namespace TiledArray::lapack +} // namespace TiledArray::local diff --git a/src/TiledArray/algebra/rank-local.h b/src/TiledArray/algebra/rank-local.h new file mode 100644 index 0000000000..a69beba83d --- /dev/null +++ b/src/TiledArray/algebra/rank-local.h @@ -0,0 +1,50 @@ +#ifndef TILEDARRAY_ALGEBRA_RANK_LOCAL_H__INCLUDED +#define TILEDARRAY_ALGEBRA_RANK_LOCAL_H__INCLUDED + +#include + +#include +#include + +#include + +namespace TiledArray::algebra::rank_local { + +template +using Matrix = ::Eigen::Matrix; + +// template +// using Vector = ::Eigen::Matrix; + +template +void cholesky(Matrix &A); + +template +void cholesky_linv(Matrix &A); + +template +void cholesky_solve(Matrix &A, Matrix &X); + +template +void cholesky_lsolve(TransposeFlag transpose, + Matrix &A, Matrix &X); + +template +void heig(Matrix &A, std::vector &W); + +template +void heig(Matrix &A, Matrix &B, std::vector &W); + +template +void svd(Matrix &A, std::vector &S, + Matrix *U, Matrix *VT); + +template +void lu_solve(Matrix &A, Matrix &B); + +template +void lu_inv(Matrix &A); + +} // namespace TiledArray::local + +#endif // TILEDARRAY_ALGEBRA_RANK_LOCAL_H__INCLUDED diff --git a/src/TiledArray/algebra/svd.h b/src/TiledArray/algebra/svd.h index 3c87fbc087..bd10a289a4 100644 --- a/src/TiledArray/algebra/svd.h +++ b/src/TiledArray/algebra/svd.h @@ -28,7 +28,7 @@ #ifdef TILEDARRAY_HAS_SCALAPACK #include #endif // TILEDARRAY_HAS_SCALAPACK -#include +#include namespace TiledArray { @@ -40,7 +40,7 @@ auto svd(const Array& A, TiledRange u_trange = TiledRange(), TiledRange vt_trang return scalapack::svd(A, u_trange, vt_trange); } #endif - return lapack::svd(A, u_trange, vt_trange); + return non_distributed::svd(A, u_trange, vt_trange); } } // namespace TiledArray diff --git a/src/TiledArray/conversions/eigen.h b/src/TiledArray/conversions/eigen.h index b43602c98e..58e9db51f0 100644 --- a/src/TiledArray/conversions/eigen.h +++ b/src/TiledArray/conversions/eigen.h @@ -28,7 +28,7 @@ #include #include -#include +#include #include #include #include diff --git a/src/TiledArray/cuda/btas_cublas.h b/src/TiledArray/cuda/btas_cublas.h index 5ff6837588..69227c8b88 100644 --- a/src/TiledArray/cuda/btas_cublas.h +++ b/src/TiledArray/cuda/btas_cublas.h @@ -25,7 +25,7 @@ #define TILEDARRAY_BTAS_CUDA_CUBLAS_H__INCLUDED #include -#include +#include #ifdef TILEDARRAY_HAS_CUDA diff --git a/src/TiledArray/math/cublas.h b/src/TiledArray/cuda/cublas.h similarity index 100% rename from src/TiledArray/math/cublas.h rename to src/TiledArray/cuda/cublas.h diff --git a/src/TiledArray/expressions/cont_engine.h b/src/TiledArray/expressions/cont_engine.h index 5237fc3883..fedce43711 100644 --- a/src/TiledArray/expressions/cont_engine.h +++ b/src/TiledArray/expressions/cont_engine.h @@ -250,14 +250,14 @@ class ContEngine : public BinaryEngine { // Initialize the tile operation in this function because it is used to // evaluate the tiled range and shape. - const madness::cblas::CBLAS_TRANSPOSE left_op = + const blas::TransposeFlag left_op = (left_outer_permtype_ == PermutationType::matrix_transpose - ? madness::cblas::Trans - : madness::cblas::NoTrans); - const madness::cblas::CBLAS_TRANSPOSE right_op = + ? blas::Transpose + : blas::NoTranspose); + const blas::TransposeFlag right_op = (right_outer_permtype_ == PermutationType::matrix_transpose - ? madness::cblas::Trans - : madness::cblas::NoTrans); + ? blas::Transpose + : blas::NoTranspose); if (outer(target_indices) != outer(indices_)) { // Initialize permuted structure @@ -397,7 +397,7 @@ class ContEngine : public BinaryEngine { /// \return The result shape shape_type make_shape() const { const TiledArray::math::GemmHelper shape_gemm_helper( - madness::cblas::NoTrans, madness::cblas::NoTrans, + blas::NoTranspose, blas::NoTranspose, op_.gemm_helper().result_rank(), op_.gemm_helper().left_rank(), op_.gemm_helper().right_rank()); return left_.shape().gemm(right_.shape(), factor_, shape_gemm_helper); @@ -409,7 +409,7 @@ class ContEngine : public BinaryEngine { /// \return The result shape shape_type make_shape(const Permutation& perm) const { const TiledArray::math::GemmHelper shape_gemm_helper( - madness::cblas::NoTrans, madness::cblas::NoTrans, + blas::NoTranspose, blas::NoTranspose, op_.gemm_helper().result_rank(), op_.gemm_helper().left_rank(), op_.gemm_helper().right_rank()); return left_.shape().gemm(right_.shape(), factor_, shape_gemm_helper, perm); diff --git a/src/TiledArray/expressions/permopt.h b/src/TiledArray/expressions/permopt.h index 01e1152f5d..ef5a58e17d 100644 --- a/src/TiledArray/expressions/permopt.h +++ b/src/TiledArray/expressions/permopt.h @@ -43,12 +43,12 @@ namespace expressions { // clang-format on enum class PermutationType { identity = 1, matrix_transpose = 2, general = 3 }; -inline madness::cblas::CBLAS_TRANSPOSE to_cblas_op(PermutationType permtype) { +inline blas::TransposeFlag to_cblas_op(PermutationType permtype) { TA_ASSERT(permtype == PermutationType::matrix_transpose || permtype == PermutationType::identity); return permtype == PermutationType::matrix_transpose - ? madness::cblas::Trans - : madness::cblas::NoTrans; + ? blas::Transpose + : blas::NoTranspose; } /// Abstract optimizer of permutations for a binary operation diff --git a/src/TiledArray/external/btas.h b/src/TiledArray/external/btas.h index acd3465c78..8a5c7497ef 100644 --- a/src/TiledArray/external/btas.h +++ b/src/TiledArray/external/btas.h @@ -632,9 +632,11 @@ inline btas::Tensor gemm( T factor_t(factor); - TiledArray::math::gemm(gemm_helper.left_op(), gemm_helper.right_op(), m, n, k, - factor_t, left.data(), lda, right.data(), ldb, T(0), - result.data(), n); + TiledArray::blas::gemm( + gemm_helper.left_op(), gemm_helper.right_op(), m, n, k, + factor_t, left.data(), lda, right.data(), ldb, T(0), + result.data(), n + ); return result; } @@ -706,9 +708,11 @@ inline void gemm(btas::Tensor& result, T factor_t(factor); - TiledArray::math::gemm(gemm_helper.left_op(), gemm_helper.right_op(), m, n, k, - factor_t, left.data(), lda, right.data(), ldb, T(1), - result.data(), n); + TiledArray::blas::gemm( + gemm_helper.left_op(), gemm_helper.right_op(), m, n, k, + factor_t, left.data(), lda, right.data(), ldb, T(1), + result.data(), n + ); } // sum of the hyperdiagonal elements diff --git a/src/TiledArray/external/eigen.h b/src/TiledArray/external/eigen.h index 90e8196c2d..8bf9dc031c 100644 --- a/src/TiledArray/external/eigen.h +++ b/src/TiledArray/external/eigen.h @@ -48,21 +48,10 @@ TILEDARRAY_PRAGMA_GCC(system_header) #define EIGEN_USE_MKL_ALL 1 #endif #else - //# ifndef EIGEN_USE_BLAS //# define EIGEN_USE_BLAS 1 //# endif - -#ifdef TILEDARRAY_EIGEN_USE_LAPACKE -#ifndef EIGEN_USE_LAPACKE -#define EIGEN_USE_LAPACKE 1 -#endif -#ifndef EIGEN_USE_LAPACKE_STRICT -#define EIGEN_USE_LAPACKE_STRICT 1 -#endif - -#endif -#endif +#endif // HAVE_INTEL_MKL ///////////////////////////////////////////////// // define lapacke types to prevent inclusion of complex.h by @@ -75,6 +64,7 @@ TILEDARRAY_PRAGMA_GCC(system_header) #if defined(EIGEN_USE_LAPACKE_STRICT) && defined(EIGEN_USE_LAPACKE) #undef EIGEN_USE_LAPACKE #endif + #include #if defined(EIGEN_USE_LAPACKE) || defined(EIGEN_USE_LAPACKE_STRICT) diff --git a/src/TiledArray/initialize.h b/src/TiledArray/initialize.h index 001d34047d..f1392aa21f 100644 --- a/src/TiledArray/initialize.h +++ b/src/TiledArray/initialize.h @@ -10,7 +10,7 @@ #include #ifdef TILEDARRAY_HAS_CUDA #include -#include +#include #include #endif #ifdef HAVE_INTEL_MKL diff --git a/src/TiledArray/math/blas.h b/src/TiledArray/math/blas.h index 7969790724..d33fb97be8 100644 --- a/src/TiledArray/math/blas.h +++ b/src/TiledArray/math/blas.h @@ -23,21 +23,31 @@ * */ -#ifndef TILEDARRAY_BLAS_H__INCLUDED -#define TILEDARRAY_BLAS_H__INCLUDED +#ifndef TILEDARRAY_MATH_BLAS_H__INCLUDED +#define TILEDARRAY_MATH_BLAS_H__INCLUDED -#include +#include #include #include -namespace TiledArray { -namespace math { +namespace TiledArray::math::blas { + +using TransposeFlag = madness::cblas::CBLAS_TRANSPOSE; +static constexpr auto NoTranspose = madness::cblas::NoTrans; +static constexpr auto Transpose = madness::cblas::Trans; +static constexpr auto ConjTranspose = madness::cblas::ConjTrans; + +template +using Matrix = ::Eigen::Matrix; + +template +using Vector = ::Eigen::Matrix; // BLAS _GEMM wrapper functions template -inline void gemm(madness::cblas::CBLAS_TRANSPOSE op_a, - madness::cblas::CBLAS_TRANSPOSE op_b, const integer m, +inline void gemm(TransposeFlag op_a, + TransposeFlag op_b, const integer m, const integer n, const integer k, const S1 alpha, const T1* a, const integer lda, const T2* b, const integer ldb, const S2 beta, T3* c, const integer ldc) { @@ -129,8 +139,8 @@ inline void gemm(madness::cblas::CBLAS_TRANSPOSE op_a, } } -inline void gemm(madness::cblas::CBLAS_TRANSPOSE op_a, - madness::cblas::CBLAS_TRANSPOSE op_b, const integer m, +inline void gemm(TransposeFlag op_a, + TransposeFlag op_b, const integer m, const integer n, const integer k, const float alpha, const float* a, const integer lda, const float* b, const integer ldb, const float beta, float* c, @@ -139,8 +149,8 @@ inline void gemm(madness::cblas::CBLAS_TRANSPOSE op_a, ldc); } -inline void gemm(madness::cblas::CBLAS_TRANSPOSE op_a, - madness::cblas::CBLAS_TRANSPOSE op_b, const integer m, +inline void gemm(TransposeFlag op_a, + TransposeFlag op_b, const integer m, const integer n, const integer k, const double alpha, const double* a, const integer lda, const double* b, const integer ldb, const double beta, double* c, @@ -149,8 +159,8 @@ inline void gemm(madness::cblas::CBLAS_TRANSPOSE op_a, ldc); } -inline void gemm(madness::cblas::CBLAS_TRANSPOSE op_a, - madness::cblas::CBLAS_TRANSPOSE op_b, const integer m, +inline void gemm(TransposeFlag op_a, + TransposeFlag op_b, const integer m, const integer n, const integer k, const std::complex alpha, const std::complex* a, const integer lda, const std::complex* b, @@ -160,8 +170,8 @@ inline void gemm(madness::cblas::CBLAS_TRANSPOSE op_a, ldc); } -inline void gemm(madness::cblas::CBLAS_TRANSPOSE op_a, - madness::cblas::CBLAS_TRANSPOSE op_b, const integer m, +inline void gemm(TransposeFlag op_a, + TransposeFlag op_b, const integer m, const integer n, const integer k, const std::complex alpha, const std::complex* a, const integer lda, @@ -177,7 +187,7 @@ inline void gemm(madness::cblas::CBLAS_TRANSPOSE op_a, template inline typename std::enable_if>::type scale( const integer n, const T alpha, U* x) { - eigen_map(x, n) *= alpha; + Vector::Map(x, n) *= alpha; } inline void scale(const integer n, const float alpha, float* x) { @@ -211,7 +221,7 @@ inline void scale(const integer n, const double alpha, template T dot(const integer n, const T* x, const U* y) { - return eigen_map(x, n).dot(eigen_map(y, n)); + return Vector::Map(x, n).dot(Vector::Map(y, n)); } inline float dot(integer n, const float* x, const float* y) { @@ -235,7 +245,10 @@ inline std::complex dot(integer n, const std::complex* x, // Import the madness dot functions into the TiledArray namespace using madness::cblas::dot; -} // namespace math -} // namespace TiledArray +} // namespace TiledArray::local + +namespace TiledArray { + namespace blas = TiledArray::math::blas; +} -#endif // TILEDARRAY_BLAS_H__INCLUDED +#endif // TILEDARRAY_MATH_BLAS_H__INCLUDED diff --git a/src/TiledArray/math/eigen.h b/src/TiledArray/math/eigen.h deleted file mode 100644 index 84e941c4df..0000000000 --- a/src/TiledArray/math/eigen.h +++ /dev/null @@ -1,112 +0,0 @@ -/* - * This file is a part of TiledArray. - * Copyright (C) 2013 Virginia Tech - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - * - * Justus Calvin - * Department of Chemistry, Virginia Tech - * - * eigen.h - * Sep 16, 2013 - * - */ - -#ifndef TILEDARRAY_MATH_EIGEN_H__INCLUDED -#define TILEDARRAY_MATH_EIGEN_H__INCLUDED - -#include -#include - -namespace TiledArray { -namespace math { - -/// Construct a const Eigen::Map object for a given Tensor object - -/// \tparam T The element type -/// \param t The buffer pointer -/// \param m The number of rows in the result matrix -/// \param n The number of columns in the result matrix -/// \return An m x n Eigen matrix map for \c tensor -/// \throw TiledArray::Exception When m * n is not equal to \c tensor size -template -inline Eigen::Map< - const Eigen::Matrix, - Eigen::AutoAlign> -eigen_map(const T* t, const std::size_t m, const std::size_t n) { - TA_ASSERT(t); - TA_ASSERT(m > 0); - TA_ASSERT(n > 0); - return Eigen::Map< - const Eigen::Matrix, - Eigen::AutoAlign>(t, m, n); -} - -/// Construct an Eigen::Map object for a given Tensor object - -/// \tparam T The tensor element type -/// \param t The tensor object -/// \param m The number of rows in the result matrix -/// \param n The number of columns in the result matrix -/// \return An m x n Eigen matrix map for \c tensor -/// \throw TiledArray::Exception When m * n is not equal to \c tensor size -template -inline Eigen::Map< - Eigen::Matrix, - Eigen::AutoAlign> -eigen_map(T* t, const std::size_t m, const std::size_t n) { - TA_ASSERT(t); - TA_ASSERT(m > 0); - TA_ASSERT(n > 0); - return Eigen::Map< - Eigen::Matrix, - Eigen::AutoAlign>(t, m, n); -} - -/// Construct a const Eigen::Map object for a given Tensor object - -/// \tparam T The element type -/// \param t The vector pointer -/// \param n The number of elements in the result matrix -/// \return An n element Eigen vector map for \c tensor -/// \throw TiledArray::Exception When n is not equal to \c tensor size -template -inline Eigen::Map, Eigen::AutoAlign> -eigen_map(const T* t, const std::size_t n) { - TA_ASSERT(t); - TA_ASSERT(n > 0); - return Eigen::Map, - Eigen::AutoAlign>(t, n); -} - -/// Construct an Eigen::Map object for a given Tensor object - -/// \tparam T The element type -/// \param t The vector pointer -/// \param n The number of elements in the result matrix -/// \return An n element Eigen vector map for \c tensor -/// \throw TiledArray::Exception When m * n is not equal to \c tensor size -template -inline Eigen::Map, Eigen::AutoAlign> -eigen_map(T* t, const std::size_t n) { - TA_ASSERT(t); - TA_ASSERT(n > 0); - return Eigen::Map, Eigen::AutoAlign>(t, - n); -} - -} // namespace math -} // namespace TiledArray - -#endif // TILEDARRAY_MATH_EIGEN_H__INCLUDED diff --git a/src/TiledArray/math/gemm_helper.h b/src/TiledArray/math/gemm_helper.h index fe947448f5..4eb3e256f0 100644 --- a/src/TiledArray/math/gemm_helper.h +++ b/src/TiledArray/math/gemm_helper.h @@ -28,9 +28,9 @@ #include #include +#include -namespace TiledArray { -namespace math { +namespace TiledArray::math { /// Contraction to *GEMM helper @@ -38,9 +38,9 @@ namespace math { /// providing information on how to fuse dimensions class GemmHelper { private: - madness::cblas::CBLAS_TRANSPOSE left_op_; + blas::TransposeFlag left_op_; ///< Transpose operation that is applied to the left-hand argument - madness::cblas::CBLAS_TRANSPOSE right_op_; + blas::TransposeFlag right_op_; ///< Transpose operation that is applied to the right-hand argument unsigned int result_rank_; ///< The rank of the result tensor @@ -57,8 +57,8 @@ class GemmHelper { right_; ///< Right-hand argument range data public: - GemmHelper(const madness::cblas::CBLAS_TRANSPOSE left_op, - const madness::cblas::CBLAS_TRANSPOSE right_op, + GemmHelper(const blas::TransposeFlag left_op, + const blas::TransposeFlag right_op, const unsigned int result_rank, const unsigned int left_rank, const unsigned int right_rank) : left_op_(left_op), @@ -74,7 +74,7 @@ class GemmHelper { const unsigned int contract_size = num_contract_ranks(); // Store the inner and outer dimension ranges for the left-hand argument. - if (left_op == madness::cblas::NoTrans) { + if (left_op == blas::NoTranspose) { left_.outer[0] = 0u; left_.outer[1] = left_.inner[0] = left_rank - contract_size; left_.inner[1] = left_rank; @@ -85,7 +85,7 @@ class GemmHelper { } // Store the inner and outer dimension ranges for the right-hand argument. - if (right_op == madness::cblas::NoTrans) { + if (right_op == blas::NoTranspose) { right_.inner[0] = 0u; right_.inner[1] = right_.outer[0] = contract_size; right_.outer[1] = right_rank; @@ -271,11 +271,10 @@ class GemmHelper { n *= right_extent[i]; } - madness::cblas::CBLAS_TRANSPOSE left_op() const { return left_op_; } - madness::cblas::CBLAS_TRANSPOSE right_op() const { return right_op_; } + blas::TransposeFlag left_op() const { return left_op_; } + blas::TransposeFlag right_op() const { return right_op_; } }; // class GemmHelper -} // namespace math -} // namespace TiledArray +} // namespace TiledArray::math #endif // TILEDARRAY_MATH_GEMM_HELPER_H__INCLUDED diff --git a/src/TiledArray/math/parallel_gemm.h b/src/TiledArray/math/parallel_gemm.h index 8d286e0ad8..fd65599356 100644 --- a/src/TiledArray/math/parallel_gemm.h +++ b/src/TiledArray/math/parallel_gemm.h @@ -129,7 +129,7 @@ class MatrixBlockTask : public tbb::task { template class GemmTask : public tbb::task { - const madness::cblas::CBLAS_TRANSPOSE op_a_, op_b_; + const blas::TransposeFlag op_a_, op_b_; const integer m_, n_, k_; const Alpha alpha_; std::shared_ptr a_; @@ -140,8 +140,8 @@ class GemmTask : public tbb::task { const integer ldc_; public: - GemmTask(madness::cblas::CBLAS_TRANSPOSE op_a, - madness::cblas::CBLAS_TRANSPOSE op_b, const integer m, + GemmTask(blas::TransposeFlag op_a, + blas::TransposeFlag op_b, const integer m, const integer n, const integer k, const Alpha alpha, const std::shared_ptr& a, const std::shared_ptr& b, const Beta beta, const std::shared_ptr& c, const integer ldc) diff --git a/src/TiledArray/proc_grid.h b/src/TiledArray/proc_grid.h index f1b150b6b7..a401e0ac1e 100644 --- a/src/TiledArray/proc_grid.h +++ b/src/TiledArray/proc_grid.h @@ -26,7 +26,6 @@ #ifndef TILEDARRAY_GRID_H__INCLUDED #define TILEDARRAY_GRID_H__INCLUDED -#include #include namespace TiledArray { diff --git a/src/TiledArray/tensor/kernels.h b/src/TiledArray/tensor/kernels.h index 7521e357d3..bfbe25f083 100644 --- a/src/TiledArray/tensor/kernels.h +++ b/src/TiledArray/tensor/kernels.h @@ -26,7 +26,6 @@ #ifndef TILEDARRAY_TENSOR_KENERLS_H__INCLUDED #define TILEDARRAY_TENSOR_KENERLS_H__INCLUDED -#include #include #include diff --git a/src/TiledArray/tensor/tensor.h b/src/TiledArray/tensor/tensor.h index cce3e42abd..0e174a5fee 100644 --- a/src/TiledArray/tensor/tensor.h +++ b/src/TiledArray/tensor/tensor.h @@ -1635,13 +1635,15 @@ class Tensor { // Get the leading dimension for left and right matrices. const integer lda = - (gemm_helper.left_op() == madness::cblas::NoTrans ? k : m); + (gemm_helper.left_op() == blas::NoTranspose ? k : m); const integer ldb = - (gemm_helper.right_op() == madness::cblas::NoTrans ? n : k); + (gemm_helper.right_op() == blas::NoTranspose ? n : k); - math::gemm(gemm_helper.left_op(), gemm_helper.right_op(), m, n, k, factor, - pimpl_->data_, lda, other.data(), ldb, numeric_type(0), - result.data(), n); + blas::gemm( + gemm_helper.left_op(), gemm_helper.right_op(), m, n, k, factor, + pimpl_->data_, lda, other.data(), ldb, numeric_type(0), + result.data(), n + ); #ifdef TA_ENABLE_TILE_OPS_LOGGING if (TiledArray::TileOpsLogger::get_instance_ptr() != nullptr && @@ -1790,9 +1792,9 @@ class Tensor { // Get the leading dimension for left and right matrices. const integer lda = - (gemm_helper.left_op() == madness::cblas::NoTrans ? k : m); + (gemm_helper.left_op() == blas::NoTranspose ? k : m); const integer ldb = - (gemm_helper.right_op() == madness::cblas::NoTrans ? n : k); + (gemm_helper.right_op() == blas::NoTranspose ? n : k); // may need to split gemm into multiply + accumulate for tracing purposes #ifdef TA_ENABLE_TILE_OPS_LOGGING @@ -1809,10 +1811,12 @@ class Tensor { std::copy(pimpl_->data_, pimpl_->data_ + tile_volume, data_copy.get()); } - math::gemm(gemm_helper.left_op(), gemm_helper.right_op(), m, n, k, - factor, left.data(), lda, right.data(), ldb, - twostep ? numeric_type(0) : numeric_type(1), pimpl_->data_, - n); + non_distributed::gemm( + gemm_helper.left_op(), gemm_helper.right_op(), m, n, k, + factor, left.data(), lda, right.data(), ldb, + twostep ? numeric_type(0) : numeric_type(1), pimpl_->data_, + n + ); if (TiledArray::TileOpsLogger::get_instance_ptr() != nullptr && TiledArray::TileOpsLogger::get_instance().gemm) { @@ -1862,9 +1866,11 @@ class Tensor { } } #else // TA_ENABLE_TILE_OPS_LOGGING - math::gemm(gemm_helper.left_op(), gemm_helper.right_op(), m, n, k, factor, - left.data(), lda, right.data(), ldb, numeric_type(1), - pimpl_->data_, n); + blas::gemm( + gemm_helper.left_op(), gemm_helper.right_op(), m, n, k, factor, + left.data(), lda, right.data(), ldb, numeric_type(1), + pimpl_->data_, n + ); #endif // TA_ENABLE_TILE_OPS_LOGGING } @@ -1928,18 +1934,18 @@ class Tensor { // Get the leading dimension for left and right matrices. const integer lda = - (gemm_helper.left_op() == madness::cblas::NoTrans ? K : M); + (gemm_helper.left_op() == blas::NoTranspose ? K : M); const integer ldb = - (gemm_helper.right_op() == madness::cblas::NoTrans ? N : K); + (gemm_helper.right_op() == blas::NoTranspose ? N : K); for (integer m = 0; m != M; ++m) { for (integer n = 0; n != N; ++n) { auto c_offset = m * N + n; for (integer k = 0; k != K; ++k) { - auto a_offset = gemm_helper.left_op() == madness::cblas::NoTrans + auto a_offset = gemm_helper.left_op() == blas::NoTranspose ? m * lda + k : k * lda + m; - auto b_offset = gemm_helper.right_op() == madness::cblas::NoTrans + auto b_offset = gemm_helper.right_op() == blas::NoTranspose ? k * ldb + n : n * ldb + k; elem_muladd_op(*(pimpl_->data_ + c_offset), *(left.data() + a_offset), diff --git a/tests/linear_algebra.cpp b/tests/linear_algebra.cpp index 07f9a58be2..87e5ae2264 100644 --- a/tests/linear_algebra.cpp +++ b/tests/linear_algebra.cpp @@ -4,10 +4,10 @@ //#include "range_fixture.h" #include "unit_test_config.h" -#include "TiledArray/algebra/lapack/cholesky.h" -#include "TiledArray/algebra/lapack/heig.h" -#include "TiledArray/algebra/lapack/lu.h" -#include "TiledArray/algebra/lapack/svd.h" +#include "TiledArray/algebra/non-distributed/cholesky.h" +#include "TiledArray/algebra/non-distributed/heig.h" +#include "TiledArray/algebra/non-distributed/lu.h" +#include "TiledArray/algebra/non-distributed/svd.h" #include "TiledArray/algebra/cholesky.h" #include "TiledArray/algebra/heig.h" @@ -15,7 +15,7 @@ #include "TiledArray/algebra/svd.h" namespace TA = TiledArray; -namespace lapack = TA::lapack; +namespace non_dist = TA::non_distributed; #if TILEDARRAY_HAS_SCALAPACK @@ -23,9 +23,9 @@ namespace lapack = TA::lapack; namespace scalapack = TA::scalapack; #define TILEDARRAY_SCALAPACK_TEST(F, E) \ GlobalFixture::world->gop.fence(); \ - compare("TiledArray::scalapack", lapack::F, scalapack::F, E); \ + compare("TiledArray::scalapack", non_dist::F, scalapack::F, E); \ GlobalFixture::world->gop.fence(); \ - compare("TiledArray", lapack::F, TiledArray::F, E); + compare("TiledArray", non_dist::F, TiledArray::F, E); #else #define TILEDARRAY_SCALAPACK_TEST(...) #endif @@ -107,11 +107,11 @@ struct LinearAlgebraFixture : ReferenceFixture { } } template - static void compare(const char* context, const A& lapack, const A& result, + static void compare(const char* context, const A& non_dist, const A& result, double e) { BOOST_TEST_CONTEXT(context); - auto diff_with_lapack = (lapack("i,j") - result("i,j")).norm().get(); - BOOST_CHECK_SMALL(diff_with_lapack, e); + auto diff_with_non_dist = (non_dist("i,j") - result("i,j")).norm().get(); + BOOST_CHECK_SMALL(diff_with_non_dist, e); } #endif }; @@ -411,15 +411,15 @@ BOOST_AUTO_TEST_CASE(heig_same_tiling) { return this->make_ta_reference(t, range); }); - auto [evals, evecs] = lapack::heig(ref_ta); - auto [evals_lapack, evecs_lapack] = lapack::heig(ref_ta); + auto [evals, evecs] = non_dist::heig(ref_ta); + auto [evals_non_dist, evecs_non_dist] = non_dist::heig(ref_ta); // auto evals = heig( ref_ta ); BOOST_CHECK(evecs.trange() == ref_ta.trange()); - // check eigenvectors against lapack only, for now ... + // check eigenvectors against non_dist only, for now ... decltype(evecs) evecs_error; - evecs_error("i,j") = evecs_lapack("i,j") - evecs("i,j"); + evecs_error("i,j") = evecs_non_dist("i,j") - evecs("i,j"); // TODO need to fix phases of the eigenvectors to be able to compare ... // BOOST_CHECK_SMALL(evecs_error("i,j").norm().get(), // N * N * std::numeric_limits::epsilon()); @@ -428,7 +428,7 @@ BOOST_AUTO_TEST_CASE(heig_same_tiling) { double tol = N * N * std::numeric_limits::epsilon(); for (int64_t i = 0; i < N; ++i) { BOOST_CHECK_SMALL(std::abs(evals[i] - exact_evals[i]), tol); - BOOST_CHECK_SMALL(std::abs(evals_lapack[i] - exact_evals[i]), tol); + BOOST_CHECK_SMALL(std::abs(evals_non_dist[i] - exact_evals[i]), tol); } GlobalFixture::world->gop.fence(); @@ -445,14 +445,14 @@ BOOST_AUTO_TEST_CASE(heig_diff_tiling) { }); auto new_trange = gen_trange(N, {64ul}); - auto [evals, evecs] = lapack::heig(ref_ta, new_trange); - auto [evals_lapack, evecs_lapack] = lapack::heig(ref_ta, new_trange); + auto [evals, evecs] = non_dist::heig(ref_ta, new_trange); + auto [evals_non_dist, evecs_non_dist] = non_dist::heig(ref_ta, new_trange); BOOST_CHECK(evecs.trange() == new_trange); - // check eigenvectors against lapack only, for now ... + // check eigenvectors against non_dist only, for now ... decltype(evecs) evecs_error; - evecs_error("i,j") = evecs_lapack("i,j") - evecs("i,j"); + evecs_error("i,j") = evecs_non_dist("i,j") - evecs("i,j"); // TODO need to fix phases of the eigenvectors to be able to compare ... // BOOST_CHECK_SMALL(evecs_error("i,j").norm().get(), // N * N * std::numeric_limits::epsilon()); @@ -461,7 +461,7 @@ BOOST_AUTO_TEST_CASE(heig_diff_tiling) { double tol = N * N * std::numeric_limits::epsilon(); for (int64_t i = 0; i < N; ++i) { BOOST_CHECK_SMALL(std::abs(evals[i] - exact_evals[i]), tol); - BOOST_CHECK_SMALL(std::abs(evals_lapack[i] - exact_evals[i]), tol); + BOOST_CHECK_SMALL(std::abs(evals_non_dist[i] - exact_evals[i]), tol); } GlobalFixture::world->gop.fence(); @@ -492,7 +492,7 @@ BOOST_AUTO_TEST_CASE(heig_generalized) { }); GlobalFixture::world->gop.fence(); - auto [evals, evecs] = lapack::heig(ref_ta, dense_iden); + auto [evals, evecs] = non_dist::heig(ref_ta, dense_iden); // auto evals = heig( ref_ta ); BOOST_CHECK(evecs.trange() == ref_ta.trange()); @@ -518,7 +518,7 @@ BOOST_AUTO_TEST_CASE(cholesky) { return this->make_ta_reference(t, range); }); - auto L = lapack::cholesky(A); + auto L = non_dist::cholesky(A); BOOST_CHECK(L.trange() == A.trange()); @@ -528,8 +528,8 @@ BOOST_AUTO_TEST_CASE(cholesky) { BOOST_CHECK_SMALL(A_minus_LLt("i,j").norm().get(), N * N * std::numeric_limits::epsilon()); - // check against LAPACK also - auto L_ref = TiledArray::lapack::cholesky(A); + // check against NON_DIST also + auto L_ref = non_dist::cholesky(A); decltype(L) L_diff; L_diff("i,j") = L("i,j") - L_ref("i,j"); @@ -551,7 +551,7 @@ BOOST_AUTO_TEST_CASE(cholesky_linv) { }); decltype(A) Acopy = A.clone(); - auto Linv = lapack::cholesky_linv(A); + auto Linv = non_dist::cholesky_linv(A); BOOST_CHECK(Linv.trange() == A.trange()); @@ -591,7 +591,7 @@ BOOST_AUTO_TEST_CASE(cholesky_linv_retl) { return this->make_ta_reference(t, range); }); - auto [L, Linv] = lapack::cholesky_linv(A); + auto [L, Linv] = non_dist::cholesky_linv(A); BOOST_CHECK(Linv.trange() == A.trange()); BOOST_CHECK(L.trange() == A.trange()); @@ -631,12 +631,12 @@ BOOST_AUTO_TEST_CASE(cholesky_solve) { return this->make_ta_reference(t, range); }); - auto iden = lapack::cholesky_solve(A, A); + auto iden = non_dist::cholesky_solve(A, A); BOOST_CHECK(iden.trange() == A.trange()); - auto iden_lapack = lapack::cholesky_solve(A, A); + auto iden_non_dist = non_dist::cholesky_solve(A, A); decltype(iden) iden_error; - iden_error("i,j") = iden("i,j") - iden_lapack("i,j"); + iden_error("i,j") = iden("i,j") - iden_non_dist("i,j"); BOOST_CHECK_SMALL(iden_error("i,j").norm().get(), N * N * std::numeric_limits::epsilon()); @@ -669,19 +669,19 @@ BOOST_AUTO_TEST_CASE(cholesky_lsolve) { }); // Should produce X = L**H - auto [L, X] = lapack::cholesky_lsolve(TA::TransposeFlag::NoTranspose, A, A); + auto [L, X] = non_dist::cholesky_lsolve(TA::TransposeFlag::NoTranspose, A, A); BOOST_CHECK(X.trange() == A.trange()); BOOST_CHECK(L.trange() == A.trange()); - // first, test against LAPACK - auto [L_lapack, X_lapack] = - lapack::cholesky_lsolve(TA::TransposeFlag::NoTranspose, A, A); + // first, test against NON_DIST + auto [L_non_dist, X_non_dist] = + non_dist::cholesky_lsolve(TA::TransposeFlag::NoTranspose, A, A); decltype(L) L_error; - L_error("i,j") = L("i,j") - L_lapack("i,j"); + L_error("i,j") = L("i,j") - L_non_dist("i,j"); BOOST_CHECK_SMALL(L_error("i,j").norm().get(), N * N * std::numeric_limits::epsilon()); decltype(X) X_error; - X_error("i,j") = X("i,j") - X_lapack("i,j"); + X_error("i,j") = X("i,j") - X_non_dist("i,j"); BOOST_CHECK_SMALL(X_error("i,j").norm().get(), N * N * std::numeric_limits::epsilon()); @@ -704,7 +704,7 @@ BOOST_AUTO_TEST_CASE(lu_solve) { return this->make_ta_reference(t, range); }); - auto iden = lapack::lu_solve(ref_ta, ref_ta); + auto iden = non_dist::lu_solve(ref_ta, ref_ta); BOOST_CHECK(iden.trange() == ref_ta.trange()); @@ -741,7 +741,7 @@ BOOST_AUTO_TEST_CASE(lu_inv) { TA::TArray iden(*GlobalFixture::world, trange); - auto Ainv = lapack::lu_inv(ref_ta); + auto Ainv = non_dist::lu_inv(ref_ta); iden("i,j") = Ainv("i,k") * ref_ta("k,j"); BOOST_CHECK(iden.trange() == ref_ta.trange()); @@ -778,7 +778,7 @@ BOOST_AUTO_TEST_CASE(svd_values_only) { return this->make_ta_reference(t, range); }); - auto S = lapack::svd(ref_ta, trange, trange); + auto S = non_dist::svd(ref_ta, trange, trange); std::vector exact_singular_values = exact_evals; std::sort(exact_singular_values.begin(), exact_singular_values.end(), @@ -802,7 +802,7 @@ BOOST_AUTO_TEST_CASE(svd_leftvectors) { return this->make_ta_reference(t, range); }); - auto [S, U] = lapack::svd(ref_ta, trange, trange); + auto [S, U] = non_dist::svd(ref_ta, trange, trange); std::vector exact_singular_values = exact_evals; std::sort(exact_singular_values.begin(), exact_singular_values.end(), @@ -826,7 +826,7 @@ BOOST_AUTO_TEST_CASE(svd_rightvectors) { return this->make_ta_reference(t, range); }); - auto [S, VT] = lapack::svd(ref_ta, trange, trange); + auto [S, VT] = non_dist::svd(ref_ta, trange, trange); std::vector exact_singular_values = exact_evals; std::sort(exact_singular_values.begin(), exact_singular_values.end(), @@ -850,7 +850,7 @@ BOOST_AUTO_TEST_CASE(svd_allvectors) { return this->make_ta_reference(t, range); }); - auto [S, U, VT] = lapack::svd(ref_ta, trange, trange); + auto [S, U, VT] = non_dist::svd(ref_ta, trange, trange); std::vector exact_singular_values = exact_evals; std::sort(exact_singular_values.begin(), exact_singular_values.end(), diff --git a/tests/math_blas.cpp b/tests/math_blas.cpp index 824cdff4d3..44fcd896f7 100644 --- a/tests/math_blas.cpp +++ b/tests/math_blas.cpp @@ -69,7 +69,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(integral_gemm, T, int_types) { // Test the gemm operation BOOST_REQUIRE_NO_THROW( - TiledArray::math::gemm(madness::cblas::NoTrans, madness::cblas::NoTrans, + TiledArray::blas::gemm(madness::cblas::NoTrans, madness::cblas::NoTrans, m, n, k, 3, a, lda, b, ldb, 0, c, ldc)); for (integer i = 0; i < m; ++i) { @@ -119,7 +119,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(integral_gemm_ld, T, int_types) { // Test the gemm operation BOOST_REQUIRE_NO_THROW( - TiledArray::math::gemm(madness::cblas::NoTrans, madness::cblas::NoTrans, + TiledArray::blas::gemm(madness::cblas::NoTrans, madness::cblas::NoTrans, m, n, k, 3, a, lda, b, ldb, 0, c, ldc)); for (integer i = 0; i < m; ++i) { @@ -168,7 +168,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(floating_point_gemm, T, floating_point_types) { // Test the gemm operation BOOST_REQUIRE_NO_THROW( - TiledArray::math::gemm(madness::cblas::NoTrans, madness::cblas::NoTrans, + TiledArray::blas::gemm(madness::cblas::NoTrans, madness::cblas::NoTrans, m, n, k, 3, a, lda, b, ldb, 0, c, ldc)); for (integer i = 0; i < m; ++i) { for (integer j = 0; j < n; ++j) { @@ -216,7 +216,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(floating_point_gemm_ld, T, floating_point_types) { // Test the gemm operation BOOST_REQUIRE_NO_THROW( - TiledArray::math::gemm(madness::cblas::NoTrans, madness::cblas::NoTrans, + TiledArray::blas::gemm(madness::cblas::NoTrans, madness::cblas::NoTrans, m, n, k, 3, a, lda, b, ldb, 0, c, ldc)); for (integer i = 0; i < m; ++i) { @@ -263,7 +263,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(complex_gemm, T, floating_point_types) { // Test the gemm operation BOOST_REQUIRE_NO_THROW( - TiledArray::math::gemm(madness::cblas::NoTrans, madness::cblas::NoTrans, + TiledArray::blas::gemm(madness::cblas::NoTrans, madness::cblas::NoTrans, m, n, k, 3, a, lda, b, ldb, 0, c, ldc)); for (integer i = 0; i < m; ++i) { @@ -315,7 +315,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(complex_gemm_ld, T, floating_point_types) { // Test the gemm operation BOOST_REQUIRE_NO_THROW( - TiledArray::math::gemm(madness::cblas::NoTrans, madness::cblas::NoTrans, + TiledArray::blas::gemm(madness::cblas::NoTrans, madness::cblas::NoTrans, m, n, k, 3, a, lda, b, ldb, 0, c, ldc)); for (integer i = 0; i < m; ++i) { From 9f2086c0a07c65c14928efcb4ce03bb9ed9fdbfc Mon Sep 17 00:00:00 2001 From: asadchev Date: Mon, 7 Dec 2020 12:43:57 -0500 Subject: [PATCH 3/7] Move TiledArray/algebra to TiledArray/math/linalg --- examples/scalapack/conversion.cpp | 9 +- examples/scalapack/evp.cpp | 5 +- src/CMakeLists.txt | 33 +++-- src/TiledArray/algebra/types.h | 66 --------- src/TiledArray/algebra/utils.h | 139 ------------------ src/TiledArray/conversions/retile.h | 4 +- src/TiledArray/dist_array.h | 40 +++++ src/TiledArray/math/linalg.h | 36 +++++ src/TiledArray/math/linalg/basic.h | 62 ++++++++ .../{algebra => math/linalg}/cholesky.h | 27 ++-- .../{algebra => math/linalg}/conjgrad.h | 53 ++++--- .../{algebra => math/linalg}/diis.h | 51 ++++--- src/TiledArray/math/linalg/forward.h | 51 +++++++ .../{algebra => math/linalg}/heig.h | 18 ++- src/TiledArray/{algebra => math/linalg}/lu.h | 10 +- .../linalg}/non-distributed/cholesky.h | 60 ++++---- .../linalg}/non-distributed/heig.h | 24 +-- .../linalg}/non-distributed/lu.h | 24 +-- .../linalg}/non-distributed/svd.h | 29 ++-- .../{algebra => math/linalg}/rank-local.cpp | 6 +- .../{algebra => math/linalg}/rank-local.h | 13 +- .../{algebra => math/linalg}/scalapack/all.h | 14 +- .../linalg/scalapack}/block_cyclic.h | 10 +- .../linalg}/scalapack/cholesky.h | 25 ++-- .../{algebra => math/linalg}/scalapack/heig.h | 15 +- .../{algebra => math/linalg}/scalapack/lu.h | 15 +- .../{algebra => math/linalg}/scalapack/svd.h | 26 ++-- .../{algebra => math/linalg}/scalapack/util.h | 17 +-- src/TiledArray/{algebra => math/linalg}/svd.h | 25 ++-- .../non-distributed => math/linalg}/util.h | 35 +++-- src/TiledArray/math/scalapack.h | 4 +- src/tiledarray.h | 9 +- tests/linear_algebra.cpp | 39 +++-- tests/solvers.cpp | 2 +- 34 files changed, 490 insertions(+), 506 deletions(-) delete mode 100644 src/TiledArray/algebra/types.h delete mode 100644 src/TiledArray/algebra/utils.h create mode 100644 src/TiledArray/math/linalg.h create mode 100644 src/TiledArray/math/linalg/basic.h rename src/TiledArray/{algebra => math/linalg}/cholesky.h (76%) rename src/TiledArray/{algebra => math/linalg}/conjgrad.h (85%) rename src/TiledArray/{algebra => math/linalg}/diis.h (93%) create mode 100644 src/TiledArray/math/linalg/forward.h rename src/TiledArray/{algebra => math/linalg}/heig.h (80%) rename src/TiledArray/{algebra => math/linalg}/lu.h (80%) rename src/TiledArray/{algebra => math/linalg}/non-distributed/cholesky.h (81%) rename src/TiledArray/{algebra => math/linalg}/non-distributed/heig.h (84%) rename src/TiledArray/{algebra => math/linalg}/non-distributed/lu.h (76%) rename src/TiledArray/{algebra => math/linalg}/non-distributed/svd.h (76%) rename src/TiledArray/{algebra => math/linalg}/rank-local.cpp (97%) rename src/TiledArray/{algebra => math/linalg}/rank-local.h (73%) rename src/TiledArray/{algebra => math/linalg}/scalapack/all.h (71%) rename src/TiledArray/{conversions => math/linalg/scalapack}/block_cyclic.h (97%) rename src/TiledArray/{algebra => math/linalg}/scalapack/cholesky.h (94%) rename src/TiledArray/{algebra => math/linalg}/scalapack/heig.h (94%) rename src/TiledArray/{algebra => math/linalg}/scalapack/lu.h (92%) rename src/TiledArray/{algebra => math/linalg}/scalapack/svd.h (87%) rename src/TiledArray/{algebra => math/linalg}/scalapack/util.h (87%) rename src/TiledArray/{algebra => math/linalg}/svd.h (67%) rename src/TiledArray/{algebra/non-distributed => math/linalg}/util.h (79%) diff --git a/examples/scalapack/conversion.cpp b/examples/scalapack/conversion.cpp index edb60fb77a..bca97c081f 100644 --- a/examples/scalapack/conversion.cpp +++ b/examples/scalapack/conversion.cpp @@ -25,8 +25,11 @@ */ #include +#include #include +namespace scalapack = TiledArray::math::linalg::scalapack; + template int64_t div_ceil(Integral1 x, Integral2 y) { int64_t x_ll = x; @@ -66,7 +69,7 @@ int main(int argc, char** argv) { // Create Test Matrix blacspp::Grid grid = blacspp::Grid::square_grid(MPI_COMM_WORLD); - TA::scalapack::BlockCyclicMatrix ref_matrix(world, grid, N, N, NB, NB); + scalapack::BlockCyclicMatrix ref_matrix(world, grid, N, N, NB, NB); for (size_t i = 0; i < N; ++i) for (size_t j = 0; j < N; ++j) @@ -125,7 +128,7 @@ int main(int argc, char** argv) { TA::make_array >(world, trange, make_ta_reference); world.gop.fence(); - TA::scalapack::BlockCyclicMatrix test_matrix(ref_ta, grid, NB, NB); + scalapack::BlockCyclicMatrix test_matrix(ref_ta, grid, NB, NB); world.gop.fence(); double local_norm_diff = @@ -172,7 +175,7 @@ int main(int argc, char** argv) { TA::make_array >(world, trange, make_ta_reference); world.gop.fence(); - TA::scalapack::BlockCyclicMatrix test_matrix(ref_ta, grid, NB, NB); + scalapack::BlockCyclicMatrix test_matrix(ref_ta, grid, NB, NB); world.gop.fence(); double local_norm_diff = diff --git a/examples/scalapack/evp.cpp b/examples/scalapack/evp.cpp index 22964ccae0..9a5bbd177e 100644 --- a/examples/scalapack/evp.cpp +++ b/examples/scalapack/evp.cpp @@ -28,10 +28,11 @@ #include #include -#include +#include using Array = TA::TArray; // using Array = TA::TSpArray; +namespace scalapack = TiledArray::math::linalg::scalapack; TA::TiledRange gen_trange(size_t N, const std::vector& TA_NBs) { assert(TA_NBs.size() > 0); @@ -93,7 +94,7 @@ int main(int argc, char** argv) { tensor_symm("i,j") = 0.5 * (tensor("i,j") + tensor("j,i")); tensor("i,j") = tensor_symm("i,j"); - auto [evals, evecs_ta] = TA::scalapack::heig(tensor); + auto [evals, evecs_ta] = scalapack::heig(tensor); //// Check EVP with TA Array tmp = diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 21bc79749b..9e8d7b01da 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -57,18 +57,20 @@ TiledArray/utility.h TiledArray/val_array.h TiledArray/version.h TiledArray/zero_tensor.h -TiledArray/algebra/conjgrad.h -TiledArray/algebra/diis.h -TiledArray/algebra/utils.h -TiledArray/algebra/cholesky.h -TiledArray/algebra/heig.h -TiledArray/algebra/lu.h -TiledArray/algebra/svd.h -TiledArray/algebra/types.h -TiledArray/algebra/scalapack/cholesky.h -TiledArray/algebra/scalapack/heig.h -TiledArray/algebra/scalapack/lu.h -TiledArray/algebra/scalapack/svd.h +TiledArray/math/linalg/forward.h +TiledArray/math/linalg/conjgrad.h +TiledArray/math/linalg/diis.h +TiledArray/math/linalg/util.h +TiledArray/math/linalg/cholesky.h +TiledArray/math/linalg/heig.h +TiledArray/math/linalg/lu.h +TiledArray/math/linalg/svd.h +TiledArray/math/linalg/scalapack/util.h +TiledArray/math/linalg/scalapack/block_cyclic.h +TiledArray/math/linalg/scalapack/cholesky.h +TiledArray/math/linalg/scalapack/heig.h +TiledArray/math/linalg/scalapack/lu.h +TiledArray/math/linalg/scalapack/svd.h TiledArray/conversions/btas.h TiledArray/conversions/clone.h TiledArray/conversions/dense_to_sparse.h @@ -77,7 +79,6 @@ TiledArray/conversions/foreach.h TiledArray/conversions/vector_of_arrays.h TiledArray/conversions/make_array.h TiledArray/conversions/sparse_to_dense.h -TiledArray/conversions/block_cyclic.h TiledArray/conversions/to_new_tile_type.h TiledArray/conversions/truncate.h TiledArray/conversions/retile.h @@ -175,8 +176,8 @@ TiledArray/util/singleton.h TiledArray/util/time.h TiledArray/util/vector.h -TiledArray/algebra/rank-local.h -TiledArray/algebra/rank-local.cpp +TiledArray/math/linalg/rank-local.h +TiledArray/math/linalg/rank-local.cpp ) @@ -209,7 +210,7 @@ TiledArray/array_impl.cpp TiledArray/dist_array.cpp TiledArray/util/backtrace.cpp TiledArray/util/bug.cpp -TiledArray/algebra/rank-local.cpp +TiledArray/math/linalg/rank-local.cpp ) # the list of libraries on which TiledArray depends on, will be cached later diff --git a/src/TiledArray/algebra/types.h b/src/TiledArray/algebra/types.h deleted file mode 100644 index 1cc0a50e57..0000000000 --- a/src/TiledArray/algebra/types.h +++ /dev/null @@ -1,66 +0,0 @@ -/* - * This file is a part of TiledArray. - * Copyright (C) 2020 Virginia Tech - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - * - * David Williams-Young - * Computational Research Division, Lawrence Berkeley National Laboratory - * - * svd.h - * Created: 12 June, 2020 - * - */ -#ifndef TILEDARRAY_ALGEBRA_SVD_UTILS_H__INCLUDED -#define TILEDARRAY_ALGEBRA_SVD_UTILS_H__INCLUDED - -#include -#include - -namespace TiledArray { - -enum TransposeFlag { NoTranspose, Transpose, ConjTranspose }; - -struct SVDReturnType {}; -struct SVDValuesOnly : public SVDReturnType {}; -struct SVDLeftVectors : public SVDReturnType {}; -struct SVDRightVectors : public SVDReturnType {}; -struct SVDAllVectors : public SVDReturnType {}; - -namespace detail { - -template -struct is_svd_return_type : public std::false_type {}; - -template -struct is_svd_return_type< - SVDType, std::enable_if_t>> - : public std::true_type {}; - -template -inline constexpr bool is_svd_return_type_v = is_svd_return_type::value; - -template -struct enable_if_svd_return_type - : public std::enable_if, U> {}; - -template -using enable_if_svd_return_type_t = - typename enable_if_svd_return_type::type; - -} // namespace detail - -} // namespace TiledArray - -#endif // TILEDARRAY_ALGEBRA_SVD_UTILS_H__INCLUDED diff --git a/src/TiledArray/algebra/utils.h b/src/TiledArray/algebra/utils.h deleted file mode 100644 index 008a3d218b..0000000000 --- a/src/TiledArray/algebra/utils.h +++ /dev/null @@ -1,139 +0,0 @@ -/* - * This file is a part of TiledArray. - * Copyright (C) 2013 Virginia Tech - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - * - * Eduard Valeyev - * Department of Chemistry, Virginia Tech - * - * utils.h - * May 20, 2013 - * - */ - -#ifndef TILEDARRAY_ALGEBRA_UTILS_H__INCLUDED -#define TILEDARRAY_ALGEBRA_UTILS_H__INCLUDED - -#include - -#include "../dist_array.h" -#include "../expressions/expr.h" -#include "TiledArray/util/annotation.h" - -namespace TiledArray { - -template -inline size_t size(const DistArray& a) { - // this is the number of tiles - if (a.size() > 0) // assuming dense shape - return a.trange().elements_range().volume(); - else - return 0; -} - -template -inline DistArray copy(const DistArray& a) { - return a; -} - -template -inline void zero(DistArray& a) { - const std::string vars = - detail::dummy_annotation(a.trange().tiles_range().rank()); - a(vars) = typename DistArray::element_type(0) * a(vars); -} - -template -typename DistArray::element_type inline minabs_value( - const DistArray& a) { - return a(detail::dummy_annotation(a.trange().tiles_range().rank())).abs_min(); -} - -template -inline typename DistArray::element_type maxabs_value( - const DistArray& a) { - return a(detail::dummy_annotation(a.trange().tiles_range().rank())).abs_max(); -} - -template -inline void vec_multiply(DistArray& a1, - const DistArray& a2) { - const std::string vars = - detail::dummy_annotation(a1.trange().tiles_range().rank()); - a1(vars) = a1(vars) * a2(vars); -} - -template -inline typename DistArray::element_type dot_product( - const DistArray& a1, const DistArray& a2) { - const std::string vars = - detail::dummy_annotation(a1.trange().tiles_range().rank()); - return a1(vars).dot(a2(vars)).get(); -} - -template -inline typename TiledArray::expressions::ExprTrait::scalar_type dot( - const TiledArray::expressions::Expr& a1, - const TiledArray::expressions::Expr& a2) { - static_assert( - TiledArray::expressions::is_aliased::value, - "no_alias() expressions are not allowed on the right-hand side of the " - "assignment operator."); - static_assert( - TiledArray::expressions::is_aliased::value, - "no_alias() expressions are not allowed on the right-hand side of the " - "assignment operator."); - return a1.dot(a2).get(); -} - -template -inline void scale( - DistArray& a, - typename DistArray::element_type scaling_factor) { - const std::string vars = - detail::dummy_annotation(a.trange().tiles_range().rank()); - a(vars) = scaling_factor * a(vars); -} - -template -inline void axpy(DistArray& y, - typename DistArray::element_type a, - const DistArray& x) { - const std::string vars = - detail::dummy_annotation(y.trange().tiles_range().rank()); - y(vars) = y(vars) + a * x(vars); -} - -template -inline void assign(DistArray& m1, - const DistArray& m2) { - m1 = m2; -} - -template -inline typename DistArray::scalar_type norm2( - const DistArray& a) { - return std::sqrt(a(detail::dummy_annotation(a.trange().tiles_range().rank())) - .squared_norm()); -} - -template -inline void print(const DistArray& a, const char* label) { - std::cout << label << ":\n" << a << "\n"; -} - -} // namespace TiledArray - -#endif // TILEDARRAY_ALGEBRA_UTILS_H__INCLUDED diff --git a/src/TiledArray/conversions/retile.h b/src/TiledArray/conversions/retile.h index cf5c858f05..26440166c4 100644 --- a/src/TiledArray/conversions/retile.h +++ b/src/TiledArray/conversions/retile.h @@ -22,8 +22,8 @@ #ifndef TILEDARRAY_RETILE_H #define TILEDARRAY_RETILE_H -#include "../algebra/utils.h" -#include "../special/diagonal_array.h" +#include "TiledArray/util/annotation.h" +#include "TiledArray/special/diagonal_array.h" /// \name Retile function /// \brief Retiles a tensor with a provided TiledRange diff --git a/src/TiledArray/dist_array.h b/src/TiledArray/dist_array.h index 49ca61972d..b2c7f74d6a 100644 --- a/src/TiledArray/dist_array.h +++ b/src/TiledArray/dist_array.h @@ -1605,6 +1605,46 @@ inline std::ostream& operator<<(std::ostream& os, return os; } + +template +auto rank(const DistArray &a) { + return a.trange().tiles_range().rank(); +} + +template +size_t volume(const DistArray& a) { + // this is the number of tiles + if (a.size() > 0) // assuming dense shape + return a.trange().elements_range().volume(); + return 0; +} + +template +auto abs_min(const DistArray& a) { + return a(detail::dummy_annotation(rank(a))).abs_min(); +} + +template +auto abs_max(const DistArray& a) { + return a(detail::dummy_annotation(rank(a))).abs_max(); +} + +template +auto dot(const DistArray& a, const DistArray& b) { + return (a(detail::dummy_annotation(rank(a))) + .dot(b(detail::dummy_annotation(rank(b))))).get(); +} + +template +auto squared_norm(const DistArray& a) { + return a(detail::dummy_annotation(rank(a))).squared_norm(); +} + +template +auto norm2(const DistArray& a) { + return std::sqrt(squared_norm(a)); +} + } // namespace TiledArray // serialization diff --git a/src/TiledArray/math/linalg.h b/src/TiledArray/math/linalg.h new file mode 100644 index 0000000000..c6dc015e7a --- /dev/null +++ b/src/TiledArray/math/linalg.h @@ -0,0 +1,36 @@ +/* + * This file is a part of TiledArray. + * Copyright (C) 2020 Virginia Tech + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + * + * David Williams-Young + * Computational Research Division, Lawrence Berkeley National Laboratory + * + * scalapack.h + * Created: 25 May, 2020 + * + */ + +#ifndef TILEDARRAY_MATH_LINALG_LINALG_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_LINALG_H__INCLUDED + +#include +#include +#include +#include +#include +#include + +#endif // TILEDARRAY_MATH_LINALG_LINALG_H__INCLUDED diff --git a/src/TiledArray/math/linalg/basic.h b/src/TiledArray/math/linalg/basic.h new file mode 100644 index 0000000000..32d7cbb94d --- /dev/null +++ b/src/TiledArray/math/linalg/basic.h @@ -0,0 +1,62 @@ +/* + * This file is a part of TiledArray. + * Copyright (C) 2013 Virginia Tech + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + * + * Eduard Valeyev + * Department of Chemistry, Virginia Tech + * + * util.h + * May 20, 2013 + * + */ + +#ifndef TILEDARRAY_MATH_LINALG_BASIC_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_BASIC_H__INCLUDED + +#include "TiledArray/dist_array.h" + +namespace TiledArray::math::linalg { + +template +inline void vec_multiply(DistArray& a1, + const DistArray& a2) { + auto vars = TiledArray::detail::dummy_annotation(rank(a1)); + a1(vars) = a1(vars) * a2(vars); +} + +template +inline void scale(DistArray& a, S scaling_factor) { + using numeric_type = typename DistArray::numeric_type; + auto vars = TiledArray::detail::dummy_annotation(rank(a)); + a(vars) = numeric_type(scaling_factor) * a(vars); +} + +template +inline void zero(DistArray& a) { + scale(a, 0); +} + +template +inline void axpy(DistArray& y, S alpha, + const DistArray& x) { + using numeric_type = typename DistArray::numeric_type; + auto vars = TiledArray::detail::dummy_annotation(rank(y)); + y(vars) = y(vars) + numeric_type(alpha) * x(vars); +} + +} // namespace TiledArray::math::linalg + +#endif // TILEDARRAY_MATH_LINALG_BASIC_H__INCLUDED diff --git a/src/TiledArray/algebra/cholesky.h b/src/TiledArray/math/linalg/cholesky.h similarity index 76% rename from src/TiledArray/algebra/cholesky.h rename to src/TiledArray/math/linalg/cholesky.h index 7f3a4e716d..3678c46d80 100644 --- a/src/TiledArray/algebra/cholesky.h +++ b/src/TiledArray/math/linalg/cholesky.h @@ -21,16 +21,16 @@ * Created: 16 October, 2020 * */ -#ifndef TILEDARRAY_ALGEBRA_CHOL_H__INCLUDED -#define TILEDARRAY_ALGEBRA_CHOL_H__INCLUDED +#ifndef TILEDARRAY_MATH_LINALG_CHOL_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_CHOL_H__INCLUDED #include #if TILEDARRAY_HAS_SCALAPACK -#include +#include #endif -#include +#include -namespace TiledArray { +namespace TiledArray::math::linalg { template auto cholesky(const Array& A, TiledRange l_trange = TiledRange()) { @@ -41,13 +41,13 @@ auto cholesky(const Array& A, TiledRange l_trange = TiledRange()) { return non_distributed::cholesky(A, l_trange); } -template +template 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(A, l_trange); + return scalapack::cholesky_linv(A, l_trange); #endif - return non_distributed::cholesky_linv(A, l_trange); + return non_distributed::cholesky_linv(A, l_trange); } template @@ -72,6 +72,13 @@ auto cholesky_lsolve(TransposeFlag transpose, const Array& A, const Array& B, return non_distributed::cholesky_lsolve(transpose, A, B, l_trange, x_trange); } -} // namespace TiledArray +} // namespace TiledArray::math::linalg + +namespace TiledArray { + using TiledArray::math::linalg::cholesky; + using TiledArray::math::linalg::cholesky_linv; + using TiledArray::math::linalg::cholesky_solve; + using TiledArray::math::linalg::cholesky_lsolve; +} -#endif // TILEDARRAY_ALGEBRA_CHOL_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_CHOL_H__INCLUDED diff --git a/src/TiledArray/algebra/conjgrad.h b/src/TiledArray/math/linalg/conjgrad.h similarity index 85% rename from src/TiledArray/algebra/conjgrad.h rename to src/TiledArray/math/linalg/conjgrad.h index 54bf2dd7f5..fdda051d9e 100644 --- a/src/TiledArray/algebra/conjgrad.h +++ b/src/TiledArray/math/linalg/conjgrad.h @@ -23,15 +23,14 @@ * */ -#ifndef TILEDARRAY_ALGEBRA_CONJGRAD_H__INCLUDED -#define TILEDARRAY_ALGEBRA_CONJGRAD_H__INCLUDED +#ifndef TILEDARRAY_MATH_LINALG_CONJGRAD_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_CONJGRAD_H__INCLUDED -#include -#include -#include -#include "../dist_array.h" +#include +#include +#include "TiledArray/dist_array.h" -namespace TiledArray { +namespace TiledArray::math::linalg { /// Solves real linear system a(x) = b , with \c a is a linear /// function of \c x , using conjugate gradient solver with a diagonal @@ -69,12 +68,8 @@ struct ConjugateGradientSolver { /// elements in the residual. value_type operator()(F& a, const D& b, D& x, const D& preconditioner, value_type convergence_target = -1.0) { -#ifdef TILEDARRAY_CXX_COMPILER_IS_ICC - std::size_t n = TiledArray::size(preconditioner); -#else - using TiledArray::size; - std::size_t n = size(preconditioner); -#endif + + std::size_t n = volume(preconditioner); const bool use_diis = false; DIIS diis; @@ -92,8 +87,8 @@ struct ConjugateGradientSolver { // approximate the conditio number as the ratio of the min and max elements // of the preconditioner assuming that preconditioner is the approximate // inverse of A in Ax - b =0 - const value_type precond_min = minabs_value(preconditioner); - const value_type precond_max = maxabs_value(preconditioner); + const value_type precond_min = abs_min(preconditioner); + const value_type precond_max = abs_max(preconditioner); const value_type cond_number = precond_max / precond_min; // std::cout << "condition number = " << precond_max << " / " << precond_min // << " = " << cond_number << std::endl; @@ -112,10 +107,10 @@ struct ConjugateGradientSolver { bool converged = false; const unsigned int max_niter = n; value_type rnorm2 = 0.0; - const std::size_t rhs_size = size(b); + const std::size_t rhs_size = volume(b); // starting guess: x_0 = D^-1 . b - XX_i = copy(b); + XX_i = b; vec_multiply(XX_i, preconditioner); // r_0 = b - a(x) @@ -126,19 +121,19 @@ struct ConjugateGradientSolver { if (use_diis) diis.extrapolate(XX_i, RR_i, true); // z_0 = D^-1 . r_0 - ZZ_i = copy(RR_i); + ZZ_i = RR_i; vec_multiply(ZZ_i, preconditioner); // p_0 = z_0 - PP_i = copy(ZZ_i); + PP_i = ZZ_i; unsigned int iter = 0; while (not converged) { // alpha_i = (r_i . z_i) / (p_i . A . p_i) - value_type rz_norm2 = dot_product(RR_i, ZZ_i); + value_type rz_norm2 = dot(RR_i, ZZ_i); a(PP_i, APP_i); - const value_type pAp_i = dot_product(PP_i, APP_i); + const value_type pAp_i = dot(PP_i, APP_i); const value_type alpha_i = rz_norm2 / pAp_i; // x_i += alpha_i p_i @@ -156,10 +151,10 @@ struct ConjugateGradientSolver { } // z_i = D^-1 . r_i - ZZ_i = copy(RR_i); + ZZ_i = RR_i; vec_multiply(ZZ_i, preconditioner); - const value_type rz_ip1_norm2 = dot_product(ZZ_i, RR_i); + const value_type rz_ip1_norm2 = dot(ZZ_i, RR_i); const value_type beta_i = rz_ip1_norm2 / rz_norm2; @@ -173,18 +168,22 @@ struct ConjugateGradientSolver { // std::cout << "iter=" << iter << " dnorm=" << r_ip1_norm << std::endl; if (converged) { - assign(x, XX_i); + x = XX_i; } else if (iter >= max_niter) throw std::domain_error( "ConjugateGradient: max # of iterations exceeded"); } // solver loop - assign(x, XX_i); + x = XX_i; return rnorm2; } }; -}; // namespace TiledArray +} // namespace TiledArray::math::linalg + +namespace TiledArray { + using TiledArray::math::linalg::ConjugateGradientSolver; +} -#endif // TILEDARRAY_ALGEBRA_CONJGRAD_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_CONJGRAD_H__INCLUDED diff --git a/src/TiledArray/algebra/diis.h b/src/TiledArray/math/linalg/diis.h similarity index 93% rename from src/TiledArray/algebra/diis.h rename to src/TiledArray/math/linalg/diis.h index 5b50ede9f3..723be23d4a 100644 --- a/src/TiledArray/algebra/diis.h +++ b/src/TiledArray/math/linalg/diis.h @@ -23,15 +23,17 @@ * */ -#ifndef TILEDARRAY_ALGEBRA_DIIS_H__INCLUDED -#define TILEDARRAY_ALGEBRA_DIIS_H__INCLUDED +#ifndef TILEDARRAY_MATH_LINALG_DIIS_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_DIIS_H__INCLUDED + +#include "TiledArray/dist_array.h" +#include +#include "TiledArray/external/eigen.h" -#include #include #include -#include "../dist_array.h" -namespace TiledArray { +namespace TiledArray::math::linalg { /// DIIS (``direct inversion of iterative subspace'') extrapolation @@ -81,11 +83,11 @@ template class DIIS { public: typedef typename D::element_type value_type; - typedef typename detail::scalar_t scalar_type; + typedef typename TiledArray::detail::scalar_t scalar_type; typedef Eigen::Matrix - EigenMatrixX; - typedef Eigen::Matrix EigenVectorX; + Matrix; + typedef Eigen::Matrix Vector; /// Constructor @@ -140,6 +142,7 @@ class DIIS { /// \param extrapolate_error whether to extrapolate the error (default = /// false). void extrapolate(D& x, D& error, bool extrapolate_error = false) { + iter++; // compute extrapolation coefficients C_ and number of skipped vectors @@ -172,8 +175,9 @@ class DIIS { /// \param nskip number of old vectors to skip (default = 0) /// \param increase_iter whether to increase the diis iteration index /// (default = false) - void extrapolate(D& x, const EigenVectorX& c, unsigned int nskip = 0, + void extrapolate(D& x, const Vector& c, unsigned int nskip = 0, bool increase_iter = false) { + if (increase_iter) { iter++; } @@ -230,6 +234,7 @@ class DIIS { /// (default = false) void compute_extrapolation_parameters(const D& error, bool increase_iter = false) { + if (increase_iter) { iter++; } @@ -238,7 +243,7 @@ class DIIS { if (errors_.size() == ndiis) { // holding max # of vectors already? drop // the least recent error errors_.pop_front(); - EigenMatrixX Bcrop = B_.bottomRightCorner(ndiis - 1, ndiis - 1); + Matrix Bcrop = B_.bottomRightCorner(ndiis - 1, ndiis - 1); Bcrop.conservativeResize(ndiis, ndiis); B_ = Bcrop; } @@ -250,8 +255,8 @@ class DIIS { // and compute the most recent elements of B, B(i,j) = for (unsigned int i = 0; i < nvec - 1; i++) B_(i, nvec - 1) = B_(nvec - 1, i) = - dot_product(errors_[i], errors_[nvec - 1]); - B_(nvec - 1, nvec - 1) = dot_product(errors_[nvec - 1], errors_[nvec - 1]); + dot(errors_[i], errors_[nvec - 1]); + B_(nvec - 1, nvec - 1) = dot(errors_[nvec - 1], errors_[nvec - 1]); using std::abs; using std::sqrt; const auto current_error_2norm = sqrt(abs(B_(nvec - 1, nvec - 1))); @@ -279,13 +284,13 @@ class DIIS { const unsigned int rank = nvec - nskip_ + 1; // size of matrix A // set up the DIIS linear system: A c = rhs - EigenMatrixX A(rank, rank); + Matrix A(rank, rank); C_.resize(rank); A.col(0).setConstant(-1.0); A.row(0).setConstant(-1.0); A(0, 0) = 0.0; - EigenVectorX rhs = EigenVectorX::Zero(rank); + Vector rhs = Vector::Zero(rank); rhs[0] = -1.0; scalar_type norm = 1.0; @@ -310,7 +315,7 @@ class DIIS { #endif // finally, solve the DIIS linear system - Eigen::ColPivHouseholderQR A_QR = A.colPivHouseholderQr(); + Eigen::ColPivHouseholderQR A_QR = A.colPivHouseholderQr(); C_ = A_QR.solve(rhs); absdetA = A_QR.absDeterminant(); @@ -350,7 +355,7 @@ class DIIS { } /// calling this function returns extrapolation coefficients - const EigenVectorX& get_coeffs() { + const Vector& get_coeffs() { TA_USER_ASSERT( parameters_computed_ && C_.size() > 0, "DIIS: empty coefficients, because they have not been computed"); @@ -379,8 +384,8 @@ class DIIS { //!< decreasing damping factor once //!< error 2-norm falls below this - EigenMatrixX B_; //!< B(i,j) = - EigenVectorX C_; //! DIIS coefficients + Matrix B_; //!< B(i,j) = + Vector C_; //! DIIS coefficients bool parameters_computed_; //! whether diis parameters C_ and nskip_ have //! been computed unsigned int nskip_; //! number of skipped vectors in extrapolation @@ -399,7 +404,7 @@ class DIIS { void init() { iter = 0; - B_ = EigenMatrixX::Zero(ndiis, ndiis); + B_ = Matrix::Zero(ndiis, ndiis); C_.resize(0); parameters_computed_ = false; nskip_ = 0; @@ -416,6 +421,10 @@ class DIIS { }; // class DIIS -} // namespace TiledArray +} // namespace TiledArray::math::linalg + +namespace TiledArray { + using TiledArray::math::linalg::DIIS; +} -#endif // TILEDARRAY_ALGEBRA_DIIS_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_DIIS_H__INCLUDED diff --git a/src/TiledArray/math/linalg/forward.h b/src/TiledArray/math/linalg/forward.h new file mode 100644 index 0000000000..2ce19e21f6 --- /dev/null +++ b/src/TiledArray/math/linalg/forward.h @@ -0,0 +1,51 @@ +/* + * This file is a part of TiledArray. + * Copyright (C) 2020 Virginia Tech + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + * + * David Williams-Young + * Computational Research Division, Lawrence Berkeley National Laboratory + * + * svd.h + * Created: 12 June, 2020 + * + */ +#ifndef TILEDARRAY_MATH_LINALG_FORWARD_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_FORWARD_H__INCLUDED + +#include +#include + +namespace TiledArray::math::linalg { + +enum TransposeFlag { NoTranspose, Transpose, ConjTranspose }; + +struct SVD { + enum Vectors { + ValuesOnly, + LeftVectors, + RightVectors, + AllVectors + }; +}; + +} // namespace TiledArray::math::linalg + +namespace TiledArray { + using TiledArray::math::linalg::TransposeFlag; + using TiledArray::math::linalg::SVD; +} + +#endif // TILEDARRAY_MATH_LINALG_FORWARD_H__INCLUDED diff --git a/src/TiledArray/algebra/heig.h b/src/TiledArray/math/linalg/heig.h similarity index 80% rename from src/TiledArray/algebra/heig.h rename to src/TiledArray/math/linalg/heig.h index 23d28329bc..139f8d0266 100644 --- a/src/TiledArray/algebra/heig.h +++ b/src/TiledArray/math/linalg/heig.h @@ -21,16 +21,16 @@ * Created: 16 October, 2020 * */ -#ifndef TILEDARRAY_ALGEBRA_HEIG_H__INCLUDED -#define TILEDARRAY_ALGEBRA_HEIG_H__INCLUDED +#ifndef TILEDARRAY_MATH_LINALG_HEIG_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_HEIG_H__INCLUDED #include #if TILEDARRAY_HAS_SCALAPACK -#include +#include #endif -#include +#include -namespace TiledArray { +namespace TiledArray::math::linalg { template auto heig(const Array& A, TiledRange evec_trange = TiledRange()) { @@ -52,6 +52,10 @@ auto heig(const ArrayA& A, const ArrayB& B, TiledRange evec_trange = TiledRange( return non_distributed::heig(A, B, evec_trange); } -} // namespace TiledArray +} // namespace TiledArray::math::linalg + +namespace TiledArray { + using TiledArray::math::linalg::heig; +} -#endif // TILEDARRAY_ALGEBRA_HEIG_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_HEIG_H__INCLUDED diff --git a/src/TiledArray/algebra/lu.h b/src/TiledArray/math/linalg/lu.h similarity index 80% rename from src/TiledArray/algebra/lu.h rename to src/TiledArray/math/linalg/lu.h index bedb2df59a..e35eaee99b 100644 --- a/src/TiledArray/algebra/lu.h +++ b/src/TiledArray/math/linalg/lu.h @@ -21,14 +21,14 @@ * Created: 16 October, 2020 * */ -#ifndef TILEDARRAY_ALGEBRA_LU_H__INCLUDED -#define TILEDARRAY_ALGEBRA_LU_H__INCLUDED +#ifndef TILEDARRAY_MATH_LINALG_LU_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_LU_H__INCLUDED #include #if TILEDARRAY_HAS_SCALAPACK -#include +#include #endif -#include +#include namespace TiledArray { @@ -40,4 +40,4 @@ using scalapack::lu_solve; } // namespace TiledArray -#endif // TILEDARRAY_ALGEBRA_LU_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_LU_H__INCLUDED diff --git a/src/TiledArray/algebra/non-distributed/cholesky.h b/src/TiledArray/math/linalg/non-distributed/cholesky.h similarity index 81% rename from src/TiledArray/algebra/non-distributed/cholesky.h rename to src/TiledArray/math/linalg/non-distributed/cholesky.h index 037685430a..b72a06e20e 100644 --- a/src/TiledArray/algebra/non-distributed/cholesky.h +++ b/src/TiledArray/math/linalg/non-distributed/cholesky.h @@ -21,21 +21,19 @@ * Created: 16 October, 2020 * */ -#ifndef TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_CHOL_H__INCLUDED -#define TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_CHOL_H__INCLUDED +#ifndef TILEDARRAY_MATH_LINALG_NON_DISTRIBUTED_CHOL_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_NON_DISTRIBUTED_CHOL_H__INCLUDED #include -#include +#include #include -#include +#include -namespace TiledArray::non_distributed { - -namespace detail { +namespace TiledArray::math::linalg::non_distributed { template -auto make_L_eig(const DistArray& A) { +auto rank_local_cholesky(const DistArray& A) { using Array = DistArray; using numeric_type = typename Array::numeric_type; static_assert(std::is_same_v, @@ -43,16 +41,14 @@ auto make_L_eig(const DistArray& A) { "scalar types"); World& world = A.world(); - auto A_eig = detail::to_eigen(A); + auto A_eig = detail::make_matrix(A); if (world.rank() == 0) { - algebra::rank_local::cholesky(A_eig); + linalg::rank_local::cholesky(A_eig); } world.gop.broadcast_serializable(A_eig, 0); return A_eig; } -} // namespace detail - /** * @brief Compute the Cholesky factorization of a HPD rank-2 tensor * @@ -74,7 +70,7 @@ auto make_L_eig(const DistArray& A) { template >> auto cholesky(const Array& A, TiledRange l_trange = TiledRange()) { - auto L_eig = detail::make_L_eig(A); + auto L_eig = rank_local_cholesky(A); detail::zero_out_upper_triangle(L_eig); if (l_trange.rank() == 0) l_trange = A.trange(); return eigen_to_array(A.world(), l_trange, L_eig); @@ -101,10 +97,10 @@ template >> auto cholesky(const ContiguousTensor& A) { - auto A_eig = detail::to_eigen(A); - algebra::rank_local::cholesky(A_eig); + auto A_eig = detail::make_matrix(A); + linalg::rank_local::cholesky(A_eig); detail::zero_out_upper_triangle(A_eig); - return detail::from_eigen(A_eig, A.range()); + return detail::make_array(A_eig, A.range()); } /** @@ -128,26 +124,26 @@ auto cholesky(const ContiguousTensor& A) { * @returns The inverse lower triangular Cholesky factor in TA format * @note this is a collective operation with respect to the world of @p A */ -template >> auto cholesky_linv(const Array& A, TiledRange l_trange = TiledRange()) { World& world = A.world(); - auto L_eig = detail::make_L_eig(A); - if constexpr (RetL) detail::zero_out_upper_triangle(L_eig); + auto L_eig = rank_local_cholesky(A); + if constexpr (Both) detail::zero_out_upper_triangle(L_eig); // if need to return L use its copy to compute inverse decltype(L_eig) L_inv_eig; if (world.rank() == 0) { - if (RetL) L_inv_eig = L_eig; - auto& L_inv_eig_ref = RetL ? L_inv_eig : L_eig; - algebra::rank_local::cholesky_linv(L_inv_eig_ref); + if (Both) L_inv_eig = L_eig; + auto& L_inv_eig_ref = Both ? L_inv_eig : L_eig; + linalg::rank_local::cholesky_linv(L_inv_eig_ref); detail::zero_out_upper_triangle(L_inv_eig_ref); } - world.gop.broadcast_serializable(RetL ? L_inv_eig : L_eig, 0); + world.gop.broadcast_serializable(Both ? L_inv_eig : L_eig, 0); if (l_trange.rank() == 0) l_trange = A.trange(); - if constexpr (RetL) + if constexpr (Both) return std::make_tuple(eigen_to_array(world, l_trange, L_eig), eigen_to_array(world, l_trange, L_inv_eig)); else @@ -163,11 +159,11 @@ auto cholesky_solve(const Array& A, const Array& B, "TA::lapack::{cholesky*} are only usable with a DistArray of " "scalar types"); - auto A_eig = detail::to_eigen(A); - auto X_eig = detail::to_eigen(B); + auto A_eig = detail::make_matrix(A); + auto X_eig = detail::make_matrix(B); World& world = A.world(); if (world.rank() == 0) { - algebra::rank_local::cholesky_solve(A_eig, X_eig); + linalg::rank_local::cholesky_solve(A_eig, X_eig); } world.gop.broadcast_serializable(X_eig, 0); if (x_trange.rank() == 0) x_trange = B.trange(); @@ -180,7 +176,7 @@ auto cholesky_lsolve(TransposeFlag transpose, const Array& A, const Array& B, TiledRange l_trange = TiledRange(), TiledRange x_trange = TiledRange()) { World& world = A.world(); - auto L_eig = detail::make_L_eig(A); + auto L_eig = rank_local_cholesky(A); detail::zero_out_upper_triangle(L_eig); using numeric_type = typename Array::numeric_type; @@ -188,9 +184,9 @@ auto cholesky_lsolve(TransposeFlag transpose, const Array& A, const Array& B, "TA::lapack::{cholesky*} are only usable with a DistArray of " "scalar types"); - auto X_eig = detail::to_eigen(B); + auto X_eig = detail::make_matrix(B); if (world.rank() == 0) { - algebra::rank_local::cholesky_lsolve(transpose, L_eig, X_eig); + linalg::rank_local::cholesky_lsolve(transpose, L_eig, X_eig); } world.gop.broadcast_serializable(X_eig, 0); if (l_trange.rank() == 0) l_trange = A.trange(); @@ -199,6 +195,6 @@ auto cholesky_lsolve(TransposeFlag transpose, const Array& A, const Array& B, eigen_to_array(world, x_trange, X_eig)); } -} // namespace TiledArray +} // namespace TiledArray::math::linalg -#endif // TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_CHOL_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_NON_DISTRIBUTED_CHOL_H__INCLUDED diff --git a/src/TiledArray/algebra/non-distributed/heig.h b/src/TiledArray/math/linalg/non-distributed/heig.h similarity index 84% rename from src/TiledArray/algebra/non-distributed/heig.h rename to src/TiledArray/math/linalg/non-distributed/heig.h index 484f0bdb20..6a7cccac3a 100644 --- a/src/TiledArray/algebra/non-distributed/heig.h +++ b/src/TiledArray/math/linalg/non-distributed/heig.h @@ -21,16 +21,16 @@ * Created: 19 October, 2020 * */ -#ifndef TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_HEIG_H__INCLUDED -#define TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_HEIG_H__INCLUDED +#ifndef TILEDARRAY_MATH_LINALG_NON_DISTRIBUTED_HEIG_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_NON_DISTRIBUTED_HEIG_H__INCLUDED #include -#include -#include +#include +#include #include -namespace TiledArray::non_distributed { +namespace TiledArray::math::linalg::non_distributed { /** * @brief Solve the standard eigenvalue problem with LAPACK @@ -54,10 +54,10 @@ template auto heig(const Array& A, TiledRange evec_trange = TiledRange()) { using numeric_type = typename detail::array_traits::numeric_type; World& world = A.world(); - auto A_eig = detail::to_eigen(A); + auto A_eig = detail::make_matrix(A); std::vector evals; if (world.rank() == 0) { - algebra::rank_local::heig(A_eig, evals); + linalg::rank_local::heig(A_eig, evals); } world.gop.broadcast_serializable(A_eig, 0); world.gop.broadcast_serializable(evals, 0); @@ -96,11 +96,11 @@ auto heig(const ArrayA& A, const ArrayB& B, TiledRange evec_trange = TiledRange( using numeric_type = typename detail::array_traits::numeric_type; (void)detail::array_traits{}; World& world = A.world(); - auto A_eig = detail::to_eigen(A); - auto B_eig = detail::to_eigen(B); + auto A_eig = detail::make_matrix(A); + auto B_eig = detail::make_matrix(B); std::vector evals; if (world.rank() == 0) { - algebra::rank_local::heig(A_eig, B_eig, evals); + linalg::rank_local::heig(A_eig, B_eig, evals); } world.gop.broadcast_serializable(A_eig, 0); world.gop.broadcast_serializable(evals, 0); @@ -111,6 +111,6 @@ auto heig(const ArrayA& A, const ArrayB& B, TiledRange evec_trange = TiledRange( ); } -} // namespace TiledArray::lapack +} // namespace TiledArray::math::linalg -#endif // TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_HEIG_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_NON_DISTRIBUTED_HEIG_H__INCLUDED diff --git a/src/TiledArray/algebra/non-distributed/lu.h b/src/TiledArray/math/linalg/non-distributed/lu.h similarity index 76% rename from src/TiledArray/algebra/non-distributed/lu.h rename to src/TiledArray/math/linalg/non-distributed/lu.h index 76beba06ac..d1b06bbb1c 100644 --- a/src/TiledArray/algebra/non-distributed/lu.h +++ b/src/TiledArray/math/linalg/non-distributed/lu.h @@ -22,16 +22,16 @@ * Created: 19 June, 2020 * */ -#ifndef TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_LU_H__INCLUDED -#define TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_LU_H__INCLUDED +#ifndef TILEDARRAY_MATH_LINALG_NON_DISTRIBUTED_LU_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_NON_DISTRIBUTED_LU_H__INCLUDED #include -#include -#include +#include +#include #include -namespace TiledArray::non_distributed { +namespace TiledArray::math::linalg::non_distributed { /** * @brief Solve a linear system via LU factorization @@ -41,10 +41,10 @@ auto lu_solve(const ArrayA& A, const ArrayB& B, TiledRange x_trange = TiledRange (void)detail::array_traits{}; (void)detail::array_traits{}; auto& world = A.world(); - auto A_eig = detail::to_eigen(A); - auto B_eig = detail::to_eigen(B); + auto A_eig = detail::make_matrix(A); + auto B_eig = detail::make_matrix(B); if (world.rank() == 0) { - algebra::rank_local::lu_solve(A_eig, B_eig); + linalg::rank_local::lu_solve(A_eig, B_eig); } world.gop.broadcast_serializable(B_eig, 0); if (x_trange.rank() == 0) x_trange = B.trange(); @@ -58,15 +58,15 @@ template auto lu_inv(const Array& A, TiledRange ainv_trange = TiledRange()) { (void)detail::array_traits{}; auto& world = A.world(); - auto A_eig = detail::to_eigen(A); + auto A_eig = detail::make_matrix(A); if (world.rank() == 0) { - algebra::rank_local::lu_inv(A_eig); + linalg::rank_local::lu_inv(A_eig); } world.gop.broadcast_serializable(A_eig, 0); if (ainv_trange.rank() == 0) ainv_trange = A.trange(); return eigen_to_array(A.world(), ainv_trange, A_eig); } -} // namespace TiledArray::lapack +} // namespace TiledArray::math::linalg::lapack -#endif // TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_LU_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_NON_DISTRIBUTED_LU_H__INCLUDED diff --git a/src/TiledArray/algebra/non-distributed/svd.h b/src/TiledArray/math/linalg/non-distributed/svd.h similarity index 76% rename from src/TiledArray/algebra/non-distributed/svd.h rename to src/TiledArray/math/linalg/non-distributed/svd.h index 01f3ef8eee..9c146784ef 100644 --- a/src/TiledArray/algebra/non-distributed/svd.h +++ b/src/TiledArray/math/linalg/non-distributed/svd.h @@ -22,16 +22,16 @@ * Created: 12 June, 2020 * */ -#ifndef TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_SVD_H__INCLUDED -#define TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_SVD_H__INCLUDED +#ifndef TILEDARRAY_MATH_LINALG_NON_DISTRIBUTED_SVD_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_NON_DISTRIBUTED_SVD_H__INCLUDED #include -#include -#include +#include +#include #include -namespace TiledArray::non_distributed { +namespace TiledArray::math::linalg::non_distributed { /** * @brief Compute the singular value decomposition (SVD) via ScaLAPACK @@ -55,19 +55,18 @@ namespace TiledArray::non_distributed { * @returns A tuple containing the eigenvalues and eigenvectors of input array * as std::vector and in TA format, respectively. */ -template > +template auto svd(const Array& A, TiledRange u_trange = TiledRange(), TiledRange vt_trange = TiledRange()) { using T = typename Array::numeric_type; - using Matrix = algebra::rank_local::Matrix; + using Matrix = linalg::rank_local::Matrix; World& world = A.world(); - auto A_eig = detail::to_eigen(A); + auto A_eig = detail::make_matrix(A); - constexpr bool svd_all_vectors = std::is_same_v; - constexpr bool need_u = std::is_same_v or svd_all_vectors; - constexpr bool need_vt = std::is_same_v or svd_all_vectors; + constexpr bool svd_all_vectors = (Vectors == SVD::AllVectors); + constexpr bool need_u = (Vectors == SVD::LeftVectors) or svd_all_vectors; + constexpr bool need_vt = (Vectors == SVD::RightVectors) or svd_all_vectors; std::vector S; std::unique_ptr U, VT; @@ -76,7 +75,7 @@ auto svd(const Array& A, TiledRange u_trange = TiledRange(), TiledRange vt_trang if constexpr (need_vt) VT = std::make_unique(); if (world.rank() == 0) { - algebra::rank_local::svd(A_eig, S, U.get(), VT.get()); + linalg::rank_local::svd(A_eig, S, U.get(), VT.get()); } world.gop.broadcast_serializable(S, 0); @@ -101,6 +100,6 @@ auto svd(const Array& A, TiledRange u_trange = TiledRange(), TiledRange vt_trang } -} // namespace TiledArray::non_distributed +} // namespace TiledArray::math::linalg::non_distributed -#endif // TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_SVD_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_NON_DISTRIBUTED_SVD_H__INCLUDED diff --git a/src/TiledArray/algebra/rank-local.cpp b/src/TiledArray/math/linalg/rank-local.cpp similarity index 97% rename from src/TiledArray/algebra/rank-local.cpp rename to src/TiledArray/math/linalg/rank-local.cpp index 20fd9dcaa0..48f286eb92 100644 --- a/src/TiledArray/algebra/rank-local.cpp +++ b/src/TiledArray/math/linalg/rank-local.cpp @@ -22,7 +22,7 @@ * */ -#include +#include #include @@ -46,7 +46,7 @@ else std::abort(); \ } -namespace TiledArray::algebra::rank_local { +namespace TiledArray::math::linalg::rank_local { template void cholesky(Matrix& A) { @@ -215,4 +215,4 @@ void lu_inv(Matrix& A) { TA_LAPACK_EXPLICIT(Matrix, std::vector); TA_LAPACK_EXPLICIT(Matrix, std::vector); -} // namespace TiledArray::local +} // namespace TiledArray::math::linalg::rank_local diff --git a/src/TiledArray/algebra/rank-local.h b/src/TiledArray/math/linalg/rank-local.h similarity index 73% rename from src/TiledArray/algebra/rank-local.h rename to src/TiledArray/math/linalg/rank-local.h index a69beba83d..220ad8461f 100644 --- a/src/TiledArray/algebra/rank-local.h +++ b/src/TiledArray/math/linalg/rank-local.h @@ -1,14 +1,14 @@ -#ifndef TILEDARRAY_ALGEBRA_RANK_LOCAL_H__INCLUDED -#define TILEDARRAY_ALGEBRA_RANK_LOCAL_H__INCLUDED +#ifndef TILEDARRAY_MATH_LINALG_RANK_LOCAL_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_RANK_LOCAL_H__INCLUDED #include -#include +#include #include #include -namespace TiledArray::algebra::rank_local { +namespace TiledArray::math::linalg::rank_local { template using Matrix = ::Eigen::Matrix; @@ -26,8 +26,7 @@ template void cholesky_solve(Matrix &A, Matrix &X); template -void cholesky_lsolve(TransposeFlag transpose, - Matrix &A, Matrix &X); +void cholesky_lsolve(TransposeFlag transpose, Matrix &A, Matrix &X); template void heig(Matrix &A, std::vector &W); @@ -47,4 +46,4 @@ void lu_inv(Matrix &A); } // namespace TiledArray::local -#endif // TILEDARRAY_ALGEBRA_RANK_LOCAL_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_RANK_LOCAL_H__INCLUDED diff --git a/src/TiledArray/algebra/scalapack/all.h b/src/TiledArray/math/linalg/scalapack/all.h similarity index 71% rename from src/TiledArray/algebra/scalapack/all.h rename to src/TiledArray/math/linalg/scalapack/all.h index fbba6826d0..3b35e7b3ff 100644 --- a/src/TiledArray/algebra/scalapack/all.h +++ b/src/TiledArray/math/linalg/scalapack/all.h @@ -23,15 +23,15 @@ * */ -#ifndef TILEDARRAY_ALGEBRA_SCALAPACK_ALL_H__INCLUDED -#define TILEDARRAY_ALGEBRA_SCALAPACK_ALL_H__INCLUDED +#ifndef TILEDARRAY_MATH_LINALG_SCALAPACK_ALL_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_SCALAPACK_ALL_H__INCLUDED #include #if TILEDARRAY_HAS_SCALAPACK -#include -#include -#include -#include +#include +#include +#include +#include #endif -#endif // TILEDARRAY_ALGEBRA_SCALAPACK_ALL_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_SCALAPACK_ALL_H__INCLUDED diff --git a/src/TiledArray/conversions/block_cyclic.h b/src/TiledArray/math/linalg/scalapack/block_cyclic.h similarity index 97% rename from src/TiledArray/conversions/block_cyclic.h rename to src/TiledArray/math/linalg/scalapack/block_cyclic.h index f723717ec6..55f687f96a 100644 --- a/src/TiledArray/conversions/block_cyclic.h +++ b/src/TiledArray/math/linalg/scalapack/block_cyclic.h @@ -24,8 +24,8 @@ * */ -#ifndef TILEDARRAY_CONVERSIONS_TO_BLOCKCYCLIC_H__INCLUDED -#define TILEDARRAY_CONVERSIONS_TO_BLOCKCYCLIC_H__INCLUDED +#ifndef TILEDARRAY_MATH_LINALG_SCALAPACK_TO_BLOCKCYCLIC_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_SCALAPACK_TO_BLOCKCYCLIC_H__INCLUDED #include #if TILEDARRAY_HAS_SCALAPACK @@ -41,8 +41,7 @@ #include #include -namespace TiledArray { -namespace scalapack { +namespace TiledArray::math::linalg::scalapack { template > @@ -298,8 +297,7 @@ std::remove_cv_t block_cyclic_to_array( return matrix.template tensor_from_matrix>(trange); } -} // namespace scalapack } // namespace TiledArray #endif // TILEDARRAY_HAS_SCALAPACK -#endif // TILEDARRAY_CONVERSIONS_TO_BLOCKCYCLIC_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_SCALAPACK_TO_BLOCKCYCLIC_H__INCLUDED diff --git a/src/TiledArray/algebra/scalapack/cholesky.h b/src/TiledArray/math/linalg/scalapack/cholesky.h similarity index 94% rename from src/TiledArray/algebra/scalapack/cholesky.h rename to src/TiledArray/math/linalg/scalapack/cholesky.h index 3556c81ecd..a11080f209 100644 --- a/src/TiledArray/algebra/scalapack/cholesky.h +++ b/src/TiledArray/math/linalg/scalapack/cholesky.h @@ -22,23 +22,21 @@ * Created: 8 June, 2020 * */ -#ifndef TILEDARRAY_ALGEBRA_SCALAPACK_CHOL_H__INCLUDED -#define TILEDARRAY_ALGEBRA_SCALAPACK_CHOL_H__INCLUDED +#ifndef TILEDARRAY_MATH_LINALG_SCALAPACK_CHOL_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_SCALAPACK_CHOL_H__INCLUDED #include #if TILEDARRAY_HAS_SCALAPACK -#include -#include -#include +#include +#include #include #include #include #include -namespace TiledArray { -namespace scalapack { +namespace TiledArray::math::linalg::scalapack { /** * @brief Compute the Cholesky factorization of a HPD rank-2 tensor @@ -103,7 +101,7 @@ auto cholesky(const Array& A, TiledRange l_trange = TiledRange(), * auto [L,Linv] = cholesky_Linv(A, ...) * * @tparam Array Input array type, must be convertible to BlockCyclicMatrix - * @tparam RetL Whether or not to return the cholesky factor + * @tparam Both Whether or not to return the cholesky factor * * @param[in] A Input array to be diagonalized. Must be rank-2 * @param[in] l_trange TiledRange for resulting inverse Cholesky factor. @@ -112,7 +110,7 @@ auto cholesky(const Array& A, TiledRange l_trange = TiledRange(), * * @returns The inverse lower triangular Cholesky factor in TA format */ -template + template auto cholesky_linv(const Array& A, TiledRange l_trange = TiledRange(), size_t NB = default_block_size()) { using value_type = typename Array::element_type; @@ -140,7 +138,7 @@ auto cholesky_linv(const Array& A, TiledRange l_trange = TiledRange(), // Copy L if needed std::shared_ptr> L_sca = nullptr; - if constexpr (RetL) { + if constexpr (Both) { L_sca = std::make_shared>( world, grid, N, N, NB, NB); L_sca->local_mat() = matrix.local_mat(); @@ -158,7 +156,7 @@ auto cholesky_linv(const Array& A, TiledRange l_trange = TiledRange(), auto Linv = scalapack::block_cyclic_to_array(matrix, l_trange); world.gop.fence(); - if constexpr (RetL) { + if constexpr (Both) { auto L = scalapack::block_cyclic_to_array(*L_sca, l_trange); world.gop.fence(); return std::tuple(L, Linv); @@ -276,8 +274,7 @@ auto cholesky_lsolve(TransposeFlag trans, const Array& A, const Array& B, return std::tuple(L, X); } -} // namespace scalapack -} // namespace TiledArray +} // namespace TiledArray::math::linalg::scalapack #endif // TILEDARRAY_HAS_SCALAPACK -#endif // TILEDARRAY_ALGEBRA_SCALAPACK_CHOL_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_SCALAPACK_CHOL_H__INCLUDED diff --git a/src/TiledArray/algebra/scalapack/heig.h b/src/TiledArray/math/linalg/scalapack/heig.h similarity index 94% rename from src/TiledArray/algebra/scalapack/heig.h rename to src/TiledArray/math/linalg/scalapack/heig.h index 488867b857..d6a6d26984 100644 --- a/src/TiledArray/algebra/scalapack/heig.h +++ b/src/TiledArray/math/linalg/scalapack/heig.h @@ -23,18 +23,18 @@ * Edited: 8 June, 2020 * */ -#ifndef TILEDARRAY_ALGEBRA_SCALAPACK_HEIG_H__INCLUDED -#define TILEDARRAY_ALGEBRA_SCALAPACK_HEIG_H__INCLUDED +#ifndef TILEDARRAY_MATH_LINALG_SCALAPACK_HEIG_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_SCALAPACK_HEIG_H__INCLUDED #include #if TILEDARRAY_HAS_SCALAPACK -#include +#include + #include #include -namespace TiledArray { -namespace scalapack { +namespace TiledArray::math::linalg::scalapack { /** * @brief Solve the standard eigenvalue problem with ScaLAPACK @@ -165,8 +165,7 @@ auto heig(const ArrayA& A, const ArrayB& B, return std::tuple(evals, evecs_ta); } -} // namespace scalapack -} // namespace TiledArray +} // namespace TiledArray::math::linalg::scalapack #endif // TILEDARRAY_HAS_SCALAPACK -#endif // TILEDARRAY_ALGEBRA_SCALAPACK_HEIG_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_SCALAPACK_HEIG_H__INCLUDED diff --git a/src/TiledArray/algebra/scalapack/lu.h b/src/TiledArray/math/linalg/scalapack/lu.h similarity index 92% rename from src/TiledArray/algebra/scalapack/lu.h rename to src/TiledArray/math/linalg/scalapack/lu.h index 8c65383846..56b06e8f38 100644 --- a/src/TiledArray/algebra/scalapack/lu.h +++ b/src/TiledArray/math/linalg/scalapack/lu.h @@ -22,21 +22,19 @@ * Created: 19 June, 2020 * */ -#ifndef TILEDARRAY_ALGEBRA_SCALAPACK_LU_H__INCLUDED -#define TILEDARRAY_ALGEBRA_SCALAPACK_LU_H__INCLUDED +#ifndef TILEDARRAY_MATH_LINALG_SCALAPACK_LU_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_SCALAPACK_LU_H__INCLUDED #include #if TILEDARRAY_HAS_SCALAPACK -#include -#include +#include #include #include #include -namespace TiledArray { -namespace scalapack { +namespace TiledArray::math::linalg::scalapack { /** * @brief Solve a linear system via LU factorization @@ -129,8 +127,7 @@ auto lu_inv(const Array& A, TiledRange ainv_trange = TiledRange(), return Ainv; } -} // namespace scalapack -} // namespace TiledArray +} // namespace TiledArray::math::linalg::scalapack #endif // TILEDARRAY_HAS_SCALAPACK -#endif // TILEDARRAY_ALGEBRA_SCALAPACK_LU_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_SCALAPACK_LU_H__INCLUDED diff --git a/src/TiledArray/algebra/scalapack/svd.h b/src/TiledArray/math/linalg/scalapack/svd.h similarity index 87% rename from src/TiledArray/algebra/scalapack/svd.h rename to src/TiledArray/math/linalg/scalapack/svd.h index d73f9e05dd..29315777d8 100644 --- a/src/TiledArray/algebra/scalapack/svd.h +++ b/src/TiledArray/math/linalg/scalapack/svd.h @@ -22,18 +22,18 @@ * Created: 12 June, 2020 * */ -#ifndef TILEDARRAY_ALGEBRA_SCALAPACK_SVD_H__INCLUDED -#define TILEDARRAY_ALGEBRA_SCALAPACK_SVD_H__INCLUDED +#ifndef TILEDARRAY_MATH_LINALG_SCALAPACK_SVD_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_SCALAPACK_SVD_H__INCLUDED #include #if TILEDARRAY_HAS_SCALAPACK -#include -#include +#include +#include + #include -namespace TiledArray { -namespace scalapack { +namespace TiledArray::math::linalg::scalapack { /** * @brief Compute the singular value decomposition (SVD) via ScaLAPACK @@ -59,8 +59,7 @@ namespace scalapack { * @returns A tuple containing the eigenvalues and eigenvectors of input array * as std::vector and in TA format, respectively. */ -template > +template auto svd(const Array& A, TiledRange u_trange, TiledRange vt_trange, size_t MB = default_block_size(), size_t NB = default_block_size()) { using value_type = typename Array::element_type; @@ -88,9 +87,9 @@ auto svd(const Array& A, TiledRange u_trange, TiledRange vt_trange, std::vector S(SVD_SIZE); - constexpr bool need_uv = std::is_same_v; - constexpr bool need_u = std::is_same_v or need_uv; - constexpr bool need_vt = std::is_same_v or need_uv; + constexpr bool need_uv = (Vectors == SVD::AllVectors); + constexpr bool need_u = (Vectors == SVD::LeftVectors) or need_uv; + constexpr bool need_vt = (Vectors == SVD::RightVectors) or need_uv; std::shared_ptr> U = nullptr, VT = nullptr; @@ -148,8 +147,7 @@ auto svd(const Array& A, TiledRange u_trange, TiledRange vt_trange, } } -} // namespace scalapack -} // namespace TiledArray +} // namespace TiledArray::math::linalg::scalapack #endif // TILEDARRAY_HAS_SCALAPACK -#endif // TILEDARRAY_ALGEBRA_SCALAPACK_SVD_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_SCALAPACK_SVD_H__INCLUDED diff --git a/src/TiledArray/algebra/scalapack/util.h b/src/TiledArray/math/linalg/scalapack/util.h similarity index 87% rename from src/TiledArray/algebra/scalapack/util.h rename to src/TiledArray/math/linalg/scalapack/util.h index 439d2fa590..ea4c76fe1e 100644 --- a/src/TiledArray/algebra/scalapack/util.h +++ b/src/TiledArray/math/linalg/scalapack/util.h @@ -22,18 +22,16 @@ * Created: 19 June, 2020 * */ -#ifndef TILEDARRAY_ALGEBRA_SCALAPACK_UTIL_H__INCLUDED -#define TILEDARRAY_ALGEBRA_SCALAPACK_UTIL_H__INCLUDED +#ifndef TILEDARRAY_MATH_LINALG_SCALAPACK_UTIL_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_SCALAPACK_UTIL_H__INCLUDED #include #if TILEDARRAY_HAS_SCALAPACK -#include -#include +#include +#include -namespace TiledArray { - -namespace scalapack { +namespace TiledArray::math::linalg::scalapack { inline scalapackpp::TransposeFlag to_scalapackpp_transposeflag( TransposeFlag t) { @@ -97,8 +95,7 @@ inline void set_default_block_size(std::size_t NB) { detail::default_block_size_accessor() = NB; } -} // namespace scalapack -} // namespace TiledArray +} // namespace TiledArray::math::linalg::scalapack #endif // TILEDARRAY_HAS_SCALAPACK -#endif // TILEDARRAY_ALGEBRA_SCALAPACK_UTIL_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_SCALAPACK_UTIL_H__INCLUDED diff --git a/src/TiledArray/algebra/svd.h b/src/TiledArray/math/linalg/svd.h similarity index 67% rename from src/TiledArray/algebra/svd.h rename to src/TiledArray/math/linalg/svd.h index bd10a289a4..3beab1def3 100644 --- a/src/TiledArray/algebra/svd.h +++ b/src/TiledArray/math/linalg/svd.h @@ -21,28 +21,31 @@ * Created: 16 October, 2020 * */ -#ifndef TILEDARRAY_ALGEBRA_SVD_H__INCLUDED -#define TILEDARRAY_ALGEBRA_SVD_H__INCLUDED +#ifndef TILEDARRAY_MATH_LINALG_SVD_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_SVD_H__INCLUDED #include #ifdef TILEDARRAY_HAS_SCALAPACK -#include +#include #endif // TILEDARRAY_HAS_SCALAPACK -#include +#include -namespace TiledArray { +namespace TiledArray::math::linalg { -template > +template auto svd(const Array& A, TiledRange u_trange = TiledRange(), TiledRange vt_trange = TiledRange()) { #if TILEDARRAY_HAS_SCALAPACK if (A.world().size() > 1 && A.range().volume() > 10000000) { - return scalapack::svd(A, u_trange, vt_trange); + return scalapack::svd(A, u_trange, vt_trange); } #endif - return non_distributed::svd(A, u_trange, vt_trange); + return non_distributed::svd(A, u_trange, vt_trange); } -} // namespace TiledArray +} // namespace TiledArray::math::linalg + +namespace TiledArray { + using TiledArray::math::linalg::svd; +} -#endif // TILEDARRAY_ALGEBRA_SVD_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_SVD_H__INCLUDED diff --git a/src/TiledArray/algebra/non-distributed/util.h b/src/TiledArray/math/linalg/util.h similarity index 79% rename from src/TiledArray/algebra/non-distributed/util.h rename to src/TiledArray/math/linalg/util.h index ec72ce963d..e86db719bd 100644 --- a/src/TiledArray/algebra/non-distributed/util.h +++ b/src/TiledArray/math/linalg/util.h @@ -1,6 +1,6 @@ /* * This file is a part of TiledArray. - * Copyright (C) 2020 Virginia Tech + * Copyright (C) 2013 Virginia Tech * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -16,18 +16,21 @@ * along with this program. If not, see . * * Eduard Valeyev + * Department of Chemistry, Virginia Tech * * util.h - * Created: 19 October, 2020 + * May 20, 2013 * */ -#ifndef TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_UTIL_H__INCLUDED -#define TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_UTIL_H__INCLUDED + +#ifndef TILEDARRAY_MATH_LINALG_UTIL_H__INCLUDED +#define TILEDARRAY_MATH_LINALG_UTIL_H__INCLUDED #include +#include "TiledArray/dist_array.h" #include -namespace TiledArray::non_distributed::detail { +namespace TiledArray::math::linalg::detail { template struct array_traits { @@ -40,7 +43,7 @@ struct array_traits { template -auto to_eigen(const DistArray& A) { +auto make_matrix(const DistArray& A) { auto A_repl = A; A_repl.make_replicated(); return array_to_eigen(A_repl); @@ -49,7 +52,7 @@ auto to_eigen(const DistArray& A) { template >> -auto to_eigen(const ContiguousTensor& A) { +auto make_matrix(const ContiguousTensor& A) { using numeric_type = TiledArray::detail::numeric_t; static_assert( std::is_same_v, @@ -79,20 +82,16 @@ auto to_eigen(const ContiguousTensor& A) { return result; } -template >> -auto from_eigen( - const Eigen::Matrix& A, - typename ContiguousTensor::range_type range = {}) { +auto make_array(const Eigen::Matrix& A, + typename ContiguousTensor::range_type range = {}) { using numeric_type = TiledArray::detail::numeric_t; static_assert( std::is_same_v, - "TA::lapack::{cholesky*} are only usable with a ContiguousTensor of " - "scalar types"); + "TA::math::linalg is only usable with a ContiguousTensor of scalar types" + ); using range_type = typename ContiguousTensor::range_type; if (range.rank() == 0) range = range_type(A.rows(), A.cols()); @@ -109,6 +108,6 @@ void zero_out_upper_triangle(Eigen::MatrixBase& A) { A.template triangularView().setZero(); } -} // namespace TiledArray +} // namespace TiledArray::math::linalg -#endif // TILEDARRAY_ALGEBRA_NON_DISTRIBUTED_UTIL_H__INCLUDED +#endif // TILEDARRAY_MATH_LINALG_UTIL_H__INCLUDED diff --git a/src/TiledArray/math/scalapack.h b/src/TiledArray/math/scalapack.h index c7ed2c3a73..7a135edce6 100644 --- a/src/TiledArray/math/scalapack.h +++ b/src/TiledArray/math/scalapack.h @@ -26,8 +26,8 @@ #define TILEDARRAY_MATH_SCALAPACK_H__INCLUDED #warning \ - "TiledArray/math/scalapack.h header is deprecated, please include TiledArray/algebra/scalapack/all.h" + "TiledArray/math/scalapack.h header is deprecated, please include TiledArray/math/linalg/scalapack/all.h" -#include +#include #endif diff --git a/src/tiledarray.h b/src/tiledarray.h index a20d1c86bc..916c228cd1 100644 --- a/src/tiledarray.h +++ b/src/tiledarray.h @@ -55,13 +55,8 @@ #include // Linear algebra -#include -#include +#include -// ScaLAPACK functions -#ifdef TILEDARRAY_HAS_SCALAPACK -#include -#include -#endif +#include #endif // TILEDARRAY_H__INCLUDED diff --git a/tests/linear_algebra.cpp b/tests/linear_algebra.cpp index 87e5ae2264..d27ff9f554 100644 --- a/tests/linear_algebra.cpp +++ b/tests/linear_algebra.cpp @@ -4,23 +4,22 @@ //#include "range_fixture.h" #include "unit_test_config.h" -#include "TiledArray/algebra/non-distributed/cholesky.h" -#include "TiledArray/algebra/non-distributed/heig.h" -#include "TiledArray/algebra/non-distributed/lu.h" -#include "TiledArray/algebra/non-distributed/svd.h" +#include "TiledArray/math/linalg/non-distributed/cholesky.h" +#include "TiledArray/math/linalg/non-distributed/heig.h" +#include "TiledArray/math/linalg/non-distributed/lu.h" +#include "TiledArray/math/linalg/non-distributed/svd.h" -#include "TiledArray/algebra/cholesky.h" -#include "TiledArray/algebra/heig.h" -#include "TiledArray/algebra/lu.h" -#include "TiledArray/algebra/svd.h" +#include "TiledArray/math/linalg/cholesky.h" +#include "TiledArray/math/linalg/heig.h" +#include "TiledArray/math/linalg/lu.h" +#include "TiledArray/math/linalg/svd.h" namespace TA = TiledArray; -namespace non_dist = TA::non_distributed; +namespace non_dist = TA::math::linalg::non_distributed; #if TILEDARRAY_HAS_SCALAPACK - -#include "TiledArray/algebra/scalapack/all.h" -namespace scalapack = TA::scalapack; +namespace scalapack = TA::math::linalg::scalapack; +#include "TiledArray/math/linalg/scalapack/all.h" #define TILEDARRAY_SCALAPACK_TEST(F, E) \ GlobalFixture::world->gop.fence(); \ compare("TiledArray::scalapack", non_dist::F, scalapack::F, E); \ @@ -551,7 +550,7 @@ BOOST_AUTO_TEST_CASE(cholesky_linv) { }); decltype(A) Acopy = A.clone(); - auto Linv = non_dist::cholesky_linv(A); + auto Linv = TA::cholesky_linv(A); BOOST_CHECK(Linv.trange() == A.trange()); @@ -575,7 +574,7 @@ BOOST_AUTO_TEST_CASE(cholesky_linv) { BOOST_CHECK_SMALL(norm, epsilon); - TILEDARRAY_SCALAPACK_TEST(cholesky_linv(Acopy), epsilon); + TILEDARRAY_SCALAPACK_TEST(cholesky_linv(Acopy), epsilon); GlobalFixture::world->gop.fence(); } @@ -591,7 +590,7 @@ BOOST_AUTO_TEST_CASE(cholesky_linv_retl) { return this->make_ta_reference(t, range); }); - auto [L, Linv] = non_dist::cholesky_linv(A); + auto [L, Linv] = TA::cholesky_linv(A); BOOST_CHECK(Linv.trange() == A.trange()); BOOST_CHECK(L.trange() == A.trange()); @@ -615,7 +614,7 @@ BOOST_AUTO_TEST_CASE(cholesky_linv_retl) { BOOST_CHECK_SMALL(norm, epsilon); - TILEDARRAY_SCALAPACK_TEST(cholesky_linv(A), epsilon); + //TILEDARRAY_SCALAPACK_TEST(cholesky_linv(A), epsilon); GlobalFixture::world->gop.fence(); } @@ -778,7 +777,7 @@ BOOST_AUTO_TEST_CASE(svd_values_only) { return this->make_ta_reference(t, range); }); - auto S = non_dist::svd(ref_ta, trange, trange); + auto S = non_dist::svd(ref_ta, trange, trange); std::vector exact_singular_values = exact_evals; std::sort(exact_singular_values.begin(), exact_singular_values.end(), @@ -802,7 +801,7 @@ BOOST_AUTO_TEST_CASE(svd_leftvectors) { return this->make_ta_reference(t, range); }); - auto [S, U] = non_dist::svd(ref_ta, trange, trange); + auto [S, U] = non_dist::svd(ref_ta, trange, trange); std::vector exact_singular_values = exact_evals; std::sort(exact_singular_values.begin(), exact_singular_values.end(), @@ -826,7 +825,7 @@ BOOST_AUTO_TEST_CASE(svd_rightvectors) { return this->make_ta_reference(t, range); }); - auto [S, VT] = non_dist::svd(ref_ta, trange, trange); + auto [S, VT] = non_dist::svd(ref_ta, trange, trange); std::vector exact_singular_values = exact_evals; std::sort(exact_singular_values.begin(), exact_singular_values.end(), @@ -850,7 +849,7 @@ BOOST_AUTO_TEST_CASE(svd_allvectors) { return this->make_ta_reference(t, range); }); - auto [S, U, VT] = non_dist::svd(ref_ta, trange, trange); + auto [S, U, VT] = non_dist::svd(ref_ta, trange, trange); std::vector exact_singular_values = exact_evals; std::sort(exact_singular_values.begin(), exact_singular_values.end(), diff --git a/tests/solvers.cpp b/tests/solvers.cpp index 2cd2480fab..af6af21847 100644 --- a/tests/solvers.cpp +++ b/tests/solvers.cpp @@ -23,7 +23,7 @@ * */ -#include +#include #include #include "unit_test_config.h" From c7236ea0f2494117b14529a548b8dcfddb182cf6 Mon Sep 17 00:00:00 2001 From: asadchev Date: Mon, 7 Dec 2020 14:19:50 -0500 Subject: [PATCH 4/7] Add missing LU implementation --- src/TiledArray/math/linalg/lu.h | 32 ++++++++++++++++++++++++++------ src/TiledArray/math/scalapack.h | 2 +- 2 files changed, 27 insertions(+), 7 deletions(-) diff --git a/src/TiledArray/math/linalg/lu.h b/src/TiledArray/math/linalg/lu.h index e35eaee99b..ac18641a09 100644 --- a/src/TiledArray/math/linalg/lu.h +++ b/src/TiledArray/math/linalg/lu.h @@ -18,7 +18,7 @@ * Eduard Valeyev * * lu.h - * Created: 16 October, 2020 + * Created: 16 October, 2020 * */ #ifndef TILEDARRAY_MATH_LINALG_LU_H__INCLUDED @@ -30,14 +30,34 @@ #endif #include -namespace TiledArray { +namespace TiledArray::math::linalg { + +template +auto lu_solve(const ArrayA& A, const ArrayB& B, + TiledRange x_trange = TiledRange()) { +#if TILEDARRAY_HAS_SCALAPACK + if (A.world().size() > 1 && A.range().volume() > 10000000) { + return scalapack::lu_solve(A, B, x_trange); + } +#endif + return non_distributed::lu_solve(A, B, x_trange); +} +template +auto lu_inv(const Array& A, TiledRange ainv_trange = TiledRange()) { #if TILEDARRAY_HAS_SCALAPACK -using scalapack::lu_inv; -using scalapack::lu_solve; -#else + if (A.world().size() > 1 && A.range().volume() > 10000000) { + return scalapack::lu_inv(A, ainv_trange); + } #endif + return non_distributed::lu_inv(A, ainv_trange); +} + +} // namespace TiledArray::math::linalg -} // namespace TiledArray +namespace TiledArray { + using TiledArray::math::linalg::lu_inv; + using TiledArray::math::linalg::lu_solve; +} #endif // TILEDARRAY_MATH_LINALG_LU_H__INCLUDED diff --git a/src/TiledArray/math/scalapack.h b/src/TiledArray/math/scalapack.h index 7a135edce6..9fb7a22709 100644 --- a/src/TiledArray/math/scalapack.h +++ b/src/TiledArray/math/scalapack.h @@ -26,7 +26,7 @@ #define TILEDARRAY_MATH_SCALAPACK_H__INCLUDED #warning \ - "TiledArray/math/scalapack.h header is deprecated, please include TiledArray/math/linalg/scalapack/all.h" + "TiledArray/math/scalapack.h header is deprecated, please use TiledArray/math/linalg.h" #include From f6de87604962d3f34bf541b34f15e2bd332c2de7 Mon Sep 17 00:00:00 2001 From: asadchev Date: Thu, 10 Dec 2020 12:09:28 -0500 Subject: [PATCH 5/7] Disable MPI C++ in scalapackpp --- CMakeLists.txt | 2 ++ external/scalapack.cmake | 4 ++++ 2 files changed, 6 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index c630921f00..6519eeec91 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -96,6 +96,8 @@ include(FeatureSummary) # Preload versions/tags of all dependencies ==================================== include(external/versions.cmake) +set(MPI_CXX_SKIP_MPICXX TRUE CACHE BOOL "MPI_CXX_SKIP_MPICXX") + # Configure options ======================================================= option(ENABLE_MPI "Enable MPI" ON) add_feature_info(MPI ENABLE_MPI "Message-Passing Interface supports distributed-memory parallel programs") diff --git a/external/scalapack.cmake b/external/scalapack.cmake index c279f81bcf..ee8d694b30 100644 --- a/external/scalapack.cmake +++ b/external/scalapack.cmake @@ -43,6 +43,10 @@ else() add_subdirectory( ${blacspp_SOURCE_DIR} ${blacspp_BINARY_DIR} ) add_subdirectory( ${scalapackpp_SOURCE_DIR} ${scalapackpp_BINARY_DIR} ) + # propagate MPI_CXX_SKIP_MPICXX=ON + target_compile_definitions( blacspp PRIVATE ${MPI_CXX_COMPILE_DEFINITIONS} ) + target_compile_definitions( scalapackpp PRIVATE ${MPI_CXX_COMPILE_DEFINITIONS} ) + install( TARGETS blacspp scalapackpp EXPORT tiledarray COMPONENT tiledarray ) # Add these dependencies to External add_dependencies(External-tiledarray scalapackpp blacspp) From 445a777afbc3e629d0592279fd49d95d5535bdb0 Mon Sep 17 00:00:00 2001 From: asadchev Date: Mon, 14 Dec 2020 02:19:45 -0500 Subject: [PATCH 6/7] Enable linalg test --- tests/CMakeLists.txt | 2 +- tests/{linear_algebra.cpp => linalg.cpp} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename tests/{linear_algebra.cpp => linalg.cpp} (100%) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 39eb0f1444..ed19572691 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -111,7 +111,7 @@ set(ta_test_src_files ta_test.cpp t_tot_tot_contract_.cpp tot_tot_tot_contract_.cpp einsum.cpp - #linear_algebra.cpp + linalg.cpp ) if(CUDA_FOUND) diff --git a/tests/linear_algebra.cpp b/tests/linalg.cpp similarity index 100% rename from tests/linear_algebra.cpp rename to tests/linalg.cpp From 11daf44acbf31ce65903f0856c1d0e26580ee03e Mon Sep 17 00:00:00 2001 From: asadchev Date: Fri, 27 Nov 2020 19:25:46 -0500 Subject: [PATCH 7/7] Enable scalapack in gitlab-ci.cmake --- .gitlab-ci.yml | 2 -- cmake/gitlab-ci.cmake | 5 ++++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index f157c16ff7..52691d237e 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -37,5 +37,3 @@ build: - CXX: [ clang++-9 ] IMAGE : [ "ubuntu:20.04" ] TA_PYTHON : [ "TA_PYTHON=OFF" ] - - diff --git a/cmake/gitlab-ci.cmake b/cmake/gitlab-ci.cmake index d8e479c207..e73ca8dec7 100644 --- a/cmake/gitlab-ci.cmake +++ b/cmake/gitlab-ci.cmake @@ -1,2 +1,5 @@ set(TA_BUILD_UNITTEST TRUE) - +set(ENABLE_SCALAPACK ON) +set(blacs_LIBRARIES "scalapack-openmpi") +set(scalapack_LIBRARIES "scalapack-openmpi") +set(lapack_LIBRARIES "lapack")