Closed mhoemmen closed 10 months ago
transpose-extents
and _transpose-extents_t
_ changes didn't get included in PR #426. I've added that commit to the new PR #428.
(I've added the following explanation to the Design section of P1673, right above the "Future work" section.)
P1673 symmetric_*
, hermitian_*
, and triangular_*
functions always take Triangle (and DiagonalStorage, if applicable) parameters so that the functions always have the same signature, regardless of mdspan layout.
For symmetric or Hermitian algorithms, mismatching the mdspan layout's Triangle and the function's Triangle has no effect.
For triangular algorithms, mismatch has an effect that users likely would not want. (It means "use the other triangle," which is zero.) Thus, it's reasonable to make mismatch an error.
In practice, users aren't likely to encounter a triangular packed matrix in isolation. Such matrices usually occur in context of symmetric or Hermitian packed matrices. A common user error might thus be mismatching the Triangles for both symmetric or Hermitian functions, and triangular functions. The first is harmless; the second is likely an error.
Therefore, we recommend
Users aren't likely to encounter a triangular packed matrix in isolation. It generally comes as an in-place transformation (e.g., factorization) of a symmetric or Hermitian packed matrix. For example, LAPACK's DSPTRF
(Double-precision Symmetric Packed TRiangular Factorization) computes a symmetric $L D L^T$ (or $U D U^T$) factorization in place, overwriting the input symmetric packed matrix $A$. LAPACK's DSPTRS
(Double-precision Symmetric Packed TRiangular Solve) then uses the result of DSPTRF
to solve a linear system. DSPTRF
overwrites $A$ with the triangle $L$ (if $A$ uses lower triangle storage, or $U$, if $A$ uses upper triangle storage). This is an idiom for which the BLAS was designed: factorizations typically overwrite their input, and thus reinterpret its "data structure" on the fly.
For a summary of the BLAS' packed storage formats, please refer to the "Packed Storage" section of the LAPACK Users' Guide, Third Edition (1999).
BLAS routines for packed storage have only a single argument, UPLO
. This describes both whether the caller is storing the upper or lower triangle, and the triangle of the matrix on which the routine will operate. (Packed BLAS formats always store the diagonal explicitly; they don't have the analog of DiagonalStorage
.) An example of a BLAS triangular packed routine is DTPMV
, double-precision (D) triangular packed (TP) matrix-vector (MV) product.
BLAS packed formats don't represent metadata explicitly; the caller is responsible for knowing whether they are storing the upper or lower triangle. Getting the UPLO
argument wrong makes the matrix wrong. For example, suppose that the matrix is 4 x 4, and the user's array input for the matrix is [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. If the user is storing the upper triangle (in column-major order, as the Fortran BLAS requires), then the matrix looks like this.
$$ \begin{matrix} 1 & 2 & 4 & 7 \ & 3 & 5 & 8 \ & & 6 & 9 \ & & & 10 \ \end{matrix} $$
Mismatching the UPLO
argument (by passing in 'Lower triangle'
instead of 'Upper triangle'
) would result in an entirely wrong matrix -- not even the transpose. Note how the diagonal elements differ, for instance.
$$ \begin{matrix} 1 & & & \ 2 & 5 & & \ 3 & 6 & 8 & \ 4 & 7 & 9 & 10 \ \end{matrix} $$
This would be incorrect for triangular, symmetric, or Hermitian matrices.
P1673 offers packed formats that encode the Triangle
. This means that the mdspan
alone conveys the data structure. P1673 retains the function's separate Triangle
parameter so that the function's signature doesn't change based on the mdspan's layout. P1673 (up to R12 at least) requires that the function's Triangle match the mdspan's Triangle.
If P1673 were to permit mismatching the two Triangles, how would the function reasonably interpret the user's intent? For triangular matrices with explicit diagonal, mismatch would mean multiplying by or solving with a zero triangle matrix. For triangular matrices with implicit unit diagonal, mismatch would mean multiplying by or solving with a diagonal matrix of ones -- the identity matrix. Users wouldn't want to do either one of those.
For symmetric matrices, mismatch has no effect; the mdspan layout's Triangle rules. For example, the lower triangle of an upper triangle storage format is just the upper triangle again. For Hermitian matrices, again, mismatch has no effect. For example, suppose that the following is the lower triangle representation of a complex-valued Hermitian matrix (where $i$ is $\sqrt{-1}$).
$$ \begin{matrix} 1+1i & & \ 2+2i & 4+4i & \ 3+3i & 5+5i & 6+6i \ \end{matrix} $$
If the user asks the function to operate on the upper triangle of this matrix, that would imply the following.
$$ \begin{matrix} 1+1i & 2-2i & 3-3i \ & 4+4i & 5-5i \ & & 6+6i \ \end{matrix} $$
(Note that the imaginary parts now have negative sign. The matrix is Hermitian, so A[j,i]
equals conj(A[i,j])
.) That's just the "other triangle" of the matrix. These are Hermitian algorithms, so they will interpret the "other other triangle" in a way that restores the original matrix. Even though the user never stores the original matrix, it would look like this mathematically.
$$ \begin{matrix} 1+1i & 2-2i & 3-3i \ 2+2i & 4+4i & 5-5i \ 3+3i & 5+5i & 6+6i \ \end{matrix} $$
PR #428, which fixed these issues, has merged.
P1673: LWG review 2023/11/07 (Kona)
[linalg.algs.blas3.trsm]
triangular_matrix_matrix_left_solve
divide
binary function object in that case. (DONE in PR #426 .)A
is the triangular matrix" --B
and/orX
could also be mathematically triangular. We're trying to say thatt
andd
apply toA
. However, linalg.general already says that the t and d apply to the preceding argument. Instead of saying "where applicable," say "that apply toA
." Make this change throughout. (DONE in PR #426 . "taking into account theTriangle
" is a key search phrase.)std::divides<void>
as the default binary divide operator, instead of a lambda. Make this change throughout. (DONE in PR #426 .)triangular_matrix_matrix_right_solve
(See above.)
[linalg.algs.blas3.inplacetrsm]
std::divides<void>
as the default binary divide operator. (DONE in PR #426 .)triangular_matrix_matrix_left_solve
andtriangular_matrix_matrix_right_solve
. (DONE in PR #426 .)B
, notX
. (DONE in PR #426 .)[algorithms.parallel.defns]
(OK)
[algorithms.parallel.user]
(OK)
[linalg.layout.packed]
operator()
to take two parameters only, not a pack. (DONE in PR #426 .)is_always_unique()
so it checks both static extents (because one extent could be static and the other could be dynamic). Ditto foris_always_strided()
. (is_unique()
andis_strided()
are fine, because there's a precondition that the matrix be square.) (DONE in PR #428 .)r1
to0
, and remove the double right parenthesis afterstatic_extent(1)
. (DONE; fixed in PR #426 .)regular
if you put in weird template arguments.) (OK, this is fine.)[linalg.layout.packed.cons]
Para 6.1 (first Precondition) is a bit pessimistic, because we only store half the size. On the other hand, worrying about that would make generic algorithms harder to implement. (See below for the action item.)
[linalg.layout.packed.obs]
required_span_size()
expression is in code font.extents_.extent(0) + 1
could overflow, and more importantly,extents_.extent(0) * (extents_.extent(0) + 1)
is bigger than the size of the multidimensional index space. Write this in math font, because the result is definitely in range. Look for other places in the proposal that have this issue. Alternately, replace Mandate above (in Para 4) with a requirement (Mandate + Precondition) thatextent(0) * (extent(0) + 1)
be representable. (Do the latter.) (Create a function to check?) (9, 6.1, 4.6 -- use the expression without the divides by 2, in math font, for Mandates and Preconditions throughout.) (DONE in PR #428 .)operator()
operator()
: Para 12.2:index_Type
->index_type
. (DONE; fixed in PR #426 .)operator()
should take exactly two parameters. Adjust 12-14 accordingly.i
andj
toindex_type
before we calculate. (DONE in PR #426 .)operator()
overflow issue? Avoiding overflow would happen if we require thatextent(0) * (extent(0) + 1)
won't overflow. (Did we resolve this?) (DONE in PR #428 .)stride
[linalg.transp]
transpose-extents-t
andtranspose-extents
transpose-extents-t
. Write it as code, not as English. Fix the Mandates, as the code wouldn't be able to express this. Also, the wording doesn't currently say that theindex_type
is the same; fix that. Instead, havetranspose-extents
returnauto
. (DONE in PR #426 .)transpose-extent-t
andtranspose-extents
theretranspose-extents-t
in terms oftranspose-extents
, rather than the other way aroundResume tomorrow with
layout_transpose
.Tasks 2023/11/08
transpose-extents
definition (so it doesn't needauto
and so the return type is obvious). (DONE in PR #428 .)triangle_t
type alias tolayout_blas_packed
, so users can write generic code. (DONE in PR #428 .)