Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Support block matrices in several BLAS1 routines #226

Merged
merged 4 commits into from
Mar 28, 2017
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
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() ||
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you noticing this crucial portion of the extension!

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