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

Use if-constexpr more in Tensor. #16468

Merged
merged 2 commits into from
Jan 15, 2024
Merged

Use if-constexpr more in Tensor. #16468

Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
1b55d98
Use if-constexpr more in Tensor.
bangerth Jan 13, 2024
7fcdb8b
Peel off first element of the sequence.
bangerth Jan 14, 2024
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
123 changes: 45 additions & 78 deletions include/deal.II/base/tensor.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1406,43 +1406,18 @@ Tensor<rank_, dim, Number>::Tensor(Tensor<rank_, dim, Number> &&other) noexcept
}
# endif

namespace internal
{
namespace TensorSubscriptor
{
template <typename ArrayElementType, int dim>
constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE ArrayElementType &
subscript(ArrayElementType *values,
const unsigned int i,
std::integral_constant<int, dim>)
{
AssertIndexRange(i, dim);
return values[i];
}

template <typename ArrayElementType>
constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE ArrayElementType &
subscript(ArrayElementType *dummy,
const unsigned int,
std::integral_constant<int, 0>)
{
Assert(
false,
ExcMessage(
"Cannot access elements of an object of type Tensor<rank,0,Number>."));
return *dummy;
}
} // namespace TensorSubscriptor
} // namespace internal


template <int rank_, int dim, typename Number>
constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE
typename Tensor<rank_, dim, Number>::value_type &
Tensor<rank_, dim, Number>::operator[](const unsigned int i)
{
return dealii::internal::TensorSubscriptor::subscript(
values, i, std::integral_constant<int, dim>());
Assert(dim != 0,
ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
AssertIndexRange(i, dim);

return values[i];
Comment on lines +1416 to +1420
Copy link
Member

Choose a reason for hiding this comment

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

I see no constexpr in this function; given that the previous code appears to have been a workaround (introduced by #1573, see there for the message it tries to fix) to silence compiler warnings, what prevents these warnings now? Or don't we support the old compilers any more that would trigger the warnings?

Copy link
Member Author

Choose a reason for hiding this comment

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

Good question. I didn't investigate the history when I wrote this patch. I'd say that if the CI checks come back clean, we should assume that compilers have learned not to warn.

(No constexpr -- I accidentally used the commit message I had used for a bigger patch that had test problems also for the case I here that I decided to break out into its own patch. I'll fix the SUBJECT line.)

Copy link
Member Author

Choose a reason for hiding this comment

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

And just to be clear, we put the workaround in #1573 in more than 8 years ago. I don't think we still support many of the compilers we used at the time.

}


Expand All @@ -1451,6 +1426,8 @@ constexpr DEAL_II_ALWAYS_INLINE
DEAL_II_HOST_DEVICE const typename Tensor<rank_, dim, Number>::value_type &
Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
{
Assert(dim != 0,
ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
AssertIndexRange(i, dim);

return values[i];
Expand Down Expand Up @@ -1646,58 +1623,30 @@ constexpr inline DEAL_II_ALWAYS_INLINE
}


namespace internal

template <int rank_, int dim, typename Number>
template <typename OtherNumber>
constexpr inline DEAL_II_ALWAYS_INLINE
DEAL_II_HOST_DEVICE Tensor<rank_, dim, Number> &
Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
{
namespace TensorImplementation
{
template <int rank,
int dim,
typename Number,
typename OtherNumber,
std::enable_if_t<
!std::is_integral<
typename ProductType<Number, OtherNumber>::type>::value &&
!std::is_same_v<Number, Differentiation::SD::Expression>,
int> = 0>
constexpr DEAL_II_HOST_DEVICE inline DEAL_II_ALWAYS_INLINE void
division_operator(Tensor<rank, dim, Number> (&t)[dim],
const OtherNumber &factor)
if constexpr (std::is_integral<
Copy link
Member Author

Choose a reason for hiding this comment

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

Oh, I do have a constexpr here -- the SUBJECT line wasn't totally wrong ;-)

typename ProductType<Number, OtherNumber>::type>::value ||
std::is_same_v<Number, Differentiation::SD::Expression>)
{
const Number inverse_factor = Number(1.) / factor;
// recurse over the base objects
for (unsigned int d = 0; d < dim; ++d)
t[d] *= inverse_factor;
values[d] /= s;
}


template <int rank,
int dim,
typename Number,
typename OtherNumber,
std::enable_if_t<
std::is_integral<
typename ProductType<Number, OtherNumber>::type>::value ||
std::is_same_v<Number, Differentiation::SD::Expression>,
int> = 0>
constexpr DEAL_II_HOST_DEVICE inline DEAL_II_ALWAYS_INLINE void
division_operator(Tensor<rank, dim, Number> (&t)[dim],
const OtherNumber &factor)
else
{
// recurse over the base objects
// If we can, avoid division by multiplying by the inverse of the given
// factor:
const Number inverse_factor = Number(1.) / s;
for (unsigned int d = 0; d < dim; ++d)
t[d] /= factor;
values[d] *= inverse_factor;
}
} // namespace TensorImplementation
} // namespace internal


template <int rank_, int dim, typename Number>
template <typename OtherNumber>
constexpr inline DEAL_II_ALWAYS_INLINE
DEAL_II_HOST_DEVICE Tensor<rank_, dim, Number> &
Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
{
internal::TensorImplementation::division_operator(values, s);
return *this;
}

Expand Down Expand Up @@ -1749,12 +1698,30 @@ constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE
typename numbers::NumberTraits<Number>::real_type
Tensor<rank_, dim, Number>::norm_square() const
{
typename numbers::NumberTraits<Number>::real_type s = internal::NumberType<
typename numbers::NumberTraits<Number>::real_type>::value(0.0);
for (unsigned int i = 0; i < dim; ++i)
s += values[i].norm_square();
if constexpr (rank_ == 1)
{
// For rank-1 tensors, the square of the norm is simply the sum of
// squares of the elements:
typename numbers::NumberTraits<Number>::real_type s =
internal::NumberType<
typename numbers::NumberTraits<Number>::real_type>::value(0.0);
for (unsigned int i = 0; i < dim; ++i)
s += numbers::NumberTraits<Number>::abs_square(values[i]);
Copy link
Member

Choose a reason for hiding this comment

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

I'm usually preferring to peel off the first loop index because it avoids adding into zero (and compilers would need -ffast-math to avoid the instruction); in this case, it wouldn't even introduce more lines of code because we would write

Suggested change
typename numbers::NumberTraits<Number>::real_type s =
internal::NumberType<
typename numbers::NumberTraits<Number>::real_type>::value(0.0);
for (unsigned int i = 0; i < dim; ++i)
s += numbers::NumberTraits<Number>::abs_square(values[i]);
typename numbers::NumberTraits<Number>::real_type s =
numbers::NumberTraits<Number>::abs_square(values[0]);
for (unsigned int i = 1; i < dim; ++i)
s += numbers::NumberTraits<Number>::abs_square(values[i]);

The question that remains for me is if would trigger a warning for dim = 0. I do not feel particularly strongly about it, so either one is fine for me.

Copy link
Member Author

Choose a reason for hiding this comment

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

This is a loop with compile-time bounds. It's almost certainly going to be unrolled. Don't you think the compiler can make the optimization itself? I think it's easier to read this way, but I also don't feel strongly about it.

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 changed my mind and did as you suggested. I handled the dim==0 case explicitly via an other if constexpr at the top.


return s;
return s;
}
else
{
// For higher-rank tensors, the square of the norm is the sum
// of squares of sub-tensors
typename numbers::NumberTraits<Number>::real_type s =
internal::NumberType<
typename numbers::NumberTraits<Number>::real_type>::value(0.0);
for (unsigned int i = 0; i < dim; ++i)
s += values[i].norm_square();

return s;
}
}


Expand Down