Skip to content

Commit e7325de

Browse files
committed
linalg : matrix_rank_1_update_cを追加 (#1233)
Signed-off-by: Yuya Asano <64895419+sukeya@users.noreply.github.com>
1 parent 2a6ccb4 commit e7325de

File tree

2 files changed

+180
-1
lines changed

2 files changed

+180
-1
lines changed

reference/linalg.md

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -87,7 +87,7 @@ BLAS 1, 2, 3のアルゴリズムでテンプレートパラメータが特に
8787
| [`triangular_matrix_vector_product`](linalg/triangular_matrix_vector_product.md) | xTRMV: 三角行列とベクトルの積を求める (function template) | C++26 |
8888
| [`triangular_matrix_vector_solve`](linalg/triangular_matrix_vector_solve.md) | xTRSV: 三角行列を係数とする行列方程式を解く (function template) | C++26 |
8989
| [`matrix_rank_1_update`](linalg/matrix_rank_1_update.md) | xGER, xGERU: 行列のRank-1更新 (function template) | C++26 |
90-
| `matrix_rank_1_update_c` | xGERC: 複素行列のRank-1更新 (function template) | C++26 |
90+
| [`matrix_rank_1_update_c`](linalg/matrix_rank_1_update_c.md) | xGERC: 複素行列のRank-1更新 (function template) | C++26 |
9191
| `symmetric_matrix_rank_1_update` | xSYR: 対称行列のRank-1更新 (function template) | C++26 |
9292
| `hermitian_matrix_rank_1_update` | xHER: ハミルトニアン行列のRank-1更新 (function template) | C++26 |
9393
| `symmetric_matrix_rank_2_update` | xSYR2: 対称行列のRank-2更新 (function template) | C++26 |
Lines changed: 179 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,179 @@
1+
# matrix_rank_1_update_c
2+
3+
4+
* [mathjax enable]
5+
* linalg[meta header]
6+
* function template[meta id-type]
7+
* std::linalg[meta namespace]
8+
* cpp26[meta cpp]
9+
10+
11+
```cpp
12+
namespace std::linalg {
13+
template<in-vector InVec1,
14+
in-vector InVec2,
15+
inout-matrix InOutMat>
16+
void matrix_rank_1_update_c(
17+
InVec1 x,
18+
InVec2 y,
19+
InOutMat A); // (1)
20+
21+
template<class ExecutionPolicy,
22+
in-vector InVec1,
23+
in-vector InVec2,
24+
inout-matrix InOutMat>
25+
void matrix_rank_1_update_c(
26+
ExecutionPolicy&& exec,
27+
InVec1 x,
28+
InVec2 y,
29+
InOutMat A); // (2)
30+
}
31+
```
32+
33+
34+
## 概要
35+
非対称で共役を取るrank-1 updateを行う。
36+
37+
- (1): $A \leftarrow A + xy^*$
38+
- (2): (1)を指定された実行ポリシーで実行する。
39+
40+
41+
## 適格要件
42+
- [`possibly-multipliable`](possibly-multipliable.md)`<InOutMat, InVec2, InVec1>() == true`
43+
44+
45+
## 事前条件
46+
- [`multipliable`](multipliable.md)`(A, y, x) == true`
47+
48+
49+
## 効果
50+
- (1): [`matrix_rank_1_update`](matrix_rank_1_update.md)`(x, conjugated(y), A)`と同じ。
51+
- (2): `matrix_rank_1_update(std::forward<ExecutionPolicy>(exec), x, conjugated(y), A)`と同じ。
52+
53+
54+
## 戻り値
55+
なし
56+
57+
58+
## 計算量
59+
$O(\verb|x.extent(0)|\times \verb|y.extent(0)|)$
60+
61+
62+
## 例
63+
**[注意] 処理系にあるコンパイラで確認していないため、間違っているかもしれません。**
64+
65+
```cpp example
66+
#include <array>
67+
#include <iostream>
68+
#include <linalg>
69+
#include <mdspan>
70+
#include <vector>
71+
72+
template <class Matrix>
73+
void print_mat(const Matrix& A) {
74+
for(int i = 0; i < A.extent(0); ++i) {
75+
for(int j = 0; j < A.extent(1) - 1; ++j) {
76+
std::cout << A[i, j] << ' ';
77+
}
78+
std::cout << A[i, A.extent(1) - 1] << '\n';
79+
}
80+
}
81+
82+
template <class Vector>
83+
void init_vec(Vector& v) {
84+
for (int i = 0; i < v.extent(0); ++i) {
85+
v[i] = std::complex<double>(i, 0);
86+
}
87+
}
88+
89+
template <class Matrix>
90+
void init_mat(Matrix& A) {
91+
for(int i = 0; i < A.extent(0); ++i) {
92+
for(int j = 0; j < A.extent(1); ++j) {
93+
A[i,j] = std::complex<double>(i, j);
94+
}
95+
}
96+
}
97+
98+
int main()
99+
{
100+
constexpr size_t N = 4;
101+
102+
std::vector<std::complex<double>> A_vec(N * N);
103+
std::vector<std::complex<double>> x_vec(N);
104+
std::array<std::complex<double>, N> y_vec;
105+
106+
std::mdspan<
107+
std::complex<double>,
108+
std::extents<size_t, N, N>> A(A_vec.data());
109+
std::mdspan x(x_vec.data(), N);
110+
std::mdspan y(y_vec.data(), N);
111+
112+
init_mat(A);
113+
init_vec(x);
114+
init_vec(y);
115+
116+
for (int i = 0; i < y.extent(0); ++i) {
117+
y[i] *= std::complex<double>(0, 1);
118+
}
119+
120+
// (1)
121+
std::cout << "(1)\n";
122+
std::linalg::matrix_rank_1_update_c(
123+
x,
124+
y,
125+
A);
126+
print_mat(A);
127+
128+
init_mat(A);
129+
130+
// (2)
131+
std::cout << "(2)\n";
132+
std::linalg::matrix_rank_1_update_c(
133+
std::execution::par,
134+
x,
135+
y,
136+
A);
137+
print_mat(A);
138+
139+
return 0;
140+
}
141+
```
142+
143+
144+
### 出力
145+
```
146+
(1)
147+
(0,0) (0,1) (0,2) (0,3)
148+
(1,0) (1,0) (1,0) (1,0)
149+
(2,0) (2,-1) (2,-2) (2,-3)
150+
(3,0) (3,-2) (3,-4) (3,-6)
151+
(2)
152+
(0,0) (0,1) (0,2) (0,3)
153+
(1,0) (1,0) (1,0) (1,0)
154+
(2,0) (2,-1) (2,-2) (2,-3)
155+
(3,0) (3,-2) (3,-4) (3,-6)
156+
```
157+
158+
159+
## バージョン
160+
### 言語
161+
- C++26
162+
163+
### 処理系
164+
- [Clang](/implementation.md#clang): ??
165+
- [GCC](/implementation.md#gcc): ??
166+
- [ICC](/implementation.md#icc): ??
167+
- [Visual C++](/implementation.md#visual_cpp): ??
168+
169+
170+
## 関連項目
171+
- [`execution`](/reference/execution.md)
172+
- [`mdspan`](/reference/mdspan.md)
173+
- [`matrix_rank_1_update`](matrix_rank_1_update.md)
174+
175+
176+
## 参照
177+
- [P1673R13 A free function linear algebra interface based on the BLAS](https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2023/p1673r13.html)
178+
- [LAPACK: ger: general matrix rank-1 update](https://netlib.org/lapack/explore-html/d8/d75/group__ger.html)
179+

0 commit comments

Comments
 (0)