36#include <Eigen/Eigenvalues>
63 auto ipiv = std::vector<int>(n);
75 for (
Index i = 0; i < n; i++) {
111 std::vector<int> ipiv(n);
112 std::vector<double> rhs(n);
115 for (
Index i = 0; i < n; i++) {
116 ipiv[i] = (int)indx[i];
121 &trans, &n_int, &one, LU.
mdata, &n_int, ipiv.data(), rhs.data(), &n_int, &info);
123 for (
Index i = 0; i < n; i++) {
175 auto ipiv = std::vector<int>(n);
184 auto work = std::vector<double>(lwork);
190 "Error inverting matrix: Matrix not of full rank.");
202 int n_int = int(n), lwork = n_int, info;
203 std::vector<int> ipiv(n);
204 std::vector<Complex> work(lwork);
208 &n_int, &n_int, Ainv.
get_c_array(), &n_int, ipiv.data(), &info);
219 "Error inverting matrix: Matrix not of full rank.");
259 int LDA, LDA_L, LDA_R, n_int, info;
264 char l_eig =
'N', r_eig =
'V';
267 int lwork = 2 * n_int + n_int * n_int;
268 auto work = std::vector<double>(lwork);
269 auto rwork = std::vector<double>(2 * n_int);
272 double* adata = A_tmp.
mdata;
273 double* rpdata = P2.
mdata;
274 double* wrdata = WR2.
mdata;
275 double* widata = WI2.
mdata;
295 for (
Index i = 0; i < n; i++)
296 for (
Index j = 0; j < n; j++) P(j, i) = P2(i, j);
331 char l_eig =
'N', r_eig =
'V';
334 int lwork = 2 * n_int + n_int * n_int;
355 for (
Index i = 0; i < n; i++)
356 for (
Index j = 0; j <= i; j++) std::swap(P(j, i), P(i, j));
383 Matrix D(n, n),
N(n, n), X(n, n), cX(n, n, 0.0), B(n, n, 0.0);
384 Vector N_col_vec(n, 0.), F_col_vec(n, 0.);
389 j = 1 + floor(1. / log(2.) * log(A_norm_inf));
392 auto j_index = (
Index)(j);
407 for (
Index k = 0; k < q; k++) {
409 c *= (q_n - k_n + 1) / ((2 * q_n - k_n + 1) * k_n);
427 for (
Index i = 0; i < n; i++) {
429 lubacksub(F_col_vec, X, N_col_vec, indx);
430 F(
joker, i) = F_col_vec;
434 for (
Index k = 0; k < j_index; k++) {
455 for (
Index i = 0; i < A.
ncols(); i++) row_sum +=
abs(A(i, j));
471 for (
Index i = 0; i < n; i++) I(i, i) = 1.;
487 return A(0, 0) * A(1, 1) * A(2, 2) + A(0, 1) * A(1, 2) * A(2, 0) +
488 A(0, 2) * A(1, 0) * A(2, 1) - A(0, 2) * A(1, 1) * A(2, 0) -
489 A(0, 1) * A(1, 0) * A(2, 2) - A(0, 0) * A(1, 2) * A(2, 1);
490 if (dim == 2)
return A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
491 if (dim == 1)
return A(0, 0);
495 for (
Index j = 0; j < dim; j++) {
496 Matrix temp(dim - 1, dim - 1);
497 for (
Index I = 1; I < dim; I++)
498 for (
Index J = 0; J < dim; J++) {
500 temp(I - 1, J) = A(I, J);
502 temp(I - 1, J - 1) = A(I, J);
507 ret_val += ((j % 2 == 0) ? -1. : 1.) * tempNum * A(0, j);
557 Numeric s1 = 0, xm = 0, s3 = 0, s4 = 0;
559 for (
Index i = 0; i < n; i++) {
563 for (
Index i = 0; i < n; i++) {
571 p[0] = s3 /
Numeric(n) - p[1] * xm;
577 bool residual)
noexcept {
579 const Index n = x.nelem();
This file contains the definition of Array.
The global header file for ARTS.
Header file for helper functions for OpenMP.
The ComplexMatrixView class.
The ComplexVectorView class.
A constant view of a ComplexMatrix.
Index ncols() const noexcept
Index nrows() const noexcept
Index get_column_extent() const
Get the extent of the underlying data.
Complex * get_c_array() const noexcept
Complex * get_c_array() const noexcept
Index nelem() const noexcept
A constant view of a Matrix.
Index nrows() const noexcept
Index ncols() const noexcept
Numeric * mdata
Pointer to the plain C array that holds the data.
Range mcr
The column range of mdata that is actually used.
A constant view of a Vector.
Numeric * mdata
Pointer to the plain C array that holds the data.
Numeric * get_c_array() const noexcept
Conversion to plain C-array, const-version.
Index nelem() const noexcept
Returns the number of elements.
constexpr Index get_extent() const noexcept
Returns the extent of the range.
constexpr Index get_stride() const noexcept
Returns the stride of the range.
void resize(Index n)
Resize function.
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_ASSERT(condition,...)
#define ARTS_USER_ERROR_IF(condition,...)
Interface for the LAPACK library.
void linreg(Vector &p, ConstVectorView x, ConstVectorView y)
void id_mat(MatrixView I)
Identity Matrix.
void solve(VectorView x, ConstMatrixView A, ConstVectorView b)
Solve a linear system.
Numeric norm_inf(ConstMatrixView A)
Maximum absolute row sum norm.
void diagonalize(MatrixView P, VectorView WR, VectorView WI, ConstMatrixView A)
Matrix Diagonalization.
void lubacksub(VectorView x, ConstMatrixView LU, ConstVectorView b, const ArrayOfIndex &indx)
LU backsubstitution.
void matrix_exp(MatrixView F, ConstMatrixView A, const Index &q)
General exponential of a Matrix.
void ludcmp(Matrix &LU, ArrayOfIndex &indx, ConstMatrixView A)
LU decomposition.
Numeric det(ConstMatrixView A)
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Numeric lsf(VectorView x, ConstMatrixView A, ConstVectorView y, bool residual) noexcept
Least squares fitting by solving x for known A and y.
Linear algebra functions.
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Header file for logic.cc.
void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
NUMERIC Numeric
The type to use for all floating point numbers.
INDEX Index
The type to use for all integer numbers and indices.
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
void zgetri_(int *n, std::complex< double > *A, int *lda, int *ipiv, std::complex< double > *work, int *lwork, int *info)
void zgeev_(char *jobvl, char *jobvr, int *n, std::complex< double > *A, int *lda, std::complex< double > *W, std::complex< double > *VL, int *ldvl, std::complex< double > *VR, int *ldvr, std::complex< double > *work, int *lwork, double *rwork, int *info)
void dgetri_(int *n, double *A, int *lda, int *ipiv, double *work, int *lwork, int *info)
Matrix inversion.
void dgeev_(char *jobvl, char *jobvr, int *n, double *A, int *lda, double *WR, double *WI, double *VL, int *ldvl, double *VR, int *ldvr, double *work, int *lwork, double *rwork, int *info)
void zgetrf_(int *m, int *n, std::complex< double > *A, int *lda, int *ipiv, int *info)
void dgetrs_(char *trans, int *n, int *nrhs, double *A, int *lda, int *ipiv, double *b, int *ldb, int *info)
Solve linear system of equations.
void dgetrf_(int *m, int *n, double *A, int *lda, int *ipiv, int *info)
LU decomposition.
Numeric pow(const Rational base, Numeric exp)
Power of.