Parallel updates to block-diagonal EtE FtF

--------------------------------------------------------------------------------
Benchmark                                                                   Time
--------------------------------------------------------------------------------
PMVUpdateBlockDiagonalFtF<final/problem-13682-4456117-pre.txt>/1   5056275941 ns
PMVUpdateBlockDiagonalFtF<final/problem-13682-4456117-pre.txt>/2   3677677097 ns
PMVUpdateBlockDiagonalFtF<final/problem-13682-4456117-pre.txt>/4   1932236015 ns
PMVUpdateBlockDiagonalFtF<final/problem-13682-4456117-pre.txt>/8    984585015 ns
PMVUpdateBlockDiagonalFtF<final/problem-13682-4456117-pre.txt>/16   614918752 ns

PMVUpdateBlockDiagonalEtE<final/problem-13682-4456117-pre.txt>/1    449324491 ns
PMVUpdateBlockDiagonalEtE<final/problem-13682-4456117-pre.txt>/2    273147462 ns
PMVUpdateBlockDiagonalEtE<final/problem-13682-4456117-pre.txt>/4    150742698 ns
PMVUpdateBlockDiagonalEtE<final/problem-13682-4456117-pre.txt>/8     81602564 ns
PMVUpdateBlockDiagonalEtE<final/problem-13682-4456117-pre.txt>/16    47010769 ns

PMVUpdateBlockDiagonalFtF<venice/problem-1778-993923-pre.txt>/1     774598200 ns
PMVUpdateBlockDiagonalFtF<venice/problem-1778-993923-pre.txt>/2     611312877 ns
PMVUpdateBlockDiagonalFtF<venice/problem-1778-993923-pre.txt>/4     326701149 ns
PMVUpdateBlockDiagonalFtF<venice/problem-1778-993923-pre.txt>/8     165634457 ns
PMVUpdateBlockDiagonalFtF<venice/problem-1778-993923-pre.txt>/16     90631068 ns

PMVUpdateBlockDiagonalEtE<venice/problem-1778-993923-pre.txt>/1      80651817 ns
PMVUpdateBlockDiagonalEtE<venice/problem-1778-993923-pre.txt>/2      49688691 ns
PMVUpdateBlockDiagonalEtE<venice/problem-1778-993923-pre.txt>/4      27199153 ns
PMVUpdateBlockDiagonalEtE<venice/problem-1778-993923-pre.txt>/8      14301768 ns
PMVUpdateBlockDiagonalEtE<venice/problem-1778-993923-pre.txt>/16      8683479 ns

PMVUpdateBlockDiagonalFtF<ladybug/problem-1723-156502-pre.txt>/1    104422529 ns
PMVUpdateBlockDiagonalFtF<ladybug/problem-1723-156502-pre.txt>/2     81555176 ns
PMVUpdateBlockDiagonalFtF<ladybug/problem-1723-156502-pre.txt>/4     43227593 ns
PMVUpdateBlockDiagonalFtF<ladybug/problem-1723-156502-pre.txt>/8     22177895 ns
PMVUpdateBlockDiagonalFtF<ladybug/problem-1723-156502-pre.txt>/16    12505813 ns

PMVUpdateBlockDiagonalEtE<ladybug/problem-1723-156502-pre.txt>/1     14205253 ns
PMVUpdateBlockDiagonalEtE<ladybug/problem-1723-156502-pre.txt>/2      7102357 ns
PMVUpdateBlockDiagonalEtE<ladybug/problem-1723-156502-pre.txt>/4      3806598 ns
PMVUpdateBlockDiagonalEtE<ladybug/problem-1723-156502-pre.txt>/8      2112615 ns
PMVUpdateBlockDiagonalEtE<ladybug/problem-1723-156502-pre.txt>/16     1245771 ns

PMVUpdateBlockDiagonalFtF<dubrovnik/problem-356-226730-pre.txt>/1   190102870 ns
PMVUpdateBlockDiagonalFtF<dubrovnik/problem-356-226730-pre.txt>/2   157359897 ns
PMVUpdateBlockDiagonalFtF<dubrovnik/problem-356-226730-pre.txt>/4    82657662 ns
PMVUpdateBlockDiagonalFtF<dubrovnik/problem-356-226730-pre.txt>/8    42746490 ns
PMVUpdateBlockDiagonalFtF<dubrovnik/problem-356-226730-pre.txt>/16   23434967 ns

