4. hermitian_matrix_rank_1_update

最終更新日時:
が更新

履歴 編集

function template
<linalg>

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

namespace std::linalg {
  template<scalar Scalar,
           in-vector InVec,
           possibly-packed-out-matrix OutMat,
           class Triangle>
  void hermitian_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 hermitian_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 hermitian_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 hermitian_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^*$
  • (2): (1)を指定された実行ポリシーで実行する。
  • (3): $A = E + \alpha xx^*$
  • (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^*$
  • (3), (4): $A = E + \alpha xx^*$

戻り値

なし

計算量

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

備考

  • overwriting版(1), (2)はAに結果を上書きする。加算更新$A = A + \alpha xx^*$を行いたい場合は、updating版(3), (4)に更新前の行列をEとして渡す。
  • エルミート性を維持するため、alphaは実部のみが使用される。
  • エルミート行列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 Vector>
void init_vec(Vector& v) {
  for (int i = 0; i < v.extent(0); ++i) {
    v[i] = std::complex<double>(0, i);
  }
}

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;

  using PackedMatrix = 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>>;

  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);

  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^*
  std::cout << "overwriting (1)\n";
  std::linalg::hermitian_matrix_rank_1_update(
    2.0,
    x,
    A,
    std::linalg::upper_triangle);
  print_mat(A);

  // (3) updating: A = E + alpha xx^*
  init_mat(E);
  std::cout << "updating (3)\n";
  std::linalg::hermitian_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,0) (0,0)
(0,0) (2,0) (4,0) (6,0)
(0,0) (4,0) (8,0) (12,0)
(0,0) (6,0) (12,0) (18,0)
updating (3)
(0,0) (0,1) (0,2) (0,3)
(0,1) (3,0) (5,2) (7,3)
(0,2) (5,2) (10,0) (14,3)
(0,3) (7,3) (14,3) (21,0)

バージョン

言語

  • C++26

処理系

関連項目

参照