P3050R1: LEWG presentation

[mhoemmen]

P3050R1

Terms

• Complex number $z = x + iy$ has a real part $x$ and an imaginary part $y$.

  • In std::linalg, a "complex number" is any number type, not necessarily std::complex, where conj is ADL-findable. (Users define their own complex number types to work around various std::complex issues.)

• Conjugate of a complex number $z = x + iy$ is $x - iy$. We write this in C++ as conj(z).

• Conjugate of a noncomplex number is just the number (as its imaginary part is zero).

  • We say "noncomplex" and not "real" because "real number" carries additional, irrelevant mathematical implications. std::linalg does not require noncomplex numbers to be real numbers. They could be some other algebraic type, like 3 x 3 matrices.

• Transpose of A is a matrix B such that B[c, r] equals A[r, c].

• Conjugate transpose of A is a matrix B such that B[c, r] equals conj(A[r, c]) (the complex conjugate, that is, $x - iy$ for a given $z = x + iy$).

• BLAS (pronounced "blahz") stands for the "Basic Linear Algebra Subroutines," a Standard Fortran and C interface providing linear algebra operations; foundation of std::linalg.

Common practice: "conjugate transpose of a noncomplex matrix is just the transpose"

The conjugate transpose of a complex matrix naturally generalizes the transpose of a noncomplex matrix. Users who develop generic algorithms for either complex or noncomplex problems write the algorithm once using the conjugate transpose. BLAS and matrix-oriented programming languages like Matlab treat both using the same notation (e.g., 'C' flag means transpose for a noncomplex matrix, conjugate transpose for a complex matrix).

mdspan views conjugated, transposed, and conjugate_transposed

A key feature of linear algebra libraries is their ability to view the transpose or conjugate transpose of a matrix "in place" without actually changing its elements. Matrices may be large and copying them may be too expensive.

• BLAS implements "view (conjugate) transpose in place" with a separate flag argument: 'N', 'T', or 'C'.

• std::linalg implements this using mdspan views conjugated, transposed, and conjugate_transposed.

How does conjugated work currently?

1. If the input mdspan has accessor type conjugated_accessor<NestedAccessor>, then the result has accessor type NestedAccessor;

2. otherwise, if the input mdspan has accessor type Accessor, then the result has accessor type conjugated_accessor<Accessor>.

conjugated_accessor's (read-only) access function conjugates the element if it's a complex number, else it just returns the number.

This is correct, but can hinder optimization

The current behavior of conjugated is correct. The issue is that conjugated(A) for an mdspan-of-noncomplex-numbers A doesn't have the same accessor type as A. It should just return A!

This is an issue because both Standard Library implementations and users may want to optimize for "known accessors" such as default_accessor. Accessors communicate optimization information, like "this is a contiguous array in memory."

template<class ElementType, class IndexType, size_t Ext0, class Layout, class Accessor>
void generic_algorithm( // fully generic
  mdspan<ElementType, extent<IndexType, Ext0>, Layout, Accessor> x);

template<class ElementType, class IndexType, size_t Ext0, class Layout>
void generic_algorithm( // specialization
  mdspan<ElementType, extent<IndexType, Ext0>, Layout, default_accessor<ElementType>> x);

Currently, conjugated of a default_accessor<ElementType> mdspan has accessor conjugated_accessor<default_accessor<ElementType>>. Calling generic_algorithm with this mdspan will thus take the "generic path," rather than the specialization.

If we want to optimize the conjugated_accessor case, we have to add another specialization. This has compile-time costs. Users either have to remember to do this, or write their generic algorithms twice (once for complex and once for noncomplex).

template<class Real, class IndexType, size_t Ext0, class Layout>
  requires(not impl::is_complex_v<Real>)
void generic_algorithm( // another specialization
  mdspan<Real, extent<IndexType, Ext0>, Layout, conjugated_accessor<default_accessor<Real>>> x)
{
  // Dispatch to default_accessor specialization
  return generic_algorithm(mdspan{x.data_handle(), x.mapping(), x.accessor().nested_accessor()});
}

Fix: conjugated(A) should return A if A is noncomplex

P3050 proposes the only reasonable fix: make conjugated(A) return A if the value_type of A is not complex.

"Not complex" just means that conj is not ADL-findable for the mdspan's value_type.

Recall that mdspan has element_type (possibly cv-qualified) and value_type = remove_cvref_t<element_type>.

[mhoemmen]

I'm putting this into P3050R2.