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 second version of TensorProductMatrixSymmetricSum::apply_inverse() #14178

Merged
merged 1 commit into from
Aug 23, 2022
Merged

Add second version of TensorProductMatrixSymmetricSum::apply_inverse() #14178

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
31 changes: 30 additions & 1 deletion include/deal.II/lac/tensor_product_matrix.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,11 +116,28 @@ class TensorProductMatrixSymmetricSumBase
* described in the main documentation of TensorProductMatrixSymmetricSum.
* This function is operating on ArrayView to allow checks of
* array bounds with respect to @p dst and @p src.
*
* @warning This function works on an internal temporal array, leading to
* increased memory consumption if many instances of this class are created,
* e.g., a different object on every cell with different underlying
* coefficients each. Furthermore, only one thread run this function at once
* (ensured internally with a mutex). If these two limitations are an issue
* for you, please consider the other version of this function.
*/
void
apply_inverse(const ArrayView<Number> & dst,
const ArrayView<const Number> &src) const;

/**
* Same as above but the user can provide a user-owned temporal array,
* resolving the two issues described above. This array is resized
* internally to the needed size.
*/
Copy link
Member

Choose a reason for hiding this comment

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

It looks like the temporary array needs to be of length 2 * n: could you document that here?

Copy link
Member Author

Choose a reason for hiding this comment

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

The array is set to the right size internally.

Copy link
Member

Choose a reason for hiding this comment

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

Got it. Could we document that instead? We typically require that users set output arrays to the right sizes so we should clarify what to do for work arrays.

Copy link
Member Author

Choose a reason for hiding this comment

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

I have added a comment!

void
apply_inverse(const ArrayView<Number> & dst,
const ArrayView<const Number> &src,
AlignedVector<Number> & tmp) const;

/**
* Return the memory consumption of the allocated memory in this class.
*/
Expand Down Expand Up @@ -559,10 +576,22 @@ inline void
TensorProductMatrixSymmetricSumBase<dim, Number, n_rows_1d>::apply_inverse(
const ArrayView<Number> & dst_view,
const ArrayView<const Number> &src_view) const
{
std::lock_guard<std::mutex> lock(this->mutex);
this->apply_inverse(dst_view, src_view, this->tmp_array);
}



template <int dim, typename Number, int n_rows_1d>
inline void
TensorProductMatrixSymmetricSumBase<dim, Number, n_rows_1d>::apply_inverse(
const ArrayView<Number> & dst_view,
const ArrayView<const Number> &src_view,
AlignedVector<Number> & tmp_array) const
{
AssertDimension(dst_view.size(), this->n());
AssertDimension(src_view.size(), this->m());
std::lock_guard<std::mutex> lock(this->mutex);
const unsigned int n = n_rows_1d > 0 ? n_rows_1d : eigenvalues[0].size();
tmp_array.resize_fast(Utilities::fixed_power<dim>(n));
constexpr int kernel_size = n_rows_1d > 0 ? n_rows_1d : 0;
Expand Down