namespace std::linalg {
template<class Scalar,
in-matrix InMat,
possibly-packed-inout-matrix InOutMat,
class Triangle>
void hermitian_matrix_rank_k_update(
Scalar alpha,
InMat A,
InOutMat C,
Triangle t); // (1)
template<class ExecutionPolicy,
class Scalar,
in-matrix InMat,
possibly-packed-inout-matrix InOutMat,
class Triangle>
void hermitian_matrix_rank_k_update(
ExecutionPolicy&& exec,
Scalar alpha,
InMat A,
InOutMat C,
Triangle t); // (2)
template<in-matrix InMat,
possibly-packed-inout-matrix InOutMat,
class Triangle>
void hermitian_matrix_rank_k_update(
InMat A,
InOutMat C,
Triangle t); // (3)
template<class ExecutionPolicy,
in-matrix InMat,
possibly-packed-inout-matrix InOutMat,
class Triangle>
void hermitian_matrix_rank_k_update(
ExecutionPolicy&& exec,
InMat A,
InOutMat C,
Triangle t); // (4)
}
概要
エルミートな(対称かつ共役を取る)rank-k updateをエルミート行列に行う。
引数t
は対称行列の成分が上三角にあるのか、それとも下三角にあるのかを示す。
- (1): $C \leftarrow C + \alpha AA^*$
- (2): (1)を指定された実行ポリシーで実行する。
- (3): $C \leftarrow C + AA^*$
- (4): (3)を指定された実行ポリシーで実行する。
適格要件
- 共通:
Triangle
はupper_triangle_t
またはlower_triangle_t
InMat
がlayout_blas_packed
を持つなら、レイアウトのTriangle
テンプレート引数とこの関数のTriangle
テンプレート引数が同じ型compatible-static-extents<decltype(A), decltype(A)>(0, 1)
がtrue
(つまりA
が正方行列であること)compatible-static-extents<decltype(C), decltype(C)>(0, 1)
がtrue
(つまりC
が正方行列であること)compatible-static-extents<decltype(A), decltype(C)>(0, 0)
がtrue
(つまりA
の次元とC
の次元が同じであること)
- (2), (4):
is_execution_policy<ExecutionPolicy>::value
がtrue
事前条件
A.extent(0) == A.extent(1)
C.extent(0) == C.extent(1)
A.extent(0) == C.extent(0)
効果
- (1), (2): $C \leftarrow C + \alpha AA^*$
- (3), (4): $C \leftarrow C + AA^*$
戻り値
なし
計算量
$O(\verb|A.extent(0)| \times \verb|A.extent(1)| \times \verb|C.extent(0)|)$
例
[注意] 処理系にあるコンパイラで確認していないため、間違っているかもしれません。
#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 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] = std::complex<double>(i, j);
}
}
}
template <class Matrix>
void init_herm_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 = 2;
using Scalar = std::complex<double>;
using Vector = std::vector<Scalar>;
using HermitianMatrix = std::mdspan<
Scalar,
std::extents<size_t, N, N>,
std::linalg::layout_blas_packed<
std::linalg::upper_triangle_t,
std::linalg::row_major_t>
>;
Vector A_vec(N * N);
Vector C_vec(N * N);
std::mdspan A(A_vec.data());
HermitianMatrix C(C_vec.data());
init_mat(A);
init_herm_mat(C);
// (1)
std::cout << "(1)\n";
std::linalg::hermitian_matrix_rank_k_update(
-1.0,
A,
C,
std::linalg::upper_triangle);
print_mat(C);
// (2)
init_herm_mat(C);
std::cout << "(2)\n";
std::linalg::hermitian_matrix_rank_k_update(
std::execution::par,
-1.0,
A,
C,
std::linalg::upper_triangle);
print_mat(C);
// (3)
init_herm_mat(C);
std::cout << "(3)\n";
std::linalg::hermitian_matrix_rank_k_update(
A,
C,
std::linalg::upper_triangle);
print_mat(C);
// (4)
init_herm_mat(C);
std::cout << "(4)\n";
std::linalg::hermitian_matrix_rank_k_update(
std::execution::par,
A,
C,
std::linalg::upper_triangle);
print_mat(C);
return 0;
}
出力
(1)
(-1,0) (-1,0)
(-1,0) (-1,0)
(2)
(-1,0) (-1,0)
(-1,0) (-1,0)
(3)
(1,0) (1,2)
(1,-2) (3,0)
(4)
(1,0) (1,2)
(1,-2) (3,0)
バージョン
言語
- C++26
処理系
- Clang: ??
- GCC: ??
- ICC: ??
- Visual C++: ??