From 90b64089722389b3ca0d76eabb0f4ab2dd56a059 Mon Sep 17 00:00:00 2001
From: Nikoli Dryden <dryden1@llnl.gov>
Date: Fri, 24 Mar 2017 07:24:05 -0700
Subject: [PATCH 1/4] Update Hadamard, Dot, ColumnTwoNorms, ColumnMaxNorms, and
 HilbertSchmidt to support block matrices, and add tests for Hadamard, Dot,
 ColumnTwoNorms, and ColumnMaxNorms.

---
 include/El/blas_like/level1/Dot.hpp      |  4 +-
 include/El/blas_like/level1/Hadamard.hpp | 12 +--
 include/El/blas_like/level1/decl.hpp     | 18 ++---
 src/blas_like/level1/ColumnNorms.cpp     | 22 ++++--
 src/blas_like/level1/HilbertSchmidt.cpp  |  4 +-
 tests/blas_like/ColumnNorms.cpp          | 97 ++++++++++++++++++++++++
 tests/blas_like/Dot.cpp                  | 63 +++++++++++++++
 tests/blas_like/Hadamard.cpp             | 64 ++++++++++++++++
 8 files changed, 257 insertions(+), 27 deletions(-)
 create mode 100644 tests/blas_like/ColumnNorms.cpp
 create mode 100644 tests/blas_like/Dot.cpp
 create mode 100644 tests/blas_like/Hadamard.cpp

diff --git a/include/El/blas_like/level1/Dot.hpp b/include/El/blas_like/level1/Dot.hpp
index beda613785..ce1dc6fdde 100644
--- a/include/El/blas_like/level1/Dot.hpp
+++ b/include/El/blas_like/level1/Dot.hpp
@@ -19,7 +19,7 @@ T Dot( const Matrix<T>& A, const Matrix<T>& B )
 }
 
 template<typename T>
