From 7d581e9477eff1ca21dc4172c7465f965b2025cb Mon Sep 17 00:00:00 2001
From: Tom Benson <30674819+benson31@users.noreply.github.com>
Date: Fri, 17 Jan 2020 09:39:50 -0800
Subject: [PATCH] Fix lbann half issues (#93)

* tick up version to 1.3.3

* fix issues with CPU half

* rework the copy interface

* work on Copy

* add a warmup run to the gemm test

* Various updates to copy dispatch; no more ETI for Copy

* fix an issue with the cuda half type's assignment operators

* make gpu_half_type assignment operators into templates

* Update include/El/blas_like/level1/CopyLocal.hpp

Co-Authored-By: Tim Moon <moon13@llnl.gov>

* remove some debugging output

* add decls for BaseDistMatrix copy and copyasync

* make gpu_half_type streamable

* be a little more clever about casting to __half

* fix some things

* add overloads of sqrt and pow for half types

* add unary minus for gpu half type

* fix an issue where NVCC decided that rvalue references have value semantics

* add an overload of Log for gpu_half_type

* add exception-throwing bitwise operators for half. This is to appease Aluminum at compile-time and should never be encounted IRL.

* add a bunch of missing library symbols

* add a bunch of transendental functions, etc, for half types

* fix some missing symbols when compiling without half support

* Fixes to instantiate Read/Write with gpu half type

* Add overload for instantiate with half types

* add a few missing symbols

* add a write impl for gpu half matrices

* temp: dispatch gemv through Gemm for __half

* Fix the GEMV as GEMM call

* temporary error-throw for unhandled case

* remove the incy != 1 gemv case

* patch around an issue in the Half library

* expose AbstractMatrix interface to Print()

Co-authored-by: Tim Moon <moon13@llnl.gov>
---
 CMakeLists.txt                                |   3 +-
 cmake/modules/HydrogenETIGeneration.cmake     |  31 +
 include/El/blas_like/level1/Copy.hpp          | 842 ++----------------
 .../blas_like/level1/CopyAsyncDistMatrix.hpp  |  97 ++
 .../El/blas_like/level1/CopyAsyncLocal.hpp    | 134 +++
 .../El/blas_like/level1/CopyDistMatrix.hpp    | 183 ++++
 include/El/blas_like/level1/CopyFromRoot.hpp  |  50 ++
 include/El/blas_like/level1/CopyLocal.hpp     | 292 ++++++
 include/El/blas_like/level1/decl.hpp          |  63 +-
 include/El/blas_like/level1/impl.hpp          |   1 -
 include/El/core.hpp                           |  11 -
 include/El/core/Element/decl.hpp              |  88 +-
 include/El/core/Element/impl.hpp              |  13 +
 include/El/core/Matrix/impl_gpu.hpp           |   2 +-
 include/El/core/random/decl.hpp               |   5 +
 include/El/core/random/impl.hpp               |  28 +-
 include/El/io.hpp                             |   3 +
 include/hydrogen/blas/GPU_BLAS_impl.hpp       |  43 +
 include/hydrogen/utils/HalfPrecision.hpp      | 184 +++-
 .../hydrogen/utils/NumericTypeConversion.hpp  |  19 +
 src/blas_like/level1/CMakeLists.txt           |   1 +
 src/blas_like/level1/Copy.cpp                 | 113 +++
 src/blas_like/level1/HilbertSchmidt.cpp       |   8 +-
 src/blas_like/level2/Gemv.cpp                 |  13 +
 src/core/DistMatrix/AbstractDistMatrix.cpp    |  55 ++
 src/core/imports/blas/Dot.hpp                 |   8 +-
 src/core/imports/blas/Gemv.hpp                |  75 +-
 src/core/imports/cuda.cpp                     |   1 -
 src/hydrogen/blas/gpu/Axpy.cu                 |   3 +-
 src/hydrogen/blas/gpu/Copy.cu                 |   4 +-
 src/hydrogen/blas/gpu/Fill.cu                 |   2 +
 src/hydrogen/blas/gpu/Scale.cu                |   3 +-
 src/hydrogen/blas/gpu/Transpose.cu            |   1 +
 src/io/Print.cpp                              |   2 +
 src/io/Read.cpp                               |  14 +
 src/io/Write.cpp                              |  23 +
 src/matrices/deterministic/classical/Ones.cpp |  10 +-
 .../deterministic/classical/Zeros.cpp         |   4 +
 src/matrices/random/independent/Bernoulli.cpp |  18 +-
 src/matrices/random/independent/Gaussian.cpp  |   4 +
 tests/blas_like/Gemm.cpp                      |   8 +
 unit_test/CMakeLists.txt                      |  13 +-
 unit_test/copy_cpu_gpu_half_test.cpp          |  82 ++
 unit_test/copy_cpu_gpu_test.cpp               |  75 ++
 44 files changed, 1734 insertions(+), 898 deletions(-)
 create mode 100644 cmake/modules/HydrogenETIGeneration.cmake
 create mode 100644 include/El/blas_like/level1/CopyAsyncDistMatrix.hpp
 create mode 100644 include/El/blas_like/level1/CopyAsyncLocal.hpp
 create mode 100644 include/El/blas_like/level1/CopyDistMatrix.hpp
 create mode 100644 include/El/blas_like/level1/CopyFromRoot.hpp
 create mode 100644 include/El/blas_like/level1/CopyLocal.hpp
 create mode 100644 src/blas_like/level1/Copy.cpp
 create mode 100644 unit_test/copy_cpu_gpu_half_test.cpp
 create mode 100644 unit_test/copy_cpu_gpu_test.cpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9ab73e1ef2..5c6760d528 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -54,7 +54,7 @@ endif (__GIT_EXECUTABLE)
 # This must be set because version tags
 set(HYDROGEN_VERSION_MAJOR 1)
 set(HYDROGEN_VERSION_MINOR 3)
-set(HYDROGEN_VERSION_PATCH 2)
+set(HYDROGEN_VERSION_PATCH 3)
 set(HYDROGEN_VERSION_MAJOR_MINOR
     "${HYDROGEN_VERSION_MAJOR}.${HYDROGEN_VERSION_MINOR}")
 set(HYDROGEN_VERSION
@@ -539,6 +539,7 @@ if (HYDROGEN_HAVE_CUDA)
   target_link_libraries(Hydrogen_CUDA PUBLIC ${NVTX_LIBRARIES})
   target_link_libraries(Hydrogen_CUDA PUBLIC cuda::toolkit)
 
+  target_link_libraries(Hydrogen_CXX PUBLIC Hydrogen_CUDA)
   list(APPEND HYDROGEN_LIBRARIES Hydrogen_CUDA)
 endif ()
 
diff --git a/cmake/modules/HydrogenETIGeneration.cmake b/cmake/modules/HydrogenETIGeneration.cmake
new file mode 100644
index 0000000000..1ec1f99d6b
--- /dev/null
+++ b/cmake/modules/HydrogenETIGeneration.cmake
@@ -0,0 +1,31 @@
+# This function sets the _ARG# variable and recurs through the parameters.
+function (h_gtpe_recur OUTVAR EXPRESSION_TEMPLATE THIS_EPARAM_ID THIS_EPARAM)
+  foreach (_VAL IN LISTS ${THIS_EPARAM})
+    set(_ARG${THIS_EPARAM_ID} "${_VAL}")
+    if (ARGN)
+      math(EXPR _NEXT_ID "${THIS_EPARAM_ID} + 1")
+      h_gtpe_recur(${OUTVAR} "${EXPRESSION_TEMPLATE}" ${_NEXT_ID} ${ARGN})
+    else ()
+      string(CONFIGURE "${EXPRESSION_TEMPLATE}" _THIS_EXPRESSION @ONLY)
+      list(APPEND ${OUTVAR} "${_THIS_EXPRESSION}")
+    endif ()
+  endforeach ()
+  set(${OUTVAR} "${${OUTVAR}}" PARENT_SCOPE)
+endfunction ()
+
+# This function calls the recursive function above to generate each
+# function instance.
+#
+# DO NOT ADD A SEMICOLON TO THE EXPRESSION_TEMPLATE PARAMETER!!!!!!!
+function (h_generate_tensor_product_expression OUTVAR EXPRESSION_TEMPLATE)
+  h_gtpe_recur(_ALL_EXPRESSIONS "${EXPRESSION_TEMPLATE}" 0 ${ARGN})
+  set(${OUTVAR} ${_ALL_EXPRESSIONS} PARENT_SCOPE)
+endfunction ()
+
+# This function adds a semicolon to each function instance in the list
+# and joins them into a string with each function instance on its own
+# line.
+function (h_func_list_to_string OUTVAR INLIST)
+  list(JOIN ${INLIST} ";\n" _TMP)
+  set(${OUTVAR} "${_TMP};" PARENT_SCOPE)
+endfunction ()
diff --git a/include/El/blas_like/level1/Copy.hpp b/include/El/blas_like/level1/Copy.hpp
index 226977cdfa..7d835c4c4f 100644
--- a/include/El/blas_like/level1/Copy.hpp
+++ b/include/El/blas_like/level1/Copy.hpp
@@ -13,802 +13,154 @@
 #include <omp.h>
 #endif
 
+#include <El/core/Grid.hpp>
 #include <El/blas_like/level1/Copy/internal_decl.hpp>
 #include <El/blas_like/level1/Copy/GeneralPurpose.hpp>
 #include <El/blas_like/level1/Copy/util.hpp>
 
-#ifdef HYDROGEN_HAVE_GPU
-#include <hydrogen/blas/GPU_BLAS.hpp>
-#endif
-
-namespace El {
-namespace details {
+#include <hydrogen/meta/MetaUtilities.hpp>
 
-template <template <typename> class X, typename... Ts>
-using Expand = TypeList<X<Ts>...>;
-
-// This is replaced by a generic multiple dispatch engine in
-// DiHydrogen; this is a one-off use-case for now, so there's no need
-// to backport a robust implementation.
-template <typename FunctorT, typename LHSList, typename RHSList>
-struct CopyDispatcher
+// Introduce some metaprogramming notions.
+//
+// TODO: Move elsewhere.
+namespace El
 {
-    static void Do(FunctorT f,
-                   BaseDistMatrix const& src, BaseDistMatrix& tgt)
-    {
-        using LHead = Head<LHSList>;
-        using LTail = Tail<LHSList>;
-        if (auto const* ptr = dynamic_cast<LHead const*>(&src))
-            return CopyDispatcher<FunctorT, LHSList, RHSList>::DoRHS(
-                f, *ptr, tgt);
-        else
-            return CopyDispatcher<FunctorT, LTail, RHSList>::Do(f, src, tgt);
-    }
 
-    template <typename LHSType>
-    static void DoRHS(FunctorT f, LHSType const& src, BaseDistMatrix& tgt)
-    {
-        using RHead = Head<RHSList>;
-        using RTail = Tail<RHSList>;
-        if (auto* ptr = dynamic_cast<RHead*>(&tgt))
-            return f(src, *ptr);
-        else
-            return CopyDispatcher<FunctorT, LHSList, RTail>::DoRHS(f, src, tgt);
-    }
-};// struct CopyDispatcher
+template <bool B>
+using BoolVT = std::integral_constant<bool, B>;
 
-template <typename FunctorT, typename RHSList>
-struct CopyDispatcher<FunctorT, TypeList<>, RHSList>
+namespace details
 {
-    static void Do(FunctorT const&,
-                   BaseDistMatrix const&, BaseDistMatrix const&)
-    {
-        LogicError("Source matrix type not found.");
-    }
-};
 
-template <typename FunctorT, typename LHSList>
-struct CopyDispatcher<FunctorT, LHSList, TypeList<>>
+/** @brief A simple metafunction for interoping bitwise-equivalent
+ *         types across device interfaces.
+ */
+template <typename T, Device D>
+struct CompatibleStorageTypeT
 {
-    static void DoRHS(FunctorT const&,
-                      BaseDistMatrix const&, BaseDistMatrix const&)
-    {
-        LogicError("Target matrix type not found.");
-    }
+    using type = T;
 };
 
-struct CopyFunctor
-{
-    template <typename T, typename U>
-    void operator()(AbstractDistMatrix<T> const& src,
-                    AbstractDistMatrix<U>& tgt) const
-    {
-        return Copy(src, tgt);
-    }
-};// CopyFunctor
-
-struct CopyAsyncFunctor
-{
-    template <typename T, typename U>
-    void operator()(AbstractDistMatrix<T> const& src,
-                    AbstractDistMatrix<U>& tgt) const
-    {
-        return CopyAsync(src, tgt);
-    }
-};// CopyAsyncFunctor
-
-}// namespace details
-
-inline void Copy(BaseDistMatrix const& Source, BaseDistMatrix& Target)
-{
-    using FunctorT = details::CopyFunctor;
-    using MatrixTs = details::Expand<AbstractDistMatrix, float, double>;
-    using Dispatcher = details::CopyDispatcher<FunctorT, MatrixTs, MatrixTs>;
-    details::CopyFunctor f;
-    return Dispatcher::Do(f, Source, Target);
-}
-
-inline void CopyAsync(BaseDistMatrix const& Source, BaseDistMatrix& Target)
-{
-    using FunctorT = details::CopyAsyncFunctor;
-    using MatrixTs = details::Expand<AbstractDistMatrix, float, double>;
-    using Dispatcher = details::CopyDispatcher<FunctorT, MatrixTs, MatrixTs>;
-    details::CopyAsyncFunctor f;
-    return Dispatcher::Do(f, Source, Target);
-}
-
-template <typename T>
-void Copy(AbstractMatrix<T> const& A, AbstractMatrix<T>& B)
-{
-    switch (A.GetDevice())
-    {
-    case Device::CPU:
-        switch (B.GetDevice())
-        {
-        case Device::CPU:
-            Copy(static_cast<Matrix<T,Device::CPU> const&>(A),
-                 static_cast<Matrix<T,Device::CPU>&>(B));
-            break;
-#ifdef HYDROGEN_HAVE_CUDA
-        case Device::GPU:
-            Copy(static_cast<Matrix<T,Device::CPU> const&>(A),
-                 static_cast<Matrix<T,Device::GPU>&>(B));
-            break;
-#endif // HYDROGEN_HAVE_CUDA
-        default:
-            LogicError("Copy: Bad device.");
-        }
-        break;
-#ifdef HYDROGEN_HAVE_CUDA
-    case Device::GPU:
-        switch (B.GetDevice())
-        {
-        case Device::CPU:
-            Copy(static_cast<Matrix<T,Device::GPU> const&>(A),
-                 static_cast<Matrix<T,Device::CPU>&>(B));
-            break;
-        case Device::GPU:
-            Copy(static_cast<Matrix<T,Device::GPU> const&>(A),
-                 static_cast<Matrix<T,Device::GPU>&>(B));
-            break;
-        default:
-            LogicError("Copy: Bad device.");
-        }
-        break;
-#endif //  HYDROGEN_HAVE_CUDA
-    default:
-        LogicError("Copy: Bad device.");
-    }
-}
-
-template <typename T, typename U>
-void Copy(AbstractMatrix<T> const& A, AbstractMatrix<U>& B)
-{
-    switch (A.GetDevice())
-    {
-    case Device::CPU:
-        switch (B.GetDevice())
-        {
-        case Device::CPU:
-            Copy(static_cast<Matrix<T,Device::CPU> const&>(A),
-                 static_cast<Matrix<U,Device::CPU>&>(B));
-            break;
-#ifdef HYDROGEN_HAVE_CUDA
-        case Device::GPU:
-            Copy(static_cast<Matrix<T,Device::CPU> const&>(A),
-                 static_cast<Matrix<U,Device::GPU>&>(B));
-            break;
-#endif // HYDROGEN_HAVE_CUDA
-        default:
-            LogicError("Copy: Bad device.");
-        }
-        break;
-#ifdef HYDROGEN_HAVE_CUDA
-    case Device::GPU:
-        switch (B.GetDevice())
-        {
-        case Device::CPU:
-            Copy(static_cast<Matrix<T,Device::GPU> const&>(A),
-                 static_cast<Matrix<U,Device::CPU>&>(B));
-            break;
-        case Device::GPU:
-            Copy(static_cast<Matrix<T,Device::GPU> const&>(A),
-                 static_cast<Matrix<U,Device::GPU>&>(B));
-            break;
-        default:
-            LogicError("Copy: Bad device.");
-        }
-        break;
-#endif //  HYDROGEN_HAVE_CUDA
-    default:
-        LogicError("Copy: Bad device.");
-    }
-}
+template <typename T, Device D>
+using CompatibleStorageType = typename CompatibleStorageTypeT<T, D>::type;
 
-template<typename T>
-void Copy( const Matrix<T>& A, Matrix<T>& B )
-{
-    EL_DEBUG_CSE
-    const Int height = A.Height();
-    const Int width = A.Width();
-    const Int size = height * width;
-    B.Resize( height, width );
-    const Int ldA = A.LDim();
-    const Int ldB = B.LDim();
-    const T* EL_RESTRICT ABuf = A.LockedBuffer();
-          T* EL_RESTRICT BBuf = B.Buffer();
+#if defined(HYDROGEN_HAVE_HALF) && defined(HYDROGEN_GPU_USE_FP16)
 
-    if( ldA == height && ldB == height )
-    {
-#ifdef _OPENMP
-#if defined(HYDROGEN_HAVE_OMP_TASKLOOP)
-        const Int numThreads = omp_get_num_threads();
-        #pragma omp taskloop default(shared)
-        for(Int thread = 0; thread < numThreads; ++thread)
-        {
-#else
-        #pragma omp parallel
-        {
-            const Int numThreads = omp_get_num_threads();
-            const Int thread = omp_get_thread_num();
-#endif
-            const Int chunk = (size + numThreads - 1) / numThreads;
-            const Int start = Min(chunk * thread, size);
-            const Int end = Min(chunk * (thread + 1), size);
-            MemCopy( &BBuf[start], &ABuf[start], end - start );
-        }
-#else
-        MemCopy( BBuf, ABuf, size );
-#endif
-    }
-    else
-    {
-        EL_PARALLEL_FOR
-        for( Int j=0; j<width; ++j )
-        {
-            MemCopy(&BBuf[j*ldB], &ABuf[j*ldA], height);
-        }
-    }
-}
-
-#ifdef HYDROGEN_HAVE_CUDA
-template <typename T, typename U>
-void Copy(Matrix<T, Device::GPU> const& A, Matrix<U, Device::GPU>& B)
+template <>
+struct CompatibleStorageTypeT<cpu_half_type, El::Device::GPU>
 {
-    EL_DEBUG_CSE;
-    Int const height = A.Height();
-    Int const width = A.Width();
-    B.Resize(height, width);
-    Int const ldA = A.LDim();
-    Int const ldB = B.LDim();
-    T const* ABuf = A.LockedBuffer();
-    U* BBuf = B.Buffer();
-
-    SyncInfo<Device::GPU> syncInfoA = SyncInfoFromMatrix(A),
-        syncInfoB = SyncInfoFromMatrix(B);
-    auto syncHelper = MakeMultiSync(syncInfoB, syncInfoA);
-
-    gpu_blas::Copy(TransposeMode::NORMAL,
-                   height, width,
-                   ABuf, ldA,
-                   BBuf, ldB,
-                   syncInfoB);
-}
-
-template<typename T>
-void Copy(const Matrix<T,Device::GPU>& A, Matrix<T,Device::GPU>& B)
-{
-    EL_DEBUG_CSE
-    const Int height = A.Height();
-    const Int width = A.Width();
-    B.Resize(height, width);
-    const Int ldA = A.LDim();
-    const Int ldB = B.LDim();
-    const T* ABuf = A.LockedBuffer();
-    T* BBuf = B.Buffer();
-
-    SyncInfo<Device::GPU> syncInfoA = SyncInfoFromMatrix(A), syncInfoB = SyncInfoFromMatrix(B);
-    auto syncHelper = MakeMultiSync(syncInfoB, syncInfoA);
-
-    // Launch the copy
-    H_CHECK_CUDA(
-        cudaMemcpy2DAsync(BBuf, ldB*sizeof(T),
-                          ABuf, ldA*sizeof(T),
-                          height*sizeof(T), width,
-                          cudaMemcpyDeviceToDevice,
-                          syncInfoB.stream_));
-}
-
-// These inter-device copy functions are SYNCHRONOUS with respect to
-// the host.
-template <typename T>
-void Copy(Matrix<T,Device::CPU> const& A, Matrix<T,Device::GPU>& B)
-{
-    EL_DEBUG_CSE
-    const Int height = A.Height();
-    const Int width = A.Width();
-    B.Resize(height, width);
-    const Int ldA = A.LDim();
-    const Int ldB = B.LDim();
-    const T* EL_RESTRICT ABuf = A.LockedBuffer();
-    T* EL_RESTRICT BBuf = B.Buffer();
-
-    SyncInfo<Device::GPU> syncInfoB = SyncInfoFromMatrix(B);
-    InterDeviceCopy<Device::CPU,Device::GPU>::MemCopy2DAsync(
-        BBuf, ldB, ABuf, ldA, height, width, syncInfoB.stream_);
-    Synchronize(syncInfoB); // Is this necessary??
-}
+    using type = gpu_half_type;
+};
 
-template <typename T>
-void Copy(Matrix<T,Device::GPU> const& A, Matrix<T,Device::CPU>& B)
+template <>
+struct CompatibleStorageTypeT<gpu_half_type, El::Device::CPU>
 {
-    EL_DEBUG_CSE
-    const Int height = A.Height();
-    const Int width = A.Width();
-    B.Resize(height, width);
-    const Int ldA = A.LDim();
-    const Int ldB = B.LDim();
-    const T* EL_RESTRICT ABuf = A.LockedBuffer();
-    T* EL_RESTRICT BBuf = B.Buffer();
+    using type = cpu_half_type;
+};
 
-    SyncInfo<Device::GPU> syncInfoA = SyncInfoFromMatrix(A);
-    InterDeviceCopy<Device::GPU,Device::CPU>::MemCopy2DAsync(
-        BBuf, ldB, ABuf, ldA, height, width, syncInfoA.stream_);
-    Synchronize(syncInfoA); // Is this necessary??
-}
+#endif // defined(HYDROGEN_HAVE_HALF) && defined(HYDROGEN_GPU_USE_FP16)
 
-// These inter-device copy functions are ASYNCHRONOUS with respect to
-// the host.
 template <typename T>
-void CopyAsync(Matrix<T,Device::CPU> const& A, Matrix<T,Device::GPU>& B)
-{
-    EL_DEBUG_CSE
-    const Int height = A.Height();
-    const Int width = A.Width();
-    B.Resize(height, width);
-    const Int ldA = A.LDim();
-    const Int ldB = B.LDim();
-    const T* EL_RESTRICT ABuf = A.LockedBuffer();
-    T* EL_RESTRICT BBuf = B.Buffer();
-
-    InterDeviceCopy<Device::CPU, Device::GPU>::MemCopy2DAsync(
-        BBuf, ldB, ABuf, ldA, height, width, B.Stream());
-}
+using CPUStorageType = CompatibleStorageType<T, Device::CPU>;
 
+#ifdef HYDROGEN_HAVE_GPU
 template <typename T>
-void CopyAsync(Matrix<T,Device::GPU> const& A, Matrix<T,Device::CPU>& B)
-{
-    EL_DEBUG_CSE
-    const Int height = A.Height();
-    const Int width = A.Width();
-    B.Resize(height, width);
-    const Int ldA = A.LDim();
-    const Int ldB = B.LDim();
-    const T* EL_RESTRICT ABuf = A.LockedBuffer();
-    T* EL_RESTRICT BBuf = B.Buffer();
-
-    InterDeviceCopy<Device::GPU, Device::CPU>::MemCopy2DAsync(
-        BBuf, ldB, ABuf, ldA, height, width, A.Stream());
-}
-
-#endif // HYDROGEN_HAVE_CUDA
-
-template<typename S,typename T,
-         typename/*=EnableIf<CanCast<S,T>>*/>
-void Copy(const Matrix<S>& A, Matrix<T>& B)
-{
-    EL_DEBUG_CSE
-    EntrywiseMap(A, B, MakeFunction(Caster<S,T>::Cast));
-}
-
-template<typename T,Dist U,Dist V,Device D,
-         typename = EnableIf<IsDeviceValidType<T,D>>>
-void Copy(const ElementalMatrix<T>& A,
-          DistMatrix<T,U,V,ELEMENT,D>& B)
-{
-    EL_DEBUG_CSE
-    B = A;
-}
-
-template<typename T,Dist U,Dist V,Device D,
-         typename = DisableIf<IsDeviceValidType<T,D>>,
-         typename = void>
-void Copy(const ElementalMatrix<T>& A,
-          DistMatrix<T,U,V,ELEMENT,D>& B)
-{
-    EL_DEBUG_CSE
-    LogicError("Copy: bad data/device combination.");
-}
-
-// Datatype conversions should not be very common, and so it is likely best to
-// avoid explicitly instantiating every combination
-template<typename S,typename T,Dist U,Dist V,Device D,
-         typename = EnableIf<IsDeviceValidType<T,D>>>
-void Copy( const ElementalMatrix<S>& A, DistMatrix<T,U,V,ELEMENT,D>& B )
-{
-    EL_DEBUG_CSE
-    if (A.Grid() == B.Grid() && A.ColDist() == U && A.RowDist() == V
-        && A.GetLocalDevice() == D)
-    {
-        if( !B.RootConstrained() )
-            B.SetRoot( A.Root() );
-        if( !B.ColConstrained() )
-            B.AlignCols( A.ColAlign() );
-        if( !B.RowConstrained() )
-            B.AlignRows( A.RowAlign() );
-        if( A.Root() == B.Root() &&
-            A.ColAlign() == B.ColAlign() && A.RowAlign() == B.RowAlign() )
-        {
-            B.Resize( A.Height(), A.Width() );
-            Copy( A.LockedMatrix(), B.Matrix() );
-            return;
-        }
-    }
-    DistMatrix<S,U,V,ELEMENT,D> BOrig(A.Grid());
-    BOrig.AlignWith( B );
-    BOrig = A;
-    B.Resize( A.Height(), A.Width() );
-    Copy( BOrig.LockedMatrix(), B.Matrix() );
-}
-
-template<typename S,typename T,Dist U,Dist V,Device D,
-         typename=DisableIf<IsDeviceValidType<T,D>>,
-         typename=void>
-void Copy( const ElementalMatrix<S>& A, DistMatrix<T,U,V,ELEMENT,D>& B )
-{
-    EL_DEBUG_CSE
-    LogicError("Copy: bad data/device combination.");
-}
-
-template<typename T,Dist U,Dist V>
-void Copy( const BlockMatrix<T>& A, DistMatrix<T,U,V,BLOCK>& B )
-{
-    EL_DEBUG_CSE
-    B = A;
-}
-
-// Datatype conversions should not be very common, and so it is likely best to
-// avoid explicitly instantiating every combination
-template<typename S,typename T,Dist U,Dist V>
-void Copy( const BlockMatrix<S>& A, DistMatrix<T,U,V,BLOCK>& B )
-{
-    EL_DEBUG_CSE
-    if( A.Grid() == B.Grid() && A.ColDist() == U && A.RowDist() == V )
-    {
-        if( !B.RootConstrained() )
-            B.SetRoot( A.Root() );
-        if( !B.ColConstrained() )
-            B.AlignColsWith( A.DistData() );
-        if( !B.RowConstrained() )
-            B.AlignRowsWith( A.DistData() );
-        if( A.Root() == B.Root() &&
-            A.ColAlign() == B.ColAlign() &&
-            A.RowAlign() == B.RowAlign() &&
-            A.ColCut() == B.ColCut() &&
-            A.RowCut() == B.RowCut() )
-        {
-            B.Resize( A.Height(), A.Width() );
-            Copy( A.LockedMatrix(), B.Matrix() );
-            return;
-        }
-    }
-    DistMatrix<S,U,V,BLOCK> BOrig(A.Grid());
-    BOrig.AlignWith( B );
-    BOrig = A;
-    B.Resize( A.Height(), A.Width() );
-    Copy( BOrig.LockedMatrix(), B.Matrix() );
-}
-
-template<typename S,typename T,
-         typename/*=EnableIf<CanCast<S,T>>*/>
-void Copy( const ElementalMatrix<S>& A, ElementalMatrix<T>& B )
-{
-    EL_DEBUG_CSE
-#define GUARD(CDIST,RDIST,WRAP,DEVICE)                                        \
-        (B.ColDist() == CDIST) && (B.RowDist() == RDIST)                \
-            && (B.Wrap() == WRAP) && (B.GetLocalDevice() == DEVICE)
-#define PAYLOAD(CDIST,RDIST,WRAP,DEVICE)                                \
-        auto& BCast =                                                   \
-            static_cast<DistMatrix<T,CDIST,RDIST,ELEMENT,DEVICE>&>(B);  \
-        Copy( A, BCast );
-    #include <El/macros/DeviceGuardAndPayload.h>
-}
-
-template<typename T>
-void Copy( const AbstractDistMatrix<T>& A, AbstractDistMatrix<T>& B )
-{
-    EL_DEBUG_CSE
-    const DistWrap wrapA=A.Wrap(), wrapB=B.Wrap();
-    if( wrapA == ELEMENT && wrapB == ELEMENT )
-    {
-        auto& ACast = static_cast<const ElementalMatrix<T>&>(A);
-        auto& BCast = static_cast<ElementalMatrix<T>&>(B);
-        Copy( ACast, BCast );
-    }
-    else if( wrapA == BLOCK && wrapB == BLOCK )
-    {
-        auto& ACast = static_cast<const BlockMatrix<T>&>(A);
-        auto& BCast = static_cast<BlockMatrix<T>&>(B);
-        Copy( ACast, BCast );
-    }
-    else
-    {
-        copy::GeneralPurpose( A, B );
-    }
-}
+using GPUStorageType = CompatibleStorageType<T, Device::GPU>;
+#endif
 
-template <typename T, Dist U, Dist V, Device D1, Device D2>
-void CopyAsync(DistMatrix<T,U,V,ELEMENT,D1> const& A,
-               DistMatrix<T,U,V,ELEMENT,D2>& B)
+// This layer of indirection checks the Tgt types and launches the
+// copy if possible.
+template <typename CopyFunctor,
+          typename T, typename U, Device D1, Device D2,
+          EnableWhen<IsStorageType<T, D1>, int> = 0>
+void LaunchCopy(Matrix<T, D1> const& src, Matrix<U, D2>& tgt,
+                CopyFunctor const& F)
 {
-    EL_DEBUG_CSE
-#ifndef EL_RELEASE
-    auto const Adata = A.DistData(), Bdata = B.DistData();
-    if (!((Adata.blockHeight == Bdata.blockHeight) &&
-          (Adata.blockWidth == Bdata.blockWidth) &&
-          (Adata.colAlign == Bdata.colAlign) &&
-          (Adata.rowAlign == Bdata.rowAlign) &&
-          (Adata.colCut == Bdata.colCut) &&
-          (Adata.rowCut == Bdata.rowCut) &&
-          (Adata.root == Bdata.root) &&
-          (Adata.grid == Bdata.grid)))
-    {
-        LogicError("CopyAsync: "
-                   "A and B must have the same DistData, except device.");
-    }
-#endif // !defined(EL_RELEASE)
-    B.Resize(A.Height(), A.Width());
-    CopyAsync(A.LockedMatrix(), B.Matrix());
+   return F(src, tgt);
 }
 
-template <typename T, Dist U, Dist V, Device D>
-void CopyAsync(DistMatrix<T,U,V,ELEMENT,D> const& A,
-               DistMatrix<T,U,V,ELEMENT,D>& B)
+template <typename CopyFunctor,
+          typename T, typename U, Device D1, Device D2,
+          EnableUnless<IsStorageType<T, D1>, int> = 0>
+void LaunchCopy(Matrix<T, D1> const&, Matrix<U, D2>&,
+                CopyFunctor const&)
 {
-    LogicError("CopyAsync: Both matrices on same device (D=",
-               DeviceName<D>(), ").");
+    LogicError("The combination U=", TypeTraits<U>::Name(), " "
+               "and D=", DeviceName<D2>(), " is not supported.");
 }
 
-template <typename T, Dist U, Dist V, Device D>
-void CopyAsync(ElementalMatrix<T> const& A, DistMatrix<T,U,V,ELEMENT,D>& B)
+// This layer of indirection checks the Src types; this overload is
+// also useful for some DistMatrix instantiations.
+template <typename CopyFunctor,
+          typename T, typename U, Device D2,
+          EnableWhen<IsStorageType<U, D2>, int> = 0>
+void LaunchCopy(AbstractMatrix<T> const& src, Matrix<U, D2>& tgt,
+                CopyFunctor const& F)
 {
-    EL_DEBUG_CSE
-    if ((A.ColDist() == U) && (A.RowDist() == V))
+    switch (src.GetDevice())
     {
-        switch (A.GetLocalDevice())
-        {
-        case Device::CPU:
-            CopyAsync(
-                static_cast<DistMatrix<T,U,V,ELEMENT,Device::CPU> const&>(A),
-                B);
-            break;
-#ifdef HYDROGEN_HAVE_CUDA
-        case Device::GPU:
-            CopyAsync(
-                static_cast<DistMatrix<T,U,V,ELEMENT,Device::GPU> const&>(A),
-                B);
-            break;
-#endif // HYDROGEN_HAVE_CUDA
-        default:
-            LogicError("CopyAsync: Unknown device type.");
-        }
+    case Device::CPU:
+        return LaunchCopy(
+            static_cast<Matrix<T, Device::CPU> const&>(src), tgt, F);
+#ifdef HYDROGEN_HAVE_GPU
+    case Device::GPU:
+        return LaunchCopy(
+            static_cast<Matrix<T, Device::GPU> const&>(src), tgt, F);
+#endif // HYDROGEN_HAVE_GPU
+    default:
+        LogicError("Copy: Bad device.");
     }
-    else
-        LogicError("CopyAsync requires A and B to have the same distribution.");
 }
 
-template <typename T>
-void CopyAsync(ElementalMatrix<T> const& A, ElementalMatrix<T>& B)
+template <typename CopyFunctor,
+          typename T, typename U, Device D2,
+          EnableUnless<IsStorageType<U, D2>, int> = 0>
+void LaunchCopy(AbstractMatrix<T> const&, Matrix<U, D2>&,
+                CopyFunctor const&)
 {
-    EL_DEBUG_CSE
-#define GUARD(CDIST,RDIST,WRAP,DEVICE)                              \
-    (B.ColDist() == CDIST) && (B.RowDist() == RDIST)                \
-        && (B.Wrap() == WRAP) && (B.GetLocalDevice() == DEVICE)
-#define PAYLOAD(CDIST,RDIST,WRAP,DEVICE)                            \
-    auto& BCast =                                                   \
-        static_cast<DistMatrix<T,CDIST,RDIST,ELEMENT,DEVICE>&>(B);  \
-    CopyAsync(A, BCast);
-    #include <El/macros/DeviceGuardAndPayload.h>
+    LogicError("The combination U=", TypeTraits<U>::Name(), " "
+               "and D=", DeviceName<D2>(), " is not supported.");
 }
 
+// The variadic templates allow these functors to be recycled across
+// sequential and distributed matrices.
 
-template <typename T>
-void CopyAsync(AbstractDistMatrix<T> const& A, AbstractDistMatrix<T>& B)
+struct CopyFunctor
 {
-    EL_DEBUG_CSE
-    const DistWrap wrapA = A.Wrap(), wrapB = B.Wrap();
-    if (wrapA == ELEMENT && wrapB == ELEMENT)
+    template <typename... Args>
+    void operator()(Args&&... args) const
     {
-        auto& ACast = static_cast<const ElementalMatrix<T>&>(A);
-        auto& BCast = static_cast<ElementalMatrix<T>&>(B);
-        CopyAsync(ACast, BCast);
+        return Copy(std::forward<Args>(args)...);
     }
-    else
-        LogicError("CopyAsync only implemented for ElementalMatrix.");
-}
+};// CopyFunctor
 
-template<typename S,typename T,
-         typename/*=EnableIf<CanCast<S,T>>*/>
-void Copy( const AbstractDistMatrix<S>& A, AbstractDistMatrix<T>& B )
+struct CopyAsyncFunctor
 {
-    EL_DEBUG_CSE
-    const DistWrap wrapA=A.Wrap(), wrapB=B.Wrap();
-    if( wrapA == ELEMENT && wrapB == ELEMENT )
-    {
-        auto& ACast = static_cast<const ElementalMatrix<S>&>(A);
-        auto& BCast = static_cast<ElementalMatrix<T>&>(B);
-        Copy( ACast, BCast );
-    }
-    else if( wrapA == BLOCK && wrapB == BLOCK )
+    template <typename... Args>
+    void operator()(Args&&... args) const
     {
-        auto& ACast = static_cast<const BlockMatrix<S>&>(A);
-        auto& BCast = static_cast<BlockMatrix<T>&>(B);
-        Copy( ACast, BCast );
+        return CopyAsync(std::forward<Args>(args)...);
     }
-    else
-    {
-        LogicError("If you see this error, please tell Tom.");
-        copy::GeneralPurpose( A, B );
-    }
-}
-
-template<typename S,typename T,
-         typename/*=EnableIf<CanCast<S,T>>*/>
-void Copy( const BlockMatrix<S>& A, BlockMatrix<T>& B )
-{
-    EL_DEBUG_CSE
-    #define GUARD(CDIST,RDIST,WRAP) \
-        B.ColDist() == CDIST && B.RowDist() == RDIST && B.Wrap() == WRAP \
-            && B.GetLocalDevice() == Device::CPU
-    #define PAYLOAD(CDIST,RDIST,WRAP) \
-      auto& BCast = static_cast<DistMatrix<T,CDIST,RDIST,BLOCK>&>(B); \
-      Copy( A, BCast );
-    #include <El/macros/GuardAndPayload.h>
-}
-
-template<typename T>
-void CopyFromRoot
-( const Matrix<T>& A, DistMatrix<T,CIRC,CIRC>& B, bool includingViewers )
-{
-    EL_DEBUG_CSE
-    if( B.CrossRank() != B.Root() )
-        LogicError("Called CopyFromRoot from non-root");
-    B.Resize( A.Height(), A.Width() );
-    B.MakeSizeConsistent( includingViewers );
-    B.Matrix() = A;
-}
+};// CopyAsyncFunctor
 
