291 KiB
[linalg]
29 Numerics library [numerics]
29.9 Basic linear algebra algorithms [linalg]
29.9.1 Overview [linalg.overview]
Subclause [linalg] defines basic linear algebra algorithms.
The algorithms that access the elements of arrays view those elements through mdspan ([views.multidim]).
29.9.2 Header synopsis [linalg.syn]
namespace std::linalg {// [linalg.tags.order], storage order tagsstruct column_major_t; inline constexpr column_major_t column_major; struct row_major_t; inline constexpr row_major_t row_major; // [linalg.tags.triangle], triangle tagsstruct upper_triangle_t; inline constexpr upper_triangle_t upper_triangle; struct lower_triangle_t; inline constexpr lower_triangle_t lower_triangle; // [linalg.tags.diagonal], diagonal tagsstruct implicit_unit_diagonal_t; inline constexpr implicit_unit_diagonal_t implicit_unit_diagonal; struct explicit_diagonal_t; inline constexpr explicit_diagonal_t explicit_diagonal; // [linalg.layout.packed], class template layout_blas_packedtemplate<class Triangle, class StorageOrder>class layout_blas_packed; // [linalg.helpers], exposition-only helpers// [linalg.helpers.concepts], linear algebra argument conceptstemplateconstexpr bool is-mdspan = see below; // exposition onlytemplateconcept in-vector = see below; // exposition onlytemplateconcept out-vector = see below; // exposition onlytemplateconcept inout-vector = see below; // exposition onlytemplateconcept in-matrix = see below; // exposition onlytemplateconcept out-matrix = see below; // exposition onlytemplateconcept inout-matrix = see below; // exposition onlytemplateconcept possibly-packed-inout-matrix = see below; // exposition onlytemplateconcept in-object = see below; // exposition onlytemplateconcept out-object = see below; // exposition onlytemplateconcept inout-object = see below; // exposition only// [linalg.scaled], scaled in-place transformation// [linalg.scaled.scaledaccessor], class template scaled_accessortemplate<class ScalingFactor, class NestedAccessor>class scaled_accessor; // [linalg.scaled.scaled], function template scaledtemplate<class ScalingFactor, class ElementType, class Extents, class Layout, class Accessor>constexpr auto scaled(ScalingFactor alpha, mdspan<ElementType, Extents, Layout, Accessor> x); // [linalg.conj], conjugated in-place transformation// [linalg.conj.conjugatedaccessor], class template conjugated_accessortemplateclass conjugated_accessor; // [linalg.conj.conjugated], function template conjugatedtemplate<class ElementType, class Extents, class Layout, class Accessor>constexpr auto conjugated(mdspan<ElementType, Extents, Layout, Accessor> a); // [linalg.transp], transpose in-place transformation// [linalg.transp.layout.transpose], class template layout_transposetemplateclass layout_transpose; // [linalg.transp.transposed], function template transposedtemplate<class ElementType, class Extents, class Layout, class Accessor>constexpr auto transposed(mdspan<ElementType, Extents, Layout, Accessor> a); // [linalg.conjtransposed], conjugated transpose in-place transformationtemplate<class ElementType, class Extents, class Layout, class Accessor>constexpr auto conjugate_transposed(mdspan<ElementType, Extents, Layout, Accessor> a); // [linalg.algs.blas1], BLAS 1 algorithms// [linalg.algs.blas1.givens], Givens rotations// [linalg.algs.blas1.givens.lartg], compute Givens rotationtemplatestruct setup_givens_rotation_result { Real c; Real s; Real r; }; templatestruct setup_givens_rotation_result<complex> { Real c; complex s; complex r; }; template setup_givens_rotation_result setup_givens_rotation(Real a, Real b) noexcept; template setup_givens_rotation_result<complex> setup_givens_rotation(complex a, complex b) noexcept; // [linalg.algs.blas1.givens.rot], apply computed Givens rotationtemplate<inout-vector InOutVec1, inout-vector InOutVec2, class Real>void apply_givens_rotation(InOutVec1 x, InOutVec2 y, Real c, Real s); template<class ExecutionPolicy, inout-vector InOutVec1, inout-vector InOutVec2, class Real>void apply_givens_rotation(ExecutionPolicy&& exec, InOutVec1 x, InOutVec2 y, Real c, Real s); template<inout-vector InOutVec1, inout-vector InOutVec2, class Real>void apply_givens_rotation(InOutVec1 x, InOutVec2 y, Real c, complex s); template<class ExecutionPolicy, inout-vector InOutVec1, inout-vector InOutVec2, class Real>void apply_givens_rotation(ExecutionPolicy&& exec, InOutVec1 x, InOutVec2 y, Real c, complex s); // [linalg.algs.blas1.swap], swap elementstemplate<inout-object InOutObj1, inout-object InOutObj2>void swap_elements(InOutObj1 x, InOutObj2 y); template<class ExecutionPolicy, inout-object InOutObj1, inout-object InOutObj2>void swap_elements(ExecutionPolicy&& exec, InOutObj1 x, InOutObj2 y); // [linalg.algs.blas1.scal], multiply elements by scalartemplate<class Scalar, inout-object InOutObj>void scale(Scalar alpha, InOutObj x); template<class ExecutionPolicy, class Scalar, inout-object InOutObj>void scale(ExecutionPolicy&& exec, Scalar alpha, InOutObj x); // [linalg.algs.blas1.copy], copy elementstemplate<in-object InObj, out-object OutObj>void copy(InObj x, OutObj y); template<class ExecutionPolicy, in-object InObj, out-object OutObj>void copy(ExecutionPolicy&& exec, InObj x, OutObj y); // [linalg.algs.blas1.add], add elementwisetemplate<in-object InObj1, in-object InObj2, out-object OutObj>void add(InObj1 x, InObj2 y, OutObj z); template<class ExecutionPolicy, in-object InObj1, in-object InObj2, out-object OutObj>void add(ExecutionPolicy&& exec, InObj1 x, InObj2 y, OutObj z); // [linalg.algs.blas1.dot], dot product of two vectorstemplate<in-vector InVec1, in-vector InVec2, class Scalar> Scalar dot(InVec1 v1, InVec2 v2, Scalar init); template<class ExecutionPolicy, in-vector InVec1, in-vector InVec2, class Scalar> Scalar dot(ExecutionPolicy&& exec, InVec1 v1, InVec2 v2, Scalar init); template<in-vector InVec1, in-vector InVec2>auto dot(InVec1 v1, InVec2 v2); template<class ExecutionPolicy, in-vector InVec1, in-vector InVec2>auto dot(ExecutionPolicy&& exec, InVec1 v1, InVec2 v2); template<in-vector InVec1, in-vector InVec2, class Scalar> Scalar dotc(InVec1 v1, InVec2 v2, Scalar init); template<class ExecutionPolicy, in-vector InVec1, in-vector InVec2, class Scalar> Scalar dotc(ExecutionPolicy&& exec, InVec1 v1, InVec2 v2, Scalar init); template<in-vector InVec1, in-vector InVec2>auto dotc(InVec1 v1, InVec2 v2); template<class ExecutionPolicy, in-vector InVec1, in-vector InVec2>auto dotc(ExecutionPolicy&& exec, InVec1 v1, InVec2 v2); // [linalg.algs.blas1.ssq], scaled sum of squares of a vector's elementstemplatestruct sum_of_squares_result { Scalar scaling_factor; Scalar scaled_sum_of_squares; }; template<in-vector InVec, class Scalar> sum_of_squares_result vector_sum_of_squares(InVec v, sum_of_squares_result init); template<class ExecutionPolicy, in-vector InVec, class Scalar> sum_of_squares_result vector_sum_of_squares(ExecutionPolicy&& exec, InVec v, sum_of_squares_result init); // [linalg.algs.blas1.nrm2], Euclidean norm of a vectortemplate<in-vector InVec, class Scalar> Scalar vector_two_norm(InVec v, Scalar init); template<class ExecutionPolicy, in-vector InVec, class Scalar> Scalar vector_two_norm(ExecutionPolicy&& exec, InVec v, Scalar init); template<in-vector InVec>auto vector_two_norm(InVec v); template<class ExecutionPolicy, in-vector InVec>auto vector_two_norm(ExecutionPolicy&& exec, InVec v); // [linalg.algs.blas1.asum], sum of absolute values of vector elementstemplate<in-vector InVec, class Scalar> Scalar vector_abs_sum(InVec v, Scalar init); template<class ExecutionPolicy, in-vector InVec, class Scalar> Scalar vector_abs_sum(ExecutionPolicy&& exec, InVec v, Scalar init); template<in-vector InVec>auto vector_abs_sum(InVec v); template<class ExecutionPolicy, in-vector InVec>auto vector_abs_sum(ExecutionPolicy&& exec, InVec v); // [linalg.algs.blas1.iamax], index of maximum absolute value of vector elementstemplate<in-vector InVec>typename InVec::extents_type vector_idx_abs_max(InVec v); template<class ExecutionPolicy, in-vector InVec>typename InVec::extents_type vector_idx_abs_max(ExecutionPolicy&& exec, InVec v); // [linalg.algs.blas1.matfrobnorm], Frobenius norm of a matrixtemplate<in-matrix InMat, class Scalar> Scalar matrix_frob_norm(InMat A, Scalar init); template<class ExecutionPolicy, in-matrix InMat, class Scalar> Scalar matrix_frob_norm(ExecutionPolicy&& exec, InMat A, Scalar init); template<in-matrix InMat>auto matrix_frob_norm(InMat A); template<class ExecutionPolicy, in-matrix InMat>auto matrix_frob_norm(ExecutionPolicy&& exec, InMat A); // [linalg.algs.blas1.matonenorm], one norm of a matrixtemplate<in-matrix InMat, class Scalar> Scalar matrix_one_norm(InMat A, Scalar init); template<class ExecutionPolicy, in-matrix InMat, class Scalar> Scalar matrix_one_norm(ExecutionPolicy&& exec, InMat A, Scalar init); template<in-matrix InMat>auto matrix_one_norm(InMat A); template<class ExecutionPolicy, in-matrix InMat>auto matrix_one_norm(ExecutionPolicy&& exec, InMat A); // [linalg.algs.blas1.matinfnorm], infinity norm of a matrixtemplate<in-matrix InMat, class Scalar> Scalar matrix_inf_norm(InMat A, Scalar init); template<class ExecutionPolicy, in-matrix InMat, class Scalar> Scalar matrix_inf_norm(ExecutionPolicy&& exec, InMat A, Scalar init); template<in-matrix InMat>auto matrix_inf_norm(InMat A); template<class ExecutionPolicy, in-matrix InMat>auto matrix_inf_norm(ExecutionPolicy&& exec, InMat A); // [linalg.algs.blas2], BLAS 2 algorithms// [linalg.algs.blas2.gemv], general matrix-vector producttemplate<in-matrix InMat, in-vector InVec, out-vector OutVec>void matrix_vector_product(InMat A, InVec x, OutVec y); template<class ExecutionPolicy, in-matrix InMat, in-vector InVec, out-vector OutVec>void matrix_vector_product(ExecutionPolicy&& exec, InMat A, InVec x, OutVec y); template<in-matrix InMat, in-vector InVec1, in-vector InVec2, out-vector OutVec>void matrix_vector_product(InMat A, InVec1 x, InVec2 y, OutVec z); template<class ExecutionPolicy, in-matrix InMat, in-vector InVec1, in-vector InVec2, out-vector OutVec>void matrix_vector_product(ExecutionPolicy&& exec, InMat A, InVec1 x, InVec2 y, OutVec z); // [linalg.algs.blas2.symv], symmetric matrix-vector producttemplate<in-matrix InMat, class Triangle, in-vector InVec, out-vector OutVec>void symmetric_matrix_vector_product(InMat A, Triangle t, InVec x, OutVec y); template<class ExecutionPolicy, in-matrix InMat, class Triangle, in-vector InVec, out-vector OutVec>void symmetric_matrix_vector_product(ExecutionPolicy&& exec, InMat A, Triangle t, InVec x, OutVec y); template<in-matrix InMat, class Triangle, in-vector InVec1, in-vector InVec2, out-vector OutVec>void symmetric_matrix_vector_product(InMat A, Triangle t, InVec1 x, InVec2 y, OutVec z); template<class ExecutionPolicy, in-matrix InMat, class Triangle, in-vector InVec1, in-vector InVec2, out-vector OutVec>void symmetric_matrix_vector_product(ExecutionPolicy&& exec, InMat A, Triangle t, InVec1 x, InVec2 y, OutVec z); // [linalg.algs.blas2.hemv], Hermitian matrix-vector producttemplate<in-matrix InMat, class Triangle, in-vector InVec, out-vector OutVec>void hermitian_matrix_vector_product(InMat A, Triangle t, InVec x, OutVec y); template<class ExecutionPolicy, in-matrix InMat, class Triangle, in-vector InVec, out-vector OutVec>void hermitian_matrix_vector_product(ExecutionPolicy&& exec, InMat A, Triangle t, InVec x, OutVec y); template<in-matrix InMat, class Triangle, in-vector InVec1, in-vector InVec2, out-vector OutVec>void hermitian_matrix_vector_product(InMat A, Triangle t, InVec1 x, InVec2 y, OutVec z); template<class ExecutionPolicy, in-matrix InMat, class Triangle, in-vector InVec1, in-vector InVec2, out-vector OutVec>void hermitian_matrix_vector_product(ExecutionPolicy&& exec, InMat A, Triangle t, InVec1 x, InVec2 y, OutVec z); // [linalg.algs.blas2.trmv], triangular matrix-vector product// Overwriting triangular matrix-vector producttemplate<in-matrix InMat, class Triangle, class DiagonalStorage, in-vector InVec, out-vector OutVec>void triangular_matrix_vector_product(InMat A, Triangle t, DiagonalStorage d, InVec x, OutVec y); template<class ExecutionPolicy, in-matrix InMat, class Triangle, class DiagonalStorage, in-vector InVec, out-vector OutVec>void triangular_matrix_vector_product(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InVec x, OutVec y); // In-place triangular matrix-vector producttemplate<in-matrix InMat, class Triangle, class DiagonalStorage, inout-vector InOutVec>void triangular_matrix_vector_product(InMat A, Triangle t, DiagonalStorage d, InOutVec y); template<class ExecutionPolicy, in-matrix InMat, class Triangle, class DiagonalStorage, inout-vector InOutVec>void triangular_matrix_vector_product(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InOutVec y); // Updating triangular matrix-vector producttemplate<in-matrix InMat, class Triangle, class DiagonalStorage, in-vector InVec1, in-vector InVec2, out-vector OutVec>void triangular_matrix_vector_product(InMat A, Triangle t, DiagonalStorage d, InVec1 x, InVec2 y, OutVec z); template<class ExecutionPolicy, in-matrix InMat, class Triangle, class DiagonalStorage, in-vector InVec1, in-vector InVec2, out-vector OutVec>void triangular_matrix_vector_product(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InVec1 x, InVec2 y, OutVec z); // [linalg.algs.blas2.trsv], solve a triangular linear system// Solve a triangular linear system, not in placetemplate<in-matrix InMat, class Triangle, class DiagonalStorage, in-vector InVec, out-vector OutVec, class BinaryDivideOp>void triangular_matrix_vector_solve(InMat A, Triangle t, DiagonalStorage d, InVec b, OutVec x, BinaryDivideOp divide); template<class ExecutionPolicy, in-matrix InMat, class Triangle, class DiagonalStorage, in-vector InVec, out-vector OutVec, class BinaryDivideOp>void triangular_matrix_vector_solve(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InVec b, OutVec x, BinaryDivideOp divide); template<in-matrix InMat, class Triangle, class DiagonalStorage, in-vector InVec, out-vector OutVec>void triangular_matrix_vector_solve(InMat A, Triangle t, DiagonalStorage d, InVec b, OutVec x); template<class ExecutionPolicy, in-matrix InMat, class Triangle, class DiagonalStorage, in-vector InVec, out-vector OutVec>void triangular_matrix_vector_solve(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InVec b, OutVec x); // Solve a triangular linear system, in placetemplate<in-matrix InMat, class Triangle, class DiagonalStorage, inout-vector InOutVec, class BinaryDivideOp>void triangular_matrix_vector_solve(InMat A, Triangle t, DiagonalStorage d, InOutVec b, BinaryDivideOp divide); template<class ExecutionPolicy, in-matrix InMat, class Triangle, class DiagonalStorage, inout-vector InOutVec, class BinaryDivideOp>void triangular_matrix_vector_solve(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InOutVec b, BinaryDivideOp divide); template<in-matrix InMat, class Triangle, class DiagonalStorage, inout-vector InOutVec>void triangular_matrix_vector_solve(InMat A, Triangle t, DiagonalStorage d, InOutVec b); template<class ExecutionPolicy, in-matrix InMat, class Triangle, class DiagonalStorage, inout-vector InOutVec>void triangular_matrix_vector_solve(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InOutVec b); // [linalg.algs.blas2.rank1], nonsymmetric rank-1 matrix updatetemplate<in-vector InVec1, in-vector InVec2, inout-matrix InOutMat>void matrix_rank_1_update(InVec1 x, InVec2 y, InOutMat A); template<class ExecutionPolicy, in-vector InVec1, in-vector InVec2, inout-matrix InOutMat>void matrix_rank_1_update(ExecutionPolicy&& exec, InVec1 x, InVec2 y, InOutMat A); template<in-vector InVec1, in-vector InVec2, inout-matrix InOutMat>void matrix_rank_1_update_c(InVec1 x, InVec2 y, InOutMat A); template<class ExecutionPolicy, in-vector InVec1, in-vector InVec2, inout-matrix InOutMat>void matrix_rank_1_update_c(ExecutionPolicy&& exec, InVec1 x, InVec2 y, InOutMat A); // [linalg.algs.blas2.symherrank1], symmetric or Hermitian rank-1 matrix updatetemplate<class Scalar, in-vector InVec, possibly-packed-inout-matrix InOutMat, class Triangle>void symmetric_matrix_rank_1_update(Scalar alpha, InVec x, InOutMat A, Triangle t); template<class ExecutionPolicy, class Scalar, in-vector InVec, possibly-packed-inout-matrix InOutMat, class Triangle>void symmetric_matrix_rank_1_update(ExecutionPolicy&& exec, Scalar alpha, InVec x, InOutMat A, Triangle t); template<in-vector InVec, possibly-packed-inout-matrix InOutMat, class Triangle>void symmetric_matrix_rank_1_update(InVec x, InOutMat A, Triangle t); template<class ExecutionPolicy, in-vector InVec, possibly-packed-inout-matrix InOutMat, class Triangle>void symmetric_matrix_rank_1_update(ExecutionPolicy&& exec, InVec x, InOutMat A, Triangle t); template<class Scalar, in-vector InVec, possibly-packed-inout-matrix InOutMat, class Triangle>void hermitian_matrix_rank_1_update(Scalar alpha, InVec x, InOutMat A, Triangle t); template<class ExecutionPolicy, class Scalar, in-vector InVec, possibly-packed-inout-matrix InOutMat, class Triangle>void hermitian_matrix_rank_1_update(ExecutionPolicy&& exec, Scalar alpha, InVec x, InOutMat A, Triangle t); template<in-vector InVec, possibly-packed-inout-matrix InOutMat, class Triangle>void hermitian_matrix_rank_1_update(InVec x, InOutMat A, Triangle t); template<class ExecutionPolicy, in-vector InVec, possibly-packed-inout-matrix InOutMat, class Triangle>void hermitian_matrix_rank_1_update(ExecutionPolicy&& exec, InVec x, InOutMat A, Triangle t); // [linalg.algs.blas2.rank2], symmetric and Hermitian rank-2 matrix updates// symmetric rank-2 matrix updatetemplate<in-vector InVec1, in-vector InVec2, possibly-packed-inout-matrix InOutMat, class Triangle>void symmetric_matrix_rank_2_update(InVec1 x, InVec2 y, InOutMat A, Triangle t); template<class ExecutionPolicy, in-vector InVec1, in-vector InVec2, possibly-packed-inout-matrix InOutMat, class Triangle>void symmetric_matrix_rank_2_update(ExecutionPolicy&& exec, InVec1 x, InVec2 y, InOutMat A, Triangle t); // Hermitian rank-2 matrix updatetemplate<in-vector InVec1, in-vector InVec2, possibly-packed-inout-matrix InOutMat, class Triangle>void hermitian_matrix_rank_2_update(InVec1 x, InVec2 y, InOutMat A, Triangle t); template<class ExecutionPolicy, in-vector InVec1, in-vector InVec2, possibly-packed-inout-matrix InOutMat, class Triangle>void hermitian_matrix_rank_2_update(ExecutionPolicy&& exec, InVec1 x, InVec2 y, InOutMat A, Triangle t); // [linalg.algs.blas3], BLAS 3 algorithms// [linalg.algs.blas3.gemm], general matrix-matrix producttemplate<in-matrix InMat1, in-matrix InMat2, out-matrix OutMat>void matrix_product(InMat1 A, InMat2 B, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, in-matrix InMat2, out-matrix OutMat>void matrix_product(ExecutionPolicy&& exec, InMat1 A, InMat2 B, OutMat C); template<in-matrix InMat1, in-matrix InMat2, in-matrix InMat3, out-matrix OutMat>void matrix_product(InMat1 A, InMat2 B, InMat3 E, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, in-matrix InMat2, in-matrix InMat3, out-matrix OutMat>void matrix_product(ExecutionPolicy&& exec, InMat1 A, InMat2 B, InMat3 E, OutMat C); // [linalg.algs.blas3.xxmm], symmetric, Hermitian, and triangular matrix-matrix producttemplate<in-matrix InMat1, class Triangle, in-matrix InMat2, out-matrix OutMat>void symmetric_matrix_product(InMat1 A, Triangle t, InMat2 B, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, class Triangle, in-matrix InMat2, out-matrix OutMat>void symmetric_matrix_product(ExecutionPolicy&& exec, InMat1 A, Triangle t, InMat2 B, OutMat C); template<in-matrix InMat1, class Triangle, in-matrix InMat2, out-matrix OutMat>void hermitian_matrix_product(InMat1 A, Triangle t, InMat2 B, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, class Triangle, in-matrix InMat2, out-matrix OutMat>void hermitian_matrix_product(ExecutionPolicy&& exec, InMat1 A, Triangle t, InMat2 B, OutMat C); template<in-matrix InMat1, class Triangle, class DiagonalStorage, in-matrix InMat2, out-matrix OutMat>void triangular_matrix_product(InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, class Triangle, class DiagonalStorage, in-matrix InMat2, out-matrix OutMat>void triangular_matrix_product(ExecutionPolicy&& exec, InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat C); template<in-matrix InMat1, in-matrix InMat2, class Triangle, out-matrix OutMat>void symmetric_matrix_product(InMat1 A, InMat2 B, Triangle t, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, in-matrix InMat2, class Triangle, out-matrix OutMat>void symmetric_matrix_product(ExecutionPolicy&& exec, InMat1 A, InMat2 B, Triangle t, OutMat C); template<in-matrix InMat1, in-matrix InMat2, class Triangle, out-matrix OutMat>void hermitian_matrix_product(InMat1 A, InMat2 B, Triangle t, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, in-matrix InMat2, class Triangle, out-matrix OutMat>void hermitian_matrix_product(ExecutionPolicy&& exec, InMat1 A, InMat2 B, Triangle t, OutMat C); template<in-matrix InMat1, in-matrix InMat2, class Triangle, class DiagonalStorage, out-matrix OutMat>void triangular_matrix_product(InMat1 A, InMat2 B, Triangle t, DiagonalStorage d, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, in-matrix InMat2, class Triangle, class DiagonalStorage, out-matrix OutMat>void triangular_matrix_product(ExecutionPolicy&& exec, InMat1 A, InMat2 B, Triangle t, DiagonalStorage d, OutMat C); template<in-matrix InMat1, class Triangle, in-matrix InMat2, in-matrix InMat3, out-matrix OutMat>void symmetric_matrix_product(InMat1 A, Triangle t, InMat2 B, InMat3 E, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, class Triangle, in-matrix InMat2, in-matrix InMat3, out-matrix OutMat>void symmetric_matrix_product(ExecutionPolicy&& exec, InMat1 A, Triangle t, InMat2 B, InMat3 E, OutMat C); template<in-matrix InMat1, class Triangle, in-matrix InMat2, in-matrix InMat3, out-matrix OutMat>void hermitian_matrix_product(InMat1 A, Triangle t, InMat2 B, InMat3 E, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, class Triangle, in-matrix InMat2, in-matrix InMat3, out-matrix OutMat>void hermitian_matrix_product(ExecutionPolicy&& exec, InMat1 A, Triangle t, InMat2 B, InMat3 E, OutMat C); template<in-matrix InMat1, class Triangle, class DiagonalStorage, in-matrix InMat2, in-matrix InMat3, out-matrix OutMat>void triangular_matrix_product(InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, InMat3 E, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, class Triangle, class DiagonalStorage, in-matrix InMat2, in-matrix InMat3, out-matrix OutMat>void triangular_matrix_product(ExecutionPolicy&& exec, InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, InMat3 E, OutMat C); template<in-matrix InMat1, in-matrix InMat2, class Triangle, in-matrix InMat3, out-matrix OutMat>void symmetric_matrix_product(InMat1 A, InMat2 B, Triangle t, InMat3 E, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, in-matrix InMat2, class Triangle, in-matrix InMat3, out-matrix OutMat>void symmetric_matrix_product(ExecutionPolicy&& exec, InMat1 A, InMat2 B, Triangle t, InMat3 E, OutMat C); template<in-matrix InMat1, in-matrix InMat2, class Triangle, in-matrix InMat3, out-matrix OutMat>void hermitian_matrix_product(InMat1 A, InMat2 B, Triangle t, InMat3 E, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, in-matrix InMat2, class Triangle, in-matrix InMat3, out-matrix OutMat>void hermitian_matrix_product(ExecutionPolicy&& exec, InMat1 A, InMat2 B, Triangle t, InMat3 E, OutMat C); template<in-matrix InMat1, in-matrix InMat2, class Triangle, class DiagonalStorage, in-matrix InMat3, out-matrix OutMat>void triangular_matrix_product(InMat1 A, InMat2 B, Triangle t, DiagonalStorage d, InMat3 E, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, in-matrix InMat2, class Triangle, class DiagonalStorage, in-matrix InMat3, out-matrix OutMat>void triangular_matrix_product(ExecutionPolicy&& exec, InMat1 A, InMat2 B, Triangle t, DiagonalStorage d, InMat3 E, OutMat C); // [linalg.algs.blas3.trmm], in-place triangular matrix-matrix producttemplate<in-matrix InMat, class Triangle, class DiagonalStorage, inout-matrix InOutMat>void triangular_matrix_left_product(InMat A, Triangle t, DiagonalStorage d, InOutMat C); template<class ExecutionPolicy, in-matrix InMat, class Triangle, class DiagonalStorage, inout-matrix InOutMat>void triangular_matrix_left_product(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InOutMat C); template<in-matrix InMat, class Triangle, class DiagonalStorage, inout-matrix InOutMat>void triangular_matrix_right_product(InMat A, Triangle t, DiagonalStorage d, InOutMat C); template<class ExecutionPolicy, in-matrix InMat, class Triangle, class DiagonalStorage, inout-matrix InOutMat>void triangular_matrix_right_product(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InOutMat C); // [linalg.algs.blas3.rankk], rank-k update of a symmetric or Hermitian matrix// rank-k symmetric matrix updatetemplate<class Scalar, in-matrix InMat, possibly-packed-inout-matrix InOutMat, class Triangle>void symmetric_matrix_rank_k_update(Scalar alpha, InMat A, InOutMat C, Triangle t); template<class ExecutionPolicy, class Scalar, in-matrix InMat, possibly-packed-inout-matrix InOutMat, class Triangle>void symmetric_matrix_rank_k_update(ExecutionPolicy&& exec, Scalar alpha, InMat A, InOutMat C, Triangle t); template<in-matrix InMat, possibly-packed-inout-matrix InOutMat, class Triangle>void symmetric_matrix_rank_k_update(InMat A, InOutMat C, Triangle t); template<class ExecutionPolicy, in-matrix InMat, possibly-packed-inout-matrix InOutMat, class Triangle>void symmetric_matrix_rank_k_update(ExecutionPolicy&& exec, InMat A, InOutMat C, Triangle t); // rank-k Hermitian matrix updatetemplate<class Scalar, in-matrix InMat, possibly-packed-inout-matrix InOutMat, class Triangle>void hermitian_matrix_rank_k_update(Scalar alpha, InMat A, InOutMat C, Triangle t); template<class ExecutionPolicy, class Scalar, in-matrix InMat, possibly-packed-inout-matrix InOutMat, class Triangle>void hermitian_matrix_rank_k_update(ExecutionPolicy&& exec, Scalar alpha, InMat A, InOutMat C, Triangle t); template<in-matrix InMat, possibly-packed-inout-matrix InOutMat, class Triangle>void hermitian_matrix_rank_k_update(InMat A, InOutMat C, Triangle t); template<class ExecutionPolicy, in-matrix InMat, possibly-packed-inout-matrix InOutMat, class Triangle>void hermitian_matrix_rank_k_update(ExecutionPolicy&& exec, InMat A, InOutMat C, Triangle t); // [linalg.algs.blas3.rank2k], rank-2k update of a symmetric or Hermitian matrix// rank-2k symmetric matrix updatetemplate<in-matrix InMat1, in-matrix InMat2, possibly-packed-inout-matrix InOutMat, class Triangle>void symmetric_matrix_rank_2k_update(InMat1 A, InMat2 B, InOutMat C, Triangle t); template<class ExecutionPolicy, in-matrix InMat1, in-matrix InMat2, possibly-packed-inout-matrix InOutMat, class Triangle>void symmetric_matrix_rank_2k_update(ExecutionPolicy&& exec, InMat1 A, InMat2 B, InOutMat C, Triangle t); // rank-2k Hermitian matrix updatetemplate<in-matrix InMat1, in-matrix InMat2, possibly-packed-inout-matrix InOutMat, class Triangle>void hermitian_matrix_rank_2k_update(InMat1 A, InMat2 B, InOutMat C, Triangle t); template<class ExecutionPolicy, in-matrix InMat1, in-matrix InMat2, possibly-packed-inout-matrix InOutMat, class Triangle>void hermitian_matrix_rank_2k_update(ExecutionPolicy&& exec, InMat1 A, InMat2 B, InOutMat C, Triangle t); // [linalg.algs.blas3.trsm], solve multiple triangular linear systems// solve multiple triangular systems on the left, not-in-placetemplate<in-matrix InMat1, class Triangle, class DiagonalStorage, in-matrix InMat2, out-matrix OutMat, class BinaryDivideOp>void triangular_matrix_matrix_left_solve(InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat X, BinaryDivideOp divide); template<class ExecutionPolicy, in-matrix InMat1, class Triangle, class DiagonalStorage, in-matrix InMat2, out-matrix OutMat, class BinaryDivideOp>void triangular_matrix_matrix_left_solve(ExecutionPolicy&& exec, InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat X, BinaryDivideOp divide); template<in-matrix InMat1, class Triangle, class DiagonalStorage, in-matrix InMat2, out-matrix OutMat>void triangular_matrix_matrix_left_solve(InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat X); template<class ExecutionPolicy, in-matrix InMat1, class Triangle, class DiagonalStorage, in-matrix InMat2, out-matrix OutMat>void triangular_matrix_matrix_left_solve(ExecutionPolicy&& exec, InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat X); // solve multiple triangular systems on the right, not-in-placetemplate<in-matrix InMat1, class Triangle, class DiagonalStorage, in-matrix InMat2, out-matrix OutMat, class BinaryDivideOp>void triangular_matrix_matrix_right_solve(InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat X, BinaryDivideOp divide); template<class ExecutionPolicy, in-matrix InMat1, class Triangle, class DiagonalStorage, in-matrix InMat2, out-matrix OutMat, class BinaryDivideOp>void triangular_matrix_matrix_right_solve(ExecutionPolicy&& exec, InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat X, BinaryDivideOp divide); template<in-matrix InMat1, class Triangle, class DiagonalStorage, in-matrix InMat2, out-matrix OutMat>void triangular_matrix_matrix_right_solve(InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat X); template<class ExecutionPolicy, in-matrix InMat1, class Triangle, class DiagonalStorage, in-matrix InMat2, out-matrix OutMat>void triangular_matrix_matrix_right_solve(ExecutionPolicy&& exec, InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat X); // solve multiple triangular systems on the left, in-placetemplate<in-matrix InMat, class Triangle, class DiagonalStorage, inout-matrix InOutMat, class BinaryDivideOp>void triangular_matrix_matrix_left_solve(InMat A, Triangle t, DiagonalStorage d, InOutMat B, BinaryDivideOp divide); template<class ExecutionPolicy, in-matrix InMat, class Triangle, class DiagonalStorage, inout-matrix InOutMat, class BinaryDivideOp>void triangular_matrix_matrix_left_solve(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InOutMat B, BinaryDivideOp divide); template<in-matrix InMat, class Triangle, class DiagonalStorage, inout-matrix InOutMat>void triangular_matrix_matrix_left_solve(InMat A, Triangle t, DiagonalStorage d, InOutMat B); template<class ExecutionPolicy, in-matrix InMat, class Triangle, class DiagonalStorage, inout-matrix InOutMat>void triangular_matrix_matrix_left_solve(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InOutMat B); // solve multiple triangular systems on the right, in-placetemplate<in-matrix InMat, class Triangle, class DiagonalStorage, inout-matrix InOutMat, class BinaryDivideOp>void triangular_matrix_matrix_right_solve(InMat A, Triangle t, DiagonalStorage d, InOutMat B, BinaryDivideOp divide); template<class ExecutionPolicy, in-matrix InMat, class Triangle, class DiagonalStorage, inout-matrix InOutMat, class BinaryDivideOp>void triangular_matrix_matrix_right_solve(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InOutMat B, BinaryDivideOp divide); template<in-matrix InMat, class Triangle, class DiagonalStorage, inout-matrix InOutMat>void triangular_matrix_matrix_right_solve(InMat A, Triangle t, DiagonalStorage d, InOutMat B); template<class ExecutionPolicy, in-matrix InMat, class Triangle, class DiagonalStorage, inout-matrix InOutMat>void triangular_matrix_matrix_right_solve(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InOutMat B);}
29.9.3 General [linalg.general]
For the effects of all functions in [linalg], when the effects are described as âcomputes R=Eâ or âcompute R=Eâ (for some R and mathematical expression E), the following apply:
-
E has the conventional mathematical meaning as written.
-
The pattern xT should be read as âthe transpose of x.â
-
The pattern xH should be read as âthe conjugate transpose of x.â
-
When R is the same name as a function parameter whose type is a template parameter with Out in its name, the intent is that the result of the computation is written to the elements of the function parameter R.
Some of the functions and types in [linalg] distinguish between the ârowsâ and the âcolumnsâ of a matrix.
For a matrix A and a multidimensional index i, j in A.extents(),
row i of A is the set of elements A[i, k1] for all k1 such that i, k1 is in A.extents(); and
column j of A is the set of elements A[k0, j] for all k0 such that k0, j is in A.extents().
Some of the functions in [linalg] distinguish between the âupper triangle,â âlower triangle,â and âdiagonalâ of a matrix.
-
The diagonal is the set of all elements of A accessed by A[i,i] for 0 ⤠i < min(A.extent(0), A.extent(1)).
-
The upper triangle of a matrix A is the set of all elements of A accessed by A[i,j] with i ⤠j. It includes the diagonal.
-
The lower triangle of A is the set of all elements of A accessed by A[i,j] with i ⥠j. It includes the diagonal.
For any function F that takes a parameter named t,t applies to accesses done through the parameter preceding t in the parameter list of F.
Let m be such an access-modified function parameter.
F will only access the triangle of m specified by t.
For accesses m[i, j] outside the triangle specified by t,F will use the value
conj-if-needed(m[j, i]) if the name of F starts with hermitian,
m[j, i] if the name of F starts with symmetric, or
the additive identity if the name of F starts with triangular.
[Example 1:
Small vector product accessing only specified triangle.
It would not be a precondition violation for the non-accessed matrix element to be non-zero.
templatevoid triangular_matrix_vector_2x2_product( mdspan<const float, extents<int, 2, 2>> m, Triangle t, mdspan<const float, extents<int, 2>> x, mdspan<float, extents<int, 2>> y) {static_assert(is_same_v<Triangle, lower_triangle_t> || is_same_v<Triangle, upper_triangle_t>); if constexpr (is_same_v<Triangle, lower_triangle_t>) { y[0] = m[0,0] * x[0]; // + 0 * x[1] y[1] = m[1,0] * x[0] + m[1,1] * x[1]; } else { // upper_triangle_t y[0] = m[0,0] * x[0] + m[0,1] * x[1]; y[1] = /* 0 * x[0] + */ m[1,1] * x[1]; }} â end example]
For any function F that takes a parameter named d,d applies to accesses done through the previous-of-the-previous parameter of d in the parameter list of F.
Let m be such an access-modified function parameter.
If d specifies that an implicit unit diagonal is to be assumed, then
F will not access the diagonal of m; and
the algorithm will interpret m as if it has a unit diagonal, that is, a diagonal each of whose elements behaves as a two-sided multiplicative identity (even if m's value type does not have a two-sided multiplicative identity).
Otherwise, if d specifies that an explicit diagonal is to be assumed, then F will access the diagonal of m.
Within all the functions in [linalg], any calls to abs, conj, imag, and real are unqualified.
Two mdspan objects x and y alias each other, if they have the same extents e, and for every pack of integers i which is a multidimensional index in e,x[i...] and y[i...] refer to the same element.
[Note 1:
This means thatx and y view the same elements in the same order.
â end note]
Two mdspan objects x and y overlap each other, if for some pack of integers i that is a multidimensional index in x.extents(), there exists a pack of integers j that is a multidimensional index in y.extents(), such that x[i...] and y[j...] refer to the same element.
[Note 2:
Aliasing is a special case of overlapping.
If x and y do not overlap, then they also do not alias each other.
â end note]
29.9.4 Requirements [linalg.reqs]
29.9.4.1 Linear algebra value types [linalg.reqs.val]
Throughout [linalg], the following types arelinear algebra value types:
the value_type type alias of any input or output mdspan parameter(s) of any function in [linalg]; and
the Scalar template parameter (if any) of any function or class in [linalg].
Linear algebra value types shall model semiregular.
A value-initialized object of linear algebra value type shall act as the additive identity.
29.9.4.2 Algorithm and class requirements [linalg.reqs.alg]
[linalg.reqs.alg] lists common requirements for all algorithms and classes in [linalg].
All of the following statements presume that the algorithm's asymptotic complexity requirements, if any, are satisfied.
-
The function may make arbitrarily many objects of any linear algebra value type, value-initializing or direct-initializing them with any existing object of that type.
-
The triangular solve algorithms in[linalg.algs.blas2.trsv],[linalg.algs.blas3.trmm],[linalg.algs.blas3.trsm], and[linalg.algs.blas3.inplacetrsm] either have a BinaryDivideOp template parameter (see [linalg.algs.reqs]) and a binary function object parameter divide of that type, or they have effects equivalent to invoking such an algorithm. Triangular solve algorithms interpret divide(a, b) asa times the multiplicative inverse of b. Each triangular solve algorithm uses a sequence of evaluations of*, *=, divide, unary +, binary +, +=, unary -, binary -, -=, and = operators that would produce the result specified by the algorithm's Effects and Remarks when operating on elements of a field with noncommutative multiplication. It is a precondition of the algorithm that any addend, any subtrahend, any partial sum of addends in any order (treating any difference as a sum with the second term negated), any factor, any partial product of factors respecting their order, any numerator (first argument of divide), any denominator (second argument of divide), and any assignment is a well-formed expression.
-
Each function in[linalg.algs.blas1], [linalg.algs.blas2], and [linalg.algs.blas3] that is not a triangular solve algorithm will use a sequence of evaluations of*, *=, +, +=, and = operators that would produce the result specified by the algorithm's Effects and Remarks when operating on elements of a semiring with noncommutative multiplication. It is a precondition of the algorithm that any addend, any partial sum of addends in any order, any factor, any partial product of factors respecting their order, and any assignment is a well-formed expression.
-
If the function has an output mdspan, then all addends, subtrahends (for the triangular solve algorithms), or results of the divide parameter on intermediate terms (if the function takes a divide parameter) are assignable and convertible to the output mdspan's value_type.
-
The function may reorder addends and partial sums arbitrarily. [Note 1: Factors in each product are not reordered; multiplication is not necessarily commutative. â end note]
[Note 2:
The above requirements do not prohibit implementation approaches and optimization techniques which are not user-observable.
In particular, if for all input and output arguments the value_type is a floating-point type, implementers are free to leverage approximations, use arithmetic operations not explicitly listed above, and compute floating point sums in any way that improves their accuracy.
â end note]
[Note 3:
For all functions in [linalg], suppose that all input and output mdspan have as value_type a floating-point type, and any Scalar template argument has a floating-point type.
Then, functions can do all of the following:
compute floating-point sums in any way that improves their accuracy for arbitrary input;
perform additional arithmetic operations (other than those specified by the function's wording and [linalg.reqs.alg]) in order to improve performance or accuracy; and
use approximations (that might not be exact even if computing with real numbers), instead of computations that would be exact if it were possible to compute without rounding error;
as long as
the function satisfies the complexity requirements; and
the function is logarithmically stable, as defined in Demmel 2007[bib]. Strassen's algorithm for matrix-matrix multiply is an example of a logarithmically stable algorithm.
â end note]
29.9.5 Tag classes [linalg.tags]
29.9.5.1 Storage order tags [linalg.tags.order]
The storage order tags describe the order of elements in an mdspan withlayout_blas_packed ([linalg.layout.packed]) layout.
`struct column_major_t { explicit column_major_t() = default; }; inline constexpr column_major_t column_major{};
struct row_major_t { explicit row_major_t() = default; }; inline constexpr row_major_t row_major{}; `
column_major_t indicates a column-major order, and row_major_t indicates a row-major order.
29.9.5.2 Triangle tags [linalg.tags.triangle]
`struct upper_triangle_t { explicit upper_triangle_t() = default; }; inline constexpr upper_triangle_t upper_triangle{};
struct lower_triangle_t { explicit lower_triangle_t() = default; }; inline constexpr lower_triangle_t lower_triangle{}; `
These tag classes specify whether algorithms and other users of a matrix (represented as an mdspan) access the upper triangle (upper_triangle_t) or lower triangle (lower_triangle_t) of the matrix (see also [linalg.general]).
This is also subject to the restrictions of implicit_unit_diagonal_t if that tag is also used as a function argument; see below.
29.9.5.3 Diagonal tags [linalg.tags.diagonal]
`struct implicit_unit_diagonal_t { explicit implicit_unit_diagonal_t() = default; }; inline constexpr implicit_unit_diagonal_t implicit_unit_diagonal{};
struct explicit_diagonal_t { explicit explicit_diagonal_t() = default; }; inline constexpr explicit_diagonal_t explicit_diagonal{}; `
These tag classes specify whether algorithms access the matrix's diagonal entries, and if not, then how algorithms interpret the matrix's implicitly represented diagonal values.
The implicit_unit_diagonal_t tag indicates that an implicit unit diagonal is to be assumed ([linalg.general]).
The explicit_diagonal_t tag indicates that an explicit diagonal is used ([linalg.general]).
29.9.6 Layouts for packed matrix types [linalg.layout.packed]
29.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; templatestruct 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], constructorsconstexpr mapping() noexcept = default; constexpr mapping(const mapping&) noexcept = default; constexpr mapping(const extents_type&) noexcept; templateconstexpr explicit(!is_convertible_v<OtherExtents, extents_type>) mapping(const mapping& other) noexcept; constexpr mapping& operator=(const mapping&) noexcept = default; // [linalg.layout.packed.obs], observersconstexpr 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; templatefriend constexpr bool operator==(const mapping&, const mapping&) noexcept; private: extents_type extents_{}; // exposition only}; };}
Mandates:
Triangle is either upper_triangle_t or lower_triangle_t,
StorageOrder is either column_major_t or row_major_t,
Extents is a specialization of std::extents,
Extents::rank() equals 2,
one ofextents_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
if Extents::rank_dynamic() == 0 is true, let Ns be equal to Extents::static_extent(0); then,NsÃ(Ns+1) is representable as a value of type index_type.
layout_blas_packed<T, SO>::mapping is a trivially copyable type that models regular for each T, SO, and E.
29.9.6.2 Constructors [linalg.layout.packed.cons]
constexpr mapping(const extents_type& e) noexcept;
Preconditions:
-
Let N be equal to e.extent(0). Then, NÃ(N+1) is representable as a value of type index_type ([basic.fundamental]).
-
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Ã(N+1) is representable as a value of type index_type ([basic.fundamental]).
Effects: Direct-non-list-initializes extents_ with other.extents().
29.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:
is_convertible_v<Index0, index_type> is true,
is_convertible_v<Index1, index_type> is true,
is_nothrow_constructible_v<index_type, Index0> is true, and
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
(*this)(j, i) if i > j is true; otherwise
i + j * (j + 1)/2 ifis_same_v<StorageOrder, column_major_t> && is_same_v<Triangle, upper_triangle_t> is true oris_same_v<StorageOrder, row_major_t> && is_same_v<Triangle, lower_triangle_t> is true; otherwise
j + N * i - i * (i + 1)/2.
constexpr index_type stride(rank_type r) const noexcept;
Preconditions:
is_strided() is true, and
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();
29.9.7 Exposition-only helpers [linalg.helpers]
29.9.7.1 abs-if-needed [linalg.helpers.abs]
The name abs-if-needed denotes an exposition-only function object.
The expression abs-if-needed(E) for a subexpression E whose type is T is expression-equivalent to:
E if T is an unsigned integer;
otherwise, std::abs(E) if T is an arithmetic type,
otherwise, abs(E), if that expression is valid, with overload resolution performed in a context that includes the declarationtemplate U abs(U) = delete; If the function selected by overload resolution does not return the absolute value of its input, the program is ill-formed, no diagnostic required.
29.9.7.2 conj-if-needed [linalg.helpers.conj]
The name conj-if-needed denotes an exposition-only function object.
The expression conj-if-needed(E) for a subexpression E whose type is T is expression-equivalent to:
conj(E), if T is not an arithmetic type and the expression conj(E) is valid, with overload resolution performed in a context that includes the declarationtemplate U conj(const U&) = delete; If the function selected by overload resolution does not return the complex conjugate of its input, the program is ill-formed, no diagnostic required;
otherwise, E.
29.9.7.3 real-if-needed [linalg.helpers.real]
The name real-if-needed denotes an exposition-only function object.
The expression real-if-needed(E) for a subexpression E whose type is T is expression-equivalent to:
real(E), if T is not an arithmetic type and the expression real(E) is valid, with overload resolution performed in a context that includes the declarationtemplate U real(const U&) = delete; If the function selected by overload resolution does not return the real part of its input, the program is ill-formed, no diagnostic required;
otherwise, E.
29.9.7.4 imag-if-needed [linalg.helpers.imag]
The name imag-if-needed denotes an exposition-only function object.
The expression imag-if-needed(E) for a subexpression E whose type is T is expression-equivalent to:
imag(E), if T is not an arithmetic type and the expression imag(E) is valid, with overload resolution performed in a context that includes the declarationtemplate U imag(const U&) = delete; If the function selected by overload resolution does not return the imaginary part of its input, the program is ill-formed, no diagnostic required;
otherwise, ((void)E, T{}).
29.9.7.5 Argument concepts [linalg.helpers.concepts]
The exposition-only concepts defined in this section constrain the algorithms in [linalg].
templateconstexpr bool is-mdspan = false;
template<class ElementType, class Extents, class Layout, class Accessor>constexpr bool is-mdspan<mdspan<ElementType, Extents, Layout, Accessor>> = true;
templateconcept in-vector =is-mdspan && T::rank() == 1;
templateconcept out-vector =is-mdspan && T::rank() == 1 && is_assignable_v<typename T::reference, typename T::element_type> && T::is_always_unique();
templateconcept inout-vector =is-mdspan && T::rank() == 1 && is_assignable_v<typename T::reference, typename T::element_type> && T::is_always_unique();
templateconcept in-matrix =is-mdspan && T::rank() == 2;
templateconcept out-matrix =is-mdspan && T::rank() == 2 && is_assignable_v<typename T::reference, typename T::element_type> && T::is_always_unique();
templateconcept inout-matrix =is-mdspan && T::rank() == 2 && is_assignable_v<typename T::reference, typename T::element_type> && T::is_always_unique();
templateconstexpr bool is-layout-blas-packed = false; // exposition onlytemplate<class Triangle, class StorageOrder>constexpr bool is-layout-blas-packed<layout_blas_packed<Triangle, StorageOrder>> = true;
templateconcept possibly-packed-inout-matrix =is-mdspan && T::rank() == 2 && is_assignable_v<typename T::reference, typename T::element_type> &&(T::is_always_unique() || is-layout-blas-packed);
templateconcept in-object =is-mdspan && (T::rank() == 1 || T::rank() == 2);
templateconcept out-object =is-mdspan && (T::rank() == 1 || T::rank() == 2) && is_assignable_v<typename T::reference, typename T::element_type> && T::is_always_unique();
templateconcept inout-object =is-mdspan && (T::rank() == 1 || T::rank() == 2) && is_assignable_v<typename T::reference, typename T::element_type> && T::is_always_unique();
If a function in [linalg] accesses the elements of a parameter constrained byin-vector,in-matrix, orin-object, those accesses will not modify the elements.
Unless explicitly permitted, anyinout-vector,inout-matrix,inout-object,out-vector,out-matrix,out-object, orpossibly-packed-inout-matrix parameter of a function in [linalg] shall not overlap any other mdspan parameter of the function.
29.9.7.6 Mandates [linalg.helpers.mandates]
[Note 1:
These exposition-only helper functions use the less constraining input concepts even for the output arguments, because the additional constraint for assignability of elements is not necessary, and they are sometimes used in a context where the third argument is an input type too.
â end note]
template<class MDS1, class MDS2>requires(is-mdspan && is-mdspan)constexprbool compatible-static-extents(size_t r1, size_t r2) { // exposition onlyreturn MDS1::static_extent(r1) == dynamic_extent || MDS2::static_extent(r2) == dynamic_extent || MDS1::static_extent(r1) == MDS2::static_extent(r2); }template<in-vector In1, in-vector In2, in-vector Out>constexpr bool possibly-addable() { // exposition onlyreturn compatible-static-extents<Out, In1>(0, 0) &&compatible-static-extents<Out, In2>(0, 0) &&compatible-static-extents<In1, In2>(0, 0); }template<in-matrix In1, in-matrix In2, in-matrix Out>constexpr bool possibly-addable() { // exposition onlyreturn compatible-static-extents<Out, In1>(0, 0) &&compatible-static-extents<Out, In1>(1, 1) &&compatible-static-extents<Out, In2>(0, 0) &&compatible-static-extents<Out, In2>(1, 1) &&compatible-static-extents<In1, In2>(0, 0) &&compatible-static-extents<In1, In2>(1, 1); }template<in-matrix InMat, in-vector InVec, in-vector OutVec>constexpr bool possibly-multipliable() { // exposition onlyreturn compatible-static-extents<OutVec, InMat>(0, 0) &&compatible-static-extents<InMat, InVec>(1, 0); }template<in-vector InVec, in-matrix InMat, in-vector OutVec>constexpr bool possibly-multipliable() { // exposition onlyreturn compatible-static-extents<OutVec, InMat>(0, 1) &&compatible-static-extents<InMat, InVec>(0, 0); }template<in-matrix InMat1, in-matrix InMat2, in-matrix OutMat>constexpr bool possibly-multipliable() { // exposition onlyreturn compatible-static-extents<OutMat, InMat1>(0, 0) &&compatible-static-extents<OutMat, InMat2>(1, 1) &&compatible-static-extents<InMat1, InMat2>(1, 0); }
29.9.7.7 Preconditions [linalg.helpers.precond]
[Note 1:
These exposition-only helper functions use the less constraining input concepts even for the output arguments, because the additional constraint for assignability of elements is not necessary, and they are sometimes used in a context where the third argument is an input type too.
â end note]
constexpr bool addable( // exposition onlyconst in-vector auto& in1, const in-vector auto& in2, const in-vector auto& out) {return out.extent(0) == in1.extent(0) && out.extent(0) == in2.extent(0);}constexpr bool addable( // exposition onlyconst in-matrix auto& in1, const in-matrix auto& in2, const in-matrix auto& out) {return out.extent(0) == in1.extent(0) && out.extent(1) == in1.extent(1) && out.extent(0) == in2.extent(0) && out.extent(1) == in2.extent(1);}constexpr bool multipliable( // exposition onlyconst in-matrix auto& in_mat, const in-vector auto& in_vec, const in-vector auto& out_vec) {return out_vec.extent(0) == in_mat.extent(0) && in_mat.extent(1) == in_vec.extent(0);}constexpr bool multipliable( // exposition onlyconst in-vector auto& in_vec, const in-matrix auto& in_mat, const in-vector auto& out_vec) {return out_vec.extent(0) == in_mat.extent(1) && in_mat.extent(0) == in_vec.extent(0);}constexpr bool multipliable( // exposition onlyconst in-matrix auto& in_mat1, const in-matrix auto& in_mat2, const in-matrix auto& out_mat) {return out_mat.extent(0) == in_mat1.extent(0) && out_mat.extent(1) == in_mat2.extent(1) && in_mat1.extent(1) == in_mat2.extent(0);}
29.9.8 Scaled in-place transformation [linalg.scaled]
29.9.8.1 Introduction [linalg.scaled.intro]
The scaled function takes a value alpha and an mdspan x, and returns a new read-only mdspan that represents the elementwise product of alpha with each element of x.
[Example 1: using Vec = mdspan<double, dextents<size_t, 1>>;
// z = alpha * x + yvoid z_equals_alpha_times_x_plus_y(double alpha, Vec x, Vec y, Vec z) { add(scaled(alpha, x), y, z);}// z = alpha * x + beta * yvoid z_equals_alpha_times_x_plus_beta_times_y(double alpha, Vec x, double beta, Vec y, Vec z) { add(scaled(alpha, x), scaled(beta, y), z);} â end example]
29.9.8.2 Class template scaled_accessor [linalg.scaled.scaledaccessor]
The class template scaled_accessor is an mdspan accessor policy which upon access produces scaled elements.
It is part of the implementation of scaled ([linalg.scaled.scaled]).
namespace std::linalg {template<class ScalingFactor, class NestedAccessor>class scaled_accessor {public:using element_type = add_const_t<decltype(declval() * declvalNestedAccessor::element_type())>; using reference = remove_const_t<element_type>; using data_handle_type = NestedAccessor::data_handle_type; using offset_policy = scaled_accessor<ScalingFactor, NestedAccessor::offset_policy>; constexpr scaled_accessor() = default; templateexplicit(!is_convertible_v<OtherNestedAccessor, NestedAccessor>)constexpr scaled_accessor(const scaled_accessor<ScalingFactor, OtherNestedAccessor>& other); constexpr scaled_accessor(const ScalingFactor& s, const NestedAccessor& a); constexpr reference access(data_handle_type p, size_t i) const; constexpr offset_policy::data_handle_type offset(data_handle_type p, size_t i) const; constexpr const ScalingFactor& scaling_factor() const noexcept { return scaling-factor; }constexpr const NestedAccessor& nested_accessor() const noexcept { return nested-accessor; }private: ScalingFactor scaling-factor{}; // exposition only NestedAccessor nested-accessor{}; // exposition only};}
Mandates:
element_type is valid and denotes a type,
is_copy_constructible_v is true,
is_reference_v<element_type> is false,
ScalingFactor models semiregular, and
NestedAccessor meets the accessor policy requirements ([mdspan.accessor.reqmts]).
template<class OtherNestedAccessor> explicit(!is_convertible_v<OtherNestedAccessor, NestedAccessor>) constexpr scaled_accessor(const scaled_accessor<ScalingFactor, OtherNestedAccessor>& other);
Constraints: is_constructible_v<NestedAccessor, const OtherNestedAccessor&> is true.
Effects:
Direct-non-list-initializes scaling-factor with other.scaling_factor(), and
direct-non-list-initializes nested-accessor with other.nested_accessor().
constexpr scaled_accessor(const ScalingFactor& s, const NestedAccessor& a);
Effects:
Direct-non-list-initializes scaling-factor with s, and
direct-non-list-initializes nested-accessor with a.
constexpr reference access(data_handle_type p, size_t i) const;
Returns: scaling_factor() * NestedAccessor::element_type(nested-accessor.access(p, i))
constexpr offset_policy::data_handle_type offset(data_handle_type p, size_t i) const;
Returns: nested-accessor.offset(p, i)
29.9.8.3 Function template scaled [linalg.scaled.scaled]
The scaled function template takes a scaling factor alpha and an mdspan x, and returns a new read-only mdspan with the same domain as x, that represents the elementwise product of alpha with each element of x.
template<class ScalingFactor, class ElementType, class Extents, class Layout, class Accessor> constexpr auto scaled(ScalingFactor alpha, mdspan<ElementType, Extents, Layout, Accessor> x);
Let SA be scaled_accessor<ScalingFactor, Accessor>.
Returns: mdspan<typename SA::element_type, Extents, Layout, SA>(x.data_handle(), x.mapping(), SA(alpha, x.accessor()))
[Example 1: void test_scaled(mdspan<double, extents<int, 10>> x){auto x_scaled = scaled(5.0, x); for (int i = 0; i < x.extent(0); ++i) { assert(x_scaled[i] == 5.0 * x[i]); }} â end example]
29.9.9 Conjugated in-place transformation [linalg.conj]
29.9.9.1 Introduction [linalg.conj.intro]
The conjugated function takes an mdspan x, and returns a new read-only mdspan y with the same domain as x, whose elements are the complex conjugates of the corresponding elements of x.
29.9.9.2 Class template conjugated_accessor [linalg.conj.conjugatedaccessor]
The class template conjugated_accessor is an mdspan accessor policy which upon access produces conjugate elements.
It is part of the implementation ofconjugated ([linalg.conj.conjugated]).
namespace std::linalg {templateclass conjugated_accessor {public:using element_type = add_const_t<decltype(conj-if-needed(declvalNestedAccessor::element_type()))>; using reference = remove_const_t<element_type>; using data_handle_type = typename NestedAccessor::data_handle_type; using offset_policy = conjugated_accessorNestedAccessor::offset_policy; constexpr conjugated_accessor() = default; templateexplicit(!is_convertible_v<OtherNestedAccessor, NestedAccessor>>)constexpr conjugated_accessor(const conjugated_accessor& other); constexpr reference access(data_handle_type p, size_t i) const; constexpr typename offset_policy::data_handle_type offset(data_handle_type p, size_t i) const; constexpr const NestedAccessor& nested_accessor() const noexcept { return nested-accessor_; }private: NestedAccessor nested-accessor_{}; // exposition only};}
Mandates:
element_type is valid and denotes a type,
is_copy_constructible_v is true,
is_reference_v<element_type> is false, and
NestedAccessor meets the accessor policy requirements ([mdspan.accessor.reqmts]).
constexpr conjugated_accessor(const NestedAccessor& acc);
Effects: Direct-non-list-initializesnested-accessor_ with acc.
template<class OtherNestedAccessor> explicit(!is_convertible_v<OtherNestedAccessor, NestedAccessor>>) constexpr conjugated_accessor(const conjugated_accessor<OtherNestedAccessor>& other);
Constraints: is_constructible_v<NestedAccessor, const OtherNestedAccessor&> is true.
Effects: Direct-non-list-initializes nested-accessor_ with other.nested_accessor().
constexpr reference access(data_handle_type p, size_t i) const;
Returns: conj-if-needed(NestedAccessor::element_type(nested-accessor_.access(p, i)))
constexpr typename offset_policy::data_handle_type offset(data_handle_type p, size_t i) const;
Returns: nested-accessor_.offset(p, i)
29.9.9.3 Function template conjugated [linalg.conj.conjugated]
template<class ElementType, class Extents, class Layout, class Accessor> constexpr auto [conjugated](#lib:conjugated "29.9.9.3 Function template conjugated [linalg.conj.conjugated]")(mdspan<ElementType, Extents, Layout, Accessor> a);
Let A be
remove_cvref_t<decltype(a.accessor().nested_accessor())> if Accessor is a specialization of conjugated_accessor;
otherwise,Accessor if remove_cvref_t is an arithmetic type;
otherwise,conjugated_accessor if the expression conj(E) is valid for any subexpression E whose type is remove_cvref_t with overload resolution performed in a context that includes the declarationtemplate U conj(const U&) = delete;;
otherwise,Accessor.
Returns: Let MD be mdspan<typename A::element_type, Extents, Layout, A>.
MD(a.data_handle(), a.mapping(), a.accessor().nested_accessor()) if Accessor is a
specialization of conjugated_accessor;
otherwise,a, if is_same_v<A, Accessor> is true;
otherwise,MD(a.data_handle(), a.mapping(), conjugated_accessor(a.accessor())).
[Example 1: void test_conjugated_complex(mdspan<complex, extents<int, 10>> a) {auto a_conj = conjugated(a); for (int i = 0; i < a.extent(0); ++i) { assert(a_conj[i] == conj(a[i]); }auto a_conj_conj = conjugated(a_conj); for (int i = 0; i < a.extent(0); ++i) { assert(a_conj_conj[i] == a[i]); }}void test_conjugated_real(mdspan<double, extents<int, 10>> a) {auto a_conj = conjugated(a); for (int i = 0; i < a.extent(0); ++i) { assert(a_conj[i] == a[i]); }auto a_conj_conj = conjugated(a_conj); for (int i = 0; i < a.extent(0); ++i) { assert(a_conj_conj[i] == a[i]); }} â end example]
29.9.10 Transpose in-place transformation [linalg.transp]
29.9.10.1 Introduction [linalg.transp.intro]
layout_transpose is an mdspan layout mapping policy that swaps the two indices, extents, and strides of any unique mdspan layout mapping policy.
The transposed function takes an mdspan representing a matrix, and returns a new mdspan representing the transpose of the input matrix.
29.9.10.2 Exposition-only helpers for layout_transpose and transposed [linalg.transp.helpers]
The exposition-only transpose-extents function takes an extents object representing the extents of a matrix, and returns a new extents object representing the extents of the transpose of the matrix.
The exposition-only alias templatetranspose-extents-t gives the type of transpose-extents(e) for a given extents object e of type InputExtents.
template<class IndexType, size_t InputExtent0, size_t InputExtent1> constexpr extents<IndexType, InputExtent1, InputExtent0> transpose-extents(const extents<IndexType, InputExtent0, InputExtent1>& in); // exposition only
Returns: extents<IndexType, InputExtent1, InputExtent0>(in.extent(1), in.extent(0))
templateusing transpose-extents-t =decltype(transpose-extents(declval())); // exposition only
29.9.10.3 Class template layout_transpose [linalg.transp.layout.transpose]
layout_transpose is an mdspan layout mapping policy that swaps the two indices, extents, and strides of any mdspan layout mapping policy.
namespace std::linalg {templateclass layout_transpose {public:using nested_layout_type = Layout; templatestruct mapping {private:using nested-mapping-type =typename Layout::template mapping<transpose-extents-t>; // exposition onlypublic: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_transpose; constexpr explicit mapping(const nested-mapping-type&); constexpr const extents_type& extents() const noexcept { return extents_; }constexpr index_type required_span_size() const{ return nested-mapping_.required_span_size(); template<class Index0, class Index1>constexpr index_type operator()(Index0 ind0, Index1 ind1) const{ return nested-mapping_(ind1, ind0); }constexpr const nested-mapping-type& nested_mapping() const noexcept{ return nested-mapping_; }static constexpr bool is_always_unique() noexcept{ return nested-mapping-type::is_always_unique(); }static constexpr bool is_always_exhaustive() noexcept{ return nested-mapping-type::is_always_exhaustive(); }static constexpr bool is_always_strided() noexcept{ return nested-mapping-type::is_always_strided(); }constexpr bool is_unique() const { return nested-mapping_.is_unique(); }constexpr bool is_exhaustive() const { return nested-mapping_.is_exhaustive(); }constexpr bool is_strided() const { return nested-mapping_.is_strided(); }constexpr index_type stride(size_t r) const; templatefriend constexpr bool operator==(const mapping& x, const mapping& y); }; private:nested-mapping-type nested-mapping_; // exposition only extents_type extents_; // exposition only};}
Layout shall meet the layout mapping policy requirements ([mdspan.layout.policy.reqmts]).
Mandates:
Extents is a specialization of std::extents, and
Extents::rank() equals 2.
constexpr explicit mapping(const nested-mapping-type& map);
Effects:
Initializes nested-mapping_ with map, and
initializes extents_ with transpose-extents(map.extents()).
constexpr index_type stride(size_t r) const;
Preconditions:
is_strided() is true, and
r < 2 is true.
Returns: nested-mapping_.stride(r == 0 ? 1 : 0)
template<class OtherExtents> friend constexpr bool operator==(const mapping& x, const mapping<OtherExtents>& y);
Constraints: The expressionx.nested-mapping_ == y.nested-mapping_ is well-formed and its result is convertible to bool.
Returns: x.nested-mapping_ == y.nested-mapping_.
29.9.10.4 Function template transposed [linalg.transp.transposed]
The transposed function takes a rank-2 mdspan representing a matrix, and returns a new mdspan representing the input matrix's transpose.
The input matrix's data are not modified, and the returned mdspan accesses the input matrix's data in place.
template<class ElementType, class Extents, class Layout, class Accessor> constexpr auto transposed(mdspan<ElementType, Extents, Layout, Accessor> a);
Mandates: Extents::rank() == 2 is true.
Let ReturnExtents betranspose-extents-t.
Let R bemdspan<ElementType, ReturnExtents, ReturnLayout, Accessor>, where ReturnLayout is:
layout_right if Layout is layout_left;
otherwise, layout_left if Layout is layout_right;
otherwise, layout_right_padded if Layout is
layout_left_padded for some size_t value PaddingValue;
otherwise, layout_left_padded if Layout is
layout_right_padded for some size_t value PaddingValue;
otherwise, layout_stride if Layout is layout_stride;
otherwise,layout_blas_packed<OppositeTriangle, OppositeStorageOrder>,
if Layout is
layout_blas_packed<Triangle, StorageOrder> for some Triangle and StorageOrder, where
OppositeTriangle isconditional_t<is_same_v<Triangle, upper_triangle_t>, lower_triangle_t, upper_triangle_t> and
OppositeStorageOrder isconditional_t<is_same_v<StorageOrder, column_major_t>, row_major_t, column_major_t>
otherwise, NestedLayout if Layout is layout_transpose for some NestedLayout;
otherwise, layout_transpose.
Returns: With ReturnMapping being the type typename ReturnLayout::template mapping:
if Layout is layout_left, layout_right, or a specialization of layout_blas_packed,R(a.data_handle(), ReturnMapping(transpose-extents(a.mapping().extents())), a.accessor())
otherwise,R(a.data_handle(), ReturnMapping(transpose-extents(a.mapping().extents()), a.mapping().stride(1)), a.accessor()) if Layout is layout_left_padded for some size_t value PaddingValue;
otherwise,R(a.data_handle(), ReturnMapping(transpose-extents(a.mapping().extents()), a.mapping().stride(0)), a.accessor()) if Layout is layout_right_padded for some size_t value PaddingValue;
otherwise, if Layout is layout_stride,R(a.data_handle(), ReturnMapping(transpose-extents(a.mapping().extents()), array{a.mapping().stride(1), a.mapping().stride(0)}), a.accessor())
otherwise, if Layout is a specialization of layout_transpose,R(a.data_handle(), a.mapping().nested_mapping(), a.accessor())
otherwise,R(a.data_handle(), ReturnMapping(a.mapping()), a.accessor())
[Example 1: void test_transposed(mdspan<double, extents<size_t, 3, 4>> a) {const auto num_rows = a.extent(0); const auto num_cols = a.extent(1); auto a_t = transposed(a); assert(num_rows == a_t.extent(1)); assert(num_cols == a_t.extent(0)); assert(a.stride(0) == a_t.stride(1)); assert(a.stride(1) == a_t.stride(0)); for (size_t row = 0; row < num_rows; ++row) {for (size_t col = 0; col < num_rows; ++col) { assert(a[row, col] == a_t[col, row]); }}auto a_t_t = transposed(a_t); assert(num_rows == a_t_t.extent(0)); assert(num_cols == a_t_t.extent(1)); assert(a.stride(0) == a_t_t.stride(0)); assert(a.stride(1) == a_t_t.stride(1)); for (size_t row = 0; row < num_rows; ++row) {for (size_t col = 0; col < num_rows; ++col) { assert(a[row, col] == a_t_t[row, col]); }}} â end example]
29.9.11 Conjugate transpose in-place transform [linalg.conjtransposed]
The conjugate_transposed function returns a conjugate transpose view of an object.
This combines the effects of transposed and conjugated.
template<class ElementType, class Extents, class Layout, class Accessor> constexpr auto conjugate_transposed(mdspan<ElementType, Extents, Layout, Accessor> a);
Effects: Equivalent to: return conjugated(transposed(a));
[Example 1: void test_conjugate_transposed(mdspan<complex, extents<size_t, 3, 4>> a) {const auto num_rows = a.extent(0); const auto num_cols = a.extent(1); auto a_ct = conjugate_transposed(a); assert(num_rows == a_ct.extent(1)); assert(num_cols == a_ct.extent(0)); assert(a.stride(0) == a_ct.stride(1)); assert(a.stride(1) == a_ct.stride(0)); for (size_t row = 0; row < num_rows; ++row) {for (size_t col = 0; col < num_rows; ++col) { assert(a[row, col] == conj(a_ct[col, row])); }}auto a_ct_ct = conjugate_transposed(a_ct); assert(num_rows == a_ct_ct.extent(0)); assert(num_cols == a_ct_ct.extent(1)); assert(a.stride(0) == a_ct_ct.stride(0)); assert(a.stride(1) == a_ct_ct.stride(1)); for (size_t row = 0; row < num_rows; ++row) {for (size_t col = 0; col < num_rows; ++col) { assert(a[row, col] == a_ct_ct[row, col]); assert(conj(a_ct[col, row]) == a_ct_ct[row, col]); }}} â end example]
29.9.12 Algorithm requirements based on template parameter name [linalg.algs.reqs]
Throughout[linalg.algs.blas1], [linalg.algs.blas2], and [linalg.algs.blas3], where the template parameters are not constrained, the names of template parameters are used to express the following constraints.
-
is_execution_policy::value is true ([execpol.type]).
-
Real is any type such that complex is specified ([complex.numbers.general]).
-
Triangle is either upper_triangle_t or lower_triangle_t.
-
DiagonalStorage is either implicit_unit_diagonal_t or explicit_diagonal_t.
[Note 1:
Function templates that have a template parameter named ExecutionPolicy are parallel algorithms ([algorithms.parallel.defns]).
â end note]
29.9.13 BLAS 1 algorithms [linalg.algs.blas1]
29.9.13.1 Complexity [linalg.algs.blas1.complexity]
Complexity: All algorithms in [linalg.algs.blas1] with mdspan parameters perform a count of mdspan array accesses and arithmetic operations that is linear in the maximum product of extents of any mdspan parameter.
29.9.13.2 Givens rotations [linalg.algs.blas1.givens]
29.9.13.2.1 Compute Givens rotation [linalg.algs.blas1.givens.lartg]
`template setup_givens_rotation_result setup_givens_rotation(Real a, Real b) noexcept;
template setup_givens_rotation_result<complex> setup_givens_rotation(complex a, complex b) noexcept; `
These functions compute the Givens plane rotation
represented by the two values c and s such that the 2 x 2 system of equations
[cs⯯¯sc]â[ab]=[r0]
holds, where c is always a real scalar, and c2+|s|2=1.
That is, c and s represent a 2 x 2 matrix, that when multiplied by the right by the input vector whose components are a and b, produces a result vector whose first component r is the Euclidean norm of the input vector, and whose second component is zero.
[Note 1:
These functions correspond to the LAPACK function xLARTG[bib].
â end note]
Returns: c, s, r, where c and s form the Givens plane rotation corresponding to the input a and b, and r is the Euclidean norm of the two-component vector formed by a and b.
29.9.13.2.2 Apply a computed Givens rotation to vectors [linalg.algs.blas1.givens.rot]
template<[inout-vector](#concept:inout-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutVec1, [inout-vector](#concept:inout-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutVec2, class Real> void apply_givens_rotation(InOutVec1 x, InOutVec2 y, Real c, Real s); template<class ExecutionPolicy, [inout-vector](#concept:inout-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutVec1, [inout-vector](#concept:inout-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutVec2, class Real> void apply_givens_rotation(ExecutionPolicy&& exec, InOutVec1 x, InOutVec2 y, Real c, Real s); template<[inout-vector](#concept:inout-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutVec1, [inout-vector](#concept:inout-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutVec2, class Real> void apply_givens_rotation(InOutVec1 x, InOutVec2 y, Real c, complex<Real> s); template<class ExecutionPolicy, [inout-vector](#concept:inout-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutVec1, [inout-vector](#concept:inout-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutVec2, class Real> void apply_givens_rotation(ExecutionPolicy&& exec, InOutVec1 x, InOutVec2 y, Real c, complex<Real> s);
[Note 1:
These functions correspond to the BLAS function xROT[bib].
â end note]
Mandates: compatible-static-extents<InOutVec1, InOutVec2>(0, 0) is true.
Preconditions: x.extent(0) equals y.extent(0).
Effects: Applies the plane rotation specified by c and s to the input vectors x and y, as if the rotation were a 2 x 2 matrix and the input vectors were successive rows of a matrix with two rows.
29.9.13.3 Swap matrix or vector elements [linalg.algs.blas1.swap]
template<[inout-object](#concept:inout-object "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutObj1, [inout-object](#concept:inout-object "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutObj2> void swap_elements(InOutObj1 x, InOutObj2 y); template<class ExecutionPolicy, [inout-object](#concept:inout-object "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutObj1, [inout-object](#concept:inout-object "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutObj2> void swap_elements(ExecutionPolicy&& exec, InOutObj1 x, InOutObj2 y);
[Note 1:
These functions correspond to the BLAS function xSWAP[bib].
â end note]
Constraints: x.rank() equals y.rank().
Mandates: For all r in the range [0, x.rank()),compatible-static-extents<InOutObj1, InOutObj2>(r, r) is true.
Preconditions: x.extents() equals y.extents().
Effects: Swaps all corresponding elements of x and y.
29.9.13.4 Multiply the elements of an object in place by a scalar [linalg.algs.blas1.scal]
template<class Scalar, [inout-object](#concept:inout-object "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutObj> void scale(Scalar alpha, InOutObj x); template<class ExecutionPolicy, class Scalar, [inout-object](#concept:inout-object "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutObj> void scale(ExecutionPolicy&& exec, Scalar alpha, InOutObj x);
[Note 1:
These functions correspond to the BLAS function xSCAL[bib].
â end note]
Effects: Overwrites x with the result of computing the elementwise multiplication αx, where the scalar α is alpha.
29.9.13.5 Copy elements of one matrix or vector into another [linalg.algs.blas1.copy]
template<[in-object](#concept:in-object "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InObj, [out-object](#concept:out-object "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutObj> void copy(InObj x, OutObj y); template<class ExecutionPolicy, [in-object](#concept:in-object "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InObj, [out-object](#concept:out-object "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutObj> void copy(ExecutionPolicy&& exec, InObj x, OutObj y);
[Note 1:
These functions correspond to the BLAS function xCOPY[bib].
â end note]
Constraints: x.rank() equals y.rank().
Mandates: For all r in the range [0,x.rank()),compatible-static-extents<InObj, OutObj>(r, r) is true.
Preconditions: x.extents() equals y.extents().
Effects: Assigns each element of x to the corresponding element of y.
29.9.13.6 Add vectors or matrices elementwise [linalg.algs.blas1.add]
template<[in-object](#concept:in-object "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InObj1, [in-object](#concept:in-object "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InObj2, [out-object](#concept:out-object "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutObj> void add(InObj1 x, InObj2 y, OutObj z); template<class ExecutionPolicy, [in-object](#concept:in-object "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InObj1, [in-object](#concept:in-object "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InObj2, [out-object](#concept:out-object "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutObj> void add(ExecutionPolicy&& exec, InObj1 x, InObj2 y, OutObj z);
[Note 1:
These functions correspond to the BLAS function xAXPY[bib].
â end note]
Constraints: x.rank(), y.rank(), and z.rank() are all equal.
Mandates: possibly-addable<InObj1, InObj2, OutObj>() is true.
Preconditions: addable(x,y,z) is true.
Effects: Computes z=x+y.
Remarks: z may alias x or y.
29.9.13.7 Dot product of two vectors [linalg.algs.blas1.dot]
[Note 1:
The functions in this section correspond to the BLAS functions xDOT, xDOTU, and xDOTC[bib].
â end note]
The following elements apply to all functions in [linalg.algs.blas1.dot].
Mandates: compatible-static-extents<InVec1, InVec2>(0, 0) is true.
Preconditions: v1.extent(0) equals v2.extent(0).
template<[in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, class Scalar> Scalar dot(InVec1 v1, InVec2 v2, Scalar init); template<class ExecutionPolicy, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, class Scalar> Scalar dot(ExecutionPolicy&& exec, InVec1 v1, InVec2 v2, Scalar init);
These functions compute a non-conjugated dot product with an explicitly specified result type.
Returns: Let N be v1.extent(0).
init if N is zero;
otherwise,GENERALIZED_SUM(plus<>(), init, v1[0]*v2[0], …, v1[N-1]*v2[N-1]).
Remarks: If InVec1::value_type, InVec2::value_type, and Scalar are all floating-point types or specializations of complex, and if Scalar has higher precision than InVec1::value_type or InVec2::value_type, then intermediate terms in the sum use Scalar's precision or greater.
template<[in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2> auto dot(InVec1 v1, InVec2 v2); template<class ExecutionPolicy, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2> auto dot(ExecutionPolicy&& exec, InVec1 v1, InVec2 v2);
These functions compute a non-conjugated dot product with a default result type.
Effects: Let T bedecltype(declval<typename InVec1::value_type>() * declval<typename InVec2::value_type>()).
Then,
the two-parameter overload is equivalent to:return dot(v1, v2, T{}); and
the three-parameter overload is equivalent to:return dot(std::forward(exec), v1, v2, T{});
template<[in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, class Scalar> Scalar dotc(InVec1 v1, InVec2 v2, Scalar init); template<class ExecutionPolicy, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, class Scalar> Scalar dotc(ExecutionPolicy&& exec, InVec1 v1, InVec2 v2, Scalar init);
These functions compute a conjugated dot product with an explicitly specified result type.
Effects:
The three-parameter overload is equivalent to:return dot(conjugated(v1), v2, init); and
the four-parameter overload is equivalent to:return dot(std::forward(exec), conjugated(v1), v2, init);
template<[in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2> auto dotc(InVec1 v1, InVec2 v2); template<class ExecutionPolicy, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2> auto dotc(ExecutionPolicy&& exec, InVec1 v1, InVec2 v2);
These functions compute a conjugated dot product with a default result type.
Effects: Let T be decltype(conj-if-needed(declval<typename InVec1::value_type>()) * declval<typename InVec2::value_type>()).
Then,
the two-parameter overload is equivalent to:return dotc(v1, v2, T{}); and
the three-parameter overload is equivalent toreturn dotc(std::forward(exec), v1, v2, T{});
29.9.13.8 Scaled sum of squares of a vector's elements [linalg.algs.blas1.ssq]
template<[in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, class Scalar> sum_of_squares_result<Scalar> vector_sum_of_squares(InVec v, sum_of_squares_result<Scalar> init); template<class ExecutionPolicy, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, class Scalar> sum_of_squares_result<Scalar> vector_sum_of_squares(ExecutionPolicy&& exec, InVec v, sum_of_squares_result<Scalar> init);
[Note 1:
These functions correspond to the LAPACK function xLASSQ[bib].
â end note]
Mandates: decltype(abs-if-needed(declval<typename InVec::value_type>())) is convertible to Scalar.
Effects: Returns a value result such that
result.scaling_factor is the maximum of init.scaling_factor andabs-if-needed(x[i]) for all i in the domain of v; and
let s2init beinit.scaling_factor * init.scaling_factor * init.scaled_sum_of_squares then result.scaling_factor * result.scaling_factor * result.scaled_sum_of_squares equals the sum of s2init and the squares of abs-if-needed(x[i]) for all i in the domain of v.
Remarks: If InVec::value_type, and Scalar are all floating-point types or specializations of complex, and if Scalar has higher precision than InVec::value_type, then intermediate terms in the sum use Scalar's precision or greater.
29.9.13.9 Euclidean norm of a vector [linalg.algs.blas1.nrm2]
template<[in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, class Scalar> Scalar vector_two_norm(InVec v, Scalar init); template<class ExecutionPolicy, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, class Scalar> Scalar vector_two_norm(ExecutionPolicy&& exec, InVec v, Scalar init);
[Note 1:
These functions correspond to the BLAS function xNRM2[bib].
â end note]
Mandates: Let a beabs-if-needed(declval<typename InVec::value_type>()).
Then, decltype(init + a * a is convertible to Scalar.
Returns: The square root of the sum of the square of init and the squares of the absolute values of the elements of v.
[Note 2:
For init equal to zero, this is the Euclidean norm (also called 2-norm) of the vector v.
â end note]
Remarks: If InVec::value_type, and Scalar are all floating-point types or specializations of complex, and if Scalar has higher precision than InVec::value_type, then intermediate terms in the sum use Scalar's precision or greater.
[Note 3:
An implementation of this function for floating-point types T can use the scaled_sum_of_squares result fromvector_sum_of_squares(x, {.scaling_factor=1.0, .scaled_sum_of_squares=init}).
â end note]
template<[in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec> auto vector_two_norm(InVec v); template<class ExecutionPolicy, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec> auto vector_two_norm(ExecutionPolicy&& exec, InVec v);
Effects: Let a beabs-if-needed(declval<typename InVec::value_type>()).
Let T be decltype(a * a).
Then,
the one-parameter overload is equivalent to:return vector_two_norm(v, T{}); and
the two-parameter overload is equivalent to:return vector_two_norm(std::forward(exec), v, T{});
29.9.13.10 Sum of absolute values of vector elements [linalg.algs.blas1.asum]
template<[in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, class Scalar> Scalar vector_abs_sum(InVec v, Scalar init); template<class ExecutionPolicy, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, class Scalar> Scalar vector_abs_sum(ExecutionPolicy&& exec, InVec v, Scalar init);
[Note 1:
These functions correspond to the BLAS functionsSASUM, DASUM, SCASUM, and DZASUM[bib].
â end note]
Mandates: decltype(init + abs-if-needed(real-if-needed(declval())) +abs-if-needed(imag-if-needed(declval()))) is convertible to Scalar.
Returns: Let N be v.extent(0).
init if N is zero;
otherwise, if InVec::value_type is an arithmetic type,GENERALIZED_SUM(plus<>(), init, abs-if-needed(v[0]), …, abs-if-needed(v[N-1]))
otherwise,GENERALIZED_SUM(plus<>(), init, abs-if-needed(real-if-needed(v[0])) + abs-if-needed(imag-if-needed(v[0])), …, abs-if-needed(real-if-needed(v[N-1])) + abs-if-needed(imag-if-needed(v[N-1])))
Remarks: If InVec::value_type and Scalar are all floating-point types or specializations of complex, and if Scalar has higher precision than InVec::value_type, then intermediate terms in the sum use Scalar's precision or greater.
template<[in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec> auto vector_abs_sum(InVec v); template<class ExecutionPolicy, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec> auto vector_abs_sum(ExecutionPolicy&& exec, InVec v);
Effects: Let T be typename InVec::value_type.
Then,
the one-parameter overload is equivalent to:return vector_abs_sum(v, T{}); and
the two-parameter overload is equivalent to:return vector_abs_sum(std::forward(exec), v, T{});
29.9.13.11 Index of maximum absolute value of vector elements [linalg.algs.blas1.iamax]
template<[in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec> typename InVec::extents_type vector_idx_abs_max(InVec v); template<class ExecutionPolicy, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec> typename InVec::extents_type vector_idx_abs_max(ExecutionPolicy&& exec, InVec v);
[Note 1:
These functions correspond to the BLAS function IxAMAX[bib].
â end note]
Let T bedecltype(abs-if-needed(real-if-needed(declval())) +abs-if-needed(imag-if-needed(declval())))
Mandates: declval() < declval() is a valid expression.
Returns:
numeric_limits<typename InVec::size_type>::max() if v has zero elements;
otherwise, the index of the first element of v having largest absolute value, if InVec::value_type is an arithmetic type;
otherwise, the index of the first element ve of v for whichabs-if-needed(real-if-needed(ve)) + abs-if-needed(imag-if-needed(ve)) has the largest value.
29.9.13.12 Frobenius norm of a matrix [linalg.algs.blas1.matfrobnorm]
[Note 1:
These functions exist in the BLAS standard[bib] but are not part of the reference implementation.
â end note]
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Scalar> Scalar matrix_frob_norm(InMat A, Scalar init); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Scalar> Scalar matrix_frob_norm(ExecutionPolicy&& exec, InMat A, Scalar init);
Mandates: Let a beabs-if-needed(declval<typename InMat::value_type>()).
Then, decltype(init + a * a) is convertible to Scalar.
Returns: The square root of the sum of squares of init and the absolute values of the elements of A.
[Note 2:
For init equal to zero, this is the Frobenius norm of the matrix A.
â end note]
Remarks: If InMat::value_type and Scalar are all floating-point types or specializations of complex, and if Scalar has higher precision than InMat::value_type, then intermediate terms in the sum use Scalar's precision or greater.
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat> auto matrix_frob_norm(InMat A); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat> auto matrix_frob_norm(ExecutionPolicy&& exec, InMat A);
Effects: Let a beabs-if-needed(declval<typename InMat::value_type>()).
Let T bedecltype(a * a).
Then,
the one-parameter overload is equivalent to:return matrix_frob_norm(A, T{}); and
the two-parameter overload is equivalent to:return matrix_frob_norm(std::forward(exec), A, T{});
29.9.13.13 One norm of a matrix [linalg.algs.blas1.matonenorm]
[Note 1:
These functions exist in the BLAS standard[bib] but are not part of the reference implementation.
â end note]
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Scalar> Scalar matrix_one_norm(InMat A, Scalar init); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Scalar> Scalar matrix_one_norm(ExecutionPolicy&& exec, InMat A, Scalar init);
Mandates: decltype(abs-if-needed(declval<typename InMat::value_type>())) is convertible to Scalar.
Returns:
init if A.extent(1) is zero;
otherwise, the sum of init and the one norm of the matrix A.
[Note 2:
The one norm of the matrix A is the maximum over all columns of A, of the sum of the absolute values of the elements of the column.
â end note]
Remarks: If InMat::value_type and Scalar are all floating-point types or specializations of complex, and if Scalar has higher precision than InMat::value_type, then intermediate terms in the sum use Scalar's precision or greater.
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat> auto matrix_one_norm(InMat A); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat> auto matrix_one_norm(ExecutionPolicy&& exec, InMat A);
Effects: Let T bedecltype(abs-if-needed(declval<typename InMat::value_type>()).
Then,
the one-parameter overload is equivalent to:return matrix_one_norm(A, T{}); and
the two-parameter overload is equivalent to:return matrix_one_norm(std::forward(exec), A, T{});
29.9.13.14 Infinity norm of a matrix [linalg.algs.blas1.matinfnorm]
[Note 1:
These functions exist in the BLAS standard[bib] but are not part of the reference implementation.
â end note]
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Scalar> Scalar matrix_inf_norm(InMat A, Scalar init); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Scalar> Scalar matrix_inf_norm(ExecutionPolicy&& exec, InMat A, Scalar init);
Mandates: decltype(abs-if-needed(declval<typename InMat::value_type>())) is convertible to Scalar.
Returns:
init if A.extent(0) is zero;
otherwise, the sum of init and the infinity norm of the matrix A.
[Note 2:
The infinity norm of the matrix A is the maximum over all rows of A, of the sum of the absolute values of the elements of the row.
â end note]
Remarks: If InMat::value_type and Scalar are all floating-point types or specializations of complex, and if Scalar has higher precision than InMat::value_type, then intermediate terms in the sum use Scalar's precision or greater.
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat> auto matrix_inf_norm(InMat A); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat> auto matrix_inf_norm(ExecutionPolicy&& exec, InMat A);
Effects: Let T bedecltype(abs-if-needed(declval<typename InMat::value_type>()).
Then,
the one-parameter overload is equivalent to:return matrix_inf_norm(A, T{}); and
the two-parameter overload is equivalent to:return matrix_inf_norm(std::forward(exec), A, T{});
29.9.14 BLAS 2 algorithms [linalg.algs.blas2]
29.9.14.1 General matrix-vector product [linalg.algs.blas2.gemv]
[Note 1:
These functions correspond to the BLAS function xGEMV.
â end note]
The following elements apply to all functions in [linalg.algs.blas2.gemv].
Mandates:
possibly-multipliable<decltype(A), decltype(x), decltype(y)>() is true, and
possibly-addable<decltype(x), decltype(y), decltype(z)>() is true for those overloads that take a z parameter.
Preconditions:
multipliable(A,x,y) is true, and
addable(x,y,z) is true for those overloads that take a z parameter.
Complexity: O(x.extent(0)ÃA.extent(1)).
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec> void matrix_vector_product(InMat A, InVec x, OutVec y); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec> void matrix_vector_product(ExecutionPolicy&& exec, InMat A, InVec x, OutVec y);
These functions perform an overwriting matrix-vector product.
Effects: Computes y=Ax.
[Example 1: constexpr size_t num_rows = 5;constexpr size_t num_cols = 6;
// y = 3.0 * A * xvoid scaled_matvec_1(mdspan<double, extents<size_t, num_rows, num_cols>> A, mdspan<double, extents<size_t, num_cols>> x, mdspan<double, extents<size_t, num_rows>> y) { matrix_vector_product(scaled(3.0, A), x, y);}// z = 7.0 times the transpose of A, times yvoid scaled_transposed_matvec(mdspan<double, extents<size_t, num_rows, num_cols>> A, mdspan<double, extents<size_t, num_rows>> y, mdspan<double, extents<size_t, num_cols>> z) { matrix_vector_product(scaled(7.0, transposed(A)), y, z);} â end example]
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec> void matrix_vector_product(InMat A, InVec1 x, InVec2 y, OutVec z); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec> void matrix_vector_product(ExecutionPolicy&& exec, InMat A, InVec1 x, InVec2 y, OutVec z);
These functions perform an updating matrix-vector product.
Effects: Computes z=y+Ax.
Remarks: z may alias y.
[Example 2: // y = 3.0 * A * x + 2.0 * yvoid scaled_matvec_2(mdspan<double, extents<size_t, num_rows, num_cols>> A, mdspan<double, extents<size_t, num_cols>> x, mdspan<double, extents<size_t, num_rows>> y) { matrix_vector_product(scaled(3.0, A), x, scaled(2.0, y), y);} â end example]
29.9.14.2 Symmetric matrix-vector product [linalg.algs.blas2.symv]
[Note 1:
These functions correspond to the BLAS functionsxSYMV and xSPMV[bib].
â end note]
The following elements apply to all functions in [linalg.algs.blas2.symv].
Mandates:
If InMat has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
compatible-static-extents<decltype(A), decltype(A)>(0, 1) is true;
possibly-multipliable<decltype(A), decltype(x), decltype(y)>() is true; and
possibly-addable<decltype(x), decltype(y), decltype(z)>() is true for those overloads that take a z parameter.
Preconditions:
A.extent(0) equals A.extent(1),
multipliable(A,x,y) is true, and
addable(x,y,z) is true for those overloads that take a z parameter.
Complexity: O(x.extent(0)ÃA.extent(1)).
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec> void symmetric_matrix_vector_product(InMat A, Triangle t, InVec x, OutVec y); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec> void symmetric_matrix_vector_product(ExecutionPolicy&& exec, InMat A, Triangle t, InVec x, OutVec y);
These functions perform an overwriting symmetric matrix-vector product, taking into account the Triangle parameter that applies to the symmetric matrix A ([linalg.general]).
Effects: Computes y=Ax.
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec> void symmetric_matrix_vector_product(InMat A, Triangle t, InVec1 x, InVec2 y, OutVec z); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec> void symmetric_matrix_vector_product(ExecutionPolicy&& exec, InMat A, Triangle t, InVec1 x, InVec2 y, OutVec z);
These functions perform an updating symmetric matrix-vector product, taking into account the Triangle parameter that applies to the symmetric matrix A ([linalg.general]).
Effects: Computes z=y+Ax.
Remarks: z may alias y.
29.9.14.3 Hermitian matrix-vector product [linalg.algs.blas2.hemv]
[Note 1:
These functions correspond to the BLAS functionsxHEMV and xHPMV[bib].
â end note]
The following elements apply to all functions in [linalg.algs.blas2.hemv].
Mandates:
If InMat has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
compatible-static-extents<decltype(A), decltype(A)>(0, 1) is true;
possibly-multipliable<decltype(A), decltype(x), decltype(y)>() is true; and
possibly-addable<decltype(x), decltype(y), decltype(z)>() is true for those overloads that take a z parameter.
Preconditions:
A.extent(0) equals A.extent(1),
multipliable(A, x, y) is true, and
addable(x, y, z) is true for those overloads that take a z parameter.
Complexity: O(x.extent(0)ÃA.extent(1)).
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec> void hermitian_matrix_vector_product(InMat A, Triangle t, InVec x, OutVec y); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec> void hermitian_matrix_vector_product(ExecutionPolicy&& exec, InMat A, Triangle t, InVec x, OutVec y);
These functions perform an overwriting Hermitian matrix-vector product, taking into account the Triangle parameter that applies to the Hermitian matrix A ([linalg.general]).
Effects: Computes y=Ax.
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec> void hermitian_matrix_vector_product(InMat A, Triangle t, InVec1 x, InVec2 y, OutVec z); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec> void hermitian_matrix_vector_product(ExecutionPolicy&& exec, InMat A, Triangle t, InVec1 x, InVec2 y, OutVec z);
These functions perform an updating Hermitian matrix-vector product, taking into account the Triangle parameter that applies to the Hermitian matrix A ([linalg.general]).
Effects: Computes z=y+Ax.
Remarks: z may alias y.
29.9.14.4 Triangular matrix-vector product [linalg.algs.blas2.trmv]
[Note 1:
These functions correspond to the BLAS functionsxTRMV and xTPMV[bib].
â end note]
The following elements apply to all functions in [linalg.algs.blas2.trmv].
Mandates:
If InMat has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
compatible-static-extents<decltype(A), decltype(A)>(0, 1) is true;
compatible-static-extents<decltype(A), decltype(y)>(0, 0) is true;
compatible-static-extents<decltype(A), decltype(x)>(0, 0) is true for those overloads that take an x parameter; and
compatible-static-extents<decltype(A), decltype(z)>(0, 0) is true for those overloads that take a z parameter.
Preconditions:
A.extent(0) equals A.extent(1),
A.extent(0) equals y.extent(0),
A.extent(0) equals x.extent(0) for those overloads that take an x parameter, and
A.extent(0) equals z.extent(0) for those overloads that take a z parameter.
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec> void triangular_matrix_vector_product(InMat A, Triangle t, DiagonalStorage d, InVec x, OutVec y); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec> void triangular_matrix_vector_product(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InVec x, OutVec y);
These functions perform an overwriting triangular matrix-vector product, taking into account the Triangle and DiagonalStorage parameters that apply to the triangular matrix A ([linalg.general]).
Effects: Computes y=Ax.
Complexity: O(x.extent(0)ÃA.extent(1)).
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [inout-vector](#concept:inout-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutVec> void triangular_matrix_vector_product(InMat A, Triangle t, DiagonalStorage d, InOutVec y); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [inout-vector](#concept:inout-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutVec> void triangular_matrix_vector_product(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InOutVec y);
These functions perform an in-place triangular matrix-vector product, taking into account the Triangle and DiagonalStorage parameters that apply to the triangular matrix A ([linalg.general]).
[Note 2:
Performing this operation in place hinders parallelization.
However, other ExecutionPolicy specific optimizations, such as vectorization, are still possible.
â end note]
Effects: Computes a vector yâ² such that yâ²=Ay, and assigns each element of yâ² to the corresponding element of y.
Complexity: O(y.extent(0)ÃA.extent(1)).
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec> void triangular_matrix_vector_product(InMat A, Triangle t, DiagonalStorage d, InVec1 x, InVec2 y, OutVec z); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec> void triangular_matrix_vector_product(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InVec1 x, InVec2 y, OutVec z);
These functions perform an updating triangular matrix-vector product, taking into account the Triangle and DiagonalStorage parameters that apply to the triangular matrix A ([linalg.general]).
Effects: Computes z=y+Ax.
Complexity: O(x.extent(0)ÃA.extent(1)).
Remarks: z may alias y.
29.9.14.5 Solve a triangular linear system [linalg.algs.blas2.trsv]
[Note 1:
These functions correspond to the BLAS functionsxTRSV and xTPSV[bib].
â end note]
The following elements apply to all functions in [linalg.algs.blas2.trsv].
Mandates:
If InMat has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
compatible-static-extents<decltype(A), decltype(A)>(0, 1) is true;
compatible-static-extents<decltype(A), decltype(b)>(0, 0) is true; and
compatible-static-extents<decltype(A), decltype(x)>(0, 0) is true for those overloads that take an x parameter.
Preconditions:
A.extent(0) equals A.extent(1),
A.extent(0) equals b.extent(0), and
A.extent(0) equals x.extent(0) for those overloads that take an x parameter.
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec, class BinaryDivideOp> void triangular_matrix_vector_solve(InMat A, Triangle t, DiagonalStorage d, InVec b, OutVec x, BinaryDivideOp divide); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec, class BinaryDivideOp> void triangular_matrix_vector_solve(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InVec b, OutVec x, BinaryDivideOp divide);
These functions perform a triangular solve, taking into account the Triangle and DiagonalStorage parameters that apply to the triangular matrix A ([linalg.general]).
Effects: Computes a vector xâ² such that b=Axâ², and assigns each element of xâ² to the corresponding element of x.
If no such xâ² exists, then the elements of x are valid but unspecified.
Complexity: O(A.extent(1)Ãb.extent(0)).
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec> void triangular_matrix_vector_solve(InMat A, Triangle t, DiagonalStorage d, InVec b, OutVec x);
Effects: Equivalent to:triangular_matrix_vector_solve(A, t, d, b, x, divides{});
template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [out-vector](#concept:out-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutVec> void triangular_matrix_vector_solve(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InVec b, OutVec x);
Effects: Equivalent to:triangular_matrix_vector_solve(std::forward(exec), A, t, d, b, x, divides{});
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [inout-vector](#concept:inout-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutVec, class BinaryDivideOp> void triangular_matrix_vector_solve(InMat A, Triangle t, DiagonalStorage d, InOutVec b, BinaryDivideOp divide); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [inout-vector](#concept:inout-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutVec, class BinaryDivideOp> void triangular_matrix_vector_solve(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InOutVec b, BinaryDivideOp divide);
These functions perform an in-place triangular solve, taking into account the Triangle and DiagonalStorage parameters that apply to the triangular matrix A ([linalg.general]).
[Note 2:
Performing triangular solve in place hinders parallelization.
However, other ExecutionPolicy specific optimizations, such as vectorization, are still possible.
â end note]
Effects: Computes a vector xâ² such that b=Axâ², and assigns each element of xâ² to the corresponding element of b.
If no such xâ² exists, then the elements of b are valid but unspecified.
Complexity: O(A.extent(1)Ãb.extent(0)).
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [inout-vector](#concept:inout-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutVec> void triangular_matrix_vector_solve(InMat A, Triangle t, DiagonalStorage d, InOutVec b);
Effects: Equivalent to:triangular_matrix_vector_solve(A, t, d, b, divides{});
template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [inout-vector](#concept:inout-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutVec> void triangular_matrix_vector_solve(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InOutVec b);
Effects: Equivalent to:triangular_matrix_vector_solve(std::forward(exec), A, t, d, b, divides{});
29.9.14.6 Rank-1 (outer product) update of a matrix [linalg.algs.blas2.rank1]
template<[in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, [inout-matrix](#concept:inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat> void matrix_rank_1_update(InVec1 x, InVec2 y, InOutMat A); template<class ExecutionPolicy, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, [inout-matrix](#concept:inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat> void matrix_rank_1_update(ExecutionPolicy&& exec, InVec1 x, InVec2 y, InOutMat A);
These functions perform a nonsymmetric nonconjugated rank-1 update.
[Note 1:
These functions correspond to the BLAS functionsxGER (for real element types) andxGERU (for complex element types)[bib].
â end note]
Mandates: possibly-multipliable<InOutMat, InVec2, InVec1>() is true.
Preconditions: multipliable(A, y, x) is true.
Effects: Computes a matrix Aâ² such that Aâ²=A+xyT, and assigns each element of Aâ² to the corresponding element of A.
Complexity: O(x.extent(0)Ãy.extent(0)).
template<[in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, [inout-matrix](#concept:inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat> void matrix_rank_1_update_c(InVec1 x, InVec2 y, InOutMat A); template<class ExecutionPolicy, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, [inout-matrix](#concept:inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat> void matrix_rank_1_update_c(ExecutionPolicy&& exec, InVec1 x, InVec2 y, InOutMat A);
These functions perform a nonsymmetric conjugated rank-1 update.
[Note 2:
These functions correspond to the BLAS functionsxGER (for real element types) andxGERC (for complex element types)[bib].
â end note]
Effects:
For the overloads without an ExecutionPolicy argument, equivalent to:matrix_rank_1_update(x, conjugated(y), A);
otherwise, equivalent to:matrix_rank_1_update(std::forward(exec), x, conjugated(y), A);
29.9.14.7 Symmetric or Hermitian Rank-1 (outer product) update of a matrix [linalg.algs.blas2.symherrank1]
[Note 1:
These functions correspond to the BLAS functionsxSYR, xSPR, xHER, and xHPR[bib].
They have overloads taking a scaling factor alpha, because it would be impossible to express the updateA=AâxxT otherwise.
â end note]
The following elements apply to all functions in [linalg.algs.blas2.symherrank1].
Mandates:
If InOutMat has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
compatible-static-extents<decltype(A), decltype(A)>(0, 1) is true; and
compatible-static-extents<decltype(A), decltype(x)>(0, 0) is true.
Preconditions:
A.extent(0) equals A.extent(1), and
A.extent(0) equals x.extent(0).
Complexity: O(x.extent(0)Ãx.extent(0)).
template<class Scalar, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void symmetric_matrix_rank_1_update(Scalar alpha, InVec x, InOutMat A, Triangle t); template<class ExecutionPolicy, class Scalar, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void symmetric_matrix_rank_1_update(ExecutionPolicy&& exec, Scalar alpha, InVec x, InOutMat A, Triangle t);
These functions perform a symmetric rank-1 update of the symmetric matrix A, taking into account the Triangle parameter that applies to A ([linalg.general]).
Effects: Computes a matrix Aâ² such thatAâ²=A+αxxT, where the scalar α is alpha, and assigns each element of Aâ² to the corresponding element of A.
template<[in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void symmetric_matrix_rank_1_update(InVec x, InOutMat A, Triangle t); template<class ExecutionPolicy, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void symmetric_matrix_rank_1_update(ExecutionPolicy&& exec, InVec x, InOutMat A, Triangle t);
These functions perform a symmetric rank-1 update of the symmetric matrix A, taking into account the Triangle parameter that applies to A ([linalg.general]).
Effects: Computes a matrix Aâ² such that Aâ²=A+xxT and assigns each element of Aâ² to the corresponding element of A.
template<class Scalar, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void hermitian_matrix_rank_1_update(Scalar alpha, InVec x, InOutMat A, Triangle t); template<class ExecutionPolicy, class Scalar, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void hermitian_matrix_rank_1_update(ExecutionPolicy&& exec, Scalar alpha, InVec x, InOutMat A, Triangle t);
These functions perform a Hermitian rank-1 update of the Hermitian matrix A, taking into account the Triangle parameter that applies to A ([linalg.general]).
Effects: Computes Aâ² such thatAâ²=A+αxxH, where the scalar α is alpha, and assigns each element of Aâ² to the corresponding element of A.
template<[in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void hermitian_matrix_rank_1_update(InVec x, InOutMat A, Triangle t); template<class ExecutionPolicy, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void hermitian_matrix_rank_1_update(ExecutionPolicy&& exec, InVec x, InOutMat A, Triangle t);
These functions perform a Hermitian rank-1 update of the Hermitian matrix A, taking into account the Triangle parameter that applies to A ([linalg.general]).
Effects: Computes a matrix Aâ² such that Aâ²=A+xxH and assigns each element of Aâ² to the corresponding element of A.
29.9.14.8 Symmetric and Hermitian rank-2 matrix updates [linalg.algs.blas2.rank2]
[Note 1:
These functions correspond to the BLAS functionsxSYR2,xSPR2, xHER2 and xHPR2[bib].
â end note]
The following elements apply to all functions in [linalg.algs.blas2.rank2].
Mandates:
If InOutMat has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
compatible-static-extents<decltype(A), decltype(A)>(0, 1) is true; and
possibly-multipliable<decltype(A), decltype(x), decltype(y)>() is true.
Preconditions:
A.extent(0) equals A.extent(1), and
multipliable(A, x, y) is true.
Complexity: O(x.extent(0)Ãy.extent(0)).
template<[in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void symmetric_matrix_rank_2_update(InVec1 x, InVec2 y, InOutMat A, Triangle t); template<class ExecutionPolicy, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void symmetric_matrix_rank_2_update(ExecutionPolicy&& exec, InVec1 x, InVec2 y, InOutMat A, Triangle t);
These functions perform a symmetric rank-2 update of the symmetric matrix A, taking into account the Triangle parameter that applies to A ([linalg.general]).
Effects: Computes Aâ² such that Aâ²=A+xyT+yxT and assigns each element of Aâ² to the corresponding element of A.
template<[in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void hermitian_matrix_rank_2_update(InVec1 x, InVec2 y, InOutMat A, Triangle t); template<class ExecutionPolicy, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec1, [in-vector](#concept:in-vector "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InVec2, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void hermitian_matrix_rank_2_update(ExecutionPolicy&& exec, InVec1 x, InVec2 y, InOutMat A, Triangle t);
These functions perform a Hermitian rank-2 update of the Hermitian matrix A, taking into account the Triangle parameter that applies to A ([linalg.general]).
Effects: Computes Aâ² such that Aâ²=A+xyH+yxH and assigns each element of Aâ² to the corresponding element of A.
29.9.15 BLAS 3 algorithms [linalg.algs.blas3]
29.9.15.1 General matrix-matrix product [linalg.algs.blas3.gemm]
[Note 1:
These functions correspond to the BLAS function xGEMM[bib].
â end note]
The following elements apply to all functions in [linalg.algs.blas3.gemm] in addition to function-specific elements.
Mandates: possibly-multipliable<decltype(A), decltype(B), decltype(C)>() is true.
Preconditions: multipliable(A, B, C) is true.
Complexity: O(A.extent(0)ÃA.extent(1)ÃB.extent(1)).
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat1, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat2, [out-matrix](#concept:out-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutMat> void matrix_product(InMat1 A, InMat2 B, OutMat C); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat1, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat2, [out-matrix](#concept:out-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutMat> void matrix_product(ExecutionPolicy&& exec, InMat1 A, InMat2 B, OutMat C);
Effects: Computes C=AB.
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat1, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat2, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat3, [out-matrix](#concept:out-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutMat> void matrix_product(InMat1 A, InMat2 B, InMat3 E, OutMat C); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat1, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat2, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat3, [out-matrix](#concept:out-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutMat> void matrix_product(ExecutionPolicy&& exec, InMat1 A, InMat2 B, InMat3 E, OutMat C);
Mandates: possibly-addable<InMat3, InMat3, OutMat>() is true.
Preconditions: addable(E, E, C) is true.
Effects: Computes C=E+AB.
Remarks: C may alias E.
29.9.15.2 Symmetric, Hermitian, and triangular matrix-matrix product [linalg.algs.blas3.xxmm]
[Note 1:
These functions correspond to the BLAS functionsxSYMM, xHEMM, and xTRMM[bib].
â end note]
The following elements apply to all functions in [linalg.algs.blas3.xxmm] in addition to function-specific elements.
Mandates:
possibly-multipliable<decltype(A), decltype(B), decltype(C)>() is true, and
possibly-addable<decltype(E), decltype(E), decltype(C)>() is true for those overloads that take an E parameter.
Preconditions:
multipliable(A, B, C) is true, and
addable(E, E, C) is true for those overloads that take an E parameter.
Complexity: O(A.extent(0)ÃA.extent(1)ÃB.extent(1)).
`template<in-matrix InMat1, class Triangle, in-matrix InMat2, out-matrix OutMat> void symmetric_matrix_product(InMat1 A, Triangle t, InMat2 B, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, class Triangle, in-matrix InMat2, out-matrix OutMat> void symmetric_matrix_product(ExecutionPolicy&& exec, InMat1 A, Triangle t, InMat2 B, OutMat C);
template<in-matrix InMat1, class Triangle, in-matrix InMat2, out-matrix OutMat> void hermitian_matrix_product(InMat1 A, Triangle t, InMat2 B, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, class Triangle, in-matrix InMat2, out-matrix OutMat> void hermitian_matrix_product(ExecutionPolicy&& exec, InMat1 A, Triangle t, InMat2 B, OutMat C);
template<in-matrix InMat1, class Triangle, class DiagonalStorage, in-matrix InMat2, out-matrix OutMat> void triangular_matrix_product(InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, class Triangle, class DiagonalStorage, in-matrix InMat2, out-matrix OutMat> void triangular_matrix_product(ExecutionPolicy&& exec, InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat C); `
These functions perform a matrix-matrix multiply, taking into account the Triangle and DiagonalStorage (if applicable) parameters that apply to the symmetric, Hermitian, or triangular (respectively) matrix A ([linalg.general]).
Mandates:
If InMat1 has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument; and
compatible-static-extents<InMat1, InMat1>(0, 1) is true.
Preconditions: A.extent(0) == A.extent(1) is true.
Effects: Computes C=AB.
` template<in-matrix InMat1, in-matrix InMat2, class Triangle, out-matrix OutMat> void symmetric_matrix_product(InMat1 A, InMat2 B, Triangle t, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, in-matrix InMat2, class Triangle, out-matrix OutMat> void symmetric_matrix_product(ExecutionPolicy&& exec, InMat1 A, InMat2 B, Triangle t, OutMat C);
template<in-matrix InMat1, in-matrix InMat2, class Triangle, out-matrix OutMat> void hermitian_matrix_product(InMat1 A, InMat2 B, Triangle t, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, in-matrix InMat2, class Triangle, out-matrix OutMat> void hermitian_matrix_product(ExecutionPolicy&& exec, InMat1 A, InMat2 B, Triangle t, OutMat C);
template<in-matrix InMat1, in-matrix InMat2, class Triangle, class DiagonalStorage, out-matrix OutMat> void triangular_matrix_product(InMat1 A, InMat2 B, Triangle t, DiagonalStorage d, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, in-matrix InMat2, class Triangle, class DiagonalStorage, out-matrix OutMat> void triangular_matrix_product(ExecutionPolicy&& exec, InMat1 A, InMat2 B, Triangle t, DiagonalStorage d, OutMat C); `
These functions perform a matrix-matrix multiply, taking into account the Triangle and DiagonalStorage (if applicable) parameters that apply to the symmetric, Hermitian, or triangular (respectively) matrix B ([linalg.general]).
Mandates:
If InMat2 has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument; and
compatible-static-extents<InMat2, InMat2>(0, 1) is true.
Preconditions: B.extent(0) == B.extent(1) is true.
Effects: Computes C=AB.
`template<in-matrix InMat1, class Triangle, in-matrix InMat2, in-matrix InMat3, out-matrix OutMat> void symmetric_matrix_product(InMat1 A, Triangle t, InMat2 B, InMat3 E, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, class Triangle, in-matrix InMat2, in-matrix InMat3, out-matrix OutMat> void symmetric_matrix_product(ExecutionPolicy&& exec, InMat1 A, Triangle t, InMat2 B, InMat3 E, OutMat C);
template<in-matrix InMat1, class Triangle, in-matrix InMat2, in-matrix InMat3, out-matrix OutMat> void hermitian_matrix_product(InMat1 A, Triangle t, InMat2 B, InMat3 E, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, class Triangle, in-matrix InMat2, in-matrix InMat3, out-matrix OutMat> void hermitian_matrix_product(ExecutionPolicy&& exec, InMat1 A, Triangle t, InMat2 B, InMat3 E, OutMat C);
template<in-matrix InMat1, class Triangle, class DiagonalStorage, in-matrix InMat2, in-matrix InMat3, out-matrix OutMat> void triangular_matrix_product(InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, InMat3 E, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, class Triangle, class DiagonalStorage, in-matrix InMat2, in-matrix InMat3, out-matrix OutMat> void triangular_matrix_product(ExecutionPolicy&& exec, InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, InMat3 E, OutMat C); `
These functions perform a potentially overwriting matrix-matrix multiply-add, taking into account the Triangle and DiagonalStorage (if applicable) parameters that apply to the symmetric, Hermitian, or triangular (respectively) matrix A ([linalg.general]).
Mandates:
If InMat1 has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument; and
compatible-static-extents<InMat1, InMat1>(0, 1) is true.
Preconditions: A.extent(0) == A.extent(1) is true.
Effects: Computes C=E+AB.
Remarks: C may alias E.
`template<in-matrix InMat1, in-matrix InMat2, class Triangle, in-matrix InMat3, out-matrix OutMat> void symmetric_matrix_product(InMat1 A, InMat2 B, Triangle t, InMat3 E, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, in-matrix InMat2, class Triangle, in-matrix InMat3, out-matrix OutMat> void symmetric_matrix_product(ExecutionPolicy&& exec, InMat1 A, InMat2 B, Triangle t, InMat3 E, OutMat C);
template<in-matrix InMat1, in-matrix InMat2, class Triangle, in-matrix InMat3, out-matrix OutMat> void hermitian_matrix_product(InMat1 A, InMat2 B, Triangle t, InMat3 E, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, in-matrix InMat2, class Triangle, in-matrix InMat3, out-matrix OutMat> void hermitian_matrix_product(ExecutionPolicy&& exec, InMat1 A, InMat2 B, Triangle t, InMat3 E, OutMat C);
template<in-matrix InMat1, in-matrix InMat2, class Triangle, class DiagonalStorage, in-matrix InMat3, out-matrix OutMat> void triangular_matrix_product(InMat1 A, InMat2 B, Triangle t, DiagonalStorage d, InMat3 E, OutMat C); template<class ExecutionPolicy, in-matrix InMat1, in-matrix InMat2, class Triangle, class DiagonalStorage, in-matrix InMat3, out-matrix OutMat> void triangular_matrix_product(ExecutionPolicy&& exec, InMat1 A, InMat2 B, Triangle t, DiagonalStorage d, InMat3 E, OutMat C); `
These functions perform a potentially overwriting matrix-matrix multiply-add, taking into account the Triangle and DiagonalStorage (if applicable) parameters that apply to the symmetric, Hermitian, or triangular (respectively) matrix B ([linalg.general]).
Mandates:
If InMat2 has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument; and
compatible-static-extents<InMat2, InMat2>(0, 1) is true.
Preconditions: B.extent(0) == B.extent(1) is true.
Effects: Computes C=E+AB.
Remarks: C may alias E.
29.9.15.3 In-place triangular matrix-matrix product [linalg.algs.blas3.trmm]
These functions perform an in-place matrix-matrix multiply, taking into account the Triangle and DiagonalStorage parameters that apply to the triangular matrix A ([linalg.general]).
[Note 1:
These functions correspond to the BLAS function xTRMM[bib].
â end note]
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [inout-matrix](#concept:inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat> void triangular_matrix_left_product(InMat A, Triangle t, DiagonalStorage d, InOutMat C); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [inout-matrix](#concept:inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat> void triangular_matrix_left_product(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InOutMat C);
Mandates:
If InMat has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
possibly-multipliable<InMat, InOutMat, InOutMat>() is true; and
compatible-static-extents<InMat, InMat>(0, 1) is true.
Preconditions:
multipliable(A, C, C) is true, and
A.extent(0) == A.extent(1) is true.
Effects: Computes a matrix Câ² such that Câ²=AC and assigns each element of Câ² to the corresponding element of C.
Complexity: O(A.extent(0)ÃA.extent(1)ÃC.extent(0)).
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [inout-matrix](#concept:inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat> void triangular_matrix_right_product(InMat A, Triangle t, DiagonalStorage d, InOutMat C); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [inout-matrix](#concept:inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat> void triangular_matrix_right_product(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InOutMat C);
Mandates:
If InMat has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
possibly-multipliable<InOutMat, InMat, InOutMat>() is true; and
compatible-static-extents<InMat, InMat>(0, 1) is true.
Preconditions:
multipliable(C, A, C) is true, and
A.extent(0) == A.extent(1) is true.
Effects: Computes a matrix Câ² such that Câ²=CA and assigns each element of Câ² to the corresponding element of C.
Complexity: O(A.extent(0)ÃA.extent(1)ÃC.extent(0)).
29.9.15.4 Rank-k update of a symmetric or Hermitian matrix [linalg.algs.blas3.rankk]
[Note 1:
These functions correspond to the BLAS functionsxSYRK and xHERK[bib].
â end note]
The following elements apply to all functions in [linalg.algs.blas3.rankk].
Mandates:
If InOutMat has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
compatible-static-extents<decltype(A), decltype(A)>(0, 1) is true;
compatible-static-extents<decltype(C), decltype(C)>(0, 1) is true; and
compatible-static-extents<decltype(A), decltype(C)>(0, 0) is true.
Preconditions:
A.extent(0) equals A.extent(1),
C.extent(0) equals C.extent(1), and
A.extent(0) equals C.extent(0).
Complexity: O(A.extent(0)ÃA.extent(1)ÃC.extent(0)).
template<class Scalar, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void symmetric_matrix_rank_k_update(Scalar alpha, InMat A, InOutMat C, Triangle t); template<class ExecutionPolicy, class Scalar, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void symmetric_matrix_rank_k_update(ExecutionPolicy&& exec, Scalar alpha, InMat A, InOutMat C, Triangle t);
Effects: Computes a matrix Câ² such that Câ²=C+αAAT, where the scalar α is alpha, and assigns each element of Câ² to the corresponding element of C.
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void symmetric_matrix_rank_k_update(InMat A, InOutMat C, Triangle t); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void symmetric_matrix_rank_k_update(ExecutionPolicy&& exec, InMat A, InOutMat C, Triangle t);
Effects: Computes a matrix Câ² such that Câ²=C+AAT, and assigns each element of Câ² to the corresponding element of C.
template<class Scalar, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void hermitian_matrix_rank_k_update(Scalar alpha, InMat A, InOutMat C, Triangle t); template<class ExecutionPolicy, class Scalar, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void hermitian_matrix_rank_k_update(ExecutionPolicy&& exec, Scalar alpha, InMat A, InOutMat C, Triangle t);
Effects: Computes a matrix Câ² such that Câ²=C+αAAH, where the scalar α is alpha, and assigns each element of Câ² to the corresponding element of C.
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void hermitian_matrix_rank_k_update(InMat A, InOutMat C, Triangle t); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void hermitian_matrix_rank_k_update(ExecutionPolicy&& exec, InMat A, InOutMat C, Triangle t);
Effects: Computes a matrix Câ² such that Câ²=C+AAH, and assigns each element of Câ² to the corresponding element of C.
29.9.15.5 Rank-2k update of a symmetric or Hermitian matrix [linalg.algs.blas3.rank2k]
[Note 1:
These functions correspond to the BLAS functionsxSYR2K and xHER2K[bib].
â end note]
The following elements apply to all functions in [linalg.algs.blas3.rank2k].
Mandates:
If InOutMat has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
possibly-addable<decltype(A), decltype(B), decltype(C)>() is true; and
compatible-static-extents<decltype(A), decltype(A)>(0, 1) is true.
Preconditions:
addable(A, B, C) is true, and
A.extent(0) equals A.extent(1).
Complexity: O(A.extent(0)ÃA.extent(1)ÃC.extent(0)).
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat1, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat2, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void symmetric_matrix_rank_2k_update(InMat1 A, InMat2 B, InOutMat C, Triangle t); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat1, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat2, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void symmetric_matrix_rank_2k_update(ExecutionPolicy&& exec, InMat1 A, InMat2 B, InOutMat C, Triangle t);
Effects: Computes a matrix Câ² such that Câ²=C+ABT+BAT, and assigns each element of Câ² to the corresponding element of C.
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat1, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat2, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void hermitian_matrix_rank_2k_update(InMat1 A, InMat2 B, InOutMat C, Triangle t); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat1, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat2, [possibly-packed-inout-matrix](#concept:possibly-packed-inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class Triangle> void hermitian_matrix_rank_2k_update(ExecutionPolicy&& exec, InMat1 A, InMat2 B, InOutMat C, Triangle t);
Effects: Computes a matrix Câ² such that Câ²=C+ABH+BAH, and assigns each element of Câ² to the corresponding element of C.
29.9.15.6 Solve multiple triangular linear systems [linalg.algs.blas3.trsm]
[Note 1:
These functions correspond to the BLAS function xTRSM[bib].
â end note]
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat1, class Triangle, class DiagonalStorage, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat2, [out-matrix](#concept:out-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutMat, class BinaryDivideOp> void triangular_matrix_matrix_left_solve(InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat X, BinaryDivideOp divide); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat1, class Triangle, class DiagonalStorage, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat2, [out-matrix](#concept:out-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutMat, class BinaryDivideOp> void triangular_matrix_matrix_left_solve(ExecutionPolicy&& exec, InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat X, BinaryDivideOp divide);
These functions perform multiple matrix solves, taking into account the Triangle and DiagonalStorage parameters that apply to the triangular matrix A ([linalg.general]).
Mandates:
If InMat1 has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
possibly-multipliable<InMat1, OutMat, InMat2>() is true; and
compatible-static-extents<InMat1, InMat1>(0, 1) is true.
Preconditions:
multipliable(A, X, B) is true, and
A.extent(0) == A.extent(1) is true.
Effects: Computes Xâ² such that AXâ²=B, and assigns each element of Xâ² to the corresponding element of X.
If no such Xâ² exists, then the elements of X are valid but unspecified.
Complexity: O(A.extent(0)ÃX.extent(1)ÃX.extent(1)).
[Note 2:
Since the triangular matrix is on the left, the desired divide implementation in the case of noncommutative multiplication is mathematically equivalent to yâ1x, where x is the first argument and y is the second argument, and yâ1 denotes the multiplicative inverse of y.
â end note]
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat1, class Triangle, class DiagonalStorage, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat2, [out-matrix](#concept:out-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutMat> void triangular_matrix_matrix_left_solve(InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat X);
Effects: Equivalent to:triangular_matrix_matrix_left_solve(A, t, d, B, X, divides{});
template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat1, class Triangle, class DiagonalStorage, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat2, [out-matrix](#concept:out-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutMat> void triangular_matrix_matrix_left_solve(ExecutionPolicy&& exec, InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat X);
Effects: Equivalent to:triangular_matrix_matrix_left_solve(std::forward(exec), A, t, d, B, X, divides{});
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat1, class Triangle, class DiagonalStorage, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat2, [out-matrix](#concept:out-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutMat, class BinaryDivideOp> void triangular_matrix_matrix_right_solve(InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat X, BinaryDivideOp divide); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat1, class Triangle, class DiagonalStorage, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat2, [out-matrix](#concept:out-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutMat, class BinaryDivideOp> void triangular_matrix_matrix_right_solve(ExecutionPolicy&& exec, InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat X, BinaryDivideOp divide);
These functions perform multiple matrix solves, taking into account the Triangle and DiagonalStorage parameters that apply to the triangular matrix A ([linalg.general]).
Mandates:
If InMat1 has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
possibly-multipliable<OutMat, InMat1, InMat2>() is true; and
compatible-static-extents<InMat1, InMat1>(0,1) is true.
Preconditions:
multipliable(X, A, B) is true, and
A.extent(0) == A.extent(1) is true.
Effects: Computes Xâ² such that Xâ²A=B, and assigns each element of Xâ² to the corresponding element of X.
If no such Xâ² exists, then the elements of X are valid but unspecified.
Complexity: O( B.extent(0) â B.extent(1) â A.extent(1) )
[Note 3:
Since the triangular matrix is on the right, the desired divide implementation in the case of noncommutative multiplication is mathematically equivalent to xyâ1, where x is the first argument and y is the second argument, and yâ1 denotes the multiplicative inverse of y.
â end note]
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat1, class Triangle, class DiagonalStorage, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat2, [out-matrix](#concept:out-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutMat> void triangular_matrix_matrix_right_solve(InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat X);
Effects: Equivalent to:triangular_matrix_matrix_right_solve(A, t, d, B, X, divides{});
template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat1, class Triangle, class DiagonalStorage, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat2, [out-matrix](#concept:out-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") OutMat> void triangular_matrix_matrix_right_solve(ExecutionPolicy&& exec, InMat1 A, Triangle t, DiagonalStorage d, InMat2 B, OutMat X);
Effects: Equivalent to:triangular_matrix_matrix_right_solve(std::forward(exec), A, t, d, B, X, divides{});
29.9.15.7 Solve multiple triangular linear systems in-place [linalg.algs.blas3.inplacetrsm]
[Note 1:
These functions correspond to the BLAS function xTRSM[bib].
â end note]
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [inout-matrix](#concept:inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class BinaryDivideOp> void triangular_matrix_matrix_left_solve(InMat A, Triangle t, DiagonalStorage d, InOutMat B, BinaryDivideOp divide); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [inout-matrix](#concept:inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class BinaryDivideOp> void triangular_matrix_matrix_left_solve(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InOutMat B, BinaryDivideOp divide);
These functions perform multiple in-place matrix solves, taking into account the Triangle and DiagonalStorage parameters that apply to the triangular matrix A ([linalg.general]).
[Note 2:
This algorithm makes it possible to compute factorizations like Cholesky and LU in place.
Performing triangular solve in place hinders parallelization.
However, other ExecutionPolicy specific optimizations, such as vectorization, are still possible.
â end note]
Mandates:
If InMat has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
possibly-multipliable<InMat, InOutMat, InOutMat>() is true; and
compatible-static-extents<InMat, InMat>(0, 1) is true.
Preconditions:
multipliable(A, B, B) is true, and
A.extent(0) == A.extent(1) is true.
Effects: Computes Xâ² such that AXâ²=B, and assigns each element of Xâ² to the corresponding element of B.
If so such Xâ² exists, then the elements of B are valid but unspecified.
Complexity: O(A.extent(0)ÃA.extent(1)ÃB.extent(1)).
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [inout-matrix](#concept:inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat> void triangular_matrix_matrix_left_solve(InMat A, Triangle t, DiagonalStorage d, InOutMat B);
Effects: Equivalent to:triangular_matrix_matrix_left_solve(A, t, d, B, divides{});
template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [inout-matrix](#concept:inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat> void triangular_matrix_matrix_left_solve(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InOutMat B);
Effects: Equivalent to:triangular_matrix_matrix_left_solve(std::forward(exec), A, t, d, B, divides{});
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [inout-matrix](#concept:inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class BinaryDivideOp> void triangular_matrix_matrix_right_solve(InMat A, Triangle t, DiagonalStorage d, InOutMat B, BinaryDivideOp divide); template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [inout-matrix](#concept:inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat, class BinaryDivideOp> void triangular_matrix_matrix_right_solve(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InOutMat B, BinaryDivideOp divide);
These functions perform multiple in-place matrix solves, taking into account the Triangle and DiagonalStorage parameters that apply to the triangular matrix A ([linalg.general]).
[Note 3:
This algorithm makes it possible to compute factorizations like Cholesky and LU in place.
Performing triangular solve in place hinders parallelization.
However, other ExecutionPolicy specific optimizations, such as vectorization, are still possible.
â end note]
Mandates:
If InMat has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
possibly-multipliable<InOutMat, InMat, InOutMat>() is true; and
compatible-static-extents<InMat, InMat>(0, 1) is true.
Preconditions:
multipliable(B, A, B) is true, and
A.extent(0) == A.extent(1) is true.
Effects: Computes Xâ² such that Xâ²A=B, and assigns each element of Xâ² to the corresponding element of B.
If so such Xâ² exists, then the elements of B are valid but unspecified.
Complexity: O(A.extent(0)ÃA.extent(1)ÃB.extent(1)).
template<[in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [inout-matrix](#concept:inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat> void triangular_matrix_matrix_right_solve(InMat A, Triangle t, DiagonalStorage d, InOutMat B);
Effects: Equivalent to:triangular_matrix_matrix_right_solve(A, t, d, B, divides{});
template<class ExecutionPolicy, [in-matrix](#concept:in-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InMat, class Triangle, class DiagonalStorage, [inout-matrix](#concept:inout-matrix "29.9.7.5 Argument concepts [linalg.helpers.concepts]") InOutMat> void triangular_matrix_matrix_right_solve(ExecutionPolicy&& exec, InMat A, Triangle t, DiagonalStorage d, InOutMat B);
Effects: Equivalent to:triangular_matrix_matrix_right_solve(std::forward(exec), A, t, d, B, divides{});