4. hermitian_matrix_rank_2_update

最終更新日時:
が更新

履歴 編集

function template
<linalg>

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

namespace std::linalg {
  template<in-vector InVec1,
           in-vector InVec2,
           possibly-packed-out-matrix OutMat,
           class Triangle>
  void hermitian_matrix_rank_2_update(
    InVec1 x,
    InVec2 y,
    OutMat A,
    Triangle t); // (1)

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

  template<in-vector InVec1,
           in-vector InVec2,
           in-matrix InMat,
           possibly-packed-out-matrix OutMat,
           class Triangle>
  void hermitian_matrix_rank_2_update(
    InVec1 x,
    InVec2 y,
    InMat E,
    OutMat A,
    Triangle t); // (3)

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

概要

エルミートな(対称かつ共役を取る)rank-2 updateをエルミート行列に行う。 引数tはエルミート行列の成分が上三角にあるのか、それとも下三角にあるのかを示す。

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

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

適格要件

事前条件

効果

  • (1), (2): $A = xy^* + yx^*$
  • (3), (4): $A = E + xy^* + yx^*$

戻り値

なし

計算量

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

備考

加算更新$A \leftarrow A + xy^ + yx^$を行いたい場合は、更新前の行列をEとして渡すupdating版(3), (4)を使用する。

  • エルミート行列Aの対角成分については、real-if-neededにより実部のみが使用される。対角成分が非ゼロの虚部を持っていても、その虚部は無視される。

例

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

#include <array>
#include <complex>
#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 Matrix>
void init_mat(Matrix& A) {
  for(int i = 0; i < A.extent(0); ++i) {
    A[i, i] = std::complex<double>(i, 0);
    for(int j = i + 1; j < A.extent(1); ++j) {
      A[i, j] = std::complex<double>(i, j);
    }
  }
}

int main()
{
  constexpr size_t N = 4;

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

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

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

  for (int i = 0; i < x.extent(0); ++i) {
    x[i] = std::complex<double>(i, 0);
  }
  for (int i = 0; i < y.extent(0); ++i) {
    y[i] = std::complex<double>(i, 1);
  }

  // (1) 上書き版： A = xy^* + yx^*
  std::cout << "(1)\n";
  std::linalg::hermitian_matrix_rank_2_update(
    x,
    y,
    A,
    std::linalg::upper_triangle);
  print_mat(A);

  // (3) 更新版： A = E + xy^* + yx^*
  init_mat(E);
  std::cout << "(3)\n";
  std::linalg::hermitian_matrix_rank_2_update(
    x,
    y,
    E,
    A,
    std::linalg::upper_triangle);
  print_mat(A);

  return 0;
}

出力

(1)
(0,0) (0,1) (0,2) (0,3)
(0,1) (2,0) (4,1) (6,2)
(0,2) (4,1) (8,0) (12,1)
(0,3) (6,2) (12,1) (18,0)
(3)
(0,0) (0,2) (0,4) (0,6)
(0,2) (3,0) (5,3) (7,5)
(0,4) (5,3) (10,0) (14,4)
(0,6) (7,5) (14,4) (21,0)

バージョン

言語

  • C++26

処理系

関連項目

参照