-template<typename T>
-void CopyFromNonRoot( DistMatrix<T,CIRC,CIRC>& B, bool includingViewers )
-{
-    EL_DEBUG_CSE
-    if( B.CrossRank() == B.Root() )
-        LogicError("Called CopyFromNonRoot from root");
-    B.MakeSizeConsistent( includingViewers );
-}
+}// namespace details
+}// namespace El
 
-template<typename T>
-void CopyFromRoot
-( const Matrix<T>& A, DistMatrix<T,CIRC,CIRC,BLOCK>& B,
-  bool includingViewers )
-{
-    EL_DEBUG_CSE
-    if( B.CrossRank() != B.Root() )
-        LogicError("Called CopyFromRoot from non-root");
-    B.Resize( A.Height(), A.Width() );
-    B.MakeSizeConsistent( includingViewers );
-    B.Matrix() = A;
-}
+//
+// Include all the definitions
+//
+#include "CopyLocal.hpp"
+#include "CopyAsyncLocal.hpp"
+#include "CopyDistMatrix.hpp"
+#include "CopyAsyncDistMatrix.hpp"
+#include "CopyFromRoot.hpp"
 
-template<typename T>
-void CopyFromNonRoot
-( DistMatrix<T,CIRC,CIRC,BLOCK>& B, bool includingViewers )
+namespace El
 {
-    EL_DEBUG_CSE
-    if( B.CrossRank() == B.Root() )
-        LogicError("Called CopyFromNonRoot from root");
-    B.MakeSizeConsistent( includingViewers );
-}
-
-#ifdef EL_INSTANTIATE_BLAS_LEVEL1
-# define EL_EXTERN
-#else
-# define EL_EXTERN extern
-#endif
-
-#define PROTO(T)                                                        \
-    EL_EXTERN template void Copy(                                       \
-        AbstractMatrix<T> const& A,                                     \
-        AbstractMatrix<T>& B);                                          \
-    EL_EXTERN template void Copy(                                       \
-        Matrix<T> const& A,                                             \
-        Matrix<T>& B);                                                  \
-    EL_EXTERN template void Copy(                                       \
-        AbstractDistMatrix<T> const& A,                                 \
-        AbstractDistMatrix<T>& B);                                      \
-    EL_EXTERN template void CopyFromRoot(                               \
-        Matrix<T> const& A,                                             \
-        DistMatrix<T,CIRC,CIRC>& B,                                     \
-        bool includingViewers);                                         \
-    EL_EXTERN template void CopyFromNonRoot(                            \
-        DistMatrix<T,CIRC,CIRC>& B,                                     \
-        bool includingViewers);                                         \
-    EL_EXTERN template void CopyFromRoot(                               \
-        Matrix<T> const& A,                                             \
-        DistMatrix<T,CIRC,CIRC,BLOCK>& B,                               \
-        bool includingViewers);                                         \
-    EL_EXTERN template void CopyFromNonRoot(                            \
-        DistMatrix<T,CIRC,CIRC,BLOCK>& B,                               \
-        bool includingViewers);                                         \
-    EL_EXTERN template void CopyAsync(                                  \
-        AbstractDistMatrix<T> const& A,                                 \
-        AbstractDistMatrix<T>& B);
-
-
-#ifdef HYDROGEN_HAVE_CUDA
-#ifdef HYDROGEN_GPU_USE_FP16
-EL_EXTERN template void Copy(
-    const AbstractMatrix<gpu_half_type>& A,
-    AbstractMatrix<gpu_half_type>& B);
-EL_EXTERN template void Copy(
-    const AbstractDistMatrix<gpu_half_type>& A,
-    AbstractDistMatrix<gpu_half_type>& B);
-EL_EXTERN template void Copy(
-    const Matrix<gpu_half_type,Device::GPU>& A,
-    Matrix<gpu_half_type,Device::GPU>& B);
-EL_EXTERN template void Copy(
-    const Matrix<gpu_half_type,Device::GPU>& A,
-    Matrix<gpu_half_type,Device::CPU>& B);
-EL_EXTERN template void Copy(
-    const Matrix<gpu_half_type,Device::CPU>& A,
-    Matrix<gpu_half_type,Device::GPU>& B);
-EL_EXTERN template void CopyAsync(
-    const Matrix<gpu_half_type,Device::GPU>& A,
-    Matrix<gpu_half_type,Device::CPU>& B);
-EL_EXTERN template void CopyAsync(
-    const Matrix<gpu_half_type,Device::CPU>& A,
-    Matrix<gpu_half_type,Device::GPU>& B);
-#endif
-EL_EXTERN template void Copy(
-    Matrix<float,Device::GPU> const& A,
-    Matrix<float,Device::GPU>& B);
-EL_EXTERN template void Copy(
-    Matrix<float,Device::GPU> const& A,
-    Matrix<float,Device::CPU>& B);
-EL_EXTERN template void Copy(
-    Matrix<float,Device::CPU> const& A,
-    Matrix<float,Device::GPU>& B);
-EL_EXTERN template void CopyAsync(
-    Matrix<float,Device::GPU> const& A,
-    Matrix<float,Device::CPU>& B);
-EL_EXTERN template void CopyAsync(
-    Matrix<float,Device::CPU> const& A,
-    Matrix<float,Device::GPU>& B);
-EL_EXTERN template void Copy(
-    Matrix<double,Device::GPU> const& A,
-    Matrix<double,Device::GPU>& B);
-EL_EXTERN template void Copy(
-    Matrix<double,Device::GPU> const& A,
-    Matrix<double,Device::CPU>& B);
-EL_EXTERN template void Copy(
-    Matrix<double,Device::CPU> const& A,
-    Matrix<double,Device::GPU>& B);
-EL_EXTERN template void CopyAsync(
-    Matrix<double,Device::GPU> const& A,
-    Matrix<double,Device::CPU>& B );
-EL_EXTERN template void CopyAsync(
-    Matrix<double,Device::CPU> const& A,
-    Matrix<double,Device::GPU>& B );
-#endif // HYDROGEN_HAVE_CUDA
-
-#ifdef HYDROGEN_HAVE_HALF
-EL_EXTERN template void Copy(
-    AbstractMatrix<cpu_half_type> const& A,
-    AbstractMatrix<cpu_half_type>& B );
-EL_EXTERN template void Copy(
-    AbstractDistMatrix<cpu_half_type> const& A,
-    AbstractDistMatrix<cpu_half_type>& B );
-EL_EXTERN template void Copy(
-    Matrix<cpu_half_type> const& A,
-    Matrix<cpu_half_type>& B );
-#endif // HYDROGEN_HAVE_HALF
-
-EL_EXTERN template void Copy(
-    AbstractMatrix<uint8_t> const& A,
-    AbstractMatrix<uint8_t>& B );
-EL_EXTERN template void Copy(
-    Matrix<uint8_t> const& A,
-    Matrix<uint8_t>& B );
-
-#define EL_ENABLE_DOUBLEDOUBLE
-#define EL_ENABLE_QUADDOUBLE
-#define EL_ENABLE_QUAD
-#define EL_ENABLE_BIGINT
-#define EL_ENABLE_BIGFLOAT
-#include <El/macros/Instantiate.h>
 
