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.

Sign up for GitHub

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

Jump to bottom

Add Lanczos solver #1050

Merged
merged 6 commits into from Oct 18, 2019
Merged

Add Lanczos solver #1050

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
fbfbd22
Update BlazeIterative beforing adding Lanczos
NanmiaoWu Oct 8, 2019
ccdda2f
Add Lanczos solver
NanmiaoWu Oct 10, 2019
db1fc01
Change to Clang-format
NanmiaoWu Oct 10, 2019
a5565e5
Add directory for int64_t
NanmiaoWu Oct 14, 2019
4015023
Revise typo
NanmiaoWu Oct 16, 2019
cf93db9
Revise the location of endif
NanmiaoWu Oct 16, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
13 changes: 7 additions & 6 deletions phylanx/plugins/solvers/linear_solver.hpp
View file
Open in desktop
Expand Up @@ -32,6 +32,7 @@ namespace phylanx { namespace execution_tree { namespace primitives

using arg_type = ir::node_data<double>;
using args_type = std::vector<arg_type, arguments_allocator<arg_type>>;
using storage0d_type = typename arg_type::storage0d_type;
using storage1d_type = typename arg_type::storage1d_type;
using storage2d_type = typename arg_type::storage2d_type;

Expand All @@ -43,23 +44,23 @@ namespace phylanx { namespace execution_tree { namespace primitives
std::string const& name, std::string const& codename);

using vector_function = arg_type(args_type&&);
using vector_function_ul = arg_type(
arg_type&&, arg_type&&, std::string);
using vector_function_uln = arg_type(
arg_type&&, arg_type&&, primitive_argument_type&&);

using vector_function_ptr = vector_function*;
using vector_function_ptr_ul = vector_function_ul*;
using vector_function_ptr_uln = vector_function_uln*;

private:
vector_function_ptr get_lin_solver_map(std::string const& name) const;
vector_function_ptr_ul get_lin_solver_map_ul(
vector_function_ptr_uln get_lin_solver_map_uln(
std::string const& name) const;

vector_function_ptr func_;
vector_function_ptr_ul func_ul_;
vector_function_ptr_uln func_uln_;

primitive_argument_type calculate_linear_solver(args_type&& args) const;
primitive_argument_type calculate_linear_solver(
arg_type&& lhs, arg_type&& rhs, std::string ul) const;
arg_type&& lhs, arg_type&& rhs, primitive_argument_type&& ul) const;
};

inline primitive create_linear_solver(hpx::id_type const& locality,
Expand Down
113 changes: 79 additions & 34 deletions src/plugins/solvers/linear_solver.cpp
View file
Open in desktop
@@ -1,6 +1,6 @@
// Copyright (c) 2018 Shahrzad Shirzad
// 2018 Patrick Diehl
//
// 2019 Nanmiao Wu
// Distributed under the Boost Software License, Version 1.0. (See accompanying
// file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

Expand All @@ -22,6 +22,7 @@
#include <string>
#include <utility>
#include <vector>
#include <cstdint>

#include <blaze/Math.h>

Expand Down Expand Up @@ -209,6 +210,20 @@ namespace phylanx { namespace execution_tree { namespace primitives
A matrix `x` such that `a x = b`, solved using the
CG solver utilizing the Symmetric Gauss Seidel as
its preconditioner.)"
),
PHYLANX_LIN_MATCH_DATA("iterative_solver_lanczos",
R"(a, b, n, x
Args:

a (matrix) : a matrix
b (vector) : a vector
n (scalar) : a scalar
x (vector) : a vector

Returns:

A vector of eigenvalues approximate the matrix `a`, whose dimention is decided
by n, solved using the Lanczos solver.)"
)
#endif
};
Expand Down Expand Up @@ -367,31 +382,68 @@ namespace phylanx { namespace execution_tree { namespace primitives
return lin_solver[name];
}

linear_solver::vector_function_ptr_ul linear_solver::get_lin_solver_map_ul(
std::string const& name) const
linear_solver::vector_function_ptr_uln
linear_solver::get_lin_solver_map_uln(std::string const& name) const
{
static std::map<std::string, vector_function_ptr_ul> lin_solver = {
static std::map<std::string, vector_function_ptr_uln> lin_solver = {
{"linear_solver_ldlt",
// Linear solver based on LDLT decomposition of a symmetric
// indefinite matrix
[](arg_type&& arg_0, arg_type&& rhs,
std::string ul) -> arg_type {
[](arg_type&& arg_0, arg_type&& arg_1,
primitive_argument_type&& arg_2) -> arg_type {
std::string ul = extract_string_value_strict(arg_2);
if (ul != "L" && ul != "U")
{
HPX_THROW_EXCEPTION(hpx::bad_parameter,
"linear_solver::eval",
util::generate_error_message(
"the linear_solver primitive requires for "
"the third argument to be either 'L' or 'U'"));
}
storage2d_type A{blaze::trans(arg_0.matrix())};
storage1d_type b{rhs.vector()};
storage1d_type b{arg_1.vector()};
const std::unique_ptr<int[]> ipiv(new int[b.size()]);
blaze::sysv(A, b, ul == "L" ? 'L' : 'U', ipiv.get());
return arg_type{std::move(b)};
}},
{"linear_solver_cholesky",
// Linear solver based on cholesky(LLH) decomposition of a
// positive definite matrix
[](arg_type&& arg_0, arg_type&& rhs,
std::string ul) -> arg_type {
[](arg_type&& arg_0, arg_type&& arg_1,
primitive_argument_type&& arg_2) -> arg_type {
std::string ul =
extract_string_value_strict(std::move(arg_2));
if (ul != "L" && ul != "U")
{
HPX_THROW_EXCEPTION(hpx::bad_parameter,
"linear_solver::eval",
util::generate_error_message(
"the linear_solver primitive requires for "
"the third argument to be either 'L' or 'U'"));
}
storage2d_type A{blaze::trans(arg_0.matrix())};
storage1d_type b{rhs.vector()};
storage1d_type b{arg_1.vector()};
blaze::posv(A, b, ul == "L" ? 'L' : 'U');
return arg_type{std::move(b)};
}}};
}},
#ifdef PHYLANX_HAVE_BLAZE_ITERATIVE
{"iterative_solver_lanczos",
// Iterative Lanczos solver for eigenvalues
// Note: Relies on BlazeIterative library and
// need to be explicitly enabled
[](arg_type&& arg_0, arg_type&& arg_1,
primitive_argument_type&& arg_2) -> arg_type {
std::int64_t n =
extract_scalar_integer_value_strict(std::move(arg_2));
storage2d_type A{blaze::trans(arg_0.matrix())};
storage1d_type b{arg_1.vector()};
storage1d_type x;
blaze::iterative::LanczosTag tag;
x = blaze::iterative::solve(A, b, tag, n);
return arg_type{std::move(x)};
}}
#endif
};
return lin_solver[name];
}

Expand All @@ -403,9 +455,9 @@ namespace phylanx { namespace execution_tree { namespace primitives
std::string func_name = extract_function_name(name);

func_ = get_lin_solver_map(func_name);
func_ul_ = get_lin_solver_map_ul(func_name);
func_uln_ = get_lin_solver_map_uln(func_name);

HPX_ASSERT(func_ != nullptr || func_ul_ != nullptr);
HPX_ASSERT(func_ != nullptr || func_uln_ != nullptr);
}

///////////////////////////////////////////////////////////////////////////
Expand All @@ -424,18 +476,19 @@ namespace phylanx { namespace execution_tree { namespace primitives
}

primitive_argument_type linear_solver::calculate_linear_solver(
arg_type && lhs, arg_type && rhs, std::string ul) const
arg_type && lhs, arg_type && rhs, primitive_argument_type&& uln) const
{
if (func_ul_ == nullptr)
if (func_uln_ == nullptr)
{
HPX_THROW_EXCEPTION(hpx::bad_parameter,
"linear_solver::eval",
util::generate_error_message(
"this linear_solver primitive requires exactly two operands",
name_, codename_));
}
return primitive_argument_type{
func_ul_(std::move(lhs), std::move(rhs), std::move(ul))};

return primitive_argument_type{
func_uln_(std::move(lhs), std::move(rhs), std::move(uln))};
}

hpx::future<primitive_argument_type> linear_solver::eval(
Expand All @@ -446,8 +499,8 @@ namespace phylanx { namespace execution_tree { namespace primitives
{
HPX_THROW_EXCEPTION(hpx::bad_parameter,
"linear_solver::eval",
generate_error_message(
"the linear_solver primitive requires exactly two operands"));
generate_error_message("the linear_solver primitive requires "
"exactly two operands"));
}

if (!valid(operands[0]) || !valid(operands[1]) ||
Expand All @@ -463,11 +516,11 @@ namespace phylanx { namespace execution_tree { namespace primitives
auto this_ = this->shared_from_this();
if (operands.size() == 3)
{
return hpx::dataflow(hpx::launch::sync, hpx::util::unwrapping(
[this_ = std::move(this_)](
arg_type&& lhs, arg_type&& rhs, std::string && ul)
-> primitive_argument_type
{
return hpx::dataflow(hpx::launch::sync,
hpx::util::unwrapping([this_ = std::move(this_)](arg_type&& lhs,
arg_type&& rhs,
primitive_argument_type&& uln)
-> primitive_argument_type {
if (lhs.num_dimensions() != 2 || rhs.num_dimensions() != 1)
{
HPX_THROW_EXCEPTION(hpx::bad_parameter,
Expand All @@ -477,21 +530,13 @@ namespace phylanx { namespace execution_tree { namespace primitives
"that first operand to be a matrix and "
"the second operand to be a vector"));
}
if (ul != "L" && ul != "U")
{
HPX_THROW_EXCEPTION(hpx::bad_parameter,
"linear_solver::eval",
this_->generate_error_message(
"the linear_solver primitive requires for "
"the third argument to be either 'L' or 'U'"));
}

return this_->calculate_linear_solver(
std::move(lhs), std::move(rhs), std::move(ul));
std::move(lhs), std::move(rhs), std::move(uln));
}),
numeric_operand(operands[0], args, name_, codename_, ctx),
numeric_operand(operands[1], args, name_, codename_, ctx),
string_operand(operands[2], args, name_, codename_, ctx));
value_operand(operands[2], args, name_, codename_, ctx));
}

return hpx::dataflow(hpx::launch::sync, hpx::util::unwrapping(
Expand Down
54 changes: 54 additions & 0 deletions tests/unit/plugins/solvers/linear_solver.cpp
View file
Open in desktop
Expand Up @@ -12,6 +12,7 @@
#include <iostream>
#include <string>
#include <utility>
#include <cstdint>
#include <vector>

///////////////////////////////////////////////////////////////////////////////
Expand Down Expand Up @@ -192,6 +193,24 @@ void test_linear_solver_cg_symmetric_gauss_seidel_PhySL()
phylanx::ir::node_data<double>(blaze::DynamicVector<double>{1, 2, 3}));
}

void test_linear_solver_lanczos_PhySL()
{
std::string const code = R"(block(
define(a, [[0, 0, 0, 0, 0, 0],[0, 1, 0, 0, 0, 0],[0, 0, 2, 0, 0, 0],
[0, 0, 0, 3, 0, 0],[0, 0, 0, 0, 4, 0],[0, 0, 0, 0, 0, 100000]]),
define(b, [1, 1, 1, 1, 1, 1]),
define(n, 6),
iterative_solver_lanczos(a, b, n))
)";

auto result =
phylanx::execution_tree::extract_numeric_value(compile_and_run(code));

HPX_TEST_EQ(result,
phylanx::ir::node_data<double>(
blaze::DynamicVector<double>{100000, 0, 4, 1, 2, 3}));
}

void test_linear_solver_lu(std::string const& func_name)
{
phylanx::execution_tree::primitive lhs =
Expand Down Expand Up @@ -300,6 +319,38 @@ void test_linear_solver_u(std::string const& func_name)
phylanx::ir::node_data<double>(blaze::DynamicVector<double>{1, 2, 3}));
}

void test_linear_solver_n(std::string const& func_name)
{
phylanx::execution_tree::primitive lhs =
phylanx::execution_tree::primitives::create_variable(hpx::find_here(),
phylanx::ir::node_data<double>{
blaze::DynamicMatrix<double>{{0, 0, 0, 0, 0, 0},
{0, 1, 0, 0, 0, 0}, {0, 0, 2, 0, 0, 0}, {0, 0, 0, 3, 0, 0},
{0, 0, 0, 0, 4, 0}, {0, 0, 0, 0, 0, 100000}}});

phylanx::execution_tree::primitive rhs =
phylanx::execution_tree::primitives::create_variable(hpx::find_here(),
phylanx::ir::node_data<double>{
blaze::DynamicVector<double>{1, 1, 1, 1, 1, 1}});

phylanx::execution_tree::primitive n =
phylanx::execution_tree::primitives::create_variable(
hpx::find_here(), phylanx::ir::node_data<std::int64_t>(6));

phylanx::execution_tree::primitive linear_solver =
phylanx::execution_tree::primitives::create_linear_solver(
hpx::find_here(),
phylanx::execution_tree::primitive_arguments_type{
std::move(lhs), std::move(rhs), std::move(n)},
func_name);

hpx::future<phylanx::execution_tree::primitive_argument_type> f =
linear_solver.eval();
HPX_TEST_EQ(f.get(),
phylanx::ir::node_data<double>(
blaze::DynamicVector<double>{100000, 0, 4, 1, 2, 3}));
}

int main()
{
test_linear_solver_lu_PhySL();
Expand All @@ -317,6 +368,8 @@ int main()
test_linear_solver_cg_incompleteCholesky_PhySL();
test_linear_solver_cg_symmetric_gauss_seidel_PhySL();

test_linear_solver_lanczos_PhySL();

test_linear_solver_lu("linear_solver_lu");

test_linear_solver("linear_solver_ldlt");
Expand All @@ -337,6 +390,7 @@ int main()
test_linear_solver("iterative_solver_cg_ssor");
test_linear_solver("iterative_solver_cg_incompleteCholesky");
test_linear_solver("iterative_solver_cg_symmetric_gauss_seidel");
test_linear_solver_n("iterative_solver_lanczos");
#endif
return hpx::util::report_errors();
}