Asadchev/refactor/lapack

.gitlab-ci.yml

@@ -11,20 +11,27 @@ variables:
   OMPI_MCA_btl_vader_single_copy_mechanism: none
 
 before_script:
-  #- apt-get install -yq ...
-  - echo 'localhost slots=8' > /etc/openmpi/openmpi-default-hostfile 
+  - echo 'localhost slots=8' > /etc/openmpi/openmpi-default-hostfile
+  - apt-get install -yq intel-mkl-core-c-2020.4-304
+  - ln -s /opt/intel/compilers_and_libraries_2020.4.304/linux/mkl/ /opt/intel/
 
 build:
   stage: build
   timeout: 3h
   image: valeevgroup/${IMAGE}
   variables:
     CXX: ${CXX}
+    TA_PYTHON : "TA_PYTHON=OFF"
+    ENABLE_MKL : "ENABLE_MKL=OFF"
+    ENABLE_SCALAPACK : "ENABLE_SCALAPACK=OFF"
   script:
     - mkdir build
     - cd build
-    - cmake .. -DCMAKE_TOOLCHAIN_FILE=../cmake/gitlab-ci.cmake 
+    - cmake .. -DCMAKE_TOOLCHAIN_FILE=../cmake/gitlab-ci.cmake
+        -DCMAKE_BUILD_TYPE=${BUILD_TYPE}
         -D${TA_PYTHON}
+        -D${ENABLE_MKL}
+        -D${ENABLE_SCALAPACK}
     - make tiledarray
     - make examples
     - make ta_test
@@ -33,7 +40,14 @@ build:
     matrix:
       - CXX: [ g++ ]
         IMAGE : [ "ubuntu:20.04" ]
+        BUILD_TYPE : [ "Release", "Debug" ]
+        ENABLE_MKL : [ "ENABLE_MKL=ON", "ENABLE_MKL=OFF" ]
+        ENABLE_SCALAPACK : [ "ENABLE_SCALAPACK=ON", "ENABLE_SCALAPACK=OFF" ]
+      - CXX: [ g++ ]
+        IMAGE : [ "ubuntu:20.04" ]
+        BUILD_TYPE : [ "Release", "Debug" ]
         TA_PYTHON : [ "TA_PYTHON=ON" ]
       - CXX: [ clang++-9 ]
         IMAGE : [ "ubuntu:20.04" ]
-        TA_PYTHON : [ "TA_PYTHON=OFF" ]
+        BUILD_TYPE : [ "Release", "Debug" ]
+        ENABLE_MKL : [ "ENABLE_MKL=ON", "ENABLE_MKL=OFF" ]

cmake/gitlab-ci.cmake

@@ -1,5 +1,4 @@
 set(TA_BUILD_UNITTEST TRUE)
-set(ENABLE_SCALAPACK ON)
 set(blacs_LIBRARIES "scalapack-openmpi")
 set(scalapack_LIBRARIES "scalapack-openmpi")
 set(lapack_LIBRARIES "lapack")

external/eigen.cmake

@@ -151,11 +151,10 @@ if (TARGET TiledArray_Eigen)
   target_compile_definitions(TiledArray_Eigen INTERFACE ${LAPACK_COMPILE_DEFINITIONS})
   target_compile_options(TiledArray_Eigen INTERFACE ${LAPACK_COMPILE_OPTIONS})
   set(TiledArray_Eigen_VERSION "${Eigen3_VERSION}" CACHE STRING "Eigen3_VERSION of the library interfaced by TiledArray_Eigen target")
+  # Eigen's prototypes for BLAS interface libraries do not match MADNESS cblas
   if (MADNESS_HAS_MKL)
-    target_compile_definitions(TiledArray_Eigen INTERFACE EIGEN_USE_MKL_ALL)
+    # target_compile_definitions(TiledArray_Eigen INTERFACE EIGEN_USE_MKL EIGEN_USE_BLAS)
   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
     # target_compile_definitions(TiledArray_Eigen INTERFACE EIGEN_USE_BLAS)
   endif(MADNESS_HAS_MKL)
   install(TARGETS TiledArray_Eigen EXPORT tiledarray COMPONENT tiledarray)

src/CMakeLists.txt

@@ -123,6 +123,7 @@ TiledArray/math/partial_reduce.h
 TiledArray/math/transpose.h
 TiledArray/math/vector_op.h
 TiledArray/math/scalapack.h
+TiledArray/math/linalg/rank-local.h
 TiledArray/pmap/blocked_pmap.h
 TiledArray/pmap/cyclic_pmap.h
 TiledArray/pmap/hash_pmap.h
@@ -176,9 +177,6 @@ TiledArray/util/singleton.h
 TiledArray/util/time.h
 TiledArray/util/vector.h
 