-#undef EL_EXTERN
+void Copy(BaseDistMatrix const&, BaseDistMatrix&);
+void CopyAsync(BaseDistMatrix const&, BaseDistMatrix&);
 
-} // namespace El
+}// namespace El
 
 #endif // ifndef EL_BLAS_COPY_HPP
diff --git a/include/El/blas_like/level1/CopyAsyncDistMatrix.hpp b/include/El/blas_like/level1/CopyAsyncDistMatrix.hpp
new file mode 100644
index 0000000000..bcc23dec7b
--- /dev/null
+++ b/include/El/blas_like/level1/CopyAsyncDistMatrix.hpp
@@ -0,0 +1,97 @@
+#ifndef EL_BLAS_LIKE_LEVEL1_COPYASYNCDISTMATRIX_HPP_
+#define EL_BLAS_LIKE_LEVEL1_COPYASYNCDISTMATRIX_HPP_
+
+#include <hydrogen/Device.hpp>
+
+#include "CopyLocal.hpp"
+#include "CopyDistMatrix.hpp"
+#include "CopyAsyncLocal.hpp"
+
+namespace El
+{
+
+template <typename S, typename T, Dist U, Dist V, Device D1, Device D2>
+void CopyAsync(DistMatrix<S,U,V,ELEMENT,D1> const& A,
+               DistMatrix<T,U,V,ELEMENT,D2>& B)
+{
+    EL_DEBUG_CSE;
+#ifndef EL_RELEASE
+    auto const Adata = A.DistData(), Bdata = B.DistData();
+    if (!((Adata.blockHeight == Bdata.blockHeight) &&
+          (Adata.blockWidth == Bdata.blockWidth) &&
+          (Adata.colAlign == Bdata.colAlign) &&
+          (Adata.rowAlign == Bdata.rowAlign) &&
+          (Adata.colCut == Bdata.colCut) &&
+          (Adata.rowCut == Bdata.rowCut) &&
+          (Adata.root == Bdata.root) &&
+          (Adata.grid == Bdata.grid)))
+    {
+        LogicError("CopyAsync: "
+                   "A and B must have the same DistData, except device.");
+    }
+#endif // !defined(EL_RELEASE)
+    B.Resize(A.Height(), A.Width());
+    CopyAsync(A.LockedMatrix(), B.Matrix());
+}
+
+template <typename S, typename T, Dist U, Dist V, Device D>
+void CopyAsync(ElementalMatrix<S> const& A, DistMatrix<T,U,V,ELEMENT,D>& B)
+{
+    EL_DEBUG_CSE;
+    if ((A.ColDist() == U) && (A.RowDist() == V))
+    {
+        switch (A.GetLocalDevice())
+        {
+        case Device::CPU:
+            CopyAsync(
+                static_cast<DistMatrix<S,U,V,ELEMENT,Device::CPU> const&>(A),
+                B);
+            break;
+#ifdef HYDROGEN_HAVE_CUDA
+        case Device::GPU:
+            CopyAsync(
+                static_cast<DistMatrix<S,U,V,ELEMENT,Device::GPU> const&>(A),
+                B);
+            break;
+#endif // HYDROGEN_HAVE_CUDA
+        default:
+            LogicError("CopyAsync: Unknown device type.");
+        }
+    }
+    else
+        LogicError("CopyAsync requires A and B to have the same distribution.");
+}
+
+
+template <typename T, typename U>
+void CopyAsync(ElementalMatrix<T> const& A, ElementalMatrix<U>& B)
+{
+    EL_DEBUG_CSE;
+#define GUARD(CDIST,RDIST,WRAP,DEVICE)                              \
+    (B.ColDist() == CDIST) && (B.RowDist() == RDIST)                \
+        && (B.Wrap() == WRAP) && (B.GetLocalDevice() == DEVICE)
+#define PAYLOAD(CDIST,RDIST,WRAP,DEVICE)                            \
+    auto& BCast =                                                   \
+        static_cast<DistMatrix<U,CDIST,RDIST,ELEMENT,DEVICE>&>(B);  \
+    CopyAsync(A, BCast);
+    #include <El/macros/DeviceGuardAndPayload.h>
+}
+
+template <typename T, typename U>
+void CopyAsync(AbstractDistMatrix<T> const& A, AbstractDistMatrix<U>& B)
+{
+    EL_DEBUG_CSE;
+    if (A.Wrap() == ELEMENT && B.Wrap() == ELEMENT)
+    {
+        auto& ACast = static_cast<const ElementalMatrix<T>&>(A);
+        auto& BCast = static_cast<ElementalMatrix<U>&>(B);
+        CopyAsync(ACast, BCast);
+    }
+    else
+    {
+        LogicError("CopyAsync is only supported for ElementalMatrices.");
+    }
+}
+
+}// namespace El
+#endif // EL_BLAS_LIKE_LEVEL1_COPYASYNCDISTMATRIX_HPP_
diff --git a/include/El/blas_like/level1/CopyAsyncLocal.hpp b/include/El/blas_like/level1/CopyAsyncLocal.hpp
new file mode 100644
index 0000000000..10e3add816
--- /dev/null
+++ b/include/El/blas_like/level1/CopyAsyncLocal.hpp
@@ -0,0 +1,134 @@
+#ifndef EL_BLAS_LIKE_LEVEL1_COPYASYNCLOCAL_HPP_
+#define EL_BLAS_LIKE_LEVEL1_COPYASYNCLOCAL_HPP_
+
+#include <hydrogen/Device.hpp>
+
+#include "CopyLocal.hpp"
+
+namespace El
+{
+
+template <typename T, typename U,
+          EnableWhen<And<IsStorageType<T, Device::CPU>,
+                         IsStorageType<U, Device::CPU>>, int> = 0>
+void CopyAsyncImpl(Matrix<T, Device::CPU> const& src,
+                   Matrix<U, Device::CPU>& tgt)
+{
+    // No optimized asynchronous copy on CPU. Perhaps at some point we
+    // will introduce a truly asynchronous Copy by running it on an
+    // internal thread. That is not the case yet.
+    return Copy(src, tgt);
+}
+
+#ifdef HYDROGEN_HAVE_GPU
+
+template <typename T, typename U,
+          EnableWhen<And<IsStorageType<T, Device::GPU>,
+                         IsStorageType<U, Device::GPU>>, int> = 0>
+void CopyAsyncImpl(Matrix<T, Device::GPU> const& src,
+                   Matrix<U, Device::GPU>& tgt)
+{
+    // The (intertype) copy on GPU is already asynchronous with
+    // respect to the host. Perhaps at some point we will introduce
+    // "emulated asynchrony" by allowing some operations to run on an
+    // internal stream; not now.
+    return Copy(src, tgt);
+}
+
+// These inter-device copy functions are as ASYNCHRONOUS with respect
+// to the host as possible. See the CUDA documentation of Memcpy
+// synchronization semantics for more information.
+template <typename T,
+          EnableWhen<And<IsStorageType<T, Device::CPU>,
+                         IsStorageType<T, Device::GPU>>, int> = 0>
+void CopyAsyncImpl(Matrix<T, Device::CPU> const& A,
+                   Matrix<T, Device::GPU>& B)
+{
+    EL_DEBUG_CSE;
+    const Int height = A.Height();
+    const Int width = A.Width();
+    B.Resize(height, width);
+    const Int ldA = A.LDim();
+    const Int ldB = B.LDim();
+    const T* EL_RESTRICT ABuf = A.LockedBuffer();
+    T* EL_RESTRICT BBuf = B.Buffer();
+
+    InterDeviceCopy<Device::CPU, Device::GPU>::MemCopy2DAsync(
+        BBuf, ldB, ABuf, ldA, height, width, B.Stream());
+}
+
+template <typename T,
+          EnableWhen<And<IsStorageType<T, Device::GPU>,
+                         IsStorageType<T, Device::CPU>>, int> = 0>
+void CopyAsyncImpl(Matrix<T, Device::GPU> const& A,
+                   Matrix<T, Device::CPU>& B)
+{
+    EL_DEBUG_CSE;
+    const Int height = A.Height();
+    const Int width = A.Width();
+    B.Resize(height, width);
+    const Int ldA = A.LDim();
+    const Int ldB = B.LDim();
+    const T* EL_RESTRICT ABuf = A.LockedBuffer();
+    T* EL_RESTRICT BBuf = B.Buffer();
+
+    InterDeviceCopy<Device::GPU, Device::CPU>::MemCopy2DAsync(
+        BBuf, ldB, ABuf, ldA, height, width, A.Stream());
+}
+#endif // HYDROGEN_HAVE_GPU
+
+// FIXME: If the source matrix is on the GPU, we should prefer
+// changing type before changing device. If the source matrix is on
+// the CPU, we should prefer changing the device before changing the
+// type. This ensures the most work is done asynchronously with
+// respect to the CPU.
+
+// Inter-device and inter-type
+template <typename T, typename U, Device D1, Device D2,
+          EnableWhen<And<BoolVT<D1!=D2>,
+                         IsStorageType<T, D1>,
+                         IsStorageType<U, D2>>, int> = 0>
+void CopyAsyncImpl(Matrix<T, D1> const& src, Matrix<U, D2>& tgt)
+{
+    using T_on_D2 = details::CompatibleStorageType<T, D2>;
+    Matrix<T_on_D2, D2> tmp;
+    CopyAsync(src, tmp); // Change device
+    CopyAsync(tmp, tgt); // Change type
+}
+
+template <typename T, typename U, Device D1, Device D2,
+          EnableWhen<Or<Not<IsStorageType<T,D1>>,
+                        Not<IsStorageType<U,D2>>>, int> = 0>
+void CopyAsyncImpl(Matrix<T, D1> const&, Matrix<U, D2>&)
+{
+    LogicError("Cannot dispatch CopyAsync.");
+}
+
+template <typename T, typename U, Device D1, Device D2>
+void CopyAsync(Matrix<T, D1> const& src, Matrix<U, D2>& tgt)
+{
+    CopyAsyncImpl(src, tgt);
+}
+
+template <typename T, typename U>
+void CopyAsync(AbstractMatrix<T> const& Source, AbstractMatrix<U>& Target)
+{
+    switch (Target.GetDevice())
+    {
+    case Device::CPU:
+        return details::LaunchCopy(
+            Source, static_cast<Matrix<U, Device::CPU>&>(Target),
+            details::CopyAsyncFunctor{});
+#ifdef HYDROGEN_HAVE_GPU
+    case Device::GPU:
+        return details::LaunchCopy(
+            Source, static_cast<Matrix<U, Device::GPU>&>(Target),
+            details::CopyAsyncFunctor{});
+#endif // HYDROGEN_HAVE_GPU
+    default:
+        LogicError("Copy: Bad device.");
+    }
+}
+
+}// namespace El
+#endif // EL_BLAS_LIKE_LEVEL1_COPYASYNCLOCAL_HPP_
diff --git a/include/El/blas_like/level1/CopyDistMatrix.hpp b/include/El/blas_like/level1/CopyDistMatrix.hpp
new file mode 100644
index 0000000000..cf7ca26b3e
--- /dev/null
+++ b/include/El/blas_like/level1/CopyDistMatrix.hpp
@@ -0,0 +1,183 @@
+#ifndef EL_BLAS_LIKE_LEVEL1_COPYDISTMATRIX_HPP_
+#define EL_BLAS_LIKE_LEVEL1_COPYDISTMATRIX_HPP_
+
+namespace El
+{
+
+template <typename T, Dist U, Dist V, Device D,
+          EnableWhen<IsStorageType<T,D>, int> = 0>
+void Copy(ElementalMatrix<T> const& A,
+          DistMatrix<T,U,V,ELEMENT,D>& B)
+{
+    EL_DEBUG_CSE;
+    B = A;
+}
+
+
+template <typename T, Dist U, Dist V, Device D,
+          EnableUnless<IsStorageType<T,D>, int> = 0>
+void Copy(ElementalMatrix<T> const&,
+          DistMatrix<T,U,V,ELEMENT,D>&)
+{
+    EL_DEBUG_CSE;
+    LogicError("Copy: bad data/device combination.");
+}
+
+// Datatype conversions should not be very common, and so it is likely best to
+// avoid explicitly instantiating every combination
+template <typename S,typename T, Dist U, Dist V, Device D,
+          EnableWhen<And<IsStorageType<S,D>,
+                         IsStorageType<T,D>>, int> = 0>
+void Copy(ElementalMatrix<S> const& A,
+          DistMatrix<T,U,V,ELEMENT,D>& B)
+{
+    EL_DEBUG_CSE;
+    if (A.Grid() == B.Grid() && A.ColDist() == U && A.RowDist() == V
+        && A.GetLocalDevice() == D)
+    {
+        if (!B.RootConstrained())
+            B.SetRoot(A.Root());
+        if (!B.ColConstrained())
+            B.AlignCols(A.ColAlign());
+        if (!B.RowConstrained())
+            B.AlignRows(A.RowAlign());
+        if (A.Root() == B.Root() &&
+            A.ColAlign() == B.ColAlign() && A.RowAlign() == B.RowAlign())
+        {
+            B.Resize(A.Height(), A.Width());
+            Copy(A.LockedMatrix(), B.Matrix());
+            return;
+        }
+    }
+    DistMatrix<S,U,V,ELEMENT,D> BOrig(A.Grid());
+    BOrig.AlignWith(B);
+    BOrig = A;
+    B.Resize(A.Height(), A.Width());
+    Copy(BOrig.LockedMatrix(), B.Matrix());
+}
+
+template <typename S,typename T, Dist U,Dist V,Device D,
+          EnableUnless<And<IsStorageType<S,D>,
+                         IsStorageType<T,D>>, int> = 0>
+void Copy(const ElementalMatrix<S>&, DistMatrix<T,U,V,ELEMENT,D>&)
+{
+    EL_DEBUG_CSE;
+    LogicError("Copy: bad data/device combination.");
+}
+
+template <typename T,Dist U,Dist V>
+void Copy(const BlockMatrix<T>& A, DistMatrix<T,U,V,BLOCK>& B)
+{
+    EL_DEBUG_CSE;
+    B = A;
+}
+
+// Datatype conversions should not be very common, and so it is likely best to
+// avoid explicitly instantiating every combination
+template <typename S, typename T, Dist U, Dist V>
+void Copy(const BlockMatrix<S>& A, DistMatrix<T,U,V,BLOCK>& B)
+{
+    EL_DEBUG_CSE;
+    if (A.Grid() == B.Grid() && A.ColDist() == U && A.RowDist() == V)
+    {
+        if (!B.RootConstrained())
+            B.SetRoot(A.Root());
+        if (!B.ColConstrained())
+            B.AlignColsWith(A.DistData());
+        if (!B.RowConstrained())
+            B.AlignRowsWith(A.DistData());
+        if (A.Root() == B.Root() &&
+            A.ColAlign() == B.ColAlign() &&
+            A.RowAlign() == B.RowAlign() &&
+            A.ColCut() == B.ColCut() &&
+            A.RowCut() == B.RowCut())
+        {
+            B.Resize(A.Height(), A.Width());
+            Copy(A.LockedMatrix(), B.Matrix());
+            return;
+        }
+    }
+    DistMatrix<S,U,V,BLOCK> BOrig(A.Grid());
+    BOrig.AlignWith(B);
+    BOrig = A;
+    B.Resize(A.Height(), A.Width());
+    Copy(BOrig.LockedMatrix(), B.Matrix());
+}
+
+template <typename T, typename U>
+void Copy(const ElementalMatrix<T>& A, ElementalMatrix<U>& B)
+{
+    EL_DEBUG_CSE;
+#define GUARD(CDIST,RDIST,WRAP,DEVICE)                                  \
+        (B.ColDist() == CDIST) && (B.RowDist() == RDIST)                \
+            && (B.Wrap() == WRAP) && (B.GetLocalDevice() == DEVICE)
+#define PAYLOAD(CDIST,RDIST,WRAP,DEVICE)                                \
+        auto& BCast =                                                   \
+            static_cast<DistMatrix<U,CDIST,RDIST,ELEMENT,DEVICE>&>(B);  \
+        Copy(A, BCast);
+    #include <El/macros/DeviceGuardAndPayload.h>
+}
+
+template<typename S,typename T,
+         typename/*=EnableIf<CanCast<S,T>>*/>
+void Copy(const BlockMatrix<S>& A, BlockMatrix<T>& B)
+{
+    EL_DEBUG_CSE;
+    #define GUARD(CDIST,RDIST,WRAP) \
+        B.ColDist() == CDIST && B.RowDist() == RDIST && B.Wrap() == WRAP \
+            && B.GetLocalDevice() == Device::CPU
+    #define PAYLOAD(CDIST,RDIST,WRAP) \
+      auto& BCast = static_cast<DistMatrix<T,CDIST,RDIST,BLOCK>&>(B); \
+      Copy(A, BCast);
+    #include <El/macros/GuardAndPayload.h>
+}
+
+template <typename T>
+void Copy(AbstractDistMatrix<T> const& A, AbstractDistMatrix<T>& B)
+{
+    EL_DEBUG_CSE;
+    DistWrap const wrapA=A.Wrap(), wrapB=B.Wrap();
+    if (wrapB == ELEMENT)
+    {
+        auto& ACast = static_cast<ElementalMatrix<T> const&>(A);
+        auto& BCast = static_cast<ElementalMatrix<T>&>(B);
+        Copy(ACast, BCast);
+    }
+    else if (wrapA == BLOCK && wrapB == BLOCK)
+    {
+        auto& ACast = static_cast<BlockMatrix<T> const&>(A);
+        auto& BCast = static_cast<BlockMatrix<T>&>(B);
+        Copy(ACast, BCast);
+    }
+    else
+    {
+        LogicError("If you see this error, please tell Tom.");
+        copy::GeneralPurpose(A, B);
+    }
+}
+
+template <typename T, typename U>
+void Copy(AbstractDistMatrix<T> const& A, AbstractDistMatrix<U>& B)
+{
+    EL_DEBUG_CSE;
+    DistWrap const wrapA=A.Wrap(), wrapB=B.Wrap();
+    if (wrapB == ELEMENT)
+    {
+        auto& ACast = static_cast<ElementalMatrix<T> const&>(A);
+        auto& BCast = static_cast<ElementalMatrix<U>&>(B);
+        Copy(ACast, BCast);
+    }
+    else if (wrapA == BLOCK && wrapB == BLOCK)
+    {
+        auto& ACast = static_cast<BlockMatrix<T> const&>(A);
+        auto& BCast = static_cast<BlockMatrix<U>&>(B);
+        Copy(ACast, BCast);
+    }
+    else
+    {
+        LogicError("If you see this error, please tell Tom.");
+    }
+}
+
+}// namespace El
+#endif // EL_BLAS_LIKE_LEVEL1_COPYDISTMATRIX_HPP_
diff --git a/include/El/blas_like/level1/CopyFromRoot.hpp b/include/El/blas_like/level1/CopyFromRoot.hpp
new file mode 100644
index 0000000000..6faddc3c98
--- /dev/null
+++ b/include/El/blas_like/level1/CopyFromRoot.hpp
@@ -0,0 +1,50 @@
+#ifndef EL_BLAS_LIKE_LEVEL1_COPYFROMROOT_HPP_
+#define EL_BLAS_LIKE_LEVEL1_COPYFROMROOT_HPP_
+
+namespace El
+{
+
+template<typename T>
+void CopyFromRoot(const Matrix<T>& A, DistMatrix<T,CIRC,CIRC>& B,
+                  bool includingViewers)
+{
+    EL_DEBUG_CSE;
+    if (B.CrossRank() != B.Root())
+        LogicError("Called CopyFromRoot from non-root");
+    B.Resize(A.Height(), A.Width());
+    B.MakeSizeConsistent(includingViewers);
+    B.Matrix() = A;
+}
+
+template<typename T>
+void CopyFromNonRoot(DistMatrix<T,CIRC,CIRC>& B, bool includingViewers)
+{
+    EL_DEBUG_CSE;
+    if (B.CrossRank() == B.Root())
+        LogicError("Called CopyFromNonRoot from root");
+    B.MakeSizeConsistent(includingViewers);
+}
+
+template<typename T>
+void CopyFromRoot (const Matrix<T>& A, DistMatrix<T,CIRC,CIRC,BLOCK>& B,
+                   bool includingViewers)
+{
+    EL_DEBUG_CSE;
+    if (B.CrossRank() != B.Root())
+        LogicError("Called CopyFromRoot from non-root");
+    B.Resize(A.Height(), A.Width());
+    B.MakeSizeConsistent(includingViewers);
+    B.Matrix() = A;
+}
+
+template<typename T>
+void CopyFromNonRoot (DistMatrix<T,CIRC,CIRC,BLOCK>& B, bool includingViewers)
+{
+    EL_DEBUG_CSE;
+    if (B.CrossRank() == B.Root())
+        LogicError("Called CopyFromNonRoot from root");
+    B.MakeSizeConsistent(includingViewers);
+}
+
+}// namespace El
+#endif // EL_BLAS_LIKE_LEVEL1_COPYFROMROOT_HPP_
diff --git a/include/El/blas_like/level1/CopyLocal.hpp b/include/El/blas_like/level1/CopyLocal.hpp
new file mode 100644
index 0000000000..e10d38e598
--- /dev/null
+++ b/include/El/blas_like/level1/CopyLocal.hpp
@@ -0,0 +1,292 @@
+#ifndef EL_BLAS_LIKE_LEVEL1_COPYLOCAL_HPP_
+#define EL_BLAS_LIKE_LEVEL1_COPYLOCAL_HPP_
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+#include "El/blas_like/level1/EntrywiseMap.hpp"
+
+#include <hydrogen/Device.hpp>
+#ifdef HYDROGEN_HAVE_GPU
+#include <hydrogen/blas/GPU_BLAS.hpp>
+#endif
+
+namespace El
+{
+
+// There are the following cases:
+//
+// 0. The source or target has a type that is not self-compatible.
+//    -> LogicError
+//
+// 1. The source and target have the same type on the same device.
+//    -> Simple copy operation
+//
+// 2. The source and target have the same type but on different devices.
+//    -> Simple interdevice copy operation (reduces to device-vendor API call)
+//
+// 3. The source and target have different types but the same binary
+//    layout on different devices.
+//    -> Reinterpret cast and simple interdevice copy operation.
+//
+// 4. The source and target have different types but on the same device.
+//    -> Casting copy.
+//
+// 5. The source and target have different types on different devices.
+//    -> Unclear if we can implement this directly by hand; may have
+//       to use a temporary to convert to the appropriate device/type
+//       and then to the appropriate type/device, resp.
+
+// LEGACY: Handle the cast-less case on CPU.
+// (Case 1, CPU)
+template <typename T,
+          EnableWhen<IsStorageType<T, Device::CPU>, int> = 0>
+void CopyImpl(Matrix<T, Device::CPU> const& A, Matrix<T, Device::CPU>& B)
+{
+    EL_DEBUG_CSE;
+    const Int height = A.Height();
+    const Int width = A.Width();
+    const Int size = height * width;
+    B.Resize(height, width);
+    const Int ldA = A.LDim();
+    const Int ldB = B.LDim();
+    const T* EL_RESTRICT ABuf = A.LockedBuffer();
+          T* EL_RESTRICT BBuf = B.Buffer();
+
+    if (ldA == height && ldB == height)
+    {
+// #ifdef _OPENMP
+// #if defined(HYDROGEN_HAVE_OMP_TASKLOOP)
+//         const Int numThreads = omp_get_num_threads();
+//         #pragma omp taskloop default(shared)
+//         for (Int thread = 0; thread < numThreads; ++thread)
+//         {
+// #else
+//         #pragma omp parallel
+//         {
+//             const Int numThreads = omp_get_num_threads();
+//             const Int thread = omp_get_thread_num();
+// #endif // defined(HYDROGEN_HAVE_OMP_TASKLOOP)
+//             const Int chunk = (size + numThreads - 1) / numThreads;
+//             const Int start = Min(chunk * thread, size);
+//             const Int end = Min(chunk * (thread + 1), size);
+//             MemCopy(&BBuf[start], &ABuf[start], end - start);
+//         }
+// #else
+        MemCopy(BBuf, ABuf, size);
+//#endif // _OPENMP
+    }
+    else
+    {
+        EL_PARALLEL_FOR
+        for (Int j=0; j<width; ++j)
+        {
+            MemCopy(&BBuf[j*ldB], &ABuf[j*ldA], height);
+        }
+    }
+}
+//}// DELETE ME
+
+// LEGACY: Generic implementation of casting on the CPU
+// (Case 4, CPU)
+template<typename S, typename T,
+         EnableWhen<And<CanCast<S,T>,
+                        IsStorageType<S, Device::CPU>,
+                        IsStorageType<T, Device::CPU>>, int> = 0>
+void CopyImpl(Matrix<S, Device::CPU> const& A, Matrix<T, Device::CPU>& B)
+{
+    EL_DEBUG_CSE;
+    EntrywiseMap(A, B, MakeFunction(Caster<S,T>::Cast));
+}
+
+#ifdef HYDROGEN_HAVE_GPU
+// Inter-type copy on the GPU. The gpu_blas can handle this via a
+// custom kernel.
+// (Case 4, GPU)
+template <typename T, typename U,
+          EnableWhen<And<IsStorageType<T, Device::GPU>,
+                         IsStorageType<U, Device::GPU>>, int> = 0>
+void CopyImpl(Matrix<T, Device::GPU> const& A, Matrix<U, Device::GPU>& B)
+{
+    EL_DEBUG_CSE;
+    Int const height = A.Height();
+    Int const width = A.Width();
+    B.Resize(height, width);
+    Int const ldA = A.LDim();
+    Int const ldB = B.LDim();
+    T const* ABuf = A.LockedBuffer();
+    U* BBuf = B.Buffer();
+
+    SyncInfo<Device::GPU> syncInfoA = SyncInfoFromMatrix(A),
+        syncInfoB = SyncInfoFromMatrix(B);
+    auto syncHelper = MakeMultiSync(syncInfoB, syncInfoA);
+
+    gpu_blas::Copy(TransposeMode::NORMAL,
+                   height, width,
+                   ABuf, ldA,
+                   BBuf, ldB,
+                   syncInfoB);
+}
+
+#ifdef HYDROGEN_HAVE_CUDA
+// If using CUDA, prefer the cudaMemcpy2D implementation. This is
+// ASYNCHRONOUS with respect to the host.
+// (Case 1, GPU)
+//
+// TODO: Profile to verify this is, indeed, faster.
+template <typename T,
+          EnableWhen<IsStorageType<T, Device::GPU>, int> = 0>
+void CopyImpl(Matrix<T, Device::GPU> const& A, Matrix<T, Device::GPU>& B)
+{
+    EL_DEBUG_CSE;
+    const Int height = A.Height();
+    const Int width = A.Width();
+    B.Resize(height, width);
+    const Int ldA = A.LDim();
+    const Int ldB = B.LDim();
+    const T* ABuf = A.LockedBuffer();
+    T* BBuf = B.Buffer();
+
+    SyncInfo<Device::GPU> syncInfoA = SyncInfoFromMatrix(A),
+        syncInfoB = SyncInfoFromMatrix(B);
+    auto syncHelper = MakeMultiSync(syncInfoB, syncInfoA);
+
+    // Launch the copy
+    H_CHECK_CUDA(
+        cudaMemcpy2DAsync(BBuf, ldB*sizeof(T),
+                          ABuf, ldA*sizeof(T),
+                          height*sizeof(T), width,
+                          cudaMemcpyDeviceToDevice,
+                          syncInfoB.stream_));
+}
+#endif // HYDROGEN_HAVE_CUDA
+
+namespace details
+{
+struct InterdeviceSync
+{
+    InterdeviceSync(SyncInfo<Device::CPU> const&,
+                    SyncInfo<Device::GPU> const& GPU_sync)
+        : gpu_sync_{GPU_sync}
+    {}
+
+    InterdeviceSync(SyncInfo<Device::GPU> const& GPU_sync,
+                    SyncInfo<Device::CPU> const&)
+        : gpu_sync_{GPU_sync}
+    {}
+
+    SyncInfo<Device::GPU> gpu_sync_;
+};
+}
+
+// These inter-device copy functions are SYNCHRONOUS with respect to
+// the host.
+template <typename T, Device D1, Device D2,
+          EnableWhen<And<BoolVT<D1 != D2>,
+                         IsStorageType<T, D1>,
+                         IsStorageType<T, D2>>, int> = 0>
+void CopyImpl(Matrix<T, D1> const& A, Matrix<T, D2>& B)
+{
+    EL_DEBUG_CSE;
+    const Int height = A.Height();
+    const Int width = A.Width();
+    B.Resize(height, width);
+    const Int ldA = A.LDim();
+    const Int ldB = B.LDim();
+    const T* EL_RESTRICT ABuf = A.LockedBuffer();
+    T* EL_RESTRICT BBuf = B.Buffer();
+
+    details::InterdeviceSync isync(SyncInfoFromMatrix(A),
+                                   SyncInfoFromMatrix(B));
+
+    InterDeviceCopy<D1, D2>::MemCopy2DAsync(
+        BBuf, ldB, ABuf, ldA, height, width, isync.gpu_sync_.stream_);
+    Synchronize(isync.gpu_sync_); // Is this necessary??
+}
+
+// These inter-device copy functions are SYNCHRONOUS with respect to
+// the host.
+template <typename T, Device D1, Device D2,
+          EnableWhen<And<BoolVT<D1 != D2>,
+                         Not<IsSame<T, details::CompatibleStorageType<T, D2>>>,
+                         IsStorageType<T, D1>>, int> = 0>
+void CopyImpl(Matrix<T, D1> const& A,
+              Matrix<details::CompatibleStorageType<T, D2>, D2>& B)
+{
+    EL_DEBUG_CSE;
+    const Int height = A.Height();
+    const Int width = A.Width();
+    B.Resize(height, width);
+    const Int ldA = A.LDim();
+    const Int ldB = B.LDim();
+    const T* EL_RESTRICT ABuf = A.LockedBuffer();
+    T* EL_RESTRICT BBuf = reinterpret_cast<T*>(B.Buffer());
+
+    details::InterdeviceSync isync(SyncInfoFromMatrix(A),
+                                   SyncInfoFromMatrix(B));
+    InterDeviceCopy<D1, D2>::MemCopy2DAsync(
+        BBuf, ldB, ABuf, ldA, height, width, isync.gpu_sync_.stream_);
+    Synchronize(isync.gpu_sync_); // Is this necessary??
+}
+
+#endif // HYDROGEN_HAVE_GPU
+
+// Inter-device and inter-type Copy.
+// (Case 5)
+template <typename T, typename U, Device D1, Device D2,
+          EnableWhen<And<BoolVT<D1!=D2>,
+                         CanCast<T,U>,
+                         IsStorageType<T, D1>,
+                         IsStorageType<U, D2>>, int> = 0,
+          EnableUnless<IsSame<U, details::CompatibleStorageType<T,D2>>, int> = 0>
+void CopyImpl(Matrix<T, D1> const& src, Matrix<U, D2>& tgt)
+{
+    // FIXME: Note that there might be instances in which it could
+    // feasibly be better to change type and then change device. In
+    // any case, none of this is ideal; aiming for correctness now.
+    using T_on_D2 = details::CompatibleStorageType<T, D2>;
+    Matrix<T_on_D2, D2> tmp;
+    CopyImpl(src, tmp); // Change device
+    CopyImpl(tmp, tgt); // Change type
+}
+
+// (Case 0)
+template <typename T, typename U, Device D1, Device D2,
+          EnableUnless<And<CanCast<T,U>,
+                           IsStorageType<T, D1>,
+                           IsStorageType<U, D2>>, int> = 0>
+void CopyImpl(Matrix<T, D1> const&, Matrix<U, D2>&)
+{
+    LogicError("Cannot dispatch Copy.");
+}
+
+// "Plain-ol' function"
+template <typename T, typename U, Device D1, Device D2>
+void Copy(Matrix<T, D1> const& src, Matrix<U, D2>& tgt)
+{
+    CopyImpl(src, tgt);
+}
+
+template <typename T, typename U>
+void Copy(AbstractMatrix<T> const& Source, AbstractMatrix<U>& Target)
+{
+    switch (Target.GetDevice())
+    {
+    case Device::CPU:
+        return details::LaunchCopy(
+            Source, static_cast<Matrix<U, Device::CPU>&>(Target),
+            details::CopyFunctor{});
+#ifdef HYDROGEN_HAVE_GPU
+    case Device::GPU:
+        return details::LaunchCopy(
+            Source, static_cast<Matrix<U, Device::GPU>&>(Target),
+            details::CopyFunctor{});
+#endif // HYDROGEN_HAVE_GPU
+    default:
+        LogicError("Copy: Bad device.");
+    }
+}
+
+}// namespace El
+#endif // EL_BLAS_LIKE_LEVEL1_COPYLOCAL_HPP_
diff --git a/include/El/blas_like/level1/decl.hpp b/include/El/blas_like/level1/decl.hpp
index 46f5230d7d..e41a5de239 100644
--- a/include/El/blas_like/level1/decl.hpp
+++ b/include/El/blas_like/level1/decl.hpp
@@ -366,61 +366,7 @@ void Contract( const BlockMatrix<T>& A, BlockMatrix<T>& B );
 
 // Copy
 // ====
