Class / Function / Type

std::
Header file

<>
Other / All

最終更新日時(UTC): 2024年12月05日 02時54分34秒
Akira Takahashi が更新

function template

std::linalg::triangular_matrix_vector_solve

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

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

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

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

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

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

  template<in-matrix InMat,
           class Triangle,
           class DiagonalStorage,
           inout-vector InOutVec>
  void triangular_matrix_vector_solve(
    InMat A,
    Triangle t,
    DiagonalStorage d,
    InOutVec b); // (7)

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

概要

三角行列に対して、連立一次方程式を解く。引数tは対称行列の成分が上三角にあるのか、それとも下三角にあるのかを示す。引数dには対称行列の対角成分を暗黙に乗法における単位元とみなすかどうかを指定する。引数divideには値の割り算を指定する。この引数は非可換な掛け算を持つ値型をサポートするためにある。

(1): 連立一次方程式 $A y = b$ を解き、yをxに代入する。もし解が存在しないなら、xは有効だが未規定。
(2): (1)を指定された実行ポリシーで実行する。
(3): 割り算にstd::divides<void>を用いて、(1)を行う。
(4): (3)を指定された実行ポリシーで実行する。
(5): xにbを使って、in-placeに(1)を行う。
(6): (5)を指定された実行ポリシーで実行する。
(7): 割り算にstd::divides<void>を用いて、(5)を行う。
(8): (7)を指定された実行ポリシーで実行する。

適格要件

共通:
- Triangleはupper_triangle_tまたはlower_triangle_t
- DiagonalStorageはimplicit_unit_diagonal_tまたはexplicit_diagonal_t
- InMatがlayout_blas_packedを持つなら、レイアウトのTriangleテンプレート引数とこの関数のTriangleテンプレート引数が同じ型
- compatible-static-extents<decltype(A), decltype(A)>(0, 1)がtrue (つまりAが正方行列であること)
- compatible-static-extents<decltype(A), decltype(b)>(0, 0)がtrue (つまりAの次元とbの次元が同じであること)
(1), (2), (3), (4): compatible-static-extents<decltype(A), decltype(x)>(0, 0)がtrue (つまりAの次元とbの次元が同じであること)
(2), (4), (6), (8): is_execution_policy<ExecutionPolicy>::valueがtrue

事前条件

共通:
- A.extent(0) == A.extent(1) (つまりAが正方行列であること)
- A.extent(0) == b.extent(0) (つまりAの次元とbの次元が同じであること)
(1), (2), (3), (4): A.extent(0) == x.extent(0) (つまりAの次元とxの次元が同じであること)

効果

対称行列の成分の位置を示すtと対角成分へアクセスするかどうかを示すdを考慮して、連立一次方程式の解を求める。

(1), (2): 連立一次方程式 $A y = b$ を解き、yをxに代入する。もし解が存在しないなら、xは有効だが未規定。
(3): triangular_matrix_vector_solve(A, t, d, b, x, divides<void>{})と同じ。
(4): triangular_matrix_vector_solve(std::forward<ExecutionPolicy>(exec), A, t, d, b, x, divides<void>{})と同じ。
(5), (6): xにbを使って、in-placeに(1)を行う。
(7): triangular_matrix_vector_solve(A, t, d, b, divides<void>{})と同じ。
(8): triangular_matrix_vector_solve(std::forward<ExecutionPolicy>(exec), A, t, d, b, divides<void>{})と同じ。

備考

(6), (8): in-placeアルゴリズムなので並列実行を妨げるが、ベクトル化といったExecutionPolicy特有の最適化はできる。

例

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

#include <array>
#include <functional>
#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;

  std::vector<double> A_vec(N * N);
  std::vector<double> x_vec(N);
  std::array<double, N> b_vec;

  std::mdspan<
    double,
    std::extents<size_t, N, N>,
    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(), N);
  std::mdspan b(b_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(b);

  // (1)
  std::cout << "(1)\n";
  std::linalg::triangular_matrix_vector_solve(
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    b,
    x,
    std::divides<void>{});
  print(x, "x");

  // (2)
  std::cout << "(2)\n";
  std::linalg::triangular_matrix_vector_solve(
    std::execution::par,
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    b,
    x,
    std::divides<void>{});
  print(x, "x");

  // (3)
  std::cout << "(3)\n";
  std::linalg::triangular_matrix_vector_solve(
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    b,
    x);
  print(x, "x");

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

  // (5)
  std::cout << "(5)\n";
  std::linalg::triangular_matrix_vector_solve(
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    b,
    std::divides<void>{});
  print(b, "b");

  init(b);

  // (6)
  std::cout << "(6)\n";
  std::linalg::triangular_matrix_vector_solve(
    std::execution::par,
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    b,
    std::divides<void>{});
  print(b, "b");

  init(b);

  // (7)
  std::cout << "(7)\n";
  std::linalg::triangular_matrix_vector_solve(
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    b);
  print(b, "b");

  init(b);

  // (8)
  std::cout << "(8)\n";
  std::linalg::triangular_matrix_vector_solve(
    std::execution::par,
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    b);
  print(b, "b");

  return 0;
}

#include <array>
#include <functional>
#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;

std::vector<double> A_vec(N * N);
  std::vector<double> x_vec(N);
  std::array<double, N> b_vec;

std::mdspan<
    double,
    std::extents<size_t, N, N>,
    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(), N);
  std::mdspan b(b_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(b);

// (1)
  std::cout << "(1)\n";
  std::linalg::triangular_matrix_vector_solve(
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    b,
    x,
    std::divides<void>{});
  print(x, "x");

// (2)
  std::cout << "(2)\n";
  std::linalg::triangular_matrix_vector_solve(
    std::execution::par,
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    b,
    x,
    std::divides<void>{});
  print(x, "x");

// (3)
  std::cout << "(3)\n";
  std::linalg::triangular_matrix_vector_solve(
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    b,
    x);
  print(x, "x");

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

// (5)
  std::cout << "(5)\n";
  std::linalg::triangular_matrix_vector_solve(
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    b,
    std::divides<void>{});
  print(b, "b");

init(b);

// (6)
  std::cout << "(6)\n";
  std::linalg::triangular_matrix_vector_solve(
    std::execution::par,
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    b,
    std::divides<void>{});
  print(b, "b");

init(b);

// (7)
  std::cout << "(7)\n";
  std::linalg::triangular_matrix_vector_solve(
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    b);
  print(b, "b");

init(b);

// (8)
  std::cout << "(8)\n";
  std::linalg::triangular_matrix_vector_solve(
    std::execution::par,
    A,
    std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal,
    b);
  print(b, "b");

return 0;
}

 
#include <array>
#include <functional>
#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;

出力

(1)
x[0] = -3
x[1] = -4
x[2] = -1
x[3] = 3
(2)
x[0] = -3
x[1] = -4
x[2] = -1
x[3] = 3
(3)
x[0] = -3
x[1] = -4
x[2] = -1
x[3] = 3
(4)
x[0] = -3
x[1] = -4
x[2] = -1
x[3] = 3
(5)
b[0] = -3
b[1] = -4
b[2] = -1
b[3] = 3
(6)
b[0] = -3
b[1] = -4
b[2] = -1
b[3] = 3
(7)
b[0] = -3
b[1] = -4
b[2] = -1
b[3] = 3
(8)
b[0] = -3
b[1] = -4
b[2] = -1
b[3] = 3

Class / Function / Type

Header file

Other / All

std::linalg::triangular_matrix_vector_solve

概要

適格要件

事前条件

効果

戻り値

計算量

備考

例

出力

バージョン

言語

処理系

関連項目

参照