Skip to content

Commit

Permalink
Merge pull request #226 from ndryden/master
Browse files Browse the repository at this point in the history
Support block matrices in several BLAS1 routines
  • Loading branch information
poulson authored Mar 28, 2017
2 parents 7261f9e + e9139f4 commit 776b805
Show file tree
Hide file tree
Showing 8 changed files with 366 additions and 27 deletions.
4 changes: 2 additions & 2 deletions include/El/blas_like/level1/Dot.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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 );
Expand Down Expand Up @@ -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 \
Expand Down
15 changes: 9 additions & 6 deletions include/El/blas_like/level1/Hadamard.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand All @@ -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() );
Expand Down Expand Up @@ -87,9 +90,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, \
Expand Down
18 changes: 9 additions & 9 deletions include/El/blas_like/level1/decl.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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 );
Expand Down Expand Up @@ -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 );
Expand Down Expand Up @@ -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 );

Expand Down Expand Up @@ -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,
Expand All @@ -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 );

Expand Down
22 changes: 14 additions & 8 deletions src/blas_like/level1/ColumnNorms.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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 );
Expand All @@ -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 );
Expand Down Expand Up @@ -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 \
Expand Down
7 changes: 5 additions & 2 deletions src/blas_like/level1/HilbertSchmidt.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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() )
Expand All @@ -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() )
Expand Down Expand Up @@ -112,7 +115,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 );

Expand Down
134 changes: 134 additions & 0 deletions tests/blas_like/ColumnNorms.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,134 @@
/*
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);
// 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);
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);
#if defined(EL_HAVE_QD)
TestColumnTwoNorms<DoubleDouble, ELEMENT>(m, n, g, print);
TestColumnTwoNorms<DoubleDouble, BLOCK>(m, n, g, print);
TestColumnTwoNorms<QuadDouble, ELEMENT>(m, n, g, print);
TestColumnTwoNorms<QuadDouble, BLOCK>(m, n, g, print);
#endif
#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)
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);
#if defined(EL_HAVE_QD)
TestColumnMaxNorms<DoubleDouble, ELEMENT>(m, n, g, print);
TestColumnMaxNorms<DoubleDouble, BLOCK>(m, n, g, print);
TestColumnMaxNorms<QuadDouble, ELEMENT>(m, n, g, print);
TestColumnMaxNorms<QuadDouble, BLOCK>(m, n, g, print);
#endif
#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)
TestColumnMaxNorms<BigFloat, ELEMENT>(m, n, g, print);
TestColumnMaxNorms<BigFloat, BLOCK>(m, n, g, print);
#endif
}
catch (exception& e)
{
ReportException(e);
}
}
Loading

0 comments on commit 776b805

Please sign in to comment.