PMVUpdateBlockDiagonalEtE<dubrovnik/problem-356-226730-pre.txt>/1    18894355 ns
PMVUpdateBlockDiagonalEtE<dubrovnik/problem-356-226730-pre.txt>/2    12138228 ns
PMVUpdateBlockDiagonalEtE<dubrovnik/problem-356-226730-pre.txt>/4     6808771 ns
PMVUpdateBlockDiagonalEtE<dubrovnik/problem-356-226730-pre.txt>/8     3829718 ns
PMVUpdateBlockDiagonalEtE<dubrovnik/problem-356-226730-pre.txt>/16    2103688 ns

PMVUpdateBlockDiagonalFtF<trafalgar/problem-257-65132-pre.txt>/1     34230036 ns
PMVUpdateBlockDiagonalFtF<trafalgar/problem-257-65132-pre.txt>/2     20121184 ns
PMVUpdateBlockDiagonalFtF<trafalgar/problem-257-65132-pre.txt>/4     10899938 ns
PMVUpdateBlockDiagonalFtF<trafalgar/problem-257-65132-pre.txt>/8      6186359 ns
PMVUpdateBlockDiagonalFtF<trafalgar/problem-257-65132-pre.txt>/16     4207103 ns

PMVUpdateBlockDiagonalEtE<trafalgar/problem-257-65132-pre.txt>/1      3077110 ns
PMVUpdateBlockDiagonalEtE<trafalgar/problem-257-65132-pre.txt>/2      2224104 ns
PMVUpdateBlockDiagonalEtE<trafalgar/problem-257-65132-pre.txt>/4      1274841 ns
PMVUpdateBlockDiagonalEtE<trafalgar/problem-257-65132-pre.txt>/8       721140 ns
PMVUpdateBlockDiagonalEtE<trafalgar/problem-257-65132-pre.txt>/16      437715 ns

