28 Numerics library [numerics]

28.9 Basic linear algebra algorithms [linalg]

28.9.6 Layouts for packed matrix types [linalg.layout.packed]

28.9.6.1 Overview [linalg.layout.packed.overview]

layout_blas_packed is an mdspan layout mapping policy that represents a square matrix that stores only the entries in one triangle, in a packed contiguous format.

Its Triangle template parameter determines whether an mdspan with this layout stores the upper or lower triangle of the matrix.

Its StorageOrder template parameter determines whether the layout packs the matrix's elements in column-major or row-major order.

A StorageOrder of column_major_t indicates column-major ordering.

This packs matrix elements starting with the leftmost (least column index) column, and proceeding column by column, from the top entry (least row index).

A StorageOrder of row_major_t indicates row-major ordering.

This packs matrix elements starting with the topmost (least row index) row, and proceeding row by row, from the leftmost (least column index) entry.

[Note 1:

layout_blas_packed describes the data layout used by the BLAS' Symmetric Packed (SP), Hermitian Packed (HP), and Triangular Packed (TP) matrix types.

— end note]

namespace std::linalg { template<class Triangle, class StorageOrder> class layout_blas_packed { public: using triangle_type = Triangle; using storage_order_type = StorageOrder; template<class Extents> struct mapping { public: using extents_type = Extents; using index_type = typename extents_type::index_type; using size_type = typename extents_type::size_type; using rank_type = typename extents_type::rank_type; using layout_type = layout_blas_packed; // [linalg.layout.packed.cons], constructors constexpr mapping() noexcept = default; constexpr mapping(const mapping&) noexcept = default; constexpr mapping(const extents_type&) noexcept; template<class OtherExtents> constexpr explicit(!is_convertible_v<OtherExtents, extents_type>) mapping(const mapping<OtherExtents>& other) noexcept; constexpr mapping& operator=(const mapping&) noexcept = default; // [linalg.layout.packed.obs], observers constexpr const extents_type& extents() const noexcept { return extents_; } constexpr index_type required_span_size() const noexcept; template<class Index0, class Index1> constexpr index_type operator() (Index0 ind0, Index1 ind1) const noexcept; static constexpr bool is_always_unique() noexcept { return (extents_type::static_extent(0) != dynamic_extent && extents_type::static_extent(0) < 2) || (extents_type::static_extent(1) != dynamic_extent && extents_type::static_extent(1) < 2); } static constexpr bool is_always_exhaustive() noexcept { return true; } static constexpr bool is_always_strided() noexcept { return is_always_unique(); } constexpr bool is_unique() const noexcept { return extents_.extent(0) < 2; } constexpr bool is_exhaustive() const noexcept { return true; } constexpr bool is_strided() const noexcept { return extents_.extent(0) < 2; } constexpr index_type stride(rank_type) const noexcept; template<class OtherExtents> friend constexpr bool operator==(const mapping&, const mapping<OtherExtents>&) noexcept; private: extents_type extents_{}; // exposition only }; }; }

Mandates:

(5.1)
Triangle is either upper_triangle_t or lower_triangle_t,
(5.2)
StorageOrder is either column_major_t or row_major_t,
(5.3)
Extents is a specialization of std::extents,
(5.4)
Extents::rank() equals 2,
(5.5)
one of extents_type::static_extent(0) == dynamic_extent, extents_type::static_extent(1) == dynamic_extent, or extents_type::static_extent(0) == extents_type::static_extent(1) is true, and
(5.6)
if Extents::rank_dynamic() == 0 is true, let $N_{s}$ be equal to Extents::static_extent(0); then, $N_{s} \times (N_{s} + 1)$ is representable as a value of type index_type.

layout_blas_packed<T, SO>::mapping<E> is a trivially copyable type that models regular for each T, SO, and E.

28.9.6.2 Constructors [linalg.layout.packed.cons]

🔗

constexpr mapping(const extents_type& e) noexcept;

Preconditions:

(1.1)
Let N be equal to e.extent(0).

Then, $N \times (N + 1)$ is representable as a value of type index_type ([basic.fundamental]).
(1.2)
e.extent(0) equals e.extent(1).

Effects: Direct-non-list-initializes extents_ with e.

🔗

template<class OtherExtents>
  explicit(!is_convertible_v<OtherExtents, extents_type>)
    constexpr mapping(const mapping<OtherExtents>& other) noexcept;

Constraints: is_constructible_v<extents_type, OtherExtents> is true.

Preconditions: Let N be other.extents().extent(0).

Then,

N \times (N + 1)

is representable as a value of type index_type ([basic.fundamental]).

Effects: Direct-non-list-initializes extents_ with other.extents().

28.9.6.3 Observers [linalg.layout.packed.obs]

🔗

constexpr index_type required_span_size() const noexcept;

Returns: extents_.extent(0) * (extents_.extent(0) + 1)/2.

[Note 1:

For example, a 5 x 5 packed matrix only stores 15 matrix elements.

— end note]

🔗

template<class Index0, class Index1>
  constexpr index_type operator() (Index0 ind0, Index1 ind1) const noexcept;

Constraints:

(2.1)
is_convertible_v<Index0, index_type> is true,
(2.2)
is_convertible_v<Index1, index_type> is true,
(2.3)
is_nothrow_constructible_v<index_type, Index0> is true, and
(2.4)
is_nothrow_constructible_v<index_type, Index1> is true.

Let i be extents_type::index-cast(ind0), and let j be extents_type::index-cast(ind1).

Preconditions: i, j is a multidimensional index in extents_ ([mdspan.overview]).

Returns: Let N be extents_.extent(0).

Then

(5.1)
(*this)(j, i) if i > j is true; otherwise
(5.2)
i + j * (j + 1)/2 if is_same_v<StorageOrder, column_major_t> && is_same_v<Triangle, upper_triangle_t> is true or is_same_v<StorageOrder, row_major_t> && is_same_v<Triangle, lower_triangle_t> is true; otherwise
(5.3)
j + N * i - i * (i + 1)/2.

🔗

constexpr index_type stride(rank_type r) const noexcept;

Preconditions:

(6.1)
is_strided() is true, and
(6.2)
r < extents_type::rank() is true.

Returns: 1.

🔗

template<class OtherExtents>
  friend constexpr bool operator==(const mapping& x, const mapping<OtherExtents>& y) noexcept;

Effects: Equivalent to: return x.extents() == y.extents();