4. triangular_matrix_vector_product

最終更新日時(UTC):
が更新

履歴 編集

function template
<linalg>

std::linalg::triangular_matrix_vector_product(C++26)

namespace std::linalg {
  template<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); // (1)

  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); // (2)

  template<in-matrix InMat,
           class Triangle,
           class DiagonalStorage,
           inout-vector InOutVec>
  void triangular_matrix_vector_product(
    InMat A,
    Triangle t,
    DiagonalStorage d,
    InOutVec y); // (3)

  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); // (4)

  template<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); // (5)

  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); // (6)
}

概要

三角行列とベクトルの積を計算する。 引数tは対称行列の成分が上三角にあるのか、それとも下三角にあるのかを示す。 引数dは対称行列の対角成分を暗黙に乗法における単位元とみなすかどうかを指定する。

  • (1): $y \leftarrow Ax$
  • (2): (1)を指定された実行ポリシーで実行する。
  • (3): $y \leftarrow Ay$
  • (4): (3)を指定された実行ポリシーで実行する。
  • (5): $z \leftarrow y + Ax$
  • (6): (5)を指定された実行ポリシーで実行する。

適格要件

事前条件

  • 共通:
  • A.extent(0) == A.extent(1) (つまりAが正方行列であること)
  • A.extent(0) == y.extent(0) (つまりyに結果を代入できること)
  • (1), (2), (5), (6): A.extent(0) == x.extent(0) (つまり$Ax$を計算できること)
  • (5), (6): A.extent(0) == z.extent(0) (つまりzに結果を代入できること)

効果

対称行列の成分の位置を示すtと対角成分へアクセスするかどうかを示すdを考慮した、対称行列とベクトルの積を計算する。

  • (1), (2): $y \leftarrow Ax$
  • (3), (4): $y \leftarrow Ay$
  • (5), (6): $z \leftarrow y + Ax$

戻り値

なし

計算量

$O(\verb|A.extent(1)|\times \verb|x.extent(0)|)$

備考

  • (4): in-placeアルゴリズムなので並列実行を妨げるが、ベクトル化といったExecutionPolicy特有の最適化はできる。
  • (5), (6): zにyを入れても良い。

例

[注意] 処理系にあるコンパイラで確認していないため、間違っているかもしれません。

#include <array>
#include <iostream>
#include <linalg>
#include <mdspan>
#include <vector>

template <class Vector>
void print(const Vector& v, const std::string& name) {
  for (int i = 0; i < v.extent(0); ++i) {
    std::cout << name << "[" << i << "]" << " = " << v[i] << '\n';
  }
}

template <class Vector>
void init(Vector& v) {
  for (int i = 0; i < v.extent(0); ++i) {
    v[i] = i;
  }
}

int main()
{
  constexpr size_t N = 4;
  constexpr size_t M = 4;

  std::vector<double> A_vec(N*M);
  std::vector<double> x_vec(M);
  std::array<double, N> y_vec, z_vec;

  std::mdspan<
    double,
    std::extents<size_t, N, M>,
    std::linalg::layout_blas_packed<
      std::linalg::upper_triangle_t,
      std::linalg::row_major_t>
  > A(A_vec.data());
  std::mdspan x(x_vec.data(), M);
  std::mdspan y(y_vec.data(), N);
  std::mdspan z(z_vec.data(), N);

  for(int i = 0; i < A.extent(0); ++i) {
    for(int j = i + 1; j < A.extent(1); ++j) {
      A[i,j] = A.extent(1) * i + j;
    }
  }

  init(x);

  // (1)
  std::cout << "(1)\n";
  std::linalg::triangular_matrix_vector_product(
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    x,
    y);
  print(y, "y");

  // (2)
  std::cout << "(2)\n";
  std::linalg::triangular_matrix_vector_product(
    std::execution::par,
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    x,
    y);
  print(y, "y");

  init(y);

  // (3)
  std::cout << "(3)\n";
  std::linalg::triangular_matrix_vector_product(
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    y);
  print(y, "y");

  init(y);

  // (4)
  std::cout << "(4)\n";
  std::linalg::triangular_matrix_vector_product(
    std::execution::par,
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    y);
  print(y, "y");

  init(y);

  // (5)
  std::cout << "(5)\n";
  std::linalg::triangular_matrix_vector_product(
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    x,
    y,
    z);
  print(z, "z");

  // (6)
  std::cout << "(6)\n";
  std::linalg::triangular_matrix_vector_product(
    std::execution::par,
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    x,
    y,
    z);
  print(z, "z");

  return 0;
}

出力

(1)
y[0] = 14
y[1] = 34
y[2] = 35
y[3] = 3
(2)
y[0] = 14
y[1] = 34
y[2] = 35
y[3] = 3
(3)
y[0] = 14
y[1] = 34
y[2] = 35
y[3] = 3
(4)
y[0] = 14
y[1] = 34
y[2] = 35
y[3] = 3
(5)
z[0] = 14
z[1] = 35
z[2] = 37
z[3] = 6
(6)
z[0] = 14
z[1] = 35
z[2] = 37
z[3] = 6

バージョン

言語

  • C++26

処理系

関連項目

参照