17
18
19
24#include "../algorithms.hpp"
25#include "../basic_array.hpp"
26#include "../declarations.hpp"
27#include "../exceptions.hpp"
28#include "../layout/policies.hpp"
29#include "../layout_transforms.hpp"
30#include "../linalg.hpp"
31#include "../mapped_functions.hpp"
32#include "../matrix_functions.hpp"
34#include <itertools/itertools.hpp>
44namespace nda::lapack {
47
48
49
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
77 matrix<T> V_x_InvS_x_UH;
83 array<
double, 1> s_vec;
87
88
89
90 int n_var()
const {
return A.extent(1); }
93
94
95
99
100
101
102
103
104
105
106
107 gelss_worker(matrix<T> A_) : M(A_.extent(0)), N(A_.extent(1)), A(std::move(A_)), s_vec(std::min(M, N)) {
108 if (N > M)
NDA_RUNTIME_ERROR <<
"Error in nda::lapack::gelss_worker: Matrix A cannot have more columns than rows";
116 gesvd(A_FL, s_vec, U, VH);
119 matrix<
double,
F_layout> S_inv(N, M);
121 for (
long i : range(std::min(M, N))) S_inv(i, i) = 1.0 / s_vec(i);
122 V_x_InvS_x_UH = dagger(VH) * S_inv * dagger(U);
125 if (N < M) UH_NULL = dagger(U)(range(N, M), range(M));
129
130
131
132
133
134 std::pair<matrix<T>,
double>
operator()(matrix_const_view<T> B, std::optional<
long> = {})
const {
138 std::vector<
double> err_vec;
139 for (
long i : range(B.shape()[1])) err_vec.push_back(frobenius_norm(UH_NULL * B(range::all, range(i, i + 1))) / sqrt(B.shape()[0]));
140 err = *std::max_element(err_vec.begin(), err_vec.end());
142 return std::make_pair(V_x_InvS_x_UH * B, err);
146
147
148
149
150
151 std::pair<vector<T>,
double>
operator()(vector_const_view<T> b, std::optional<
long> = {})
const {
154 if (M != N) { err = norm(UH_NULL * b) / sqrt(b.size()); }
155 return std::make_pair(V_x_InvS_x_UH * b, err);
160
161
162
163
164
168 using dcomplex = std::complex<
double>;
181
182
183
187
188
189
193
194
195
199
200
201
202
203
204 std::pair<matrix<dcomplex>,
double>
operator()(matrix_const_view<dcomplex> B, std::optional<
long> inner_matrix_dim = {})
const {
205 if (
not inner_matrix_dim.has_value())
206 NDA_RUNTIME_ERROR <<
"Error in nda::lapack::gelss_worker_hermitian: Inner matrix dimension required for hermitian least square fitting";
207 long d = *inner_matrix_dim;
213 auto inner_adjoint = [d](
auto &M) {
214 auto idx_map = M.indexmap();
215 auto l = idx_map.lengths();
218 NDA_ASSERT2(l[1] % (d * d) == 0,
"Error in nda::lapack::gelss_worker_hermitian: Data shape incompatible with given dimension");
219 long N = l[1] / (d * d);
231 array<dcomplex, 4> arr_dag = conj(permuted_indices_view<
encode(std::array{0, 1, 3, 2}
)>(reshape(M, std::array{l[0], N, d, d})));
233 return matrix<dcomplex>{reshape(std::move(arr_dag), l)};
237 matrix<dcomplex> B_dag = inner_adjoint(B);
239 auto [x, err] = _lss_matrix
(B_stack
);
242 return {0.5
* (x
+ inner_adjoint(x)), err};
long extent(int i) const noexcept
Get the extent of the ith dimension.
basic_array(basic_array &&)=default
Default move constructor moves the memory handle and layout.
Worker class for solving linear least square problems.
array< double, 1 > const & S_vec() const
Get the singular value array.
gelss_worker(matrix< T > A_)
Construct a new worker object for a given matrix .
std::pair< vector< T >, double > operator()(vector_const_view< T > b, std::optional< long >={}) const
Solve the least-square problem for a given right hand side vector .
int n_var() const
Get the number of variables of the given problem.
std::pair< matrix< T >, double > operator()(matrix_const_view< T > B, std::optional< long >={}) const
Solve the least-square problem for a given right hand side matrix .
#define NDA_RUNTIME_ERROR
#define NDA_ASSERT2(X,...)
matrix< get_value_t< A > > vstack(A const &a, B const &b)
Stack two 2-dimensional arrays/views vertically.
decltype(auto) conj(A &&a)
Lazy, coefficient-wise complex conjugate function for nda::Array types.
Array auto operator+(L &&l, R &&r)
Addition operator for two nda::Array types.
auto operator*(L &&l, R &&r)
Multiplication operator for two nda::Array types.
constexpr uint64_t encode(std::array< int, N > const &a)
Encode a std::array<int, N> in a uint64_t.
Contiguous layout policy with Fortran-order (column-major order).
Worker class for solving linear least square problems for hermitian tail-fitting.
int n_var() const
Get the number of variables of the given problem.
array< double, 1 > const & S_vec() const
Get the singular value array.
std::pair< matrix< dcomplex >, double > operator()(matrix_const_view< dcomplex > B, std::optional< long > inner_matrix_dim={}) const
Solve the least-square problem for a given right hand side matrix .
gelss_worker_hermitian(matrix< dcomplex > A_)
Construct a new worker object for a given matrix .