-
Notifications
You must be signed in to change notification settings - Fork 110
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
Changes from 3 commits
90b6408
4afe36a
be1c42a
e9139f4
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,135 @@ | ||
/* | ||
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; | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Explicitly accumulating the squares and then square-rooting should be fine for random data but can lead to very low accuracy in extreme cases. There is the routine There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. My main goal there was to provide a simple baseline that ensures there's no major issues (e.g. block vs element distributions causing problems) while not just reimplementing the method that's being tested. If you think it would be better to switch to |
||
} | ||
expected = mpi::AllReduce(expected, g.ColComm()); | ||
expected = Sqrt(expected); | ||
if (Abs(got - expected) > 10 * limits::Epsilon<El::Base<T>>()) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. The accumulated error should be a function of the number of entries (with said function being logarithmic if a tree-based scheme is used, but linear with the current single-process implementation): would you mind multiplying the right-hand side by the number of entries? Also, it would be better to use a relative bound and to diving the left-hand side by the maximum of the norm and 1 (so that it still behaves well for zero norms). There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Just want to clarify: The new check should be |
||
{ | ||
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) && 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); | ||
#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); | ||
} | ||
} |
There was a problem hiding this comment.
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!