4. matrix_rank_1_update_c

最終更新日時:
が更新

履歴 編集

function template
<linalg>

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

namespace std::linalg {
  template<in-vector InVec1,
           in-vector InVec2,
           out-matrix OutMat>
  void matrix_rank_1_update_c(
    InVec1 x,
    InVec2 y,
    OutMat A); // (1)

  template<class ExecutionPolicy,
           in-vector InVec1,
           in-vector InVec2,
           out-matrix OutMat>
  void matrix_rank_1_update_c(
    ExecutionPolicy&& exec,
    InVec1 x,
    InVec2 y,
    OutMat A); // (2)

  template<in-vector InVec1,
           in-vector InVec2,
           in-matrix InMat,
           out-matrix OutMat>
  void matrix_rank_1_update_c(
    InVec1 x,
    InVec2 y,
    InMat E,
    OutMat A); // (3)

  template<class ExecutionPolicy,
           in-vector InVec1,
           in-vector InVec2,
           in-matrix InMat,
           out-matrix OutMat>
  void matrix_rank_1_update_c(
    ExecutionPolicy&& exec,
    InVec1 x,
    InVec2 y,
    InMat E,
    OutMat A); // (4)
}

概要

非対称で共役を取るrank-1 updateを行う。

  • (1), (2): 上書き版。$A = xy^*$ を計算してAに書き込む。
  • (3), (4): 更新版。入力行列Eに更新項を加えた $A = E + xy^*$ を計算してAに書き込む。
  • (2), (4): それぞれ(1), (3)を指定された実行ポリシーで実行する。

適格要件

事前条件

効果

  • (1): matrix_rank_1_update(x, conjugated(y), A)と同じ。
  • (2): matrix_rank_1_update(std::forward<ExecutionPolicy>(exec), x, conjugated(y), A)と同じ。
  • (3): matrix_rank_1_update(x, conjugated(y), E, A)と同じ。
  • (4): matrix_rank_1_update(std::forward<ExecutionPolicy>(exec), x, conjugated(y), E, A)と同じ。

戻り値

なし

計算量

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

備考

  • 上書き版 (1), (2) は、出力行列Aの元の値を参照せずに更新項で上書きする。
  • 加算更新 $A = A + xy^*$ を行いたい場合は、更新版 (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 < 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 = 0; j < A.extent(1); ++j) {
      A[i,j] = A.extent(1) * i + j;
    }
  }
}

int main()
{
  constexpr size_t N = 4;

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

  std::mdspan<
    double,
    std::extents<size_t, N, N>> A(A_vec.data());
  std::mdspan<
    double,
    std::extents<size_t, N, N>> E(E_vec.data());
  std::mdspan x(x_vec.data(), N);
  std::mdspan y(y_vec.data(), N);

  init_vec(x);
  init_vec(y);

  // (1) 上書き版 : A = x conj(y)^T
  std::cout << "(1)\n";
  std::linalg::matrix_rank_1_update_c(
    x,
    y,
    A);
  print_mat(A);

  // (3) 更新版 : A = E + x conj(y)^T
  init_mat(E);
  std::cout << "(3)\n";
  std::linalg::matrix_rank_1_update_c(
    x,
    y,
    E,
    A);
  print_mat(A);

  return 0;
}

この例はxとyが実数のため、共役を取っても値は変化せず、結果はmatrix_rank_1_updateと同じになる。

出力

(1)
0 0 0 0
0 1 2 3
0 2 4 6
0 3 6 9
(3)
0 1 2 3
4 6 8 10
8 11 14 17
12 16 20 24

バージョン

言語

  • C++26

処理系

関連項目

参照