ARTS 2.5.4 (git: 4c0d3b4d)
|
Linear algebra functions. More...
#include "lin_alg.h"
#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include <cmath>
#include <new>
#include <stdexcept>
#include "array.h"
#include "arts.h"
#include "arts_omp.h"
#include "lapack.h"
#include "logic.h"
#include "matpackIII.h"
Go to the source code of this file.
Linear algebra functions.
This file contains mathematical tools to solve the vector radiative transfer equation.
Definition in file lin_alg.cc.
void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__eigen | ( | MatrixView | F, |
ConstMatrixView | A | ||
) |
Special method developed by Philippe Baron to solve cases similar to Zeeman matrices,
| a b c d | A = - | b a u v | r | c -u a w |
F = exp(A)
The method is to use the Cayley-Hamilton theorem that states that you can get at the function of a matrix by fitting coefficients as found from the eigenvalues to an expansion in the matrix...
In our case, F = (C0*I + C1*A + C2*A^2 + C3*A^3) * exp(-ar), where the Cs are from solving the problem
c0 + c1*x + c2*x^2 + c3*x^3, x := exp(x) c0 - c1*x + c2*x^2 - c3*x^3, x := exp(-x) c0 + c1*y + c2*y^2 + c3*y^3, y := exp(y) c0 - c1*y + c2*y^2 - c3*y^3, y := exp(-y),
for the pairs of x, y eigenvalues that can be analytically found (using preferably symbolic expressions or via handwork) from the off-diagonal elements of A. The variables Const1 and Const2 below denotes the analytical solution for getting x and y
F | Output: The matrix exponential of A (Has to be initialized before calling the function. |
dF | Output: The derivative of the matrix exponential of A (Has to be initialized before calling the function. Page dimension is for each derivative. |
A | Input: arbitrary square matrix. |
dA | Input: derivative of A. Page dimension is for each derivative. |
Definition at line 958 of file lin_alg.cc.
References a, b, b2, c, d, MapToEigen4x4(), sqrt(), u, v, and w.
Referenced by ext2trans(), and test_matrix_exp_propmat().
void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__eigen | ( | MatrixView | F, |
Tensor3View | dF_upp, | ||
Tensor3View | dF_low, | ||
ConstMatrixView | A, | ||
ConstTensor3View | dA_upp, | ||
ConstTensor3View | dA_low | ||
) |
Special method developed by Philippe Baron to solve cases similar to Zeeman matrices,
| a b c d | A = - | b a u v | r | c -u a w |
F = exp(A)
The method is to use the Cayley-Hamilton theorem that states that you can get at the function of a matrix by fitting coefficients as found from the eigenvalues to an expansion in the matrix...
In our case, F = (C0*I + C1*A + C2*A^2 + C3*A^3) * exp(-ar), where the Cs are from solving the problem
c0 + c1*x + c2*x^2 + c3*x^3, x := exp(x) c0 - c1*x + c2*x^2 - c3*x^3, x := exp(-x) c0 + c1*y + c2*y^2 + c3*y^3, y := exp(y) c0 - c1*y + c2*y^2 - c3*y^3, y := exp(-y),
for the pairs of x, y eigenvalues that can be analytically found (using preferably symbolic expressions or via handwork) from the off-diagonal elements of A. The variables Const1 and Const2 below denotes the analytical solution for getting x and y
Derivatives then simply follows from the analytical expression such that
dF = (dC0*I + dC1*A + C1*dA + dC2*A^2 + C2*(A*dA+dA*A) + dC3*A^3 + C3*(A*A*dA+A*dA*A+A*A*dA)) * exp(-ar) - F * da * r
where
| da db dc dd | dA = - | db da du dv | r, | dc -du da dw |
and the dCs can be found analytically as well by simply taking the derivative of the problem, regarding x and y as functions of some variable
F | Output: The matrix exponential of A (Has to be initialized before calling the function. |
dF | Output: The derivative of the matrix exponential of A (Has to be initialized before calling the function. Page dimension is for each derivative. |
A | Input: arbitrary square matrix. |
dA | Input: derivative of A. Page dimension is for each derivative. |
Definition at line 1545 of file lin_alg.cc.
References a, b, b2, c, d, MapToEigen4x4(), ConstTensor3View::npages(), sqrt(), u, v, and w.
void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__explicit | ( | MatrixView | F, |
ConstMatrixView | A | ||
) |
Special method developed by Philippe Baron to solve cases similar to Zeeman matrices,
| a b c d | A = - | b a u v | r | c -u a w |
F = exp(A)
The method is to use the Cayley-Hamilton theorem that states that you can get at the function of a matrix by fitting coefficients as found from the eigenvalues to an expansion in the matrix...
In our case, F = (C0*I + C1*A + C2*A^2 + C3*A^3) * exp(-ar), where the Cs are from solving the problem
c0 + c1*x + c2*x^2 + c3*x^3, x := exp(x) c0 - c1*x + c2*x^2 - c3*x^3, x := exp(-x) c0 + c1*y + c2*y^2 + c3*y^3, y := exp(y) c0 - c1*y + c2*y^2 - c3*y^3, y := exp(-y),
for the pairs of x, y eigenvalues that can be analytically found (using preferably symbolic expressions or via handwork) from the off-diagonal elements of A. The variables Const1 and Const2 below denotes the analytical solution for getting x and y
F | Output: The matrix exponential of A (Has to be initialized before calling the function. |
dF | Output: The derivative of the matrix exponential of A (Has to be initialized before calling the function. Page dimension is for each derivative. |
A | Input: arbitrary square matrix. |
dA | Input: derivative of A. Page dimension is for each derivative. |
Definition at line 784 of file lin_alg.cc.
References a, b, b2, c, d, sqrt(), u, v, and w.
Referenced by test_matrix_exp_propmat().
void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__explicit | ( | MatrixView | F, |
Tensor3View | dF_upp, | ||
Tensor3View | dF_low, | ||
ConstMatrixView | A, | ||
ConstTensor3View | dA_upp, | ||
ConstTensor3View | dA_low | ||
) |
Special method developed by Philippe Baron to solve cases similar to Zeeman matrices,
| a b c d | A = - | b a u v | r | c -u a w |
F = exp(A)
The method is to use the Cayley-Hamilton theorem that states that you can get at the function of a matrix by fitting coefficients as found from the eigenvalues to an expansion in the matrix...
In our case,
F = (C0*I + C1*A + C2*A^2 + C3*A^3) * exp(-ar),
where the Cs are from solving the problem
c0 + c1*x + c2*x^2 + c3*x^3, x := exp(x) c0 - c1*x + c2*x^2 - c3*x^3, x := exp(-x) c0 + c1*y + c2*y^2 + c3*y^3, y := exp(y) c0 - c1*y + c2*y^2 - c3*y^3, y := exp(-y),
for the pairs of x, y eigenvalues that can be analytically found (using preferably symbolic expressions or via handwork) from the off-diagonal elements of A. The variables Const1 and Const2 below denotes the analytical solution for getting x and y
Derivatives then simply follows from the analytical expression such that
dF = (dC0*I + dC1*A + C1*dA + dC2*A^2 + C2*(A*dA+dA*A) + dC3*A^3 + C3*(A*A*dA+A*dA*A+A*A*dA)) * exp(-ar) - F * da * r
where
| da db dc dd | dA = - | db da du dv | r, | dc -du da dw |
and the dCs can be found analytically as well by simply taking the derivative of the problem, regarding x and y as functions of some variable
F | Output: The matrix exponential of A (Has to be initialized before calling the function. |
dF | Output: The derivative of the matrix exponential of A (Has to be initialized before calling the function. Page dimension is for each derivative. |
A | Input: arbitrary square matrix. |
dA | Input: derivative of A. Page dimension is for each derivative. |
Definition at line 1106 of file lin_alg.cc.
References a, b, b2, c, d, ConstTensor3View::npages(), sqrt(), u, v, and w.
Numeric det | ( | ConstMatrixView | A | ) |
Determinant of N by N matrix. Simple recursive method.
A | In: Matrix of size NxN. |
Definition at line 2251 of file lin_alg.cc.
References ARTS_ASSERT, det(), ConstMatrixView::ncols(), ConstMatrixView::nrows(), and temp.
Referenced by det().
void diagonalize | ( | ComplexMatrixView | P, |
ComplexVectorView | W, | ||
const ConstComplexMatrixView | A | ||
) |
Matrix Diagonalization.
Return P and W from A in the statement diag(P^-1*A*P)-W == 0.
The function makes many copies and is thereby not fast. There are no tests on the condition of the returned matrix, so nan and inf can occur.
[out] | P | The right eigenvectors. |
[out] | W | The eigenvalues. |
[in] | A | The matrix to diagonalize. |
Definition at line 326 of file lin_alg.cc.
References ARTS_ASSERT, ComplexVectorView::get_c_array(), ComplexMatrixView::get_c_array(), VectorView::get_c_array(), ConstComplexMatrixView::get_column_extent(), ConstComplexMatrixView::ncols(), ConstComplexVectorView::nelem(), ConstComplexMatrixView::nrows(), swap(), and lapack::zgeev_().
void diagonalize | ( | MatrixView | P, |
VectorView | WR, | ||
VectorView | WI, | ||
ConstMatrixView | A | ||
) |
Matrix Diagonalization.
Return P and W from A in the statement diag(P^-1*A*P)-W == 0. The real function will require some manipulation if the eigenvalues are imaginary.
The real version returns WR and WI as returned by dgeev. The complex version just returns W.
The function makes many copies and is thereby not fast. There are no tests on the condition of the returned matrix, so nan and inf can occur.
[out] | P | The right eigenvectors. |
[out] | WR | The real eigenvalues. |
[out] | WI | The imaginary eigenvalues. |
[in] | A | The matrix to diagonalize. |
Definition at line 241 of file lin_alg.cc.
Referenced by Absorption::LineMixing::EquivalentLines::EquivalentLines(), lm_hitran_2017::eqvlines(), matrix_exp2(), test_complex_diagonalize(), and test_real_diagonalize().
Eigen::ComplexEigenSolver< Eigen::MatrixXcd > eig | ( | const Eigen::Ref< Eigen::MatrixXcd > | A | ) |
Return the Eigen decomposition of the eigen matrix.
[in] | A | a matrix to eigen value decompose |
Definition at line 2368 of file lin_alg.cc.
Referenced by Absorption::LineMixing::eigenvalue_adaptation_of_relmat().
void id_mat | ( | MatrixView | I | ) |
Identity Matrix.
I | Output: identity matrix |
Definition at line 2235 of file lin_alg.cc.
References ARTS_ASSERT, ConstMatrixView::ncols(), and ConstMatrixView::nrows().
Referenced by get_ppath_transmat(), matrix_exp(), MatrixIdentity(), MCGeneral(), mcPathTraceRadar(), rte_step_doit_replacement(), sensor_responseInit(), SparseIdentity(), test47(), test_identity(), test_inv(), and test_oem_gauss_newton_sparse().
void inv | ( | ComplexMatrixView | Ainv, |
const ConstComplexMatrixView | A | ||
) |
Definition at line 202 of file lin_alg.cc.
References ARTS_ASSERT, ARTS_USER_ERROR_IF, ComplexVectorView::get_c_array(), ComplexMatrixView::get_c_array(), ConstComplexMatrixView::ncols(), ConstComplexMatrixView::nrows(), lapack::zgetrf_(), and lapack::zgetri_().
void inv | ( | MatrixView | Ainv, |
ConstMatrixView | A | ||
) |
Matrix Inverse.
Compute the inverse of a matrix such that I = Ainv*A = A*Ainv. Both MatrixViews must be square and have the same size n. During the inversion one additional n times n Matrix is allocated and work space memory for faster inversion is allocated and freed.
[out] | Ainv | The MatrixView to contain the inverse of A. |
[in] | A | The matrix to be inverted. |
Definition at line 167 of file lin_alg.cc.
Referenced by benchmark_inv(), benchmark_oem_linear(), lm_hitran_2017::eqvlines(), LineShape::ComputeData::interp_add_even(), LineShape::ComputeData::interp_add_triplequad(), CovarianceMatrix::invert_correlation_block(), matrix_exp2(), test_complex_diagonalize(), test_inv(), test_oem_gauss_newton(), test_oem_levenberg_marquardt(), and test_real_diagonalize().
void linreg | ( | Vector & | p, |
ConstVectorView | x, | ||
ConstVectorView | y | ||
) |
Determines coefficients for linear regression
Performs a least squares estimation of the model
y = p[0] + p[1] * x
p | Out: Fitted coefficients. |
x | In: x-value of data points |
y | In: y-value of data points |
Definition at line 2297 of file lin_alg.cc.
References ARTS_ASSERT, ConstVectorView::nelem(), and Vector::resize().
Referenced by derive_scat_species_a_and_b(), and ppathFromRtePos2().
|
noexcept |
Least squares fitting by solving x for known A and y.
(A^T A)x = A^T y
Returns the squared residual, i.e., <(A^T A)x-A^T y|(A^T A)x-A^T y>.
[in] | x | As equation |
[in] | A | As equation |
[in] | y | As equation |
[in] | residual | (optional) Returns the residual if true |
Definition at line 2345 of file lin_alg.cc.
References mult(), solve(), and transpose().
void lubacksub | ( | VectorView | x, |
ConstMatrixView | LU, | ||
ConstVectorView | b, | ||
const ArrayOfIndex & | indx | ||
) |
LU backsubstitution.
Solves a set of linear equations Ax=b. It is neccessairy to do a L decomposition using the function ludcp before using this function. The backsubstitution is in-place, i.e. x and b may be the same vector.
x | Output: Solution vector of the equation system. |
LU | Input: LU decomposition of the matrix (output of function ludcp). |
b | Input: Right-hand-side vector of equation system. |
indx | Input: Pivoting information (output of function ludcp). |
Definition at line 91 of file lin_alg.cc.
Referenced by solve(), test_lusolve1D(), test_lusolve4D(), and test_solve_linear_system().
void ludcmp | ( | Matrix & | LU, |
ArrayOfIndex & | indx, | ||
ConstMatrixView | A | ||
) |
LU decomposition.
This function performes a LU Decomposition of the matrix A. (Compare Numerical Recipies in C, pages 36-48.)
LU | Output: returns L and U in one matrix |
indx | Output: Vector that records the row permutation. |
A | Input: Matrix for which the LU decomposition is performed |
Definition at line 56 of file lin_alg.cc.
Referenced by solve(), test_lusolve1D(), test_lusolve4D(), and test_solve_linear_system().
MatrixViewMap MapToEigen | ( | Tensor3View & | A, |
const Index & | i | ||
) |
Definition at line 745 of file lin_alg.cc.
References joker, and MapToEigen().
Referenced by Zeeman::dsum(), MapToEigen(), Zeeman::sum(), test03(), test48(), and test_matrix_exp_propmat().
Matrix4x4ViewMap MapToEigen4x4 | ( | Tensor3View & | A, |
const Index & | i | ||
) |
Definition at line 740 of file lin_alg.cc.
References joker, and MapToEigen4x4().
Referenced by cayley_hamilton_fitted_method_4x4_propmat_to_transmat__eigen(), MapToEigen4x4(), and matrix_exp2_4x4().
void matrix_exp | ( | MatrixView | F, |
ConstMatrixView | A, | ||
const Index & | q | ||
) |
General exponential of a Matrix.
The exponential of a matrix is computed using the Pade-Approximation. The method is decribed in: Golub, G. H. and C. F. Van Loan, Matrix Computation, p. 384, Johns Hopkins University Press, 1983.
The Pade-approximation is applied on all cases. If a faster option can be applied has to be checked before calling the function.
F | Output: The matrix exponential of A (Has to be initialized before calling the function. |
A | Input: arbitrary square matrix |
q | Input: Parameter for the accuracy of the computation |
Definition at line 387 of file lin_alg.cc.
References ARTS_ASSERT, c, id_mat(), is_size(), N, ConstMatrixView::ncols(), norm_inf(), pow(), and q.
Referenced by ext2trans(), test_matrix_exp1D(), test_matrix_exp3D(), test_matrix_exp4D(), test_matrix_exp_propmat(), and test_real_diagonalize().
void matrix_exp2 | ( | MatrixView | F, |
ConstMatrixView | A | ||
) |
Definition at line 454 of file lin_alg.cc.
References ARTS_ASSERT, diagonalize(), inv(), is_size(), and ConstMatrixView::nrows().
Referenced by test_real_diagonalize().
void matrix_exp2_4x4 | ( | MatrixView | arts_f, |
ConstMatrixView | A | ||
) |
Definition at line 516 of file lin_alg.cc.
References ARTS_ASSERT, is_size(), MapToEigen4x4(), and q.
void matrix_exp_4x4 | ( | MatrixView | arts_f, |
ConstMatrixView | A, | ||
const Index & | q | ||
) |
Definition at line 476 of file lin_alg.cc.
References ARTS_ASSERT, is_size(), and norm_inf().
void matrix_exp_dmatrix_exp | ( | MatrixView | F, |
MatrixView | dF, | ||
ConstMatrixView | A, | ||
ConstMatrixView | dA, | ||
const Index & | q | ||
) |
General exponential of a Matrix with their derivatives.
The exponential of a matrix is computed using the Pade-Approximation. The method is decribed in: Golub, G. H. and C. F. Van Loan, Matrix Computation, p. 384, Johns Hopkins University Press, 1983.
The extension of Pade-Approximation to the derivative is explained by Lubomír Brančík, MATLAB PROGRAMS FOR MATRIX EXPONENTIAL FUNCTION DERIVATIVE EVALUATION, 2008
The Pade-approximation is applied on all cases. If a faster option can be applied has to be checked before calling the function.
F | Output: The matrix exponential of A (Has to be initialized before calling the function. |
dF | Output: The derivative of the matrix exponential of A (Has to be initialized before calling the function. |
A | Input: arbitrary square matrix |
dA | Input: derivative of arbitrary square matrix |
q | Input: Parameter for the accuracy of the computation |
Definition at line 2082 of file lin_alg.cc.
References ARTS_ASSERT, is_size(), and ConstMatrixView::ncols().
void matrix_exp_dmatrix_exp | ( | MatrixView | F, |
Tensor3View | dF, | ||
ConstMatrixView | A, | ||
ConstTensor3View | dA, | ||
const Index & | q | ||
) |
General exponential of a Matrix with their derivatives.
The exponential of a matrix is computed using the Pade-Approximation. The method is decribed in: Golub, G. H. and C. F. Van Loan, Matrix Computation, p. 384, Johns Hopkins University Press, 1983.
The extension of Pade-Approximation to the derivative is explained by Lubomír Brančík, MATLAB PROGRAMS FOR MATRIX EXPONENTIAL FUNCTION DERIVATIVE EVALUATION, 2008
The Pade-approximation is applied on all cases. If a faster option can be applied has to be checked before calling the function.
F | Output: The matrix exponential of A (Has to be initialized before calling the function. |
dF | Output: The derivative of the matrix exponential of A (Has to be initialized before calling the function. Page dimension is for each derivative. |
A | Input: arbitrary square matrix. |
dA | Input: derivative of A. Page dimension is for each derivative. |
q | Input: Parameter for the accuracy of the computation, Matlab default is 6. |
Definition at line 1921 of file lin_alg.cc.
References ARTS_ASSERT, is_size(), joker, ConstMatrixView::ncols(), and ConstTensor3View::npages().
Numeric norm_inf | ( | ConstMatrixView | A | ) |
Maximum absolute row sum norm.
This function returns the maximum absolute row sum norm of a matrix A (see user guide for the definition).
A | Input: arbitrary matrix |
Definition at line 2218 of file lin_alg.cc.
References abs, ConstMatrixView::ncols(), norm_inf(), and ConstMatrixView::nrows().
Referenced by matrix_exp(), matrix_exp_4x4(), norm_inf(), and propmat4x4_to_transmat4x4().
void propmat4x4_to_transmat4x4 | ( | MatrixView | F, |
Tensor3View | dF_upp, | ||
Tensor3View | dF_low, | ||
ConstMatrixView | A, | ||
ConstTensor3View | dA_upp, | ||
ConstTensor3View | dA_low, | ||
const Index & | q | ||
) |
Definition at line 1778 of file lin_alg.cc.
References ARTS_ASSERT, is_size(), joker, norm_inf(), and ConstTensor3View::npages().
void solve | ( | VectorView | x, |
ConstMatrixView | A, | ||
ConstVectorView | b | ||
) |
Solve a linear system.
Solves the linear system A*x = b for a general matrix A. For the solution of the system an additional n-times-n matrix and a size-n index vector are allocated.
x | The solution vector x. |
A | The matrix A defining the system. |
b | The vector b. |
Definition at line 138 of file lin_alg.cc.
References ARTS_ASSERT, b, lubacksub(), ludcmp(), ConstMatrixView::ncols(), ConstVectorView::nelem(), and ConstMatrixView::nrows().
Referenced by lsf().
void special_matrix_exp_and_dmatrix_exp_dx_for_rt | ( | MatrixView | F, |
Tensor3View | dF_upp, | ||
Tensor3View | dF_low, | ||
ConstMatrixView | A, | ||
ConstTensor3View | dA_upp, | ||
ConstTensor3View | dA_low, | ||
const Index & | q | ||
) |
Special exponential of a Matrix with their derivatives.
The exponential of a matrix is computed using the Pade-Approximation. The method is decribed in: Golub, G. H. and C. F. Van Loan, Matrix Computation, p. 384, Johns Hopkins University Press, 1983.
The extension of Pade-Approximation to the derivative is explained by Lubomír Brančík, MATLAB PROGRAMS FOR MATRIX EXPONENTIAL FUNCTION DERIVATIVE EVALUATION, 2008
The Pade-approximation is applied on all cases. If a faster option can be applied has to be checked before calling the function.
F | Output: The matrix exponential of A (Has to be initialized before calling the function. |
dF | Output: The derivative of the matrix exponential of A (Has to be initialized before calling the function. Page dimension is for each derivative. |
A | Input: arbitrary square matrix. |
dA | Input: derivative of A. Page dimension is for each derivative. |
q | Input: Parameter for the accuracy of the computation, Matlab default is 6. |
Definition at line 554 of file lin_alg.cc.
References ARTS_ASSERT, is_size(), joker, ConstMatrixView::ncols(), and ConstTensor3View::npages().
Referenced by test_matrix_exp_propmat().