-TiledArray/math/linalg/rank-local.h
-TiledArray/math/linalg/rank-local.cpp
-
 )
 
 if(CUDA_FOUND)
@@ -203,6 +201,7 @@ if(CUDA_FOUND)
 endif(CUDA_FOUND)
 
 set(TILEDARRAY_SOURCE_FILES
+TiledArray/tiledarray.cpp
 TiledArray/tensor/tensor.cpp
 TiledArray/sparse_shape.cpp
 TiledArray/tensor_impl.cpp

src/TiledArray/external/eigen.h

@@ -37,27 +37,6 @@
 TILEDARRAY_PRAGMA_GCC(diagnostic push)
 TILEDARRAY_PRAGMA_GCC(system_header)
 
-////////////////////////////////////////////////
-// this duplicates TiledArray_Eigen definitions
-#if HAVE_INTEL_MKL
-#ifndef EIGEN_USE_MKL  // strangely, defining EIGEN_USE_MKL_ALL does not imply
-                       // EIGEN_USE_MKL
-#define EIGEN_USE_MKL 1
-#endif
-#ifndef EIGEN_USE_MKL_ALL
-#define EIGEN_USE_MKL_ALL 1
-#endif
-#else
-//# ifndef EIGEN_USE_BLAS
-//#  define EIGEN_USE_BLAS 1
-//# endif
-#endif // HAVE_INTEL_MKL
-
-/////////////////////////////////////////////////
-// define lapacke types to prevent inclusion of complex.h by
-// Eigen/src/misc/lapacke.h
-#include <madness/tensor/lapacke_types.h>
-
 // If EIGEN_USE_LAPACKE_STRICT is defined, Eigen doesn't check if
 // EIGEN_USE_LAPACKE is defined before setting it, leading to a warning when it
 // is already set, so we unset here to avoid that warning
@@ -70,9 +49,9 @@ TILEDARRAY_PRAGMA_GCC(system_header)
 #if defined(EIGEN_USE_LAPACKE) || defined(EIGEN_USE_LAPACKE_STRICT)
 #if !EIGEN_VERSION_AT_LEAST(3, 3, 7)
 #error "Eigen3 < 3.3.7 with LAPACKE enabled may give wrong eigenvalue results"
-#error "Either turn off TILEDARRAY_EIGEN_USE_LAPACKE or use Eigen3 3.3.7"
-#endif
+#error "Either turn off EIGEN_USE_LAPACKE/EIGEN_USE_LAPACKE_STRICT or use Eigen3 3.3.7"
 #endif
+#endif // EIGEN_USE_LAPACKE || EIGEN_USE_LAPACKE_STRICT
 
 TILEDARRAY_PRAGMA_GCC(diagnostic pop)
 

src/TiledArray/initialize.h

@@ -6,68 +6,16 @@
 #define TILEDARRAY_INITIALIZE_H__INCLUDED
 
 #include <TiledArray/config.h>
-
 #include <TiledArray/external/madness.h>
-#ifdef TILEDARRAY_HAS_CUDA
-#include <TiledArray/external/cuda.h>
-#include <TiledArray/cuda/cublas.h>
-#include <cutt.h>
-#endif
-#ifdef HAVE_INTEL_MKL
-#include <mkl.h>
-#endif
 
 namespace TiledArray {
 
-#ifdef TILEDARRAY_HAS_CUDA
-/// initialize cuda environment
-inline void cuda_initialize() {
-  /// initialize cudaGlobal
-  cudaEnv::instance();
-  //
-  cuBLASHandlePool::handle();
-  // initialize cuTT
-  cuttInitialize();
-}
-
-/// finalize cuda environment
-inline void cuda_finalize() {
-  CudaSafeCall(cudaDeviceSynchronize());
-  cuttFinalize();
-  cublasDestroy(cuBLASHandlePool::handle());
-  delete &cuBLASHandlePool::handle();
-  cudaEnv::instance().reset(nullptr);
-}
-#endif
-
-namespace detail {
-inline bool& initialized_madworld_accessor() {
-  static bool flag = false;
-  return flag;
-}
-inline bool initialized_madworld() { return initialized_madworld_accessor(); }
-inline bool& initialized_accessor() {
-  static bool flag = false;
-  return flag;
-}
-inline bool& finalized_accessor() {
-  static bool flag = false;
-  return flag;
-}
-#ifdef HAVE_INTEL_MKL
-inline int& mklnumthreads_accessor() {
-  static int value = -1;
-  return value;
-}
-#endif
-}  // namespace detail
-
 /// @return true if TiledArray (and, necessarily, MADWorld runtime) is in an
 /// initialized state
-inline bool initialized() { return detail::initialized_accessor(); }
+bool initialized();
 
 /// @return true if TiledArray has been finalized at least once
-inline bool finalized() { return detail::finalized_accessor(); }
+bool finalized();
 
 /// @name TiledArray initialization.
 ///       These functions initialize TiledArray and (if needed) MADWorld
@@ -81,40 +29,11 @@ inline bool finalized() { return detail::finalized_accessor(); }
 
 /// @throw TiledArray::Exception if TiledArray initialized MADWorld and
 /// TiledArray::finalize() had been called
-inline World& initialize(int& argc, char**& argv,
-                         const SafeMPI::Intracomm& comm, bool quiet = true) {
-  if (detail::initialized_madworld() && finalized())
-    throw Exception(
-        "TiledArray finalized MADWorld already, cannot re-initialize MADWorld "
-        "again");
-  if (!initialized()) {
-    if (!madness::initialized())
-      detail::initialized_madworld_accessor() = true;
-    else {  // if MADWorld initialized, we must assume that comm is its default
-            // World.
-      if (madness::World::is_default(comm))
-        throw Exception(
-            "MADWorld initialized before TiledArray::initialize(argc, argv, "
-            "comm), but not initialized with comm");
-    }
-    auto& default_world = detail::initialized_madworld()
-                              ? madness::initialize(argc, argv, comm, quiet)
-                              : *madness::World::find_instance(comm);
-    TiledArray::set_default_world(default_world);
-#ifdef TILEDARRAY_HAS_CUDA
-    TiledArray::cuda_initialize();
-#endif
-#ifdef HAVE_INTEL_MKL
-    // record number of MKL threads and set to 1
-    detail::mklnumthreads_accessor() = mkl_get_max_threads();
-    mkl_set_num_threads(1);
-#endif
-    madness::print_meminfo_disable();
-    detail::initialized_accessor() = true;
-    return default_world;
-  } else
-    throw Exception("TiledArray already initialized");
-}
+World& initialize(
+  int& argc, char**& argv,
+  const SafeMPI::Intracomm& comm,
+  bool quiet = true
+);
 
 inline World& initialize(int& argc, char**& argv, bool quiet = true) {
   return TiledArray::initialize(argc, argv, SafeMPI::COMM_WORLD, quiet);
@@ -129,23 +48,7 @@ inline World& initialize(int& argc, char**& argv, const MPI_Comm& comm,
 
 /// Finalizes TiledArray (and MADWorld runtime, if it had not been initialized
 /// when TiledArray::initialize was called).
-inline void finalize() {
-#ifdef HAVE_INTEL_MKL
-  // reset number of MKL threads
-  mkl_set_num_threads(detail::mklnumthreads_accessor());
-#endif
-#ifdef TILEDARRAY_HAS_CUDA
-  TiledArray::cuda_finalize();
-#endif
-  TiledArray::get_default_world()
-      .gop.fence();  // TODO remove when madness::finalize() fences
-  if (detail::initialized_madworld()) {
-    madness::finalize();
-  }
-  TiledArray::reset_default_world();
-  detail::initialized_accessor() = false;
-  detail::finalized_accessor() = true;
-}
+void finalize();
 
 }  // namespace TiledArray
 

src/TiledArray/math/lapack.h

@@ -0,0 +1,43 @@
+#ifndef TILEDARRAY_MATH_LAPACK_H__INCLUDED
+#define TILEDARRAY_MATH_LAPACK_H__INCLUDED
+
+#include <TiledArray/config.h>
+
+#if defined(BTAS_HAS_INTEL_MKL)
+
+#include <mkl_lapack.h>
+#include <mkl_lapacke.h> // lapack_int
+
+#elif defined(BTAS_HAS_LAPACKE)
+
+    // see https://github.com/xianyi/OpenBLAS/issues/1992 why this is needed to prevent lapacke.h #define'ing I
+#   include <complex>
+#   ifndef lapack_complex_float
+#     define lapack_complex_float std::complex<float>
+#   else // lapack_complex_float
+      static_assert(sizeof(std::complex<float>)==sizeof(lapack_complex_float), "sizes of lapack_complex_float and std::complex<float> do not match");
+#   endif // lapack_complex_float
+#   ifndef lapack_complex_double
+#     define lapack_complex_double std::complex<double>
+#   else // lapack_complex_double
+      static_assert(sizeof(std::complex<double>)==sizeof(lapack_complex_double), "sizes of lapack_complex_double and std::complex<double> do not match");
+#   endif // lapack_complex_double
+
+#if defined(BTAS_LAPACKE_HEADER)
+#include BTAS_LAPACKE_HEADER
+#else
+#include <lapacke.h>
+#endif
+
+#elif defined(__APPLE__)
+
+#include <Accelerate/Accelerate.h>
+using lapack_int = __CLPK_integer;
+
+#else
+
+#error "Could not find Lapack/e"
+
+#endif
+
+#endif // TILEDARRAY_MATH_LAPACK_H__INCLUDED