4. symmetric_matrix_rank_1_update

最終更新日時:
が更新

履歴 編集

function template
<linalg>

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

namespace std::linalg {
  template<scalar Scalar,
           in-vector InVec,
           possibly-packed-out-matrix OutMat,
           class Triangle>
  void symmetric_matrix_rank_1_update(
    Scalar alpha,
    InVec x,
    OutMat A,
    Triangle t); // (1)

  template<class ExecutionPolicy,
           scalar Scalar,
           in-vector InVec,
           possibly-packed-out-matrix OutMat,
           class Triangle>
  void symmetric_matrix_rank_1_update(
    ExecutionPolicy&& exec,
    Scalar alpha,
    InVec x,
    OutMat A,
    Triangle t); // (2)

  template<scalar Scalar,
           in-vector InVec,
           in-matrix InMat,
           possibly-packed-out-matrix OutMat,
           class Triangle>
  void symmetric_matrix_rank_1_update(
    Scalar alpha,
    InVec x,
    InMat E,
    OutMat A,
    Triangle t); // (3)

  template<class ExecutionPolicy,
           scalar Scalar,
           in-vector InVec,
           in-matrix InMat,
           possibly-packed-out-matrix OutMat,
           class Triangle>
  void symmetric_matrix_rank_1_update(
    ExecutionPolicy&& exec,
    Scalar alpha,
    InVec x,
    InMat E,
    OutMat A,
    Triangle t); // (4)
}

概要

対称かつ共役を取らないrank-1 updateを対称行列に行う。 引数tは対称行列の成分が上三角にあるのか、それとも下三角にあるのかを示す。

(1), (2)は結果をAに上書きするoverwriting版、(3), (4)は入力行列Eに更新項を加えてAに書き込むupdating版である。

  • (1): $A = \alpha xx^T$
  • (2): (1)を指定された実行ポリシーで実行する。
  • (3): $A = E + \alpha xx^T$
  • (4): (3)を指定された実行ポリシーで実行する。

適格要件

事前条件

  • A.extent(0) == A.extent(1)
  • A.extent(0) == x.extent(0)
  • (3), (4): addable(A, E, A) == true

効果

  • (1), (2): $A = \alpha xx^T$
  • (3), (4): $A = E + \alpha xx^T$

戻り値

なし

計算量

$O((\verb|x.extent(0)|)^2)$

備考

  • overwriting版(1), (2)はAに結果を上書きする。加算更新$A = A + \alpha xx^T$を行いたい場合は、updating版(3), (4)に更新前の行列をEとして渡す。

例

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

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

template <class Matrix>
void print_mat(const Matrix& A) {
  for(int i = 0; i < A.extent(0); ++i) {
    for(int j = 0; j < i; ++j) {
      std::cout << A[j, i] << ' ';
    }
    for(int j = i; j < A.extent(1) - 1; ++j) {
      std::cout << A[i, j] << ' ';
    }
    std::cout << A[i, A.extent(1) - 1] << '\n';
  }
}

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

template <class Matrix>
void init_mat(Matrix& A) {
  for(int i = 0; i < A.extent(0); ++i) {
    for(int j = i; j < A.extent(1); ++j) {
      A[i,j] = A.extent(1) * i + j;
    }
  }
}

int main()
{
  constexpr size_t N = 4;

  using PackedMatrix = std::mdspan<
    double,
    std::extents<size_t, N, N>,
    std::linalg::layout_blas_packed<
      std::linalg::upper_triangle_t,
      std::linalg::row_major_t>>;

  std::vector<double> A_vec(N * N);
  std::vector<double> E_vec(N * N);
  std::vector<double> x_vec(N);

  PackedMatrix A(A_vec.data());
  PackedMatrix E(E_vec.data());
  std::mdspan  x(x_vec.data(), N);

  init_vec(x);

  // (1) overwriting: A = alpha xx^T
  std::cout << "overwriting (1)\n";
  std::linalg::symmetric_matrix_rank_1_update(
    2.0,
    x,
    A,
    std::linalg::upper_triangle);
  print_mat(A);

  // (3) updating: A = E + alpha xx^T
  init_mat(E);
  std::cout << "updating (3)\n";
  std::linalg::symmetric_matrix_rank_1_update(
    2.0,
    x,
    E,
    A,
    std::linalg::upper_triangle);
  print_mat(A);

  return 0;
}

出力

overwriting (1)
0 0 0 0
0 2 4 6
0 4 8 12
0 6 12 18
updating (3)
0 1 2 3
1 7 10 13
2 10 18 23
3 13 23 33

バージョン

言語

  • C++26

処理系

関連項目

参照