Change-Id: If5a342a063869bd9c0505bf96b6d957da5169c1d
This commit is contained in:
Dmitriy Korchemkin
2022-11-17 19:29:56 +03:00
parent 946fa50de7
commit 77a54dd3d2
4 changed files with 330 additions and 81 deletions
+75 -11
View File
@@ -147,19 +147,28 @@ struct BALData {
const PartitionedView* PartitionedMatrixViewJacobian(
const LinearSolver::Options& options) {
auto block_sparse = BlockSparseJacobian(options.context);
auto block_sparse = BlockSparseJacobianWithTranspose(options.context);
partitioned_view_jacobian =
std::make_unique<PartitionedView>(options, *block_sparse);
return partitioned_view_jacobian.get();
}
const PartitionedView* PartitionedMatrixViewJacobianWithTranspose(
const LinearSolver::Options& options) {
auto block_sparse_transpose =
BlockSparseJacobianWithTranspose(options.context);
partitioned_view_jacobian_with_transpose =
std::make_unique<PartitionedView>(options, *block_sparse_transpose);
return partitioned_view_jacobian_with_transpose.get();
BlockSparseMatrix* BlockDiagonalEtE(const LinearSolver::Options& options) {
if (!block_diagonal_ete) {
auto partitioned_view =
PartitionedMatrixViewJacobian(options);
block_diagonal_ete = partitioned_view->CreateBlockDiagonalEtE();
}
return block_diagonal_ete.get();
}
BlockSparseMatrix* BlockDiagonalFtF(const LinearSolver::Options& options) {
if (!block_diagonal_ftf) {
auto partitioned_view =
PartitionedMatrixViewJacobian(options);
block_diagonal_ftf = partitioned_view->CreateBlockDiagonalFtF();
}
return block_diagonal_ftf.get();
}
Vector parameters;
@@ -169,7 +178,8 @@ struct BALData {
std::unique_ptr<BlockSparseMatrix> block_sparse_jacobian_with_transpose;
std::unique_ptr<CompressedRowSparseMatrix> crs_jacobian;
std::unique_ptr<PartitionedView> partitioned_view_jacobian;
std::unique_ptr<PartitionedView> partitioned_view_jacobian_with_transpose;
std::unique_ptr<BlockSparseMatrix> block_diagonal_ete;
std::unique_ptr<BlockSparseMatrix> block_diagonal_ftf;
};
static void Residuals(benchmark::State& state,
@@ -262,7 +272,7 @@ static void PMVLeftMultiplyAndAccumulateF(benchmark::State& state,
options.num_threads = state.range(0);
options.elimination_groups.push_back(data->bal_problem->num_points());
options.context = context;
auto jacobian = data->PartitionedMatrixViewJacobianWithTranspose(options);
auto jacobian = data->PartitionedMatrixViewJacobian(options);
Vector y = Vector::Zero(jacobian->num_cols_f());
Vector x = Vector::Random(jacobian->num_rows());
@@ -298,7 +308,7 @@ static void PMVLeftMultiplyAndAccumulateE(benchmark::State& state,
options.num_threads = state.range(0);
options.elimination_groups.push_back(data->bal_problem->num_points());
options.context = context;
auto jacobian = data->PartitionedMatrixViewJacobianWithTranspose(options);
auto jacobian = data->PartitionedMatrixViewJacobian(options);
Vector y = Vector::Zero(jacobian->num_cols_e());
Vector x = Vector::Random(jacobian->num_rows());
@@ -309,6 +319,36 @@ static void PMVLeftMultiplyAndAccumulateE(benchmark::State& state,
CHECK_GT(y.squaredNorm(), 0.);
}
static void PMVUpdateBlockDiagonalEtE(benchmark::State& state,
BALData* data,
ContextImpl* context) {
LinearSolver::Options options;
options.num_threads = state.range(0);
options.elimination_groups.push_back(data->bal_problem->num_points());
options.context = context;
auto jacobian = data->PartitionedMatrixViewJacobian(options);
auto block_diagonal_ete = data->BlockDiagonalEtE(options);
for (auto _ : state) {
jacobian->UpdateBlockDiagonalEtE(block_diagonal_ete);
}
}
static void PMVUpdateBlockDiagonalFtF(benchmark::State& state,
BALData* data,
ContextImpl* context) {
LinearSolver::Options options;
options.num_threads = state.range(0);
options.elimination_groups.push_back(data->bal_problem->num_points());
options.context = context;
auto jacobian = data->PartitionedMatrixViewJacobian(options);
auto block_diagonal_ftf = data->BlockDiagonalFtF(options);
for (auto _ : state) {
jacobian->UpdateBlockDiagonalFtF(block_diagonal_ftf);
}
}
static void JacobianRightMultiplyAndAccumulate(benchmark::State& state,
BALData* data,
ContextImpl* context) {
@@ -485,6 +525,30 @@ int main(int argc, char** argv) {
->Arg(8)
->Arg(16);
const std::string name_update_block_diagonal_ftf =
"PMVUpdateBlockDiagonalFtF<" + path + ">";
::benchmark::RegisterBenchmark(name_update_block_diagonal_ftf.c_str(),
ceres::internal::PMVUpdateBlockDiagonalFtF,
data,
&context)
->Arg(1)
->Arg(2)
->Arg(4)
->Arg(8)
->Arg(16);
const std::string name_update_block_diagonal_ete =
"PMVUpdateBlockDiagonalEtE<" + path + ">";
::benchmark::RegisterBenchmark(name_update_block_diagonal_ete.c_str(),
ceres::internal::PMVUpdateBlockDiagonalEtE,
data,
&context)
->Arg(1)
->Arg(2)
->Arg(4)
->Arg(8)
->Arg(16);
#ifndef CERES_NO_CUDA
const std::string name_right_product_cuda =
"JacobianRightMultiplyAndAccumulateCuda<" + path + ">";
+8
View File
@@ -166,7 +166,15 @@ class CERES_NO_EXPORT PartitionedMatrixView final
std::unique_ptr<BlockSparseMatrix> CreateBlockDiagonalEtE() const final;
std::unique_ptr<BlockSparseMatrix> CreateBlockDiagonalFtF() const final;
void UpdateBlockDiagonalEtE(BlockSparseMatrix* block_diagonal) const final;
void UpdateBlockDiagonalEtESingleThreaded(
BlockSparseMatrix* block_diagonal) const;
void UpdateBlockDiagonalEtEMultiThreaded(
BlockSparseMatrix* block_diagonal) const;
void UpdateBlockDiagonalFtF(BlockSparseMatrix* block_diagonal) const final;
void UpdateBlockDiagonalFtFSingleThreaded(
BlockSparseMatrix* block_diagonal) const;
void UpdateBlockDiagonalFtFMultiThreaded(
BlockSparseMatrix* block_diagonal) const;
// clang-format off
int num_col_blocks_e() const final { return num_col_blocks_e_; }
int num_col_blocks_f() const final { return num_col_blocks_f_; }
+137 -8
View File
@@ -444,10 +444,10 @@ PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
//
template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
UpdateBlockDiagonalEtE(BlockSparseMatrix* block_diagonal) const {
const CompressedRowBlockStructure* bs = matrix_.block_structure();
const CompressedRowBlockStructure* block_diagonal_structure =
block_diagonal->block_structure();
UpdateBlockDiagonalEtESingleThreaded(
BlockSparseMatrix* block_diagonal) const {
auto bs = matrix_.block_structure();
auto block_diagonal_structure = block_diagonal->block_structure();
block_diagonal->SetZero();
const double* values = matrix_.values();
@@ -470,6 +470,57 @@ void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
}
}
template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
UpdateBlockDiagonalEtEMultiThreaded(
BlockSparseMatrix* block_diagonal) const {
auto transpose_block_structure = matrix_.transpose_block_structure();
CHECK(transpose_block_structure != nullptr);
auto block_diagonal_structure = block_diagonal->block_structure();
const double* values = matrix_.values();
double* values_diagonal = block_diagonal->mutable_values();
ParallelFor(
options_.context,
0,
num_col_blocks_e_,
options_.num_threads,
[values,
transpose_block_structure,
values_diagonal,
block_diagonal_structure](int col_block_id) {
int cell_position =
block_diagonal_structure->rows[col_block_id].cells[0].position;
double* cell_values = values_diagonal + cell_position;
int col_block_size =
transpose_block_structure->rows[col_block_id].block.size;
auto& cells = transpose_block_structure->rows[col_block_id].cells;
MatrixRef(cell_values, col_block_size, col_block_size).setZero();
for (auto& c : cells) {
int row_block_size = transpose_block_structure->cols[c.block_id].size;
// clang-format off
MatrixTransposeMatrixMultiply<kRowBlockSize, kEBlockSize, kRowBlockSize, kEBlockSize, 1>(
values + c.position, row_block_size, col_block_size,
values + c.position, row_block_size, col_block_size,
cell_values, 0, 0, col_block_size, col_block_size);
// clang-format on
}
},
e_cols_partition_);
}
template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
UpdateBlockDiagonalEtE(BlockSparseMatrix* block_diagonal) const {
if (options_.num_threads == 1) {
UpdateBlockDiagonalEtESingleThreaded(block_diagonal);
} else {
CHECK(options_.context != nullptr);
UpdateBlockDiagonalEtEMultiThreaded(block_diagonal);
}
}
// Similar to the code in RightMultiplyAndAccumulateF, except instead of the
// matrix vector multiply its an outer product.
//
@@ -477,10 +528,10 @@ void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
//
template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
UpdateBlockDiagonalFtF(BlockSparseMatrix* block_diagonal) const {
const CompressedRowBlockStructure* bs = matrix_.block_structure();
const CompressedRowBlockStructure* block_diagonal_structure =
block_diagonal->block_structure();
UpdateBlockDiagonalFtFSingleThreaded(
BlockSparseMatrix* block_diagonal) const {
auto bs = matrix_.block_structure();
auto block_diagonal_structure = block_diagonal->block_structure();
block_diagonal->SetZero();
const double* values = matrix_.values();
@@ -527,4 +578,82 @@ void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
}
}
template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
UpdateBlockDiagonalFtFMultiThreaded(
BlockSparseMatrix* block_diagonal) const {
auto transpose_block_structure = matrix_.transpose_block_structure();
CHECK(transpose_block_structure != nullptr);
auto block_diagonal_structure = block_diagonal->block_structure();
const double* values = matrix_.values();
double* values_diagonal = block_diagonal->mutable_values();
const int num_col_blocks_e = num_col_blocks_e_;
const int num_row_blocks_e = num_row_blocks_e_;
ParallelFor(
options_.context,
num_col_blocks_e_,
num_col_blocks_e + num_col_blocks_f_,
options_.num_threads,
[transpose_block_structure,
block_diagonal_structure,
num_col_blocks_e,
num_row_blocks_e,
values,
values_diagonal](int col_block_id) {
const int col_block_size =
transpose_block_structure->rows[col_block_id].block.size;
const int diagonal_block_id = col_block_id - num_col_blocks_e;
const int cell_position =
block_diagonal_structure->rows[diagonal_block_id].cells[0].position;
double* cell_values = values_diagonal + cell_position;
MatrixRef(cell_values, col_block_size, col_block_size).setZero();
auto& cells = transpose_block_structure->rows[col_block_id].cells;
const int num_cells = cells.size();
int i = 0;
for (; i < num_cells; ++i) {
auto& cell = cells[i];
const int row_block_id = cell.block_id;
if (row_block_id >= num_row_blocks_e) break;
const int row_block_size =
transpose_block_structure->cols[row_block_id].size;
// clang-format off
MatrixTransposeMatrixMultiply
<kRowBlockSize, kFBlockSize, kRowBlockSize, kFBlockSize, 1>(
values + cell.position, row_block_size, col_block_size,
values + cell.position, row_block_size, col_block_size,
cell_values, 0, 0, col_block_size, col_block_size);
// clang-format on
}
for (; i < num_cells; ++i) {
auto& cell = cells[i];
const int row_block_id = cell.block_id;
const int row_block_size =
transpose_block_structure->cols[row_block_id].size;
// clang-format off
MatrixTransposeMatrixMultiply
<Eigen::Dynamic, Eigen::Dynamic, Eigen::Dynamic, Eigen::Dynamic, 1>(
values + cell.position, row_block_size, col_block_size,
values + cell.position, row_block_size, col_block_size,
cell_values, 0, 0, col_block_size, col_block_size);
// clang-format on
}
},
f_cols_partition_);
}
template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
UpdateBlockDiagonalFtF(BlockSparseMatrix* block_diagonal) const {
if (options_.num_threads == 1) {
UpdateBlockDiagonalFtFSingleThreaded(block_diagonal);
} else {
CHECK(options_.context != nullptr);
UpdateBlockDiagonalFtFMultiThreaded(block_diagonal);
}
}
} // namespace ceres::internal
+110 -62
View File
@@ -49,63 +49,6 @@ namespace internal {
const double kEpsilon = 1e-14;
class PartitionedMatrixViewTest : public ::testing::Test {
protected:
void SetUp() final {
std::unique_ptr<LinearLeastSquaresProblem> problem =
CreateLinearLeastSquaresProblemFromId(2);
CHECK(problem != nullptr);
A_ = std::move(problem->A);
num_cols_ = A_->num_cols();
num_rows_ = A_->num_rows();
num_eliminate_blocks_ = problem->num_eliminate_blocks;
LinearSolver::Options options;
options.elimination_groups.push_back(num_eliminate_blocks_);
pmv_ = PartitionedMatrixViewBase::Create(
options, *down_cast<BlockSparseMatrix*>(A_.get()));
}
double RandDouble() { return distribution_(prng_); }
int num_rows_;
int num_cols_;
int num_eliminate_blocks_;
std::unique_ptr<SparseMatrix> A_;
std::unique_ptr<PartitionedMatrixViewBase> pmv_;
std::mt19937 prng_;
std::uniform_real_distribution<double> distribution_ =
std::uniform_real_distribution<double>(0.0, 1.0);
};
TEST_F(PartitionedMatrixViewTest, BlockDiagonalEtE) {
std::unique_ptr<BlockSparseMatrix> block_diagonal_ee(
pmv_->CreateBlockDiagonalEtE());
const CompressedRowBlockStructure* bs = block_diagonal_ee->block_structure();
EXPECT_EQ(block_diagonal_ee->num_rows(), 2);
EXPECT_EQ(block_diagonal_ee->num_cols(), 2);
EXPECT_EQ(bs->cols.size(), 2);
EXPECT_EQ(bs->rows.size(), 2);
EXPECT_NEAR(block_diagonal_ee->values()[0], 10.0, kEpsilon);
EXPECT_NEAR(block_diagonal_ee->values()[1], 155.0, kEpsilon);
}
TEST_F(PartitionedMatrixViewTest, BlockDiagonalFtF) {
std::unique_ptr<BlockSparseMatrix> block_diagonal_ff(
pmv_->CreateBlockDiagonalFtF());
const CompressedRowBlockStructure* bs = block_diagonal_ff->block_structure();
EXPECT_EQ(block_diagonal_ff->num_rows(), 3);
EXPECT_EQ(block_diagonal_ff->num_cols(), 3);
EXPECT_EQ(bs->cols.size(), 3);
EXPECT_EQ(bs->rows.size(), 3);
EXPECT_NEAR(block_diagonal_ff->values()[0], 70.0, kEpsilon);
EXPECT_NEAR(block_diagonal_ff->values()[1], 17.0, kEpsilon);
EXPECT_NEAR(block_diagonal_ff->values()[2], 37.0, kEpsilon);
}
// Param = <problem_id, num_threads>
using Param = ::testing::tuple<int, int>;
@@ -116,7 +59,7 @@ static std::string ParamInfoToString(testing::TestParamInfo<Param> info) {
return ss.str();
}
class PartitionedMatrixViewSpMVTest : public ::testing::TestWithParam<Param> {
class PartitionedMatrixViewTest : public ::testing::TestWithParam<Param> {
protected:
void SetUp() final {
const int problem_id = ::testing::get<0>(GetParam());
@@ -132,6 +75,10 @@ class PartitionedMatrixViewSpMVTest : public ::testing::TestWithParam<Param> {
options_.elimination_groups.push_back(problem->num_eliminate_blocks);
pmv_ = PartitionedMatrixViewBase::Create(options_, *block_sparse);
LinearSolver::Options options_single_threaded = options_;
options_single_threaded.num_threads = 1;
pmv_single_threaded_ = PartitionedMatrixViewBase::Create(options_, *block_sparse);
EXPECT_EQ(pmv_->num_col_blocks_e(), problem->num_eliminate_blocks);
EXPECT_EQ(pmv_->num_col_blocks_f(),
block_sparse->block_structure()->cols.size() -
@@ -147,12 +94,13 @@ class PartitionedMatrixViewSpMVTest : public ::testing::TestWithParam<Param> {
std::unique_ptr<LinearLeastSquaresProblem> problem_;
std::unique_ptr<SparseMatrix> A_;
std::unique_ptr<PartitionedMatrixViewBase> pmv_;
std::unique_ptr<PartitionedMatrixViewBase> pmv_single_threaded_;
std::mt19937 prng_;
std::uniform_real_distribution<double> distribution_ =
std::uniform_real_distribution<double>(0.0, 1.0);
};
TEST_P(PartitionedMatrixViewSpMVTest, RightMultiplyAndAccumulateE) {
TEST_P(PartitionedMatrixViewTest, RightMultiplyAndAccumulateE) {
Vector x1(pmv_->num_cols_e());
Vector x2(pmv_->num_cols());
x2.setZero();
@@ -172,7 +120,7 @@ TEST_P(PartitionedMatrixViewSpMVTest, RightMultiplyAndAccumulateE) {
}
}
TEST_P(PartitionedMatrixViewSpMVTest, RightMultiplyAndAccumulateF) {
TEST_P(PartitionedMatrixViewTest, RightMultiplyAndAccumulateF) {
Vector x1(pmv_->num_cols_f());
Vector x2(pmv_->num_cols());
x2.setZero();
@@ -192,7 +140,7 @@ TEST_P(PartitionedMatrixViewSpMVTest, RightMultiplyAndAccumulateF) {
}
}
TEST_P(PartitionedMatrixViewSpMVTest, LeftMultiplyAndAccumulate) {
TEST_P(PartitionedMatrixViewTest, LeftMultiplyAndAccumulate) {
Vector x = Vector::Zero(pmv_->num_rows());
for (int i = 0; i < pmv_->num_rows(); ++i) {
x(i) = RandDouble();
@@ -215,9 +163,109 @@ TEST_P(PartitionedMatrixViewSpMVTest, LeftMultiplyAndAccumulate) {
}
}
TEST_P(PartitionedMatrixViewTest, BlockDiagonalFtF) {
std::unique_ptr<BlockSparseMatrix> block_diagonal_ff(
pmv_->CreateBlockDiagonalFtF());
const auto bs_diagonal = block_diagonal_ff->block_structure();
const int num_rows = pmv_->num_rows();
const int num_cols_f = pmv_->num_cols_f();
const int num_cols_e = pmv_->num_cols_e();
const int num_col_blocks_f = pmv_->num_col_blocks_f();
const int num_col_blocks_e = pmv_->num_col_blocks_e();
CHECK_EQ(block_diagonal_ff->num_rows(), num_cols_f);
CHECK_EQ(block_diagonal_ff->num_cols(), num_cols_f);
EXPECT_EQ(bs_diagonal->cols.size(), num_col_blocks_f);
EXPECT_EQ(bs_diagonal->rows.size(), num_col_blocks_f);
Matrix EF;
A_->ToDenseMatrix(&EF);
const auto F = EF.topRightCorner(num_rows, num_cols_f);
Matrix expected_FtF = F.transpose() * F;
Matrix actual_FtF;
block_diagonal_ff->ToDenseMatrix(&actual_FtF);
// FtF might be not block-diagonal
auto bs = down_cast<BlockSparseMatrix*>(A_.get())->block_structure();
for (int i = 0; i < num_col_blocks_f; ++i) {
const auto col_block_f = bs->cols[num_col_blocks_e + i];
const int block_size = col_block_f.size;
const int block_pos = col_block_f.position - num_cols_e;
const auto cell_expected =
expected_FtF.block(block_pos, block_pos, block_size, block_size);
auto cell_actual =
actual_FtF.block(block_pos, block_pos, block_size, block_size);
cell_actual -= cell_expected;
EXPECT_NEAR(cell_actual.norm(), 0., kEpsilon);
}
// There should be nothing remaining outside block-diagonal
EXPECT_NEAR(actual_FtF.norm(), 0., kEpsilon);
}
TEST_P(PartitionedMatrixViewTest, BlockDiagonalEtE) {
std::unique_ptr<BlockSparseMatrix> block_diagonal_ee(
pmv_->CreateBlockDiagonalEtE());
const CompressedRowBlockStructure* bs = block_diagonal_ee->block_structure();
const int num_rows = pmv_->num_rows();
const int num_cols_e = pmv_->num_cols_e();
const int num_col_blocks_e = pmv_->num_col_blocks_e();
CHECK_EQ(block_diagonal_ee->num_rows(), num_cols_e);
CHECK_EQ(block_diagonal_ee->num_cols(), num_cols_e);
EXPECT_EQ(bs->cols.size(), num_col_blocks_e);
EXPECT_EQ(bs->rows.size(), num_col_blocks_e);
Matrix EF;
A_->ToDenseMatrix(&EF);
const auto E = EF.topLeftCorner(num_rows, num_cols_e);
Matrix expected_EtE = E.transpose() * E;
Matrix actual_EtE;
block_diagonal_ee->ToDenseMatrix(&actual_EtE);
EXPECT_NEAR((expected_EtE - actual_EtE).norm(), 0., kEpsilon);
}
TEST_P(PartitionedMatrixViewTest, UpdateBlockDiagonalEtE) {
std::unique_ptr<BlockSparseMatrix> block_diagonal_ete(
pmv_->CreateBlockDiagonalEtE());
const CompressedRowBlockStructure* bs = block_diagonal_ete->block_structure();
const int num_cols = pmv_->num_cols_e();
Matrix multi_threaded(num_cols, num_cols);
pmv_->UpdateBlockDiagonalEtE(block_diagonal_ete.get());
block_diagonal_ete->ToDenseMatrix(&multi_threaded);
Matrix single_threaded(num_cols, num_cols);
pmv_single_threaded_->UpdateBlockDiagonalEtE(block_diagonal_ete.get());
block_diagonal_ete->ToDenseMatrix(&single_threaded);
EXPECT_NEAR((multi_threaded - single_threaded).norm(), 0., kEpsilon);
}
TEST_P(PartitionedMatrixViewTest, UpdateBlockDiagonalFtF) {
std::unique_ptr<BlockSparseMatrix> block_diagonal_ftf(
pmv_->CreateBlockDiagonalFtF());
const CompressedRowBlockStructure* bs = block_diagonal_ftf->block_structure();
const int num_cols = pmv_->num_cols_f();
Matrix multi_threaded(num_cols, num_cols);
pmv_->UpdateBlockDiagonalFtF(block_diagonal_ftf.get());
block_diagonal_ftf->ToDenseMatrix(&multi_threaded);
Matrix single_threaded(num_cols, num_cols);
pmv_single_threaded_->UpdateBlockDiagonalFtF(block_diagonal_ftf.get());
block_diagonal_ftf->ToDenseMatrix(&single_threaded);
EXPECT_NEAR((multi_threaded - single_threaded).norm(), 0., kEpsilon);
}
INSTANTIATE_TEST_SUITE_P(
ParallelProducts,
PartitionedMatrixViewSpMVTest,
PartitionedMatrixViewTest,
::testing::Combine(::testing::Values(2, 4, 6),
::testing::Values(1, 2, 3, 4, 5, 6, 7, 8)),
ParamInfoToString);