-
-template<typename T>
-void Copy( const AbstractMatrix<T>& A, AbstractMatrix<T>& B );
-template<typename T>
-void Copy( const Matrix<T>& A, Matrix<T>& B );
-template<typename S,typename T,
-         typename=EnableIf<And<CanCast<S,T>,Not<IsSame<S,T>>>>>
-void Copy( const Matrix<S>& A, Matrix<T>& B );
-
-#ifdef HYDROGEN_HAVE_CUDA
-template <typename T>
-void Copy(Matrix<T,Device::CPU> const&, Matrix<T,Device::GPU>&);
-template <typename T>
-void Copy(Matrix<T,Device::GPU> const&, Matrix<T,Device::CPU>&);
-template <typename T>
-void CopyAsync(Matrix<T,Device::CPU> const&, Matrix<T,Device::GPU>&);
-template <typename T>
-void CopyAsync(Matrix<T,Device::GPU> const&, Matrix<T,Device::CPU>&);
-#endif // HYDROGEN_HAVE_CUDA
-
-template <typename T>
-void CopyAsync(ElementalMatrix<T> const& A, ElementalMatrix<T>& B);
-template <typename T>
-void CopyAsync(AbstractDistMatrix<T> const& A, AbstractDistMatrix<T>& B);
-
-template<typename S,typename T,
-         typename=EnableIf<CanCast<S,T>>>
-void Copy( const ElementalMatrix<S>& A, ElementalMatrix<T>& B );
-
-template<typename S,typename T,
-         typename=EnableIf<CanCast<S,T>>>
-void Copy( const BlockMatrix<S>& A, BlockMatrix<T>& B );
-
-template<typename T>
-void Copy( const AbstractDistMatrix<T>& A, AbstractDistMatrix<T>& B );
-template<typename S,typename T,
-         typename=EnableIf<And<CanCast<S,T>,Not<IsSame<S,T>>>>>
-void Copy( const AbstractDistMatrix<S>& A, AbstractDistMatrix<T>& B );
-
-template<typename T>
-void CopyFromRoot
-( const Matrix<T>& A, DistMatrix<T,CIRC,CIRC>& B,
-  bool includingViewers=false );
-template<typename T>
-void CopyFromNonRoot( DistMatrix<T,CIRC,CIRC>& B,
-  bool includingViewers=false );
-
-template<typename T>
-void CopyFromRoot
-( const Matrix<T>& A, DistMatrix<T,CIRC,CIRC,BLOCK>& B,
-  bool includingViewers=false );
-template<typename T>
-void CopyFromNonRoot( DistMatrix<T,CIRC,CIRC,BLOCK>& B,
-  bool includingViewers=false );
-
+class BaseDistMatrix;
 
 namespace copy {
 namespace util {
@@ -667,6 +613,13 @@ void StridedUnpack(
 } // namespace util
 } // namespace copy
 
+} // namespace El
+
+#include "Copy.hpp"
+
+namespace El
+{
+
 // DiagonalScale
 // =============
 #ifdef HYDROGEN_HAVE_CUDA
diff --git a/include/El/blas_like/level1/impl.hpp b/include/El/blas_like/level1/impl.hpp
index 5f6002ebc1..911b697520 100644
--- a/include/El/blas_like/level1/impl.hpp
+++ b/include/El/blas_like/level1/impl.hpp
@@ -19,7 +19,6 @@
 #include <El/blas_like/level1/ConjugateDiagonal.hpp>
 #include <El/blas_like/level1/ConjugateSubmatrix.hpp>
 #include <El/blas_like/level1/Contract.hpp>
-#include <El/blas_like/level1/Copy.hpp>
 #include <El/blas_like/level1/DiagonalScale.hpp>
 #include <El/blas_like/level1/DiagonalScaleTrapezoid.hpp>
 #include <El/blas_like/level1/DiagonalSolve.hpp>
diff --git a/include/El/core.hpp b/include/El/core.hpp
index 0c175b95c3..746cc0756e 100644
--- a/include/El/core.hpp
+++ b/include/El/core.hpp
@@ -111,17 +111,6 @@ using hydrogen::gpu_half_type;
 #include <quadmath.h>
 #endif
 
-#ifdef HYDROGEN_HAVE_HALF
-template <>
-struct std::is_floating_point<hydrogen::cpu_half_type> : std::true_type {};
-
-template <>
-struct std::is_integral<hydrogen::cpu_half_type> : std::false_type {};
-
-template <>
-struct std::is_arithmetic<hydrogen::cpu_half_type> : std::true_type {};
-#endif
-
 namespace El
 {
 
diff --git a/include/El/core/Element/decl.hpp b/include/El/core/Element/decl.hpp
index 383f0bb13d..72d30b19c2 100644
--- a/include/El/core/Element/decl.hpp
+++ b/include/El/core/Element/decl.hpp
@@ -684,10 +684,96 @@ template<typename Real,
 Complex<Real> Sqrt( const Complex<Real>& alpha );
 
 template<> Int Sqrt( const Int& alpha );
+#ifdef HYDROGEN_HAVE_HALF
+inline cpu_half_type Exp(cpu_half_type const& alpha)
+{
+    return half_float::exp(alpha);
+}
+inline cpu_half_type Log(cpu_half_type const& alpha)
+{
+    return half_float::log(alpha);
+}
+inline cpu_half_type Pow(cpu_half_type const& base, cpu_half_type const& exp)
+{
+    return half_float::pow(base, exp);
+}
+inline cpu_half_type Sqrt(cpu_half_type const& x)
+{
+    return half_float::sqrt(x);
+}
+inline cpu_half_type Expm1(cpu_half_type const& alpha)
+{
+    return half_float::expm1(alpha);
+}
+inline cpu_half_type Cos(cpu_half_type const& alpha)
+{
+    return half_float::cos(alpha);
+}
+inline cpu_half_type Sin(cpu_half_type const& alpha)
+{
+    return half_float::sin(alpha);
+}
+inline cpu_half_type Tan(cpu_half_type const& alpha)
+{
+    return half_float::tan(alpha);
+}
+inline cpu_half_type Acos(cpu_half_type const& alpha)
+{
+    return half_float::acos(alpha);
+}
+inline cpu_half_type Asin(cpu_half_type const& alpha)
+{
+    return half_float::asin(alpha);
+}
+inline cpu_half_type Atan(cpu_half_type const& alpha)
+{
+    return half_float::atan(alpha);
+}
+inline cpu_half_type Cosh(cpu_half_type const& alpha)
+{
+    return half_float::cosh(alpha);
+}
+inline cpu_half_type Sinh(cpu_half_type const& alpha)
+{
+    return half_float::sinh(alpha);
+}
+inline cpu_half_type Tanh(cpu_half_type const& alpha)
+{
+    return half_float::tanh(alpha);
+}
+inline cpu_half_type Acosh(cpu_half_type const& alpha)
+{
+    return half_float::acosh(alpha);
+}
+inline cpu_half_type Asinh(cpu_half_type const& alpha)
+{
+    return half_float::asinh(alpha);
+}
+inline cpu_half_type Atanh(cpu_half_type const& alpha)
+{
+    return half_float::atanh(alpha);
+}
+#endif // HYDROGEN_HAVE_HALF
 #ifdef HYDROGEN_GPU_USE_FP16
+inline gpu_half_type Exp(gpu_half_type const& x)
+{
+    return std::exp(float(x));
+}
+inline gpu_half_type Log(gpu_half_type const& x)
+{
+    return std::log(float(x));
+}
+inline gpu_half_type Pow(gpu_half_type const& base, gpu_half_type const& exp)
+{
+    return std::pow(float(base), float(exp));
+}
 inline gpu_half_type Sqrt(gpu_half_type const& x)
 {
-    return Sqrt(float(x));
+    return std::sqrt(float(x));
+}
+inline gpu_half_type Expm1(gpu_half_type const& x)
+{
+    return std::expm1(float(x));
 }
 #endif // HYDROGEN_GPU_USE_FP16
 #ifdef HYDROGEN_HAVE_QD
diff --git a/include/El/core/Element/impl.hpp b/include/El/core/Element/impl.hpp
index 4996bf6d37..4680bd2955 100644
--- a/include/El/core/Element/impl.hpp
+++ b/include/El/core/Element/impl.hpp
@@ -354,6 +354,19 @@ template<typename Real,
 Complex<Real> Exp( const Complex<Real>& alpha ) EL_NO_EXCEPT
 { return ComplexFromPolar( Exp(RealPart(alpha)), ImagPart(alpha) ); }
 
+template<typename F,
+         typename/*=EnableIf<IsField<F>>*/,
+         typename/*=EnableIf<IsStdScalar<F>>*/>
+F Expm1( const F& alpha ) EL_NO_EXCEPT
+{ return std::expm1(alpha); }
+
+template<typename Real,
+         typename/*=EnableIf<IsField<Real>>*/,
+         typename/*=DisableIf<IsStdScalar<Real>>*/,
+         typename/*=void*/>
+Complex<Real> Expm1( const Complex<Real>& alpha ) EL_NO_EXCEPT
+{ return ComplexFromPolar( Expm1(RealPart(alpha)), ImagPart(alpha) ); }
+
 template<typename FBase,typename TExp,
          typename/*=EnableIf<IsStdField<FBase>>*/,
          typename/*=EnableIf<IsStdField<TExp>>*/>
diff --git a/include/El/core/Matrix/impl_gpu.hpp b/include/El/core/Matrix/impl_gpu.hpp
index bde23abef8..cc11bcbdfb 100644
--- a/include/El/core/Matrix/impl_gpu.hpp
+++ b/include/El/core/Matrix/impl_gpu.hpp
@@ -81,7 +81,7 @@ Matrix<T, Device::GPU>::operator=(Matrix<T, Device::CPU> const& A)
 
 template <typename T>
 Matrix<T, Device::GPU>::Matrix(Matrix<T, Device::GPU>&& A) EL_NO_EXCEPT
-    : AbstractMatrix<T>{std::move(A)},
+    : AbstractMatrix<T>(std::move(A)),
       memory_{std::move(A.memory_)}, data_{A.data_}
 {
     A.data_ = nullptr;
diff --git a/include/El/core/random/decl.hpp b/include/El/core/random/decl.hpp
index a0e6f65a2e..57e8a45fe7 100644
--- a/include/El/core/random/decl.hpp
+++ b/include/El/core/random/decl.hpp
@@ -34,6 +34,11 @@ Int CoinFlip();
 template<typename T>
 T UnitCell();
 
+#ifdef HYDROGEN_HAVE_HALF
+cpu_half_type SampleUniform(
+    cpu_half_type const& a=0._h, cpu_half_type const& b=1._h);
+#endif
+
 template<typename Real=double,typename=EnableIf<IsReal<Real>>>
 Real SampleUniformNaive
 ( const Real& a=Real(0), const Real& b=UnitCell<Real>() );
diff --git a/include/El/core/random/impl.hpp b/include/El/core/random/impl.hpp
index f8c24343b7..00d17777a2 100644
--- a/include/El/core/random/impl.hpp
+++ b/include/El/core/random/impl.hpp
@@ -103,12 +103,22 @@ T UnitCell()
 {
     typedef Base<T> Real;
     T cell;
-    SetRealPart( cell, Real(1) );
+    SetRealPart( cell, TypeTraits<Real>::One() );
     if( IsComplex<T>::value )
-        SetImagPart( cell, Real(1) );
+        SetImagPart( cell, TypeTraits<Real>::One() );
     return cell;
 }
 
+#ifdef HYDROGEN_HAVE_HALF
+inline cpu_half_type SampleUniform(
+    cpu_half_type const& a, cpu_half_type const& b)
+{
+    std::mt19937& gen = Generator();
+    std::uniform_real_distribution<float> uni((float)a,(float)b);
+    return cpu_half_type{uni(gen)};
+}
+#endif // HYDROGEN_HAVE_HALF
+
 template<typename Real,typename>
 Real SampleUniformNaive( const Real& a, const Real& b )
 {
@@ -151,7 +161,7 @@ F SampleNormalMarsiglia( const F& mean, const Base<F>& stddev )
 
     Real stddevAdj = stddev;
     if( IsComplex<F>::value )
-        stddevAdj /= Sqrt(Real(2));
+        stddevAdj /= Sqrt(Real(2.));
 
     // Run Marsiglia's polar method
     // ============================
@@ -159,12 +169,12 @@ F SampleNormalMarsiglia( const F& mean, const Base<F>& stddev )
     //       F is real.
     while( true )
     {
-        const Real U = SampleUniform(Real(-1),Real(1));
-        const Real V = SampleUniform(Real(-1),Real(1));
+        const Real U = SampleUniform(Real(-1.),TypeTraits<Real>::One());
+        const Real V = SampleUniform(Real(-1.),TypeTraits<Real>::One());
         const Real S = Sqrt(U*U+V*V);
-        if( S > Real(0) && S < Real(1) )
+        if( S > TypeTraits<Real>::Zero() && S < TypeTraits<Real>::One() )
         {
-            const Real W = Sqrt(-2*Log(S)/S);
+            const Real W = Sqrt(Real(-2.)*Log(S)/S);
             SetRealPart( sample, RealPart(mean) + stddevAdj*U*W );
             if( IsComplex<F>::value )
                 SetImagPart( sample, ImagPart(mean) + stddevAdj*V*W );
@@ -211,8 +221,8 @@ template<typename F>
 F SampleBall( const F& center, const Base<F>& radius )
 {
     typedef Base<F> Real;
-    const Real r = SampleUniform(Real(0),radius);
-    const Real angle = SampleUniform(Real(0),2*Pi<Real>());
+    const Real r = SampleUniform(TypeTraits<Real>::Zero(),radius);
+    const Real angle = SampleUniform(TypeTraits<Real>::Zero(),2*Pi<Real>());
     return center + F(r*Cos(angle),r*Sin(angle));
 }
 
diff --git a/include/El/io.hpp b/include/El/io.hpp
index 97808146b8..2dad5e8b26 100644
--- a/include/El/io.hpp
+++ b/include/El/io.hpp
@@ -56,6 +56,9 @@ void Display( const AbstractDistMatrix<T>& A, string title="DistMatrix" );
 // Dense
 // -----
 template<typename T>
+void Print(const AbstractMatrix<T>& A, string title="Matrix",
+           ostream& os=cout);
+template<typename T>
 void Print( const Matrix<T>& A, string title="Matrix", ostream& os=cout );
 template<typename T>
 void Print
diff --git a/include/hydrogen/blas/GPU_BLAS_impl.hpp b/include/hydrogen/blas/GPU_BLAS_impl.hpp
index 24fb1289e5..82e7caf2a5 100644
--- a/include/hydrogen/blas/GPU_BLAS_impl.hpp
+++ b/include/hydrogen/blas/GPU_BLAS_impl.hpp
@@ -474,8 +474,51 @@ void AxpyImpl(SizeT const&, T const&,
     throw std::logic_error(oss.str());
 }
 
+template <typename T, typename SizeT,
+          typename=EnableUnless<IsSupportedType<T, BLAS_Op::GEMV>>,
+          typename=EnableWhen<IsSupportedType<T, BLAS_Op::GEMM>>>
+void GemvImpl(
+    TransposeMode transA,
+    SizeT nrows, SizeT ncols,
+    T const& alpha,
+    T const* A, SizeT lda,
+    T const* x, SizeT incx,
+    T const& beta,
+    T* y, SizeT incy,
+    SyncInfo<Device::GPU> const& si)
+{
+    using NTP = MakePointer<NativeType<T>>;
+    using CNTP = MakePointerToConst<NativeType<T>>;
+
+    if (incy != SizeT(1))
+        throw std::runtime_error("incy must be 1 right now. "
+                                 "Let Tom know you've hit this case.");
+
+    auto const ATrans = transA;
+    auto const BTrans = (incx == SizeT(1)
+                         ? TransposeMode::NORMAL
+                         : TransposeMode::TRANSPOSE);
+    auto const m = (ATrans == TransposeMode::NORMAL ? nrows : ncols);
+    auto const k = (ATrans == TransposeMode::NORMAL ? ncols : nrows);
+    auto const n = SizeT(1);
+    auto const LDB = (incx == SizeT(1) ? ncols : incx);
+    auto const LDC = nrows;
+
+    SyncManager mgr(GetLibraryHandle(), si);
+    gpu_blas_impl::Gemm(
+        GetLibraryHandle(),
+        ToNativeTransposeMode(ATrans), ToNativeTransposeMode(BTrans),
+        ToSizeT(m), ToSizeT(n), ToSizeT(k),
+        alpha,
+        reinterpret_cast<CNTP>(A), ToSizeT(lda),
+        reinterpret_cast<CNTP>(x), ToSizeT(LDB),
+        beta,
+        reinterpret_cast<NTP>(y), ToSizeT(LDC));
+}
+
 template <typename T, typename SizeT,
           typename=EnableUnless<IsSupportedType<T,BLAS_Op::GEMV>>,
+          typename=EnableUnless<IsSupportedType<T,BLAS_Op::GEMM>>,
           typename=void>
 void GemvImpl(
     TransposeMode const&,
diff --git a/include/hydrogen/utils/HalfPrecision.hpp b/include/hydrogen/utils/HalfPrecision.hpp
index 55593fd534..117e6795df 100644
--- a/include/hydrogen/utils/HalfPrecision.hpp
+++ b/include/hydrogen/utils/HalfPrecision.hpp
@@ -6,64 +6,202 @@
  *  Here we handle a few high-level metaprogramming tasks for 16-bit
  *  floating point support. This currently just includes IEEE-754
  *  support but will likely soon expand to include "bfloat16" as well.
+ *
+ *  When compiling using v2.1.0 of the Half library on OSX (10.14 with
+ *  Apple LLVM version 10.0.1 (clang-1001.0.46.4)), there were many
+ *  interesting issues that showed up. Hence the strange layout of
+ *  this file.
  */
 
+#ifdef HYDROGEN_HAVE_HALF
+
 #include <El/hydrogen_config.h>
 #include <hydrogen/meta/TypeTraits.hpp>
 
-#ifdef HYDROGEN_HAVE_HALF
+// Forward-declare things so I can start specializing templates.
+namespace half_float
+{
+class half;
+}// namespace half_float
+
+// Specialize some STL-y things
+#include <type_traits>
+namespace std
+{
+template <>
+struct is_floating_point<half_float::half> : true_type {};
+
+template <>
+struct is_integral<half_float::half> : false_type {};
+
+template <>
+struct is_arithmetic<half_float::half> : true_type {};
+}// namespace std
+
+// Fix an issue inside the Half library
+#ifndef HALF_ENABLE_F16C_INTRINSICS
+#define HALF_ENABLE_F16C_INTRINSICS __F16C__
+#endif
+
+// Now include the actual Half library.
 #include <half.hpp>
+
+// Declare the hydrogen typedef
+namespace hydrogen
+{
+using namespace half_float::literal;
+
+using cpu_half_type = half_float::half;
+
+template <>
+struct TypeTraits<cpu_half_type>
+{
+    static cpu_half_type One() noexcept { return 1._h; }
+    static cpu_half_type Zero() noexcept { return 0._h; }
+    static std::string Name() { return typeid(cpu_half_type).name(); }
+};// struct TypeTraits<cpu_half_type>
+
+}// namespace hydrogen
+
+/** @name Bitwise operators, for MPI reduction. */
+///@{
+inline hydrogen::cpu_half_type operator~(hydrogen::cpu_half_type const&)
+{
+    throw std::logic_error(
+        "Bitwise operations not supported for floating-point types.");
+}
+inline hydrogen::cpu_half_type operator|(hydrogen::cpu_half_type const&,
+                                         hydrogen::cpu_half_type const&)
+{
+    throw std::logic_error(
+        "Bitwise operations not supported for floating-point types.");
+}
+inline hydrogen::cpu_half_type operator&(hydrogen::cpu_half_type const&,
+                                         hydrogen::cpu_half_type const&)
+{
+    throw std::logic_error(
+        "Bitwise operations not supported for floating-point types.");
+}
+inline hydrogen::cpu_half_type operator^(hydrogen::cpu_half_type const&,
+                                         hydrogen::cpu_half_type const&)
+{
+    throw std::logic_error(
+        "Bitwise operations not supported for floating-point types.");
+}
+///@}
 #endif // HYDROGEN_HAVE_HALF
 
+// Finally, do the GPU stuff
 #ifdef HYDROGEN_GPU_USE_FP16
+
+// Grab the right header
 #if defined(HYDROGEN_HAVE_CUDA)
 #include <cuda_fp16.h>
 #elif defined(HYDROGEN_HAVE_AMDGPU)
 #include <rocblas-types.h>
 #endif // HYDROGEN_HAVE_CUDA
-#endif // HYDROGEN_GPU_USE_FP16
 
 namespace hydrogen
 {
 
-#ifdef HYDROGEN_HAVE_HALF
-/** @brief Unified name for the FP16 type on CPU */
-using cpu_half_type = half_float::half;
-#endif // HYDROGEN_HAVE_HALF
-
-
-#ifdef HYDROGEN_GPU_USE_FP16
 #if defined(HYDROGEN_HAVE_CUDA)
 /** @brief Unified name for the FP16 type on GPU */
 using gpu_half_type = __half;
 
+template <>
+struct TypeTraits<gpu_half_type>
+{
+    static gpu_half_type One() noexcept { return 1.f; }
+    static gpu_half_type Zero() noexcept { return 0.f; }
+    static std::string Name() { return typeid(gpu_half_type).name(); }
+};// struct TypeTraits<gpu_half_type>
+
+#elif defined(HYDROGEN_HAVE_AMDGPU)
+/** @brief Unified name for the FP16 type on GPU */
+using gpu_half_type = rocblas_half;
+#endif // HYDROGEN_HAVE_CUDA
+
+}// namespace hydrogen
+
+#if defined(HYDROGEN_HAVE_CUDA) && !defined(__CUDACC__)
+
 /** @brief Enable "update" functionality for __half. */
-inline gpu_half_type& operator+=(gpu_half_type& val, gpu_half_type const& rhs)
+template <typename T>
+hydrogen::gpu_half_type& operator+=(
+    hydrogen::gpu_half_type& val, T const& rhs)
 {
     val = float(val) + float(rhs);
     return val;
 }
 
+/** @brief Enable subtract-equal functionality for __half. */
+template <typename T>
+hydrogen::gpu_half_type& operator-=(
+    hydrogen::gpu_half_type& val, T const& rhs)
+{
+    val = float(val) - float(rhs);
+    return val;
+}
+
 /** @brief Enable "scale" functionality for __half. */
-inline gpu_half_type& operator*=(gpu_half_type& val, gpu_half_type const& rhs)
+template <typename T>
+hydrogen::gpu_half_type& operator*=(
+    hydrogen::gpu_half_type& val, T const& rhs)
 {
     val = float(val) * float(rhs);
     return val;
 }
 
-template <>
-struct TypeTraits<gpu_half_type>
+/** @brief Enable divide-equal functionality for __half. */
+template <typename T>
+hydrogen::gpu_half_type& operator/=(
+    hydrogen::gpu_half_type& val, T const& rhs)
 {
-    static gpu_half_type One() noexcept { return 1.f; }
-    static gpu_half_type Zero() noexcept { return 0.f; }
-    static std::string Name() { return typeid(gpu_half_type).name(); }
-};// struct TypeTraits<gpu_half_type>
+    val = float(val) / float(rhs);
+    return val;
+}
 
-#elif defined(HYDROGEN_HAVE_AMDGPU)
-/** @brief Unified name for the FP16 type on GPU */
-using gpu_half_type = rocblas_half;
-#endif // HYDROGEN_HAVE_CUDA
-#endif // HYDROGEN_GPU_USE_FP16
+/** @brief Enable add functionality for __half. */
+inline hydrogen::gpu_half_type operator+(
+    hydrogen::gpu_half_type const& val, hydrogen::gpu_half_type const& rhs)
+{
+    return float(val) + float(rhs);
+}
 
-}// namespace hydrogen
+/** @brief Enable subtract functionality for __half. */
+inline hydrogen::gpu_half_type operator-(
+    hydrogen::gpu_half_type const& val, hydrogen::gpu_half_type const& rhs)
+{
+    return float(val) - float(rhs);
+}
+
+/** @brief Enable multiply functionality for __half. */
+inline hydrogen::gpu_half_type operator*(
+    hydrogen::gpu_half_type const& val, hydrogen::gpu_half_type const& rhs)
+{
+    return float(val) * float(rhs);
+}
+
+/** @brief Enable divide functionality for __half. */
+inline hydrogen::gpu_half_type operator/(
+    hydrogen::gpu_half_type const& val, hydrogen::gpu_half_type const& rhs)
+{
+    return float(val) / float(rhs);
+}
+
+/** @brief Enable unary minus functionality for __half. */
+inline hydrogen::gpu_half_type operator-(
+    hydrogen::gpu_half_type const& val)
+{
+    return -float(val);
+}
+
+#endif // defined(HYDROGEN_HAVE_CUDA) && !defined(__CUDACC__)
+
+inline std::ostream& operator<<(std::ostream& os, hydrogen::gpu_half_type const& x)
+{
+    return os << float(x) << "_h";
+}
+
+#endif // HYDROGEN_GPU_USE_FP16
 #endif // HYDROGEN_UTILS_HALFPRECISION_HPP_
diff --git a/include/hydrogen/utils/NumericTypeConversion.hpp b/include/hydrogen/utils/NumericTypeConversion.hpp
index bd17e43517..b30f4c3e32 100644
--- a/include/hydrogen/utils/NumericTypeConversion.hpp
+++ b/include/hydrogen/utils/NumericTypeConversion.hpp
@@ -41,6 +41,25 @@ struct Caster<F, half_float::half>
     }
 };
 
+#ifdef HYDROGEN_GPU_USE_FP16
+template <typename F>
+struct Caster<F, __half>
+{
+    template <typename DF,
+              EnableWhen<std::is_integral<DF>, int> = 0>
+    static __half Cast(DF const& x)
+    {
+        return __half(float(x));
+    }
+
+    template <typename DF,
+              EnableUnless<std::is_integral<DF>, int> = 0>
+    static __half Cast(DF const& x)
+    {
+        return static_cast<__half>(x);
+    }
+};
+#endif // HYDROGEN_GPU_USE_FP16
 
 template <typename T, typename F>
 T To(F const& x)
diff --git a/src/blas_like/level1/CMakeLists.txt b/src/blas_like/level1/CMakeLists.txt
index b2a137af98..3703543526 100644
--- a/src/blas_like/level1/CMakeLists.txt
+++ b/src/blas_like/level1/CMakeLists.txt
@@ -2,6 +2,7 @@
 set_full_path(THIS_DIR_SOURCES
   ColumnMinAbs.cpp
   ColumnNorms.cpp
+  Copy.cpp
   HilbertSchmidt.cpp
   Instantiate.cpp
   Max.cpp
diff --git a/src/blas_like/level1/Copy.cpp b/src/blas_like/level1/Copy.cpp
new file mode 100644
index 0000000000..ed0f85fe1c
--- /dev/null
+++ b/src/blas_like/level1/Copy.cpp
@@ -0,0 +1,113 @@
+#include "El/core.hpp"
+#include "El/blas_like/level1/Copy.hpp"
+
+static_assert(std::is_integral<El::Int>::value,
+              "El::Int should be integral!");
+
+namespace El
+{
+namespace
+{
+
+// For now, I just want to generate the tensor-product of
+// {MatrixTypes}^2 for Copy. This will cover most of the usecases in
+// LBANN. AFAIK, the integer matrices are always converted to a
+// floating-point type before they interact "virtually"; all
+// operations on those matrices are dispatched statically.
+
+using MatrixTypes = TypeList<
+    float, double, Complex<float>, Complex<double>
+#ifdef HYDROGEN_HAVE_HALF
+    , cpu_half_type
+#endif // HYDROGEN_HAVE_HALF
+#ifdef HYDROGEN_GPU_USE_FP16
+    , gpu_half_type
+#endif // HYDROGEN_GPU_USE_FP16
+    >;
+
+template <template <typename> class X, typename... Ts>
+using Expand = TypeList<X<Ts>...>;
+
+template <template <typename> class X, typename List>
+struct ExpandTLT {};
+
+template <template <typename> class X, typename... Ts>
+struct ExpandTLT<X, TypeList<Ts...>>
+{
+    using type = Expand<X, Ts...>;
+};
+
+template <template <typename> class X, typename List>
+using ExpandTL = typename ExpandTLT<X, List>::type;
+
+// This is replaced by a generic multiple dispatch engine in
+// DiHydrogen; this is a one-off use-case for now, so there's no need
+// to backport a robust implementation.
+template <typename FunctorT, typename LHSList, typename RHSList>
+struct CopyDispatcher
+{
+    static void Do(FunctorT f,
+                   BaseDistMatrix const& src, BaseDistMatrix& tgt)
+    {
+        using LHead = Head<LHSList>;
+        using LTail = Tail<LHSList>;
+        if (auto const* ptr = dynamic_cast<LHead const*>(&src))
+            return CopyDispatcher<FunctorT, LHSList, RHSList>::DoRHS(
+                f, *ptr, tgt);
+        else
+            return CopyDispatcher<FunctorT, LTail, RHSList>::Do(f, src, tgt);
+    }
+
+    template <typename LHSType>
+    static void DoRHS(FunctorT f, LHSType const& src, BaseDistMatrix& tgt)
+    {
+        using RHead = Head<RHSList>;
+        using RTail = Tail<RHSList>;
+        if (auto* ptr = dynamic_cast<RHead*>(&tgt))
+            return f(src, *ptr);
+        else
+            return CopyDispatcher<FunctorT, LHSList, RTail>::DoRHS(f, src, tgt);
+    }
+};// struct CopyDispatcher
+
+template <typename FunctorT, typename RHSList>
+struct CopyDispatcher<FunctorT, TypeList<>, RHSList>
+{
+    static void Do(FunctorT const&,
+                   BaseDistMatrix const&, BaseDistMatrix const&)
+    {
+        LogicError("Source matrix type not found.");
+    }
+};
+
+template <typename FunctorT, typename LHSList>
+struct CopyDispatcher<FunctorT, LHSList, TypeList<>>
+{
+    static void DoRHS(FunctorT const&,
+                      BaseDistMatrix const&, BaseDistMatrix const&)
+    {
+        LogicError("Target matrix type not found.");
+    }
+};
+
+}// namespace <anon>
+
+void Copy(BaseDistMatrix const& Source, BaseDistMatrix& Target)
+{
+    using FunctorT = details::CopyFunctor;
+    using MatrixTs = ExpandTL<AbstractDistMatrix, MatrixTypes>;
+    using Dispatcher = CopyDispatcher<FunctorT, MatrixTs, MatrixTs>;
+    FunctorT f;
+    return Dispatcher::Do(f, Source, Target);
+}
+
+void CopyAsync(BaseDistMatrix const& Source, BaseDistMatrix& Target)
+{
+    using FunctorT = details::CopyAsyncFunctor;
+    using MatrixTs = ExpandTL<AbstractDistMatrix, MatrixTypes>;
+    using Dispatcher = CopyDispatcher<FunctorT, MatrixTs, MatrixTs>;
+    FunctorT f;
+    return Dispatcher::Do(f, Source, Target);
+}
+
+}// namespace El
diff --git a/src/blas_like/level1/HilbertSchmidt.cpp b/src/blas_like/level1/HilbertSchmidt.cpp
index 530471670a..c24d79b564 100644
--- a/src/blas_like/level1/HilbertSchmidt.cpp
+++ b/src/blas_like/level1/HilbertSchmidt.cpp
@@ -23,7 +23,7 @@ Ring HilbertSchmidt(
     if (A.GetDevice() != Device::CPU || A.GetDevice() != B.GetDevice())
         LogicError("HilbertSchmidt not supported for this device.");
 
-    Ring innerProd(0);
+    Ring innerProd(0.);
     const Int width = A.Width();
     const Int height = A.Height();
     const Ring* ABuf = A.LockedBuffer();
@@ -72,7 +72,7 @@ Ring HilbertSchmidt
     Ring innerProd;
     if( A.Participating() )
     {
-        Ring localInnerProd(0);
+        Ring localInnerProd(0.);
         const Int localHeight = A.LocalHeight();
         const Int localWidth = A.LocalWidth();
         const Ring* ABuf = A.LockedBuffer();
@@ -104,6 +104,10 @@ Ring HilbertSchmidt
   template Ring HilbertSchmidt \
   ( const AbstractDistMatrix<Ring>& A, const AbstractDistMatrix<Ring>& B );
 
+#ifdef HYDROGEN_GPU_USE_FP16
+PROTO(gpu_half_type)
+#endif
+
 #define EL_ENABLE_DOUBLEDOUBLE
 #define EL_ENABLE_QUADDOUBLE
 #define EL_ENABLE_QUAD
diff --git a/src/blas_like/level2/Gemv.cpp b/src/blas_like/level2/Gemv.cpp
index 75385e340f..5556885af2 100644
--- a/src/blas_like/level2/Gemv.cpp
+++ b/src/blas_like/level2/Gemv.cpp
@@ -322,6 +322,19 @@ template void Gemv(Orientation orientA,
                    Matrix<double,Device::GPU> const& B,
                    double beta,
                    Matrix<double,Device::GPU>& C);
+#ifdef HYDROGEN_GPU_USE_FP16
+template void Gemv(Orientation orientA,
+                   gpu_half_type alpha,
+                   Matrix<gpu_half_type,Device::GPU> const& A,
+                   Matrix<gpu_half_type,Device::GPU> const& B,
+                   gpu_half_type beta,
+                   Matrix<gpu_half_type,Device::GPU>& C);
+template void Gemv(Orientation, gpu_half_type,
+                   AbstractMatrix<gpu_half_type> const&,
+                   AbstractMatrix<gpu_half_type> const&,
+                   gpu_half_type,
+                   AbstractMatrix<gpu_half_type>&);
+#endif // HYDROGEN_GPU_USE_FP16
 #endif // HYDROGEN_HAVE_CUDA
 
 #define PROTO(T) \
diff --git a/src/core/DistMatrix/AbstractDistMatrix.cpp b/src/core/DistMatrix/AbstractDistMatrix.cpp
index 284974254d..91e3a2ca3c 100644
--- a/src/core/DistMatrix/AbstractDistMatrix.cpp
+++ b/src/core/DistMatrix/AbstractDistMatrix.cpp
@@ -541,6 +541,61 @@ AbstractDistMatrix<double>::Instantiate
     return Instantiate_<double>(grid, root, colDist, rowDist, wrap, dev);
 }
 
+#ifdef HYDROGEN_HAVE_HALF
+
+
+template<>
+AbstractDistMatrix<cpu_half_type>*
+AbstractDistMatrix<cpu_half_type>::Instantiate
+(const El::Grid& grid, int root,
+ Dist colDist, Dist rowDist, DistWrap wrap, Device dev)
+{
+    if (dev != Device::CPU)
+    {
+        LogicError("cpu_half_type only supported on CPU");
+    }
+
+#define GUARD(CDIST,RDIST,WRAP)                                         \
+    (colDist == CDIST) && (rowDist == RDIST) && (wrap == WRAP)
+#define PAYLOAD(CDIST,RDIST,WRAP)                                       \
+    return new DistMatrix<cpu_half_type,CDIST,RDIST,WRAP,Device::CPU>(grid, root);
+#include "El/macros/GuardAndPayload.h"
+    LogicError
+        ("Invalid template arguments for DistMatrix "
+         "(colDist=",Int(colDist),", rowDist=",Int(rowDist),", "
+         "wrap=",Int(wrap),", dev=",Int(dev),")");
+    return nullptr;
+}
+
+#ifdef HYDROGEN_GPU_USE_FP16
+template<>
+AbstractDistMatrix<gpu_half_type>*
+AbstractDistMatrix<gpu_half_type>::Instantiate
+(const El::Grid& grid, int root,
+ Dist colDist, Dist rowDist, DistWrap wrap, Device dev)
+{
+    if (dev != Device::GPU)
+    {
+        LogicError("gpu_half_type only supported on GPU");
+    }
+    if (wrap != ELEMENT)
+    {
+        LogicError("GPU matrices only supported with ELEMENT distributions");
+    }
+#define GUARD(CDIST,RDIST,WRAP)                                         \
+    (colDist == CDIST) && (rowDist == RDIST) && (wrap == WRAP)
+#define PAYLOAD(CDIST,RDIST,WRAP)                                       \
+    return new DistMatrix<gpu_half_type,CDIST,RDIST,ELEMENT,Device::GPU>(grid, root);
+#include "El/macros/GuardAndPayload.h"
+    LogicError
+        ("Invalid template arguments for DistMatrix "
+         "(colDist=",Int(colDist),", rowDist=",Int(rowDist),", "
+         "wrap=",Int(wrap),", dev=",Int(dev),")");
+    return nullptr;
+}
+#endif // HYDROGEN_GPU_USE_FP16
+#endif // HYDROGEN_HAVE_HALF
+
 template<typename T>
 AbstractDistMatrix<T>*
 AbstractDistMatrix<T>::Instantiate(const El::DistData& data)
diff --git a/src/core/imports/blas/Dot.hpp b/src/core/imports/blas/Dot.hpp
index f5a8c4e7ef..b232de1d8c 100644
--- a/src/core/imports/blas/Dot.hpp
+++ b/src/core/imports/blas/Dot.hpp
@@ -32,7 +32,7 @@ T Dot
     // NOTE: Temporaries are avoided since constructing a BigInt/BigFloat
     //       involves a memory allocation
     T gamma;
-    T alpha{0};
+    T alpha(0.);
     for( BlasInt i=0; i<n; ++i )
     {
         Conj( x[i*incx], gamma );
@@ -91,6 +91,12 @@ template cpu_half_type Dot
   const cpu_half_type* x, BlasInt incx,
   const cpu_half_type* y, BlasInt incy );
 #endif
+#ifdef HYDROGEN_GPU_USE_FP16
+template gpu_half_type Dot
+( BlasInt n,
+  const gpu_half_type* x, BlasInt incx,
+  const gpu_half_type* y, BlasInt incy );
+#endif
 #ifdef HYDROGEN_HAVE_MPC
 template BigInt Dot
 ( BlasInt n,
diff --git a/src/core/imports/blas/Gemv.hpp b/src/core/imports/blas/Gemv.hpp
index 5d8d627355..804d12e970 100644
--- a/src/core/imports/blas/Gemv.hpp
+++ b/src/core/imports/blas/Gemv.hpp
@@ -2,8 +2,8 @@
    Copyright (c) 2009-2016, Jack Poulson
    All rights reserved.
 
-   This file is part of Elemental and is under the BSD 2-Clause License, 
-   which can be found in the LICENSE file in the root directory, or at 
+   This file is part of Elemental and is under the BSD 2-Clause License,
+   which can be found in the LICENSE file in the root directory, or at
    http://opensource.org/licenses/BSD-2-Clause
 */
 
@@ -54,10 +54,10 @@ void Gemv
     T gamma, delta;
     if( std::toupper(trans) == 'N' )
     {
-        if( m > 0 && n == 0 && beta == T(0) )
+        if( m > 0 && n == 0 && beta == TypeTraits<T>::Zero() )
         {
             for( BlasInt i=0; i<m; ++i )
-                y[i*incy] = 0;
+                y[i*incy] = TypeTraits<T>::Zero();
             return;
         }
 
@@ -77,10 +77,10 @@ void Gemv
     }
     else if( std::toupper(trans) == 'T' )
     {
-        if( n > 0 && m == 0 && beta == T(0) )
+        if( n > 0 && m == 0 && beta == TypeTraits<T>::Zero() )
         {
             for( BlasInt i=0; i<n; ++i )
-                y[i*incy] = 0;
+                y[i*incy] = TypeTraits<T>::Zero();
             return;
         }
         Scal( n, beta, y, incy );
@@ -106,10 +106,10 @@ void Gemv
     }
     else
     {
-        if( n > 0 && m == 0 && beta == T(0) )
+        if( n > 0 && m == 0 && beta == TypeTraits<T>::Zero() )
         {
             for( BlasInt i=0; i<n; ++i )
-                y[i*incy] = 0;
+                y[i*incy] = TypeTraits<T>::Zero();
             return;
         }
         Scal( n, beta, y, incy );
@@ -135,87 +135,96 @@ void Gemv
     }
 }
 template void Gemv
-( char trans, BlasInt m, BlasInt n, 
+( char trans, BlasInt m, BlasInt n,
   const Int& alpha,
   const Int* A, BlasInt ALDim,
-  const Int* x, BlasInt incx, 
+  const Int* x, BlasInt incx,
   const Int& beta,
         Int* y, BlasInt incy );
 #ifdef HYDROGEN_HAVE_QD
 template void Gemv
-( char trans, BlasInt m, BlasInt n, 
+( char trans, BlasInt m, BlasInt n,
   const DoubleDouble& alpha,
   const DoubleDouble* A, BlasInt ALDim,
-  const DoubleDouble* x, BlasInt incx, 
+  const DoubleDouble* x, BlasInt incx,
   const DoubleDouble& beta,
         DoubleDouble* y, BlasInt incy );
 template void Gemv
-( char trans, BlasInt m, BlasInt n, 
+( char trans, BlasInt m, BlasInt n,
   const QuadDouble& alpha,
   const QuadDouble* A, BlasInt ALDim,
-  const QuadDouble* x, BlasInt incx, 
+  const QuadDouble* x, BlasInt incx,
   const QuadDouble& beta,
         QuadDouble* y, BlasInt incy );
 template void Gemv
-( char trans, BlasInt m, BlasInt n, 
+( char trans, BlasInt m, BlasInt n,
   const Complex<DoubleDouble>& alpha,
   const Complex<DoubleDouble>* A, BlasInt ALDim,
-  const Complex<DoubleDouble>* x, BlasInt incx, 
+  const Complex<DoubleDouble>* x, BlasInt incx,
   const Complex<DoubleDouble>& beta,
         Complex<DoubleDouble>* y, BlasInt incy );
 template void Gemv
-( char trans, BlasInt m, BlasInt n, 
+( char trans, BlasInt m, BlasInt n,
   const Complex<QuadDouble>& alpha,
   const Complex<QuadDouble>* A, BlasInt ALDim,
-  const Complex<QuadDouble>* x, BlasInt incx, 
+  const Complex<QuadDouble>* x, BlasInt incx,
   const Complex<QuadDouble>& beta,
         Complex<QuadDouble>* y, BlasInt incy );
 #endif
 #ifdef HYDROGEN_HAVE_HALF
 template void Gemv
-( char trans, BlasInt m, BlasInt n, 
+( char trans, BlasInt m, BlasInt n,
   const cpu_half_type& alpha,
   const cpu_half_type* A, BlasInt ALDim,
-  const cpu_half_type* x, BlasInt incx, 
+  const cpu_half_type* x, BlasInt incx,
   const cpu_half_type& beta,
         cpu_half_type* y, BlasInt incy );
 #endif
+#ifdef HYDROGEN_GPU_USE_FP16
+template void Gemv
+( char trans, BlasInt m, BlasInt n,
+  const gpu_half_type& alpha,
+  const gpu_half_type* A, BlasInt ALDim,
+  const gpu_half_type* x, BlasInt incx,
+  const gpu_half_type& beta,
+        gpu_half_type* y, BlasInt incy );
+#endif
 #ifdef HYDROGEN_HAVE_QUADMATH
 template void Gemv
-( char trans, BlasInt m, BlasInt n, 
+( char trans, BlasInt m, BlasInt n,
   const Quad& alpha,
   const Quad* A, BlasInt ALDim,
-  const Quad* x, BlasInt incx, 
+  const Quad* x, BlasInt incx,
   const Quad& beta,
         Quad* y, BlasInt incy );
 template void Gemv
-( char trans, BlasInt m, BlasInt n, 
+( char trans, BlasInt m, BlasInt n,
   const Complex<Quad>& alpha,
-  const Complex<Quad>* A, BlasInt ALDim, 
-  const Complex<Quad>* x, BlasInt incx, 
+  const Complex<Quad>* A, BlasInt ALDim,
+  const Complex<Quad>* x, BlasInt incx,
   const Complex<Quad>& beta,
         Complex<Quad>* y, BlasInt incy );
 #endif
 #ifdef HYDROGEN_HAVE_MPC
 template void Gemv
-( char trans, BlasInt m, BlasInt n, 
+( char trans, BlasInt m, BlasInt n,
   const BigInt& alpha,
   const BigInt* A, BlasInt ALDim,
-  const BigInt* x, BlasInt incx, 
+  const BigInt* x, BlasInt incx,
   const BigInt& beta,
         BigInt* y, BlasInt incy );
 template void Gemv
-( char trans, BlasInt m, BlasInt n, 
+( char trans, BlasInt m, BlasInt n,
   const BigFloat& alpha,
   const BigFloat* A, BlasInt ALDim,
-  const BigFloat* x, BlasInt incx, 
+  const BigFloat* x, BlasInt incx,
   const BigFloat& beta,
         BigFloat* y, BlasInt incy );
 template void Gemv
-( char trans, BlasInt m, BlasInt n, 
+( char trans, BlasInt m, BlasInt n,
   const Complex<BigFloat>& alpha,
   const Complex<BigFloat>* A, BlasInt ALDim,
-  const Complex<BigFloat>* x, BlasInt incx, 
+  const Complex<BigFloat>* x, BlasInt incx,
   const Complex<BigFloat>& beta,
         Complex<BigFloat>* y, BlasInt incy );
 #endif
@@ -249,7 +258,7 @@ void Gemv
 void Gemv
 ( char trans, BlasInt m, BlasInt n,
   const scomplex& alpha,
-  const scomplex* A, BlasInt ALDim, 
+  const scomplex* A, BlasInt ALDim,
   const scomplex* x, BlasInt incx,
   const scomplex& beta,
         scomplex* y, BlasInt incy )
@@ -259,7 +268,7 @@ void Gemv
 void Gemv
 ( char trans, BlasInt m, BlasInt n,
   const dcomplex& alpha,
-  const dcomplex* A, BlasInt ALDim, 
+  const dcomplex* A, BlasInt ALDim,
   const dcomplex* x, BlasInt incx,
   const dcomplex& beta,
         dcomplex* y, BlasInt incy )
diff --git a/src/core/imports/cuda.cpp b/src/core/imports/cuda.cpp
index a393dc57a9..24ab734c87 100644
--- a/src/core/imports/cuda.cpp
+++ b/src/core/imports/cuda.cpp
@@ -28,7 +28,6 @@ void InitializeCUDA(int argc, char* argv[])
 
     unsigned int numDevices = 0;
     int device = 0;
-
     nvmlReturn_t r = nvmlInit();
     if (r != NVML_SUCCESS) { throw std::runtime_error("NVML error"); }
     r = nvmlDeviceGetCount(&numDevices);
diff --git a/src/hydrogen/blas/gpu/Axpy.cu b/src/hydrogen/blas/gpu/Axpy.cu
index dc5d6b0536..1773463df3 100644
--- a/src/hydrogen/blas/gpu/Axpy.cu
+++ b/src/hydrogen/blas/gpu/Axpy.cu
@@ -1,6 +1,7 @@
 #include <hydrogen/blas/gpu/Axpy.hpp>
 
 #include <El/hydrogen_config.h>
+#include <hydrogen/meta/TypeTraits.hpp>
 #include <hydrogen/device/gpu/CUDA.hpp>
 
 #include <cuda_runtime.h>
@@ -115,7 +116,7 @@ void Axpy_GPU_impl(
     {
         return;
     }
-    
+
     constexpr int TILE_SIZE = 32;
     constexpr int BLK_COLS = 8;
 
diff --git a/src/hydrogen/blas/gpu/Copy.cu b/src/hydrogen/blas/gpu/Copy.cu
index 6d31cd636d..ee03c7e401 100644
--- a/src/hydrogen/blas/gpu/Copy.cu
+++ b/src/hydrogen/blas/gpu/Copy.cu
@@ -1,7 +1,9 @@
 #include <hydrogen/blas/gpu/Copy.hpp>
 
 #include <El/hydrogen_config.h>
+#include <hydrogen/meta/TypeTraits.hpp>
 #include <hydrogen/device/gpu/CUDA.hpp>
+
 #include <cuda_runtime.h>
 
 namespace
@@ -82,7 +84,7 @@ void Copy_GPU_impl(
 {
     if (num_entries <= TypeTraits<SizeT>::Zero())
         return;
-    
+
 #ifdef HYDROGEN_DO_BOUNDS_CHECKING
     // The kernel parameters are __restrict__-ed. This helps ensure
     // that's not a lie.
diff --git a/src/hydrogen/blas/gpu/Fill.cu b/src/hydrogen/blas/gpu/Fill.cu
index 32c4e178f4..456d8608bd 100644
--- a/src/hydrogen/blas/gpu/Fill.cu
+++ b/src/hydrogen/blas/gpu/Fill.cu
@@ -1,7 +1,9 @@
 #include <hydrogen/blas/gpu/Fill.hpp>
 
 #include <El/hydrogen_config.h>
+#include <hydrogen/meta/TypeTraits.hpp>
 #include <hydrogen/device/gpu/CUDA.hpp>
+
 #include <cuda_runtime.h>
 
 namespace hydrogen
diff --git a/src/hydrogen/blas/gpu/Scale.cu b/src/hydrogen/blas/gpu/Scale.cu
index 1e8507a3eb..f55bcde28f 100644
--- a/src/hydrogen/blas/gpu/Scale.cu
+++ b/src/hydrogen/blas/gpu/Scale.cu
@@ -1,6 +1,7 @@
 #include <hydrogen/blas/gpu/Scale.hpp>
 
 #include <El/hydrogen_config.h>
+#include <hydrogen/meta/TypeTraits.hpp>
 #include <hydrogen/device/gpu/CUDA.hpp>
 
 #include <cuda_runtime.h>
@@ -66,7 +67,7 @@ void Scale_GPU_impl(
     {
         return;
     }
-    
+
     constexpr int TILE_DIM = 32;
     constexpr int BLK_COLS = 8;
 
diff --git a/src/hydrogen/blas/gpu/Transpose.cu b/src/hydrogen/blas/gpu/Transpose.cu
index ad3304ae05..e48bf5497d 100644
--- a/src/hydrogen/blas/gpu/Transpose.cu
+++ b/src/hydrogen/blas/gpu/Transpose.cu
@@ -1,6 +1,7 @@
 #include <hydrogen/blas/gpu/Transpose.hpp>
 
 #include <El/hydrogen_config.h>
+#include <hydrogen/meta/TypeTraits.hpp>
 #include <hydrogen/device/gpu/CUDA.hpp>
 
 #include <cuda_runtime.h>
diff --git a/src/io/Print.cpp b/src/io/Print.cpp
index 3aad9276f6..60e8d7a8de 100644
--- a/src/io/Print.cpp
+++ b/src/io/Print.cpp
@@ -115,6 +115,8 @@ void Print( const vector<T>& x, string title, ostream& os )
   template void Print \
   ( const Matrix<T>& A, string title, ostream& os ); \
   template void Print \
+  ( const AbstractMatrix<T>& A, string title, ostream& os ); \
+  template void Print \
   ( const AbstractDistMatrix<T>& A, string title, ostream& os );
 
 #ifdef HYDROGEN_GPU_USE_FP16
diff --git a/src/io/Read.cpp b/src/io/Read.cpp
index d19da92870..8762f6d1b9 100644
--- a/src/io/Read.cpp
+++ b/src/io/Read.cpp
@@ -125,6 +125,20 @@ void Read
     }
 }
 
+#ifdef HYDROGEN_GPU_USE_FP16
+template <>
+void Read<gpu_half_type>(
+    AbstractDistMatrix<gpu_half_type>& A,
+    const string filename, FileFormat format,
+    bool sequential)
+{
+    std::unique_ptr<AbstractDistMatrix<float>> A_tmp(
+        AbstractDistMatrix<float>::Instantiate(A.DistData()));
+    Read(*A_tmp, filename, format, sequential);
+    Copy(*A_tmp, A);
+}
+#endif // HYDROGEN_GPU_USE_FP16
+
 #define PROTO(T) \
   template void Read \
   ( Matrix<T>& A, const string filename, FileFormat format ); \
diff --git a/src/io/Write.cpp b/src/io/Write.cpp
index 5549c13ec6..7399e25746 100644
--- a/src/io/Write.cpp
+++ b/src/io/Write.cpp
@@ -86,6 +86,29 @@ void Write
     }
 }
 
+#ifdef HYDROGEN_GPU_USE_FP16
+template <>
+void Write<gpu_half_type>(AbstractMatrix<gpu_half_type> const& A,
+                          string basename,
+                          FileFormat format, string title)
+{
+    Matrix<float> A_tmp;
+    Copy(A, A_tmp);
+    return Write(A_tmp, basename, format, title);
+}
+
+template <>
+void Write<gpu_half_type>(
+    const AbstractDistMatrix<gpu_half_type>& A,
+    string basename, FileFormat format, string title)
+{
+    std::unique_ptr<AbstractDistMatrix<float>> A_tmp(
+        AbstractDistMatrix<float>::Instantiate(A.DistData()));
+    Copy(A, *A_tmp);
+    Write(*A_tmp, basename, format, title);
+}
+#endif // HYDROGEN_GPU_USE_FP16
+
 #define PROTO(T) \
     template void Write                         \
     ( const AbstractMatrix<T>&, string, FileFormat, string );   \
diff --git a/src/matrices/deterministic/classical/Ones.cpp b/src/matrices/deterministic/classical/Ones.cpp
index 821ae2518c..3c527887fc 100644
--- a/src/matrices/deterministic/classical/Ones.cpp
+++ b/src/matrices/deterministic/classical/Ones.cpp
@@ -17,7 +17,7 @@ void Ones( Matrix<T>& A, Int m, Int n )
 {
     EL_DEBUG_CSE
     A.Resize( m, n );
-    Fill( A, T(1) );
+    Fill( A, TypeTraits<T>::One() );
 }
 
 template<typename T>
@@ -25,7 +25,7 @@ void Ones( AbstractMatrix<T>& A, Int m, Int n )
 {
     EL_DEBUG_CSE
     A.Resize( m, n );
-    Fill( A, T(1) );
+    Fill( A, TypeTraits<T>::One() );
 }
 
 template<typename T>
@@ -33,7 +33,7 @@ void Ones( AbstractDistMatrix<T>& A, Int m, Int n )
 {
     EL_DEBUG_CSE
     A.Resize( m, n );
-    Fill( A, T(1) );
+    Fill( A, TypeTraits<T>::One() );
 }
 
 
@@ -42,6 +42,10 @@ void Ones( AbstractDistMatrix<T>& A, Int m, Int n )
   template void Ones( AbstractMatrix<T>& A, Int m, Int n ); \
   template void Ones( AbstractDistMatrix<T>& A, Int m, Int n );
 
+#ifdef HYDROGEN_GPU_USE_FP16
+PROTO(gpu_half_type)
+#endif
+
 #define EL_ENABLE_DOUBLEDOUBLE
 #define EL_ENABLE_QUADDOUBLE
 #define EL_ENABLE_QUAD
diff --git a/src/matrices/deterministic/classical/Zeros.cpp b/src/matrices/deterministic/classical/Zeros.cpp
index b2e738c22d..978532caf3 100644
--- a/src/matrices/deterministic/classical/Zeros.cpp
+++ b/src/matrices/deterministic/classical/Zeros.cpp
@@ -69,6 +69,10 @@ void Zeros( AbstractDistMatrix<T>& A, Int m, Int n )
   template void Zeros( AbstractMatrix<T>& A, Int m, Int n ); \
   template void Zeros( AbstractDistMatrix<T>& A, Int m, Int n );
 
+#ifdef HYDROGEN_GPU_USE_FP16
+PROTO(gpu_half_type)
+#endif
+
 #define EL_ENABLE_DOUBLEDOUBLE
 #define EL_ENABLE_QUADDOUBLE
 #define EL_ENABLE_QUAD
diff --git a/src/matrices/random/independent/Bernoulli.cpp b/src/matrices/random/independent/Bernoulli.cpp
index c25b3e90ad..a45af4d6cb 100644
--- a/src/matrices/random/independent/Bernoulli.cpp
+++ b/src/matrices/random/independent/Bernoulli.cpp
@@ -2,8 +2,8 @@
    Copyright (c) 2009-2016, Jack Poulson
    All rights reserved.
 
-   This file is part of Elemental and is under the BSD 2-Clause License, 
-   which can be found in the LICENSE file in the root directory, or at 
+   This file is part of Elemental and is under the BSD 2-Clause License,
+   which can be found in the LICENSE file in the root directory, or at
    http://opensource.org/licenses/BSD-2-Clause
 */
 #include <El-lite.hpp>
@@ -14,7 +14,7 @@ namespace El {
 
 template<typename T>
 void Bernoulli( Matrix<T>& A, Int m, Int n, double p )
-{ 
+{
     EL_DEBUG_CSE
     if( p < 0. || p > 1. )
         LogicError
@@ -24,8 +24,8 @@ void Bernoulli( Matrix<T>& A, Int m, Int n, double p )
     auto doubleCoin = [=]() -> T
     {
         const double alpha = SampleUniform<double>(0,1);
-        if( alpha <= q ) return T(0); 
-        else             return T(1);
+        if( alpha <= q ) return TypeTraits<T>::Zero();
+        else             return TypeTraits<T>::One();
     };
     EntrywiseFill( A, function<T()>(doubleCoin) );
 }
@@ -42,8 +42,8 @@ void Bernoulli( AbstractDistMatrix<T>& A, Int m, Int n, double p )
     auto doubleCoin = [=]() -> T
     {
         const double alpha = SampleUniform<double>(0,1);
-        if( alpha <= q ) return T(0); 
-        else             return T(1);
+        if( alpha <= q ) return TypeTraits<T>::Zero();
+        else             return TypeTraits<T>::One();
     };
     EntrywiseFill( A, function<T()>(doubleCoin) );
 }
@@ -52,6 +52,10 @@ void Bernoulli( AbstractDistMatrix<T>& A, Int m, Int n, double p )
   template void Bernoulli( Matrix<T>& A, Int m, Int n, double p ); \
   template void Bernoulli( AbstractDistMatrix<T>& A, Int m, Int n, double p );
 
+#ifdef HYDROGEN_GPU_USE_FP16
+PROTO(gpu_half_type)
+#endif
+
 #define EL_ENABLE_DOUBLEDOUBLE
 #define EL_ENABLE_QUADDOUBLE
 #define EL_ENABLE_QUAD
diff --git a/src/matrices/random/independent/Gaussian.cpp b/src/matrices/random/independent/Gaussian.cpp
index 8c4bd864c2..6584744108 100644
--- a/src/matrices/random/independent/Gaussian.cpp
+++ b/src/matrices/random/independent/Gaussian.cpp
@@ -83,6 +83,10 @@ void Gaussian
   template void Gaussian \
   ( AbstractDistMatrix<F>& A, Int m, Int n, F mean, Base<F> stddev );
 
+#ifdef HYDROGEN_GPU_USE_FP16
+PROTO(gpu_half_type)
+#endif
+
 #define EL_NO_INT_PROTO
 #define EL_ENABLE_DOUBLEDOUBLE
 #define EL_ENABLE_QUADDOUBLE
diff --git a/tests/blas_like/Gemm.cpp b/tests/blas_like/Gemm.cpp
index 2025d46e86..5b74dcad44 100644
--- a/tests/blas_like/Gemm.cpp
+++ b/tests/blas_like/Gemm.cpp
@@ -118,6 +118,14 @@ void TestGemm
     cudaEventCreate(&start);
     cudaEventCreate(&stop);
     float cudaTime;
+
+    // Warmup run -- doesn't matter in CPU land
+    if (D == Device::GPU)
+    {
+        C = COrig;
+        Gemm(orientA, orientB, alpha, A, B, beta, C, GEMM_SUMMA_A);
+        mpi::Barrier(g.Comm());
+    }
 #endif
 
     // Test the variant of Gemm that keeps A stationary
diff --git a/unit_test/CMakeLists.txt b/unit_test/CMakeLists.txt
index 1cc6048a00..0318babe07 100644
--- a/unit_test/CMakeLists.txt
+++ b/unit_test/CMakeLists.txt
@@ -1,5 +1,16 @@
 list(APPEND HYDROGEN_CATCH2_TEST_FILES
-  matrix_test.cpp)
+  matrix_test.cpp
+  )
+if (HYDROGEN_HAVE_GPU)
+  list(APPEND HYDROGEN_CATCH2_TEST_FILES
+    copy_cpu_gpu_test.cpp
+    )
+  if (HYDROGEN_HAVE_HALF AND HYDROGEN_GPU_USE_FP16)
+    list(APPEND HYDROGEN_CATCH2_TEST_FILES
+      copy_cpu_gpu_half_test.cpp
+      )
+  endif ()
+endif (HYDROGEN_HAVE_GPU)
 
 # Add the sequential test main() function
 add_executable(seq-catch-tests
diff --git a/unit_test/copy_cpu_gpu_half_test.cpp b/unit_test/copy_cpu_gpu_half_test.cpp
new file mode 100644
index 0000000000..84755d9217
--- /dev/null
+++ b/unit_test/copy_cpu_gpu_half_test.cpp
@@ -0,0 +1,82 @@
+#include <catch2/catch.hpp>
+
+// Other includes
+#include <El.hpp>
+#include <hydrogen/meta/TypeTraits.hpp>
+
+TEST_CASE(
+  "Testing Interdevice Matrix Copy -- half types",
+  "[seq][matrix][blas][copy][half]")
+{
+    using cpu_matrix_type = El::Matrix<El::cpu_half_type, El::Device::CPU>;
+    using gpu_matrix_type = El::Matrix<El::gpu_half_type, El::Device::GPU>;
+    using size_type = typename cpu_matrix_type::size_type;
+
+    auto const zero_size = hydrogen::TypeTraits<size_type>::Zero();
+
+    SECTION ("Copy host to device.")
+    {
+        cpu_matrix_type cpu_src(7, 11);
+        gpu_matrix_type gpu_tgt;
+
+        CHECK(cpu_src.Height() == 7);
+        CHECK(cpu_src.Width() == 11);
+        CHECK(gpu_tgt.Height() == zero_size);
+        CHECK(gpu_tgt.Width() == zero_size);
+
+        // Set some values
+        MakeUniform(cpu_src, El::cpu_half_type{0.f}, El::cpu_half_type{1.f});
+
+        // Do the copy
+        Copy(cpu_src, gpu_tgt);
+
+        CHECK(gpu_tgt.Height() == cpu_src.Height());
+        CHECK(gpu_tgt.Width() == cpu_src.Width());
+
+        // This takes a different codepath, for better or worse. In
+        // this case, it is SUPER hacky.
+        El::Matrix<El::gpu_half_type, El::Device::CPU> cpu_tgt_tmp(gpu_tgt);
+        cpu_matrix_type cpu_tgt(
+            cpu_tgt_tmp.Height(), cpu_tgt_tmp.Width(),
+            reinterpret_cast<El::cpu_half_type*>(cpu_tgt_tmp.Buffer()),
+            cpu_tgt_tmp.LDim());
+        for (size_type col = zero_size; col < cpu_src.Width(); ++col)
+            for (size_type row = zero_size; row < cpu_src.Height(); ++row)
+            {
+                INFO("Row = " << row << "; Col = " << col);
+                CHECK(cpu_src(row, col) == cpu_tgt(row, col));
+            }
+    }
+
+    SECTION ("Copy device to host.")
+    {
+        gpu_matrix_type gpu_src(13, 11, 17);
+        cpu_matrix_type cpu_tgt;
+
+        CHECK(gpu_src.Height() == 13);
+        CHECK(gpu_src.Width() == 11);
+        CHECK(cpu_tgt.Height() == zero_size);
+        CHECK(cpu_tgt.Width() == zero_size);
+
+        // Set some values
+        MakeUniform(gpu_src, El::gpu_half_type{0.f}, El::gpu_half_type{1.f});
+
+        // Do the copy
+        Copy(gpu_src, cpu_tgt);
+
+        CHECK(cpu_tgt.Height() == gpu_src.Height());
+        CHECK(cpu_tgt.Width() == gpu_src.Width());
+
+        El::Matrix<El::gpu_half_type, El::Device::CPU> cpu_src_tmp(gpu_src);
+        cpu_matrix_type cpu_src(
+            cpu_src_tmp.Height(), cpu_src_tmp.Width(),
+            reinterpret_cast<El::cpu_half_type*>(cpu_src_tmp.Buffer()),
+            cpu_src_tmp.LDim());
+        for (size_type col = zero_size; col < cpu_src.Width(); ++col)
+            for (size_type row = zero_size; row < cpu_src.Height(); ++row)
+            {
+                INFO("Row = " << row << "; Col = " << col);
+                CHECK(cpu_src(row, col) == cpu_tgt(row, col));
+            }
+    }
+}
diff --git a/unit_test/copy_cpu_gpu_test.cpp b/unit_test/copy_cpu_gpu_test.cpp
new file mode 100644
index 0000000000..17a4753f9d
--- /dev/null
+++ b/unit_test/copy_cpu_gpu_test.cpp
@@ -0,0 +1,75 @@
+#include <catch2/catch.hpp>
+
+// Other includes
+#include <El.hpp>
+#include <hydrogen/meta/TypeTraits.hpp>
+
+TEMPLATE_TEST_CASE(
+  "Testing Interdevice Matrix Copy",
+  "[seq][matrix][blas][copy]",
+  float, double)
+{
+    using value_type = TestType;
+    using cpu_matrix_type = El::Matrix<value_type, El::Device::CPU>;
+    using gpu_matrix_type = El::Matrix<value_type, El::Device::GPU>;
+    using size_type = typename cpu_matrix_type::size_type;
+
+    auto const zero_size = hydrogen::TypeTraits<size_type>::Zero();
+
+    SECTION ("Copy host to device.")
+    {
+        cpu_matrix_type cpu_src(7, 11);
+        gpu_matrix_type gpu_tgt;
+
+        CHECK(cpu_src.Height() == 7);
+        CHECK(cpu_src.Width() == 11);
+        CHECK(gpu_tgt.Height() == zero_size);
+        CHECK(gpu_tgt.Width() == zero_size);
+
+        // Set some values
+        MakeUniform(cpu_src);
+
+        // Do the copy
+        Copy(cpu_src, gpu_tgt);
+
+        CHECK(gpu_tgt.Height() == cpu_src.Height());
+        CHECK(gpu_tgt.Width() == cpu_src.Width());
+
+        // This takes a different codepath, for better or worse
+        cpu_matrix_type cpu_tgt(gpu_tgt);
+        for (size_type col = zero_size; col < cpu_src.Width(); ++col)
+            for (size_type row = zero_size; row < cpu_src.Height(); ++row)
+            {
+                INFO("Row = " << row << "; Col = " << col);
+                CHECK(cpu_src(row, col) == cpu_tgt(row, col));
+            }
+    }
+
+    SECTION ("Copy device to host.")
+    {
+        gpu_matrix_type gpu_src(13, 11, 17);
+        cpu_matrix_type cpu_tgt;
+
+        CHECK(gpu_src.Height() == 13);
+        CHECK(gpu_src.Width() == 11);
+        CHECK(cpu_tgt.Height() == zero_size);
+        CHECK(cpu_tgt.Width() == zero_size);
+
+        // Set some values
+        MakeUniform(gpu_src);
+
+        // Do the copy
+        Copy(gpu_src, cpu_tgt);
+
+        CHECK(cpu_tgt.Height() == gpu_src.Height());
+        CHECK(cpu_tgt.Width() == gpu_src.Width());
+
+        cpu_matrix_type cpu_src(gpu_src);
+        for (size_type col = zero_size; col < cpu_src.Width(); ++col)
+            for (size_type row = zero_size; row < cpu_src.Height(); ++row)
+            {
+                INFO("Row = " << row << "; Col = " << col);
+                CHECK(cpu_src(row, col) == cpu_tgt(row, col));
+            }
+    }
+}