Commit 370f8cca authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Use a block-diagonal matrix for evaluation matrix blocks

The blocks of the evaluation matrix are sparse.  In particular,
they have a block-diagonal structure.  Previously, the implementation
used a FieldMatrix for these blocks.  Which, in theory, is wasteful
because of the matrix sparsity.

This commit introduces a block-diagonal matrix with compile-time size
for this.

Unfortunately, initial tests do not show any speedup.
parent e3519fc9
Pipeline #8458 passed with stage
in 3 minutes and 50 seconds
......@@ -24,4 +24,5 @@ install(FILES
damagedelasticenergydensity_nondecomposed.hh
ductilefracturejet.hh
smallstrainbrittlefracturejetevaluator.hh
unblockeddiagonalmatrix.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/fracture-phasefields)
......@@ -3,6 +3,8 @@
#ifndef DUNE_FRACTUREPHASEFIELDS_GRADIENTTOSTRAINVALUEMAP_HH
#define DUNE_FRACTUREPHASEFIELDS_GRADIENTTOSTRAINVALUEMAP_HH
#include <dune/fracture-phasefields/unblockeddiagonalmatrix.hh>
namespace Dune::FracturePhasefields
{
......@@ -55,29 +57,30 @@ namespace Dune::FracturePhasefields
// Type of entries for transformed matrix
// A strain is a vector of four entries, and there are three shape functions
using block_type = std::conditional_t<transposed,
FieldMatrix<field_type, dim+1, dim*(dim+1)/2 + 1>,
FieldMatrix<field_type, dim*(dim+1)/2 + 1, dim+1> >;
using StrainMatrix = std::conditional_t<transposed,
FieldMatrix<field_type, dim, dim*(dim+1)/2>,
FieldMatrix<field_type, dim*(dim+1)/2, dim> >;
using block_type = UnblockedDiagonalMatrix<StrainMatrix,FieldMatrix<field_type,1,1> >;
/** \brief Transform (virtualRow,virtualCol)-th copy of (row,col)-th entry A original matrix
/** \brief Transform (virtualRow,virtualCol)-th copy of (row,col)-th entry A of original matrix
*/
block_type apply(const typename Matrix::block_type& A, int row, int col, int virtualRow, int virtualCol) const
{
block_type result(0);
StrainMatrix strain;
// Compute strain in i-th direction
for (unsigned int i=0; i<dim; i++)
{
for (unsigned int j=0; j<dim; j++)
if (transposed)
result[i][linearIndex(i,j)] = ((i==j) ? 1.0 : 0.5) * A[0][j];
strain[i][linearIndex(i,j)] = ((i==j) ? 1.0 : 0.5) * A[0][j];
else
result[linearIndex(i,j)][i] = ((i==j) ? 1.0 : 0.5) * A[j][0];
strain[linearIndex(i,j)][i] = ((i==j) ? 1.0 : 0.5) * A[j][0];
}
result[result.rows-1][result.cols-1] = (transposed) ? A[0][A.cols-1] : A[A.rows-1][0];
FieldMatrix<field_type,1,1> value = (transposed) ? A[0][A.cols-1] : A[A.rows-1][0];
return result;
return block_type(strain, value);
}
};
......
......@@ -8,11 +8,19 @@
#include <dune/common/exceptions.hh>
#include <dune/common/typeutilities.hh>
#include <dune/solvers/common/interval.hh>
#include <dune/solvers/common/resize.hh>
namespace Dune {
namespace FracturePhasefields {
namespace Impl {
// Helper to test whether T is of type UnblockedDiagonal
template<typename T>
using isUnblockedDiagonalMatrix_t = decltype( &T::isUnblockedDiagonalMatrix );
}
/** \brief Restriction of an integral functional to a low-dimensional axis-aligned subspace
*
* In other words, if the functional is block-separable we restrict to a single block.
......@@ -113,8 +121,6 @@ public:
BT_ij.umtv(v, jet_j);
auto g_j = densityGradient(jet_j, block);
// Multiply with the entire gradient, even though we know that some entries are zero.
// I think exploiting the zeros here would complicate the code too much.
BT_ij.usmv(weights[j], g_j, gradient);
}
return gradient;
......@@ -137,11 +143,16 @@ public:
Hj *= weights[j];
// compute BT_ij * Hj * (BT_ij)^T
for (auto&& [BT_ij_k, k] : sparseRange(BT_ij))
for (auto&& [BT_ij_l, l] : sparseRange(BT_ij))
for (auto&& [Hj_r, r] : sparseRange(Hj))
for (auto&& [Hj_rs, s] : sparseRange(Hj_r))
result[k][l] += BT_ij_k[r] * Hj_rs * BT_ij_l[s];
if constexpr (Std::is_detected_v<Impl::isUnblockedDiagonalMatrix_t, decltype(BT_ij)>)
BT_ij.leftRightMultiply(Hj, result);
else
{
for (auto&& [BT_ij_k, k] : sparseRange(BT_ij))
for (auto&& [BT_ij_l, l] : sparseRange(BT_ij))
for (auto&& [Hj_r, r] : sparseRange(Hj))
for (auto&& [Hj_rs, s] : sparseRange(Hj_r))
result[k][l] += BT_ij_k[r] * Hj_rs * BT_ij_l[s];
}
}
}
......@@ -166,22 +177,28 @@ public:
auto jet_j = f_->jet()[j];
BT_ij.umtv(v, jet_j);
auto rowRange = densityHessian.getBlockRange(row);
auto colRange = densityHessian.getBlockRange(col);
auto Hj = densityHessian(jet_j, row, col);
Hj *= weights[j];
// compute BT_ij * Hj * (BT_ij)^T
for (size_t k : rowRange)
if constexpr (Std::is_detected_v<Impl::isUnblockedDiagonalMatrix_t, decltype(BT_ij)>)
BT_ij.leftRightMultiply(Hj, result, row, col);
else
{
for (size_t l : colRange)
{
auto&& BT_ij_k = BT_ij[k];
auto&& BT_ij_l = BT_ij[l];
auto rowRange = densityHessian.getBlockRange(row);
auto colRange = densityHessian.getBlockRange(col);
for (auto&& [Hj_r, r] : sparseRange(Hj))
for (auto&& [Hj_rs, s] : sparseRange(Hj_r))
result[k][l] += BT_ij_k[r] * Hj_rs * BT_ij_l[s];
for (size_t k : rowRange)
{
for (size_t l : colRange)
{
auto&& BT_ij_k = BT_ij[k];
auto&& BT_ij_l = BT_ij[l];
for (auto&& [Hj_r, r] : sparseRange(Hj))
for (auto&& [Hj_rs, s] : sparseRange(Hj_r))
result[k][l] += BT_ij_k[r] * Hj_rs * BT_ij_l[s];
}
}
}
}
......
......@@ -5,6 +5,8 @@
#include <dune/common/rangeutilities.hh>
#include <dune/fracture-phasefields/integralfunctional.hh>
namespace Dune::FracturePhasefields {
/**
......@@ -174,10 +176,15 @@ namespace Dune::FracturePhasefields {
// Compute matrix entry in upper right triangle part of Y
if (i<=j)
{
for (int ii=0; ii<Y_ij.rows; ++ii)
for (int jj=0; jj<Y_ij.cols; ++jj)
for (int l=0; l<M; ++l)
Y_ij[ii][jj] += precomp[l][ii] * rowB_j[l][jj];
if constexpr (Std::is_detected_v<Impl::isUnblockedDiagonalMatrix_t, decltype(rowB_j)>)
rowB_j.leftMultiplyTransposed(precomp, Y_ij);
else
{
for (int ii=0; ii<Y_ij.rows; ++ii)
for (int jj=0; jj<Y_ij.cols; ++jj)
for (int l=0; l<M; ++l)
Y_ij[ii][jj] += precomp[l][ii] * rowB_j[l][jj];
}
}
else
// Do not compute matrix entry in left triangle part of Y.
......
......@@ -7,3 +7,5 @@ dune_add_test(SOURCES damagedductileelasticenergydensitytest.cc)
dune_add_test(SOURCES integralfunctionaltest.cc)
dune_add_test(SOURCES largestraindamagedelasticenergydensitytest.cc)
dune_add_test(SOURCES unblockeddiagonalmatrixtest.cc)
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include "config.h"
#include <iostream>
#include <dune/common/exceptions.hh>
#include <dune/common/fmatrix.hh>
#include <dune/matrix-vector/transpose.hh>
#include <dune/fracture-phasefields/unblockeddiagonalmatrix.hh>
using namespace Dune;
template <class TestMatrix, class ReferenceMatrix>
bool compareBehavior(TestMatrix testMatrix,
ReferenceMatrix referenceMatrix,
const std::vector<size_t>& rowBlocks,
const std::vector<size_t>& colBlocks)
{
bool passed = true;
// Test whether the matrix dimensions match
if (testMatrix.rows != referenceMatrix.rows)
DUNE_THROW(Exception, "testMatrix does not have the correct number of rows!");
if (testMatrix.cols != referenceMatrix.cols)
DUNE_THROW(Exception, "testMatrix does not have the correct number of columns!");
// Test whether method umv computes the right thing
{
FieldVector<double, testMatrix.cols> v;
for (size_t i=0; i<v.size(); ++i)
{
v = 0;
v[i] = 1.0;
FieldVector<double, testMatrix.rows> testResult(0);
FieldVector<double, referenceMatrix.rows> referenceResult(0);
testMatrix.umv(v,testResult);
referenceMatrix.umv(v,referenceResult);
if ( (testResult - referenceResult).infinity_norm() > 1e-6 )
{
std::cout << "Method umv did not compute the correct result!" << std::endl;
std::cout << "Result: " << testResult << std::endl;
std::cout << "Expected: " << referenceResult << std::endl;
passed = false;
}
}
}
// Test whether method usmv computes the right thing
{
FieldVector<double, testMatrix.cols> v;
for (size_t i=0; i<v.size(); ++i)
{
v = 0;
v[i] = 1.0;
FieldVector<double, testMatrix.rows> testResult(0);
FieldVector<double, referenceMatrix.rows> referenceResult(0);
testMatrix.usmv(0.5, v,testResult);
referenceMatrix.usmv(0.5, v,referenceResult);
if ( (testResult - referenceResult).infinity_norm() > 1e-6 )
{
std::cout << "Method usmv did not compute the correct result!" << std::endl;
std::cout << "Result: " << testResult << std::endl;
std::cout << "Expected: " << referenceResult << std::endl;
passed = false;
}
}
}
// Test whether method umtv computes the right thing
{
FieldVector<double, testMatrix.rows> v;
for (size_t i=0; i<v.size(); ++i)
{
v = 0;
v[i] = 1.0;
FieldVector<double, testMatrix.cols> testResult(0);
FieldVector<double, referenceMatrix.cols> referenceResult(0);
testMatrix.umtv(v,testResult);
referenceMatrix.umtv(v,referenceResult);
if ( (testResult - referenceResult).infinity_norm() > 1e-6 )
{
std::cout << "Method umtv did not compute the correct result!" << std::endl;
std::cout << "Result: " << testResult << std::endl;
std::cout << "Expected: " << referenceResult << std::endl;
passed = false;
}
}
}
// Test whether operator* computes the right thing
{
FieldMatrix<double, 2, testMatrix.rows> v;
for (size_t i=0; i<v.rows; ++i)
for (size_t j=0; j<v.cols; ++j)
{
v = 0;
v[i][j] = 1.0;
auto testResult = v * testMatrix;
auto referenceResult = v * referenceMatrix;
if ( (testResult - referenceResult).infinity_norm() > 1e-6 )
{
std::cout << "Method operator* did not compute the correct result!" << std::endl;
std::cout << "Result: " << testResult << std::endl;
std::cout << "Expected: " << referenceResult << std::endl;
passed = false;
}
}
}
// Test whether leftMultiplyTransposed computes the right thing
{
FieldMatrix<double, testMatrix.rows, 2> v;
for (size_t i=0; i<v.rows; ++i)
{
for (size_t j=0; j<v.cols; ++j)
{
v = 0;
v[i][j] = 1.0;
// Yes, init testResult with 1!
// The method leftMultiplyTransposed is supposed to add its result,
// and we want to test whether it really does that.
FieldMatrix<double,2, testMatrix.cols> testResult(1);
testMatrix.leftMultiplyTransposed(v, testResult);
FieldMatrix<double, 2, testMatrix.rows> vTransposed;
MatrixVector::transpose(v, vTransposed);
auto referenceResult = vTransposed * referenceMatrix;
referenceResult += decltype(referenceResult)(1.0);
if ( (testResult - referenceResult).infinity_norm() > 1e-6 )
{
std::cout << "Method leftMultiplyTransposed did not compute the correct result!" << std::endl;
std::cout << "Result:\n" << testResult << std::endl;
std::cout << "Expected:\n" << referenceResult << std::endl;
passed = false;
}
}
}
}
// Test whether leftRightMultiply computes the right thing
{
FieldMatrix<double, referenceMatrix.cols, referenceMatrix.rows> referenceMatrixTransposed;
MatrixVector::transpose(referenceMatrix, referenceMatrixTransposed);
FieldMatrix<double, testMatrix.cols, testMatrix.cols> v;
for (size_t i=0; i<v.rows; ++i)
{
for (size_t j=0; j<v.cols; ++j)
{
v = 0;
v[i][j] = 1.0;
// Yes, init testResult with 1!
// The method leftRightMultiply is supposed to add its result,
// and we want to test whether it really does that.
FieldMatrix<double,testMatrix.rows, testMatrix.rows> testResult(1);
testMatrix.leftRightMultiply(v, testResult);
auto referenceResult = referenceMatrix * v * referenceMatrixTransposed;
referenceResult += decltype(referenceResult)(1.0);
if ( (testResult - referenceResult).infinity_norm() > 1e-6 )
{
std::cout << "Method leftRightMultiply did not compute the correct result!" << std::endl;
std::cout << "Result:\n" << testResult << std::endl;
std::cout << "Expected:\n" << referenceResult << std::endl;
passed = false;
}
}
}
}
// Test whether leftRightMultiply for individual blocks computes the right thing
{
FieldMatrix<double, referenceMatrix.cols, referenceMatrix.rows> referenceMatrixTransposed;
MatrixVector::transpose(referenceMatrix, referenceMatrixTransposed);
FieldMatrix<double, testMatrix.cols, testMatrix.cols> v;
// Loop over all matrix blocks that leftRightMultiply can compute.
for (size_t blockRow=0; blockRow<rowBlocks.size()-1; ++blockRow)
{
for (size_t blockCol=0; blockCol<colBlocks.size()-1; ++blockCol)
{
// Create test matrices
for (size_t i=0; i<v.rows; ++i)
{
for (size_t j=0; j<v.cols; ++j)
{
v = 0;
v[i][j] = 1.0;
// Yes, init testResult with 1!
// The method leftRightMultiply is supposed to add its result,
// and we want to test whether it really does that.
FieldMatrix<double,testMatrix.rows, testMatrix.rows> testResult(1);
testMatrix.leftRightMultiply(v, testResult, blockRow, blockCol);
auto referenceResult = referenceMatrix * v * referenceMatrixTransposed;
referenceResult += decltype(referenceResult)(1.0);
// Check whether result is correct
for (size_t ii=0; ii<referenceResult.rows; ++ii)
{
for (size_t jj=0; jj<referenceResult.cols; ++jj)
{
double expectedValue = (ii>=rowBlocks[blockRow] && ii<rowBlocks[blockRow+1] && jj>=rowBlocks[blockCol] && jj<rowBlocks[blockCol+1]) ? referenceResult[ii][jj] : 1.0;
if ( std::abs(testResult[ii][jj] - expectedValue) > 1e-6 )
{
std::cout << "Method leftRightMultiply for block " << blockRow << ", " << blockCol
<< " did not compute the correct result!" << std::endl;
std::cout << "Result:\n" << testResult << std::endl;
std::cout << "Expected:\n" << referenceResult << std::endl;
passed = false;
}
}
}
}
}
}
}
}
return passed;
}
int main(int argc, char* argv[]) try
{
bool passed = true;
// Create an UnblockedDiagonalMatrix for testing
FieldMatrix<double,1,2> block1 = {{2.0,1.0}};
FieldMatrix<double,1,1> block2 = {{5}};
FieldMatrix<double,2,2> block3 = {{3.0, 4.0},
{7.0, 8.0}};
FracturePhasefields::UnblockedDiagonalMatrix diagonalMatrix(block1,block2,block3);
// Create the same matrix as a FieldMatrix for comparison
FieldMatrix<double,4,5> referenceMatrix = {{2.0, 1.0, 0.0, 0.0, 0.0},
{0.0, 0.0, 5.0, 0.0, 0.0},
{0.0, 0.0, 0.0, 3.0, 4.0},
{0.0, 0.0, 0.0, 7.0, 8.0}};
passed = passed && compareBehavior(diagonalMatrix, referenceMatrix, {0,1,2,4}, {0,2,3,5});
return (passed) ? 0 : 1;
} catch (Exception& e)
{
std::cout << e.what() << std::endl;
return 1;
}
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_FRACTUREPHASEFIELDS_UNBLOCKEDDIAGONALMATRIX_HH
#define DUNE_FRACTUREPHASEFIELDS_UNBLOCKEDDIAGONALMATRIX_HH
#include <dune/common/hybridutilities.hh>
#include <dune/common/std/type_traits.hh>
namespace Dune::FracturePhasefields
{
/** \brief A block-diagonal matrix that hides the blocking structure
*
* \tparam Args The types of the diagonal blocks
*/
template <typename... Args>
class UnblockedDiagonalMatrix
{
std::tuple<Args...> blocks_;
template <unsigned int i>
static constexpr unsigned int subRows()
{
return std::tuple_element_t<i, std::tuple<Args...> >::rows;
}
template <unsigned int i>
static constexpr unsigned int subCols()
{
return std::tuple_element_t<i, std::tuple<Args...> >::cols;
}
public:
/** \brief Number of matrix rows */
static constexpr unsigned int rows = (Args::rows + ...);
/** \brief Number of matrix columns */
static constexpr unsigned int cols = (Args::cols + ...);
/** \brief The type used for scalars
*
* Use the `std::common_type_t` of the `Args`' field_type and use `nonesuch` if no implementation of
* `std::common_type` is provided for the given `field_type` arguments.
*/
using field_type = Std::detected_t<std::common_type_t, typename FieldTraits< std::decay_t<Args> >::field_type...>;
using size_type = std::size_t;
/** \brief Constructor taking a shape function gradient and a value
*/
UnblockedDiagonalMatrix(Args... args)
: blocks_(args...)
{}
/** \brief Dummy method that allows to test whether a class is an UnblockedDiagonalMatrix
*
* The UnblockedDiagonalMatrix class has more methods than a standard FieldMatrix.
* The calling code should use them if they exist. An ideal code would test for precisely
* the methods it needs, and call them if they exist.
*
* My C++ knowledge and time is limited, and therefore I instead test for this tag method.
* Feel free to replace this by something better.
*/
static bool isUnblockedDiagonalMatrix()
{
return true;
}
// Matrix-vector products
template<class X, class Y>
void umv(const X& x, Y& y) const
{
auto yPtr = y.begin();
auto xPtr = x.begin();
// Loop over all blocks
Hybrid::forEach(blocks_,
[&](auto block)
{
for (size_type i=0; i<block.rows; ++i)
for (size_type j=0; j<block.cols; j++)
yPtr[i] += block[i][j] * xPtr[j];
yPtr += block.rows;
xPtr += block.cols;
});
}
template<class X, class Y>
void umtv(const X& x, Y& y) const
{
auto yPtr = y.begin();
auto xPtr = x.begin();
// Loop over all blocks
Hybrid::forEach(blocks_,
[&](auto block)
{
for (size_type i=0; i<block.rows; ++i)
for (size_type j=0; j<block.cols; j++)
yPtr[j] += block[i][j] * xPtr[i];
xPtr += block.rows;
yPtr += block.cols;
});
}
//! y += alpha A x
template<class X, class Y, class F>
void usmv (F&& alpha, const X& x, Y& y) const
{
auto yPtr = y.begin();
auto xPtr = x.begin();
// Loop over all blocks
Hybrid::forEach(blocks_,
[&](auto block)
{
for (size_type i=0; i<block.rows; ++i)
for (size_type j=0; j<block.cols; j++)
yPtr[i] += alpha * block[i][j] * xPtr[j];