-T Dot( const ElementalMatrix<T>& A, const ElementalMatrix<T>& B )
+T Dot( const AbstractDistMatrix<T>& A, const AbstractDistMatrix<T>& B )
 {
     EL_DEBUG_CSE
     return HilbertSchmidt( A, B );
@@ -114,7 +114,7 @@ T Dotu( const DistMultiVec<T>& A, const DistMultiVec<T>& B )
   EL_EXTERN template T Dot \
   ( const Matrix<T>& A, const Matrix<T>& B ); \
   EL_EXTERN template T Dot \
-  ( const ElementalMatrix<T>& A, const ElementalMatrix<T>& B ); \
+  ( const AbstractDistMatrix<T>& A, const AbstractDistMatrix<T>& B ); \
   EL_EXTERN template T Dot \
   ( const DistMultiVec<T>& A, const DistMultiVec<T>& B ); \
   EL_EXTERN template T Dotu \
diff --git a/include/El/blas_like/level1/Hadamard.hpp b/include/El/blas_like/level1/Hadamard.hpp
index e1a5b062d2..4d741b1298 100644
--- a/include/El/blas_like/level1/Hadamard.hpp
+++ b/include/El/blas_like/level1/Hadamard.hpp
@@ -42,9 +42,9 @@ void Hadamard( const Matrix<T>& A, const Matrix<T>& B, Matrix<T>& C )
 
 template<typename T>
 void Hadamard
-( const ElementalMatrix<T>& A,
-  const ElementalMatrix<T>& B,
-        ElementalMatrix<T>& C )
+( const AbstractDistMatrix<T>& A,
+  const AbstractDistMatrix<T>& B,
+        AbstractDistMatrix<T>& C )
 {
     EL_DEBUG_CSE
     const DistData& ADistData = A.DistData();
@@ -87,9 +87,9 @@ void Hadamard
   EL_EXTERN template void Hadamard \
   ( const Matrix<T>& A, const Matrix<T>& B, Matrix<T>& C ); \
   EL_EXTERN template void Hadamard \
-  ( const ElementalMatrix<T>& A, \
-    const ElementalMatrix<T>& B, \
-          ElementalMatrix<T>& C ); \
+  ( const AbstractDistMatrix<T>& A, \
+    const AbstractDistMatrix<T>& B, \
+          AbstractDistMatrix<T>& C ); \
   EL_EXTERN template void Hadamard \
   ( const DistMultiVec<T>& A, \
     const DistMultiVec<T>& B, \
diff --git a/include/El/blas_like/level1/decl.hpp b/include/El/blas_like/level1/decl.hpp
index 48cd5f269e..e129c0968b 100644
--- a/include/El/blas_like/level1/decl.hpp
+++ b/include/El/blas_like/level1/decl.hpp
@@ -196,9 +196,9 @@ void Recv( Matrix<T>& A, mpi::Comm comm, int source );
 template<typename Field>
 void ColumnTwoNorms
 ( const Matrix<Field>& X, Matrix<Base<Field>>& norms );
-template<typename Field,Dist U,Dist V>
+template<typename Field,Dist U,Dist V,DistWrap W>
 void ColumnTwoNorms
-( const DistMatrix<Field,U,V>& X, DistMatrix<Base<Field>,V,STAR>& norms );
+( const DistMatrix<Field,U,V,W>& X, DistMatrix<Base<Field>,V,STAR,W>& norms );
 template<typename Field>
 void ColumnTwoNorms
 ( const DistMultiVec<Field>& X, Matrix<Base<Field>>& norms );
@@ -235,9 +235,9 @@ void ColumnTwoNorms
 template<typename Ring>
 void ColumnMaxNorms
 ( const Matrix<Ring>& X, Matrix<Base<Ring>>& norms );
-template<typename Ring,Dist U,Dist V>
+template<typename Ring,Dist U,Dist V,DistWrap W>
 void ColumnMaxNorms
-( const DistMatrix<Ring,U,V>& X, DistMatrix<Base<Ring>,V,STAR>& norms );
+( const DistMatrix<Ring,U,V,W>& X, DistMatrix<Base<Ring>,V,STAR,W>& norms );
 template<typename Ring>
 void ColumnMaxNorms
 ( const DistMultiVec<Ring>& X, Matrix<Base<Ring>>& norms );
@@ -802,7 +802,7 @@ void DiagonalSolve
 template<typename T>
 T Dot( const Matrix<T>& A, const Matrix<T>& B );
 template<typename T>
-T Dot( const ElementalMatrix<T>& A, const ElementalMatrix<T>& B );
+T Dot( const AbstractDistMatrix<T>& A, const AbstractDistMatrix<T>& B );
 template<typename T>
 T Dot( const DistMultiVec<T>& A, const DistMultiVec<T>& B );
 
@@ -1155,9 +1155,9 @@ template<typename T>
 void Hadamard( const Matrix<T>& A, const Matrix<T>& B, Matrix<T>& C );
 template<typename T>
 void Hadamard
-( const ElementalMatrix<T>& A,
-  const ElementalMatrix<T>& B,
-        ElementalMatrix<T>& C );
+( const AbstractDistMatrix<T>& A,
+  const AbstractDistMatrix<T>& B,
+        AbstractDistMatrix<T>& C );
 template<typename T>
 void Hadamard
 ( const DistMultiVec<T>& A,
@@ -1170,7 +1170,7 @@ template<typename T>
 T HilbertSchmidt( const Matrix<T>& A, const Matrix<T>& B );
 template<typename T>
 T HilbertSchmidt
-( const ElementalMatrix<T>& A, const ElementalMatrix<T>& C );
+( const AbstractDistMatrix<T>& A, const AbstractDistMatrix<T>& C );
 template<typename T>
 T HilbertSchmidt( const DistMultiVec<T>& A, const DistMultiVec<T>& B );
 
diff --git a/src/blas_like/level1/ColumnNorms.cpp b/src/blas_like/level1/ColumnNorms.cpp
index 05bfc21b6a..9ac43012f6 100644
--- a/src/blas_like/level1/ColumnNorms.cpp
+++ b/src/blas_like/level1/ColumnNorms.cpp
@@ -103,9 +103,9 @@ void ColumnMaxNorms( const Matrix<Field>& X, Matrix<Base<Field>>& norms )
     }
 }
 
-template<typename Field,Dist U,Dist V>
+template<typename Field,Dist U,Dist V,DistWrap W>
 void ColumnTwoNorms
-( const DistMatrix<Field,U,V>& A, DistMatrix<Base<Field>,V,STAR>& norms )
+( const DistMatrix<Field,U,V,W>& A, DistMatrix<Base<Field>,V,STAR,W>& norms )
 {
     EL_DEBUG_CSE
     norms.AlignWith( A );
@@ -118,9 +118,9 @@ void ColumnTwoNorms
     ColumnTwoNormsHelper( A.LockedMatrix(), norms.Matrix(), A.ColComm() );
 }
 
-template<typename Field,Dist U,Dist V>
+template<typename Field,Dist U,Dist V,DistWrap W>
 void ColumnMaxNorms
-( const DistMatrix<Field,U,V>& A, DistMatrix<Base<Field>,V,STAR>& norms )
+( const DistMatrix<Field,U,V,W>& A, DistMatrix<Base<Field>,V,STAR,W>& norms )
 {
     EL_DEBUG_CSE
     norms.AlignWith( A );
@@ -363,11 +363,17 @@ void ColumnTwoNorms
 
 #define PROTO_DIST(Field,U,V) \
   template void ColumnTwoNorms \
-  ( const DistMatrix<Field,U,V>& X, \
-          DistMatrix<Base<Field>,V,STAR>& norms ); \
+  ( const DistMatrix<Field,U,V,ELEMENT>& X, \
+          DistMatrix<Base<Field>,V,STAR,ELEMENT>& norms ); \
   template void ColumnMaxNorms \
-  ( const DistMatrix<Field,U,V>& X, \
-          DistMatrix<Base<Field>,V,STAR>& norms );
+  ( const DistMatrix<Field,U,V,ELEMENT>& X, \
+          DistMatrix<Base<Field>,V,STAR,ELEMENT>& norms ); \
+  template void ColumnTwoNorms \
+  ( const DistMatrix<Field,U,V,BLOCK>& X, \
+          DistMatrix<Base<Field>,V,STAR,BLOCK>& norms ); \
+  template void ColumnMaxNorms \
+  ( const DistMatrix<Field,U,V,BLOCK>& X, \
+          DistMatrix<Base<Field>,V,STAR,BLOCK>& norms );
 
 #define PROTO(Field) \
   template void ColumnTwoNorms \
diff --git a/src/blas_like/level1/HilbertSchmidt.cpp b/src/blas_like/level1/HilbertSchmidt.cpp
index 65dd1e7be3..c38f19556b 100644
--- a/src/blas_like/level1/HilbertSchmidt.cpp
+++ b/src/blas_like/level1/HilbertSchmidt.cpp
@@ -41,7 +41,7 @@ Ring HilbertSchmidt( const Matrix<Ring>& A, const Matrix<Ring>& B )
 
 template<typename Ring>
 Ring HilbertSchmidt
-( const ElementalMatrix<Ring>& A, const ElementalMatrix<Ring>& B )
+( const AbstractDistMatrix<Ring>& A, const AbstractDistMatrix<Ring>& B )
 {
     EL_DEBUG_CSE
     if( A.Height() != B.Height() || A.Width() != B.Width() )
@@ -112,7 +112,7 @@ Ring HilbertSchmidt( const DistMultiVec<Ring>& A, const DistMultiVec<Ring>& B )
   template Ring HilbertSchmidt \
   ( const Matrix<Ring>& A, const Matrix<Ring>& B ); \
   template Ring HilbertSchmidt \
-  ( const ElementalMatrix<Ring>& A, const ElementalMatrix<Ring>& B ); \
+  ( const AbstractDistMatrix<Ring>& A, const AbstractDistMatrix<Ring>& B ); \
   template Ring HilbertSchmidt \
   ( const DistMultiVec<Ring>& A, const DistMultiVec<Ring>& B );
 
diff --git a/tests/blas_like/ColumnNorms.cpp b/tests/blas_like/ColumnNorms.cpp
new file mode 100644
index 0000000000..ef7f3d8251
--- /dev/null
+++ b/tests/blas_like/ColumnNorms.cpp
@@ -0,0 +1,97 @@
+/*
+   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
+   http://opensource.org/licenses/BSD-2-Clause
+*/
+
+#include <El.hpp>
+using namespace El;
+
+template <typename T, DistWrap W>
+void TestColumnTwoNorms(Int m, Int n, const Grid& g, bool print) {
+  // Generate random matrix to test.
+  DistMatrix<T, MC, MR, W> A(g);
+  Uniform(A, m, n);
+  if (print) {
+    Print(A, "A");
+  }
+  DistMatrix<T, MR, STAR, W> norms(g);
+  ColumnTwoNorms(A, norms);
+  if (print) {
+    Print(norms, "norms");
+  }
+  for (Int j = 0; j < A.LocalWidth(); ++j) {
+    T got = norms.GetLocal(j, 0);
+    T expected = 0;
+    for (Int i = 0; i < A.LocalHeight(); ++i) {
+      T val = A.GetLocal(i, j);
+      expected += val * val;
+    }
+    expected = mpi::AllReduce(expected, g.ColComm());
+    expected = Sqrt(expected);
+    if (Abs(got - expected) > 1e-5) {
+      Output("Results do not match, norms(", j, ")=", got,
+             " instead of ", expected);
+      RuntimeError("got != expected");
+    }
+  }
+}
+
+template <typename T, DistWrap W>
+void TestColumnMaxNorms(Int m, Int n, const Grid& g, bool print) {
+  // Generate random matrix to test.
+  DistMatrix<T, MC, MR, W> A(g);
+  Uniform(A, m, n);
+  if (print) {
+    Print(A, "A");
+  }
+  DistMatrix<T, MR, STAR, W> norms(g);
+  ColumnMaxNorms(A, norms);
+  if (print) {
+    Print(norms, "norms");
+  }
+  for (Int j = 0; j < A.LocalWidth(); ++j) {
+    T got = norms.GetLocal(j, 0);
+    T expected = 0;
+    for (Int i = 0; i < A.LocalHeight(); ++i) {
+      expected = Max(expected, Abs(A.GetLocal(i, j)));
+    }
+    T r;
+    mpi::AllReduce(&expected, &r, 1, mpi::MAX, g.ColComm());
+    expected = r;
+    if (got != expected) {
+      Output("Results do not match, norms(", j, ")=", got,
+             " instead of ", expected);
+      RuntimeError("got != expected");
+    }
+  }
+}
+
+int main(int argc, char** argv) {
+  Environment env(argc, argv);
+  mpi::Comm comm = mpi::COMM_WORLD;
+  try {
+    const Int m = Input("--m", "height", 100);
+    const Int n = Input("--n", "width", 100);
+    const bool print = Input("--print", "print matrices?", false);
+    ProcessInput();
+    PrintInputReport();
+
+    const Grid g(comm);
+    OutputFromRoot(comm, "Testing ColumnTwoNorms");
+    TestColumnTwoNorms<float, ELEMENT>(m, n, g, print);
+    TestColumnTwoNorms<float, BLOCK>(m, n, g, print);
+    TestColumnTwoNorms<double, ELEMENT>(m, n, g, print);
+    TestColumnTwoNorms<double, BLOCK>(m, n, g, print);
+    OutputFromRoot(comm, "Testing ColumnMaxNorms");
+    TestColumnMaxNorms<float, ELEMENT>(m, n, g, print);
+    TestColumnMaxNorms<float, BLOCK>(m, n, g, print);
+    TestColumnMaxNorms<double, ELEMENT>(m, n, g, print);
+    TestColumnMaxNorms<double, BLOCK>(m, n, g, print);
+  } catch (exception& e) {
+    ReportException(e);
+  }
+}
diff --git a/tests/blas_like/Dot.cpp b/tests/blas_like/Dot.cpp
new file mode 100644
index 0000000000..2d9e3899a0
--- /dev/null
+++ b/tests/blas_like/Dot.cpp
@@ -0,0 +1,63 @@
+/*
+   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
+   http://opensource.org/licenses/BSD-2-Clause
+*/
+
+#include <El.hpp>
+using namespace El;
+
+template <typename T, DistWrap W>
+void TestDot(Int m, Int n, const Grid& g, bool print) {
+  // Generate random matrices to test.
+  DistMatrix<T, MC, MR, W> A(g);
+  Uniform(A, m, n);
+  DistMatrix<T, MC, MR, W> B(g);
+  Uniform(B, m, n);
+  if (print) {
+    Print(A, "A");
+    Print(B, "B");
+  }
+  // Do Dot.
+  T got = Dot(A, B);
+  if (print) {
+    OutputFromRoot(g.Comm(), "result=", got);
+  }
+  // Manually check results.
+  T expected = 0;
+  for (Int j = 0; j < A.LocalWidth(); ++j) {
+    for (Int i = 0; i < B.LocalHeight(); ++i) {
+      expected += A.GetLocal(i, j) * B.GetLocal(i, j);
+    }
+  }
+  expected = mpi::AllReduce(expected, g.Comm());
+  if (Abs(got - expected) > 1e-6) {
+    Output("Results do not match, got=", got,
+           " instead of ", expected);
+    RuntimeError("got != expected");
+  }
+}
+
+int main(int argc, char** argv) {
+  Environment env(argc, argv);
+  mpi::Comm comm = mpi::COMM_WORLD;
+  try {
+    const Int m = Input("--m", "height", 100);
+    const Int n = Input("--n", "width", 100);
+    const bool print = Input("--print", "print matrices?", false);
+    ProcessInput();
+    PrintInputReport();
+
+    const Grid g(comm);
+    OutputFromRoot(comm, "Testing Dot");
+    TestDot<float, ELEMENT>(m, n, g, print);
+    TestDot<float, BLOCK>(m, n, g, print);
+    TestDot<double, ELEMENT>(m, n, g, print);
+    TestDot<double, BLOCK>(m, n, g, print);
+  } catch (exception& e) {
+    ReportException(e);
+  }
+}
diff --git a/tests/blas_like/Hadamard.cpp b/tests/blas_like/Hadamard.cpp
new file mode 100644
index 0000000000..c1a1c1ceda
--- /dev/null
+++ b/tests/blas_like/Hadamard.cpp
@@ -0,0 +1,64 @@
+/*
+   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
+   http://opensource.org/licenses/BSD-2-Clause
+*/
+
+#include <El.hpp>
+using namespace El;
+
+template <typename T, DistWrap W>
+void TestHadamard(Int m, Int n, const Grid& g, bool print) {
+  // Generate random matrices to test.
+  DistMatrix<T, MC, MR, W> A(g);
+  Uniform(A, m, n);
+  DistMatrix<T, MC, MR, W> B(g);
+  Uniform(B, m, n);
+  if (print) {
+    Print(A, "A");
+    Print(B, "B");
+  }
+  // Apply Hadamard.
+  DistMatrix<T, MC, MR, W> C(g);
+  Hadamard(A, B, C);
+  mpi::Barrier(g.Comm());
+  if (print) {
+    Print(C, "C");
+  }
+  // Manually check results.
+  for (Int j = 0; j < A.LocalWidth(); ++j) {
+    for (Int i = 0; i < A.LocalHeight(); ++i) {
+      T got = C.GetLocal(i, j);
+      T expected = A.GetLocal(i, j) * B.GetLocal(i, j);
+      if (Abs(got - expected) > 1e-6) {
+        Output("Results do not match, C(", i, ",", j, ")=", got,
+               " instead of ", expected);
+        RuntimeError("got != expected");
+      }
+    }
+  }
+}
+
+int main(int argc, char** argv) {
+  Environment env(argc, argv);
+  mpi::Comm comm = mpi::COMM_WORLD;
+  try {
+    const Int m = Input("--m", "height", 100);
+    const Int n = Input("--n", "width", 100);
+    const bool print = Input("--print", "print matrices?", false);
+    ProcessInput();
+    PrintInputReport();
+
+    const Grid g(comm);
+    OutputFromRoot(comm, "Testing Hadamard");
+    TestHadamard<float, ELEMENT>(m, n, g, print);
+    TestHadamard<float, BLOCK>(m, n, g, print);
+    TestHadamard<double, ELEMENT>(m, n, g, print);
+    TestHadamard<double, BLOCK>(m, n, g, print);
+  } catch (exception& e) {
+    ReportException(e);
+  }
+}

From 4afe36a110f7adac87d7971d0f95b8952a144190 Mon Sep 17 00:00:00 2001
From: Nikoli Dryden <dryden1@llnl.gov>
Date: Fri, 24 Mar 2017 13:16:03 -0700
Subject: [PATCH 2/4] Add check for whether the same block size is used in
 Hadamard and HilbertSchmidt.

---
 include/El/blas_like/level1/Hadamard.hpp | 3 +++
 src/blas_like/level1/HilbertSchmidt.cpp  | 3 +++
 2 files changed, 6 insertions(+)

diff --git a/include/El/blas_like/level1/Hadamard.hpp b/include/El/blas_like/level1/Hadamard.hpp
index 4d741b1298..89ff1c52b7 100644
--- a/include/El/blas_like/level1/Hadamard.hpp
+++ b/include/El/blas_like/level1/Hadamard.hpp
@@ -60,6 +60,9 @@ void Hadamard
         LogicError("A, B, and C must share the same distribution");
     if( A.ColAlign() != B.ColAlign() || A.RowAlign() != B.RowAlign() )
         LogicError("A and B must be aligned");
+    if ( A.BlockHeight() != B.BlockHeight() ||
+         A.BlockWidth() != B.BlockWidth())
+      LogicError("A and B must have the same block size");
     C.AlignWith( A.DistData() );
     C.Resize( A.Height(), A.Width() );
     Hadamard( A.LockedMatrix(), B.LockedMatrix(), C.Matrix() );
diff --git a/src/blas_like/level1/HilbertSchmidt.cpp b/src/blas_like/level1/HilbertSchmidt.cpp
index c38f19556b..ca7f2c3685 100644
--- a/src/blas_like/level1/HilbertSchmidt.cpp
+++ b/src/blas_like/level1/HilbertSchmidt.cpp
@@ -53,6 +53,9 @@ Ring HilbertSchmidt
         LogicError("A and B must have the same distribution");
     if( A.ColAlign() != B.ColAlign() || A.RowAlign() != B.RowAlign() )
         LogicError("Matrices must be aligned");
+    if ( A.BlockHeight() != B.BlockHeight() ||
+         A.BlockWidth() != B.BlockWidth())
+      LogicError("A and B must have the same block size");
 
     Ring innerProd;
     if( A.Participating() )

From be1c42a2b648dcbc88471ef278a82045bb09bd69 Mon Sep 17 00:00:00 2001
From: Nikoli Dryden <dryden1@llnl.gov>
Date: Sun, 26 Mar 2017 18:25:02 -0700
Subject: [PATCH 3/4] Test Hadamard, Dot, ColumnTwoNorms, and ColumnMaxNorms at
 higher precision and in complex variants. Clean up the syntax and use
 limits::Epsilon for checking error.

---
 tests/blas_like/ColumnNorms.cpp | 78 ++++++++++++++++++++++++---------
 tests/blas_like/Dot.cpp         | 60 +++++++++++++++++++------
 tests/blas_like/Hadamard.cpp    | 56 ++++++++++++++++++-----
 3 files changed, 151 insertions(+), 43 deletions(-)

diff --git a/tests/blas_like/ColumnNorms.cpp b/tests/blas_like/ColumnNorms.cpp
index ef7f3d8251..e8411add93 100644
--- a/tests/blas_like/ColumnNorms.cpp
+++ b/tests/blas_like/ColumnNorms.cpp
@@ -11,28 +11,30 @@
 using namespace El;
 
 template <typename T, DistWrap W>
-void TestColumnTwoNorms(Int m, Int n, const Grid& g, bool print) {
+void TestColumnTwoNorms(Int m, Int n, const Grid& g, bool print)
+{
   // Generate random matrix to test.
   DistMatrix<T, MC, MR, W> A(g);
   Uniform(A, m, n);
-  if (print) {
+  if (print)
     Print(A, "A");
-  }
   DistMatrix<T, MR, STAR, W> norms(g);
   ColumnTwoNorms(A, norms);
-  if (print) {
+  if (print)
     Print(norms, "norms");
-  }
-  for (Int j = 0; j < A.LocalWidth(); ++j) {
+  for (Int j = 0; j < A.LocalWidth(); ++j)
+  {
     T got = norms.GetLocal(j, 0);
     T expected = 0;
-    for (Int i = 0; i < A.LocalHeight(); ++i) {
+    for (Int i = 0; i < A.LocalHeight(); ++i)
+    {
       T val = A.GetLocal(i, j);
       expected += val * val;
     }
     expected = mpi::AllReduce(expected, g.ColComm());
     expected = Sqrt(expected);
-    if (Abs(got - expected) > 1e-5) {
+    if (Abs(got - expected) > 10 * limits::Epsilon<El::Base<T>>())
+    {
       Output("Results do not match, norms(", j, ")=", got,
              " instead of ", expected);
       RuntimeError("got != expected");
@@ -41,28 +43,28 @@ void TestColumnTwoNorms(Int m, Int n, const Grid& g, bool print) {
 }
 
 template <typename T, DistWrap W>
-void TestColumnMaxNorms(Int m, Int n, const Grid& g, bool print) {
+void TestColumnMaxNorms(Int m, Int n, const Grid& g, bool print)
+{
   // Generate random matrix to test.
   DistMatrix<T, MC, MR, W> A(g);
   Uniform(A, m, n);
-  if (print) {
+  if (print)
     Print(A, "A");
-  }
   DistMatrix<T, MR, STAR, W> norms(g);
   ColumnMaxNorms(A, norms);
-  if (print) {
+  if (print)
     Print(norms, "norms");
-  }
-  for (Int j = 0; j < A.LocalWidth(); ++j) {
+  for (Int j = 0; j < A.LocalWidth(); ++j)
+  {
     T got = norms.GetLocal(j, 0);
     T expected = 0;
-    for (Int i = 0; i < A.LocalHeight(); ++i) {
+    for (Int i = 0; i < A.LocalHeight(); ++i)
       expected = Max(expected, Abs(A.GetLocal(i, j)));
-    }
     T r;
     mpi::AllReduce(&expected, &r, 1, mpi::MAX, g.ColComm());
     expected = r;
-    if (got != expected) {
+    if (got != expected)
+    {
       Output("Results do not match, norms(", j, ")=", got,
              " instead of ", expected);
       RuntimeError("got != expected");
@@ -70,10 +72,12 @@ void TestColumnMaxNorms(Int m, Int n, const Grid& g, bool print) {
   }
 }
 
-int main(int argc, char** argv) {
+int main(int argc, char** argv)
+{
   Environment env(argc, argv);
   mpi::Comm comm = mpi::COMM_WORLD;
-  try {
+  try
+  {
     const Int m = Input("--m", "height", 100);
     const Int n = Input("--n", "width", 100);
     const bool print = Input("--print", "print matrices?", false);
@@ -86,12 +90,46 @@ int main(int argc, char** argv) {
     TestColumnTwoNorms<float, BLOCK>(m, n, g, print);
     TestColumnTwoNorms<double, ELEMENT>(m, n, g, print);
     TestColumnTwoNorms<double, BLOCK>(m, n, g, print);
+#if defined(EL_HAVE_QD) && defined(EL_ENABLE_DOUBLEDOUBLE)
+    TestColumnTwoNorms<DoubleDouble, ELEMENT>(m, n, g, print);
+    TestColumnTwoNorms<DoubleDouble, BLOCK>(m, n, g, print);
+#endif
+#if defined(EL_HAVE_QD) && defined(EL_ENABLE_QUADDOUBLE)
+    TestColumnTwoNorms<QuadDouble, ELEMENT>(m, n, g, print);
+    TestColumnTwoNorms<QuadDouble, BLOCK>(m, n, g, print);
+#endif
+#if defined(EL_HAVE_QUAD) && defined(EL_ENABLE_QUAD)
+    TestColumnTwoNorms<Quad, ELEMENT>(m, n, g, print);
+    TestColumnTwoNorms<Quad, BLOCK>(m, n, g, print);
+#endif
+#if defined(EL_HAVE_MPC) && defined(EL_ENABLE_BIGFLOAT)
+    TestColumnTwoNorms<BigFloat, ELEMENT>(m, n, g, print);
+    TestColumnTwoNorms<BigFloat, BLOCK>(m, n, g, print);
+#endif
     OutputFromRoot(comm, "Testing ColumnMaxNorms");
     TestColumnMaxNorms<float, ELEMENT>(m, n, g, print);
     TestColumnMaxNorms<float, BLOCK>(m, n, g, print);
     TestColumnMaxNorms<double, ELEMENT>(m, n, g, print);
     TestColumnMaxNorms<double, BLOCK>(m, n, g, print);
-  } catch (exception& e) {
+#if defined(EL_HAVE_QD) && defined(EL_ENABLE_DOUBLEDOUBLE)
+    TestColumnMaxNorms<DoubleDouble, ELEMENT>(m, n, g, print);
+    TestColumnMaxNorms<DoubleDouble, BLOCK>(m, n, g, print);
+#endif
+#if defined(EL_HAVE_QD) && defined(EL_ENABLE_QUADDOUBLE)
+    TestColumnMaxNorms<QuadDouble, ELEMENT>(m, n, g, print);
+    TestColumnMaxNorms<QuadDouble, BLOCK>(m, n, g, print);
+#endif
+#if defined(EL_HAVE_QUAD) && defined(EL_ENABLE_QUAD)
+    TestColumnMaxNorms<Quad, ELEMENT>(m, n, g, print);
+    TestColumnMaxNorms<Quad, BLOCK>(m, n, g, print);
+#endif
+#if defined(EL_HAVE_MPC) && defined(EL_ENABLE_BIGFLOAT)
+    TestColumnMaxNorms<BigFloat, ELEMENT>(m, n, g, print);
+    TestColumnMaxNorms<BigFloat, BLOCK>(m, n, g, print);
+#endif
+  }
+  catch (exception& e)
+  {
     ReportException(e);
   }
 }
diff --git a/tests/blas_like/Dot.cpp b/tests/blas_like/Dot.cpp
index 2d9e3899a0..ce54b53b79 100644
--- a/tests/blas_like/Dot.cpp
+++ b/tests/blas_like/Dot.cpp
@@ -11,40 +11,44 @@
 using namespace El;
 
 template <typename T, DistWrap W>
-void TestDot(Int m, Int n, const Grid& g, bool print) {
+void TestDot(Int m, Int n, const Grid& g, bool print)
+{
   // Generate random matrices to test.
   DistMatrix<T, MC, MR, W> A(g);
   Uniform(A, m, n);
   DistMatrix<T, MC, MR, W> B(g);
   Uniform(B, m, n);
-  if (print) {
+  if (print)
+  {
     Print(A, "A");
     Print(B, "B");
   }
   // Do Dot.
   T got = Dot(A, B);
-  if (print) {
+  if (print)
     OutputFromRoot(g.Comm(), "result=", got);
-  }
   // Manually check results.
   T expected = 0;
-  for (Int j = 0; j < A.LocalWidth(); ++j) {
-    for (Int i = 0; i < B.LocalHeight(); ++i) {
-      expected += A.GetLocal(i, j) * B.GetLocal(i, j);
-    }
-  }
+  for (Int j = 0; j < A.LocalWidth(); ++j)
+    for (Int i = 0; i < B.LocalHeight(); ++i)
+      expected += Conj(A.GetLocal(i, j)) * B.GetLocal(i, j);
   expected = mpi::AllReduce(expected, g.Comm());
-  if (Abs(got - expected) > 1e-6) {
+  // The constant here is large because this is not an especially stable way
+  // to compute the dot product, but it provides a dumb implementation baseline.
+  if (Abs(got - expected) > 700 * limits::Epsilon<El::Base<T>>())
+  {
     Output("Results do not match, got=", got,
            " instead of ", expected);
     RuntimeError("got != expected");
   }
 }
 
-int main(int argc, char** argv) {
+int main(int argc, char** argv)
+{
   Environment env(argc, argv);
   mpi::Comm comm = mpi::COMM_WORLD;
-  try {
+  try
+  {
     const Int m = Input("--m", "height", 100);
     const Int n = Input("--n", "width", 100);
     const bool print = Input("--print", "print matrices?", false);
@@ -55,9 +59,39 @@ int main(int argc, char** argv) {
     OutputFromRoot(comm, "Testing Dot");
     TestDot<float, ELEMENT>(m, n, g, print);
     TestDot<float, BLOCK>(m, n, g, print);
+    TestDot<Complex<float>, ELEMENT>(m, n, g, print);
+    TestDot<Complex<float>, BLOCK>(m, n, g, print);
     TestDot<double, ELEMENT>(m, n, g, print);
     TestDot<double, BLOCK>(m, n, g, print);
-  } catch (exception& e) {
+    TestDot<Complex<double>, ELEMENT>(m, n, g, print);
+    TestDot<Complex<double>, BLOCK>(m, n, g, print);
+#if defined(EL_HAVE_QD) && defined(EL_ENABLE_DOUBLEDOUBLE)
+    TestDot<DoubleDouble, ELEMENT>(m, n, g, print);
+    TestDot<DoubleDouble, BLOCK>(m, n, g, print);
+    TestDot<Complex<DoubleDouble>, ELEMENT>(m, n, g, print);
+    TestDot<Complex<DoubleDouble>, BLOCK>(m, n, g, print);
+#endif
+#if defined(EL_HAVE_QD) && defined(EL_ENABLE_QUADDOUBLE)
+    TestDot<QuadDouble, ELEMENT>(m, n, g, print);
+    TestDot<QuadDouble, BLOCK>(m, n, g, print);
+    TestDot<Complex<QuadDouble>, ELEMENT>(m, n, g, print);
+    TestDot<Complex<QuadDouble>, BLOCK>(m, n, g, print);
+#endif
+#if defined(EL_HAVE_QUAD) && defined(EL_ENABLE_QUAD)
+    TestDot<Quad, ELEMENT>(m, n, g, print);
+    TestDot<Quad, BLOCK>(m, n, g, print);
+    TestDot<Complex<Quad>, ELEMENT>(m, n, g, print);
+    TestDot<Complex<Quad>, BLOCK>(m, n, g, print);
+#endif
+#if defined(EL_HAVE_MPC) && defined(EL_ENABLE_BIGFLOAT)
+    TestDot<BigFloat, ELEMENT>(m, n, g, print);
+    TestDot<BigFloat, BLOCK>(m, n, g, print);
+    TestDot<Complex<BigFloat>, ELEMENT>(m, n, g, print);
+    TestDot<Complex<BigFloat>, BLOCK>(m, n, g, print);
+#endif
+  }
+  catch (exception& e)
+  {
     ReportException(e);
   }
 }
diff --git a/tests/blas_like/Hadamard.cpp b/tests/blas_like/Hadamard.cpp
index c1a1c1ceda..73cc5fd83d 100644
--- a/tests/blas_like/Hadamard.cpp
+++ b/tests/blas_like/Hadamard.cpp
@@ -11,13 +11,15 @@
 using namespace El;
 
 template <typename T, DistWrap W>
-void TestHadamard(Int m, Int n, const Grid& g, bool print) {
+void TestHadamard(Int m, Int n, const Grid& g, bool print)
+{
   // Generate random matrices to test.
   DistMatrix<T, MC, MR, W> A(g);
   Uniform(A, m, n);
   DistMatrix<T, MC, MR, W> B(g);
   Uniform(B, m, n);
-  if (print) {
+  if (print)
+  {
     Print(A, "A");
     Print(B, "B");
   }
@@ -25,15 +27,17 @@ void TestHadamard(Int m, Int n, const Grid& g, bool print) {
   DistMatrix<T, MC, MR, W> C(g);
   Hadamard(A, B, C);
   mpi::Barrier(g.Comm());
-  if (print) {
+  if (print)
     Print(C, "C");
-  }
   // Manually check results.
-  for (Int j = 0; j < A.LocalWidth(); ++j) {
-    for (Int i = 0; i < A.LocalHeight(); ++i) {
+  for (Int j = 0; j < A.LocalWidth(); ++j) 
+  {
+    for (Int i = 0; i < A.LocalHeight(); ++i)
+    {
       T got = C.GetLocal(i, j);
       T expected = A.GetLocal(i, j) * B.GetLocal(i, j);
-      if (Abs(got - expected) > 1e-6) {
+      if (Abs(got - expected) > 10 * limits::Epsilon<El::Base<T>>())
+      {
         Output("Results do not match, C(", i, ",", j, ")=", got,
                " instead of ", expected);
         RuntimeError("got != expected");
@@ -42,10 +46,12 @@ void TestHadamard(Int m, Int n, const Grid& g, bool print) {
   }
 }
 
-int main(int argc, char** argv) {
+int main(int argc, char** argv)
+{
   Environment env(argc, argv);
   mpi::Comm comm = mpi::COMM_WORLD;
-  try {
+  try
+  {
     const Int m = Input("--m", "height", 100);
     const Int n = Input("--n", "width", 100);
     const bool print = Input("--print", "print matrices?", false);
@@ -56,9 +62,39 @@ int main(int argc, char** argv) {
     OutputFromRoot(comm, "Testing Hadamard");
     TestHadamard<float, ELEMENT>(m, n, g, print);
     TestHadamard<float, BLOCK>(m, n, g, print);
+    TestHadamard<Complex<float>, ELEMENT>(m, n, g, print);
+    TestHadamard<Complex<float>, BLOCK>(m, n, g, print);
     TestHadamard<double, ELEMENT>(m, n, g, print);
     TestHadamard<double, BLOCK>(m, n, g, print);
-  } catch (exception& e) {
+    TestHadamard<Complex<double>, ELEMENT>(m, n, g, print);
+    TestHadamard<Complex<double>, BLOCK>(m, n, g, print);
+#if defined(EL_HAVE_QD) && defined(EL_ENABLE_DOUBLEDOUBLE)
+    TestHadamard<DoubleDouble, ELEMENT>(m, n, g, print);
+    TestHadamard<DoubleDouble, BLOCK>(m, n, g, print);
+    TestHadamard<Complex<DoubleDouble>, ELEMENT>(m, n, g, print);
+    TestHadamard<Complex<DoubleDouble>, BLOCK>(m, n, g, print);
+#endif
+#if defined(EL_HAVE_QD) && defined(EL_ENABLE_QUADDOUBLE)
+    TestHadamard<QuadDouble, ELEMENT>(m, n, g, print);
+    TestHadamard<QuadDouble, BLOCK>(m, n, g, print);
+    TestHadamard<Complex<QuadDouble>, ELEMENT>(m, n, g, print);
+    TestHadamard<Complex<QuadDouble>, BLOCK>(m, n, g, print);
+#endif
+#if defined(EL_HAVE_QUAD) && defined(EL_ENABLE_QUAD)
+    TestHadamard<Quad, ELEMENT>(m, n, g, print);
+    TestHadamard<Quad, BLOCK>(m, n, g, print);
+    TestHadamard<Complex<Quad>, ELEMENT>(m, n, g, print);
+    TestHadamard<Complex<Quad>, BLOCK>(m, n, g, print);
+#endif
+#if defined(EL_HAVE_MPC) && defined(EL_ENABLE_BIGFLOAT)
+    TestHadamard<BigFloat, ELEMENT>(m, n, g, print);
+    TestHadamard<BigFloat, BLOCK>(m, n, g, print);
+    TestHadamard<Complex<BigFloat>, ELEMENT>(m, n, g, print);
+    TestHadamard<Complex<BigFloat>, BLOCK>(m, n, g, print);
+#endif
+  }
+  catch (exception& e)
+  {
     ReportException(e);
   }
 }

From e9139f4b8a57eaef37feb4bf255f459cd36eb701 Mon Sep 17 00:00:00 2001
From: Nikoli Dryden <dryden1@llnl.gov>
Date: Mon, 27 Mar 2017 16:39:15 -0700
Subject: [PATCH 4/4] Remove unneeded macro checks and use a more rigorous
 bound in the ColumnTwoNorms test.

---
 tests/blas_like/ColumnNorms.cpp | 21 ++++++++++-----------
 tests/blas_like/Dot.cpp         |  8 +++-----
 tests/blas_like/Hadamard.cpp    |  8 +++-----
 3 files changed, 16 insertions(+), 21 deletions(-)

diff --git a/tests/blas_like/ColumnNorms.cpp b/tests/blas_like/ColumnNorms.cpp
index e8411add93..6d91fa0e97 100644
--- a/tests/blas_like/ColumnNorms.cpp
+++ b/tests/blas_like/ColumnNorms.cpp
@@ -33,7 +33,10 @@ void TestColumnTwoNorms(Int m, Int n, const Grid& g, bool print)
     }
     expected = mpi::AllReduce(expected, g.ColComm());
     expected = Sqrt(expected);
-    if (Abs(got - expected) > 10 * limits::Epsilon<El::Base<T>>())
+    // Compute max(expected, 1) to use relative bound.
+    // (std::max and El::Max don't support BigFloat.
+    T div = expected > 1 ? expected : 1;
+    if (Abs(got - expected) / div > m * n * 10 * limits::Epsilon<El::Base<T>>())
     {
       Output("Results do not match, norms(", j, ")=", got,
              " instead of ", expected);
@@ -90,19 +93,17 @@ int main(int argc, char** argv)
     TestColumnTwoNorms<float, BLOCK>(m, n, g, print);
     TestColumnTwoNorms<double, ELEMENT>(m, n, g, print);
     TestColumnTwoNorms<double, BLOCK>(m, n, g, print);
-#if defined(EL_HAVE_QD) && defined(EL_ENABLE_DOUBLEDOUBLE)
+#if defined(EL_HAVE_QD)
     TestColumnTwoNorms<DoubleDouble, ELEMENT>(m, n, g, print);
     TestColumnTwoNorms<DoubleDouble, BLOCK>(m, n, g, print);
-#endif
-#if defined(EL_HAVE_QD) && defined(EL_ENABLE_QUADDOUBLE)
     TestColumnTwoNorms<QuadDouble, ELEMENT>(m, n, g, print);
     TestColumnTwoNorms<QuadDouble, BLOCK>(m, n, g, print);
 #endif
-#if defined(EL_HAVE_QUAD) && defined(EL_ENABLE_QUAD)
+#if defined(EL_HAVE_QUAD)
     TestColumnTwoNorms<Quad, ELEMENT>(m, n, g, print);
     TestColumnTwoNorms<Quad, BLOCK>(m, n, g, print);
 #endif
-#if defined(EL_HAVE_MPC) && defined(EL_ENABLE_BIGFLOAT)
+#if defined(EL_HAVE_MPC)
     TestColumnTwoNorms<BigFloat, ELEMENT>(m, n, g, print);
     TestColumnTwoNorms<BigFloat, BLOCK>(m, n, g, print);
 #endif
@@ -111,19 +112,17 @@ int main(int argc, char** argv)
     TestColumnMaxNorms<float, BLOCK>(m, n, g, print);
     TestColumnMaxNorms<double, ELEMENT>(m, n, g, print);
     TestColumnMaxNorms<double, BLOCK>(m, n, g, print);
-#if defined(EL_HAVE_QD) && defined(EL_ENABLE_DOUBLEDOUBLE)
+#if defined(EL_HAVE_QD)
     TestColumnMaxNorms<DoubleDouble, ELEMENT>(m, n, g, print);
     TestColumnMaxNorms<DoubleDouble, BLOCK>(m, n, g, print);
-#endif
-#if defined(EL_HAVE_QD) && defined(EL_ENABLE_QUADDOUBLE)
     TestColumnMaxNorms<QuadDouble, ELEMENT>(m, n, g, print);
     TestColumnMaxNorms<QuadDouble, BLOCK>(m, n, g, print);
 #endif
-#if defined(EL_HAVE_QUAD) && defined(EL_ENABLE_QUAD)
+#if defined(EL_HAVE_QUAD)
     TestColumnMaxNorms<Quad, ELEMENT>(m, n, g, print);
     TestColumnMaxNorms<Quad, BLOCK>(m, n, g, print);
 #endif
-#if defined(EL_HAVE_MPC) && defined(EL_ENABLE_BIGFLOAT)
+#if defined(EL_HAVE_MPC)
     TestColumnMaxNorms<BigFloat, ELEMENT>(m, n, g, print);
     TestColumnMaxNorms<BigFloat, BLOCK>(m, n, g, print);
 #endif
diff --git a/tests/blas_like/Dot.cpp b/tests/blas_like/Dot.cpp
index ce54b53b79..6ba62857ad 100644
--- a/tests/blas_like/Dot.cpp
+++ b/tests/blas_like/Dot.cpp
@@ -65,25 +65,23 @@ int main(int argc, char** argv)
     TestDot<double, BLOCK>(m, n, g, print);
     TestDot<Complex<double>, ELEMENT>(m, n, g, print);
     TestDot<Complex<double>, BLOCK>(m, n, g, print);
-#if defined(EL_HAVE_QD) && defined(EL_ENABLE_DOUBLEDOUBLE)
+#if defined(EL_HAVE_QD)
     TestDot<DoubleDouble, ELEMENT>(m, n, g, print);
     TestDot<DoubleDouble, BLOCK>(m, n, g, print);
     TestDot<Complex<DoubleDouble>, ELEMENT>(m, n, g, print);
     TestDot<Complex<DoubleDouble>, BLOCK>(m, n, g, print);
-#endif
-#if defined(EL_HAVE_QD) && defined(EL_ENABLE_QUADDOUBLE)
     TestDot<QuadDouble, ELEMENT>(m, n, g, print);
     TestDot<QuadDouble, BLOCK>(m, n, g, print);
     TestDot<Complex<QuadDouble>, ELEMENT>(m, n, g, print);
     TestDot<Complex<QuadDouble>, BLOCK>(m, n, g, print);
 #endif
-#if defined(EL_HAVE_QUAD) && defined(EL_ENABLE_QUAD)
+#if defined(EL_HAVE_QUAD)
     TestDot<Quad, ELEMENT>(m, n, g, print);
     TestDot<Quad, BLOCK>(m, n, g, print);
     TestDot<Complex<Quad>, ELEMENT>(m, n, g, print);
     TestDot<Complex<Quad>, BLOCK>(m, n, g, print);
 #endif
-#if defined(EL_HAVE_MPC) && defined(EL_ENABLE_BIGFLOAT)
+#if defined(EL_HAVE_MPC)
     TestDot<BigFloat, ELEMENT>(m, n, g, print);
     TestDot<BigFloat, BLOCK>(m, n, g, print);
     TestDot<Complex<BigFloat>, ELEMENT>(m, n, g, print);
diff --git a/tests/blas_like/Hadamard.cpp b/tests/blas_like/Hadamard.cpp
index 73cc5fd83d..0af0fa6468 100644
--- a/tests/blas_like/Hadamard.cpp
+++ b/tests/blas_like/Hadamard.cpp
@@ -68,25 +68,23 @@ int main(int argc, char** argv)
     TestHadamard<double, BLOCK>(m, n, g, print);
     TestHadamard<Complex<double>, ELEMENT>(m, n, g, print);
     TestHadamard<Complex<double>, BLOCK>(m, n, g, print);
-#if defined(EL_HAVE_QD) && defined(EL_ENABLE_DOUBLEDOUBLE)
+#if defined(EL_HAVE_QD)
     TestHadamard<DoubleDouble, ELEMENT>(m, n, g, print);
     TestHadamard<DoubleDouble, BLOCK>(m, n, g, print);
     TestHadamard<Complex<DoubleDouble>, ELEMENT>(m, n, g, print);
     TestHadamard<Complex<DoubleDouble>, BLOCK>(m, n, g, print);
-#endif
-#if defined(EL_HAVE_QD) && defined(EL_ENABLE_QUADDOUBLE)
     TestHadamard<QuadDouble, ELEMENT>(m, n, g, print);
     TestHadamard<QuadDouble, BLOCK>(m, n, g, print);
     TestHadamard<Complex<QuadDouble>, ELEMENT>(m, n, g, print);
     TestHadamard<Complex<QuadDouble>, BLOCK>(m, n, g, print);
 #endif
-#if defined(EL_HAVE_QUAD) && defined(EL_ENABLE_QUAD)
+#if defined(EL_HAVE_QUAD)
     TestHadamard<Quad, ELEMENT>(m, n, g, print);
     TestHadamard<Quad, BLOCK>(m, n, g, print);
     TestHadamard<Complex<Quad>, ELEMENT>(m, n, g, print);
     TestHadamard<Complex<Quad>, BLOCK>(m, n, g, print);
 #endif
-#if defined(EL_HAVE_MPC) && defined(EL_ENABLE_BIGFLOAT)
+#if defined(EL_HAVE_MPC)
     TestHadamard<BigFloat, ELEMENT>(m, n, g, print);
     TestHadamard<BigFloat, BLOCK>(m, n, g, print);
     TestHadamard<Complex<BigFloat>, ELEMENT>(m, n, g, print);