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"

Functions
void	ludcmp (Matrix &LU, ArrayOfIndex &indx, ConstMatrixView A)
	LU decomposition. More...

void	lubacksub (VectorView x, ConstMatrixView LU, ConstVectorView b, const ArrayOfIndex &indx)
	LU backsubstitution. More...

void	solve (VectorView x, ConstMatrixView A, ConstVectorView b)
	Solve a linear system. More...

void	inv (MatrixView Ainv, ConstMatrixView A)
	Matrix Inverse. More...

void	inv (ComplexMatrixView Ainv, const ConstComplexMatrixView A)

void	diagonalize (MatrixView P, VectorView WR, VectorView WI, ConstMatrixView A)
	Matrix Diagonalization. More...

void	diagonalize (ComplexMatrixView P, ComplexVectorView W, const ConstComplexMatrixView A)
	Matrix Diagonalization. More...

void	matrix_exp (MatrixView F, ConstMatrixView A, const Index &q)
	General exponential of a Matrix. More...

Numeric	norm_inf (ConstMatrixView A)
	Maximum absolute row sum norm. More...

void	id_mat (MatrixView I)
	Identity Matrix. More...

Numeric	det (ConstMatrixView A)

void	linreg (Vector &p, ConstVectorView x, ConstVectorView y)

Numeric	lsf (VectorView x, ConstMatrixView A, ConstVectorView y, bool residual) noexcept
	Least squares fitting by solving x for known A and y. More...

Detailed Description

Linear algebra functions.

Author: Claudia Emde claud.nosp@m.ia.e.nosp@m.mde@d.nosp@m.lr.d.nosp@m.e

Date: Thu May 2 10:59:55 2002

This file contains mathematical tools to solve the vector radiative transfer equation.

Definition in file lin_alg.cc.

Function Documentation

◆ det()

Numeric det ( ConstMatrixView A )

Determinant of N by N matrix. Simple recursive method.

Parameters

A	In: Matrix of size NxN.

Author: Richard Larsson

Date: 2012-08-03

Definition at line 504 of file lin_alg.cc.

References ARTS_ASSERT, det(), ConstMatrixView::ncols(), ConstMatrixView::nrows(), and temp.

Referenced by det().

◆ diagonalize() [1/2]

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.

Parameters

[out]	P	The right eigenvectors.
[out]	W	The eigenvalues.
[in]	A	The matrix to diagonalize.

Definition at line 335 of file lin_alg.cc.

References ARTS_ASSERT, ConstComplexVectorView::get_c_array(), ConstComplexMatrixView::get_c_array(), ConstVectorView::get_c_array(), ConstComplexMatrixView::get_column_extent(), ConstComplexMatrixView::ncols(), ConstComplexVectorView::nelem(), ConstComplexMatrixView::nrows(), swap(), and lapack::zgeev_().

◆ diagonalize() [2/2]

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.

Parameters

[out]	P	The right eigenvectors.
[out]	WR	The real eigenvalues.
[out]	WI	The imaginary eigenvalues.
[in]	A	The matrix to diagonalize.

Definition at line 250 of file lin_alg.cc.

Referenced by Absorption::LineMixing::EquivalentLines::EquivalentLines(), lm_hitran_2017::eqvlines(), test_complex_diagonalize(), and test_real_diagonalize().

◆ id_mat()

void id_mat ( MatrixView I )

Identity Matrix.

Parameters

I	Output: identity matrix

Definition at line 488 of file lin_alg.cc.

References ARTS_ASSERT, ConstMatrixView::ncols(), and ConstMatrixView::nrows().

Referenced by get_ppath_transmat(), matrix_exp(), MatrixIdentity(), MCGeneral(), mcPathTraceGeneral(), mcPathTraceRadar(), rte_step_doit_replacement(), sensor_responseInit(), SparseIdentity(), test47(), test_identity(), test_inv(), test_inverse(), test_io(), and test_oem_gauss_newton_sparse().

◆ inv() [1/2]

void inv	(	ComplexMatrixView	Ainv,
		const ConstComplexMatrixView	A
	)

Definition at line 204 of file lin_alg.cc.

References ARTS_ASSERT, ARTS_USER_ERROR_IF, ConstComplexVectorView::get_c_array(), ConstComplexMatrixView::get_c_array(), ConstComplexMatrixView::ncols(), ConstComplexMatrixView::nrows(), lapack::zgetrf_(), and lapack::zgetri_().

◆ inv() [2/2]

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.

Parameters

[out]	Ainv	The MatrixView to contain the inverse of A.
[in]	A	The matrix to be inverted.

Definition at line 168 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(), test_addition(), test_complex_diagonalize(), test_inv(), test_oem_gauss_newton(), test_oem_levenberg_marquardt(), and test_real_diagonalize().

◆ linreg()

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

Parameters

p	Out: Fitted coefficients.
x	In: x-value of data points
y	In: y-value of data points

Author: Patrick Eriksson

Date: 2013-01-25

Definition at line 548 of file lin_alg.cc.

References ARTS_ASSERT, ConstVectorView::nelem(), and Vector::resize().

Referenced by derive_scat_species_a_and_b(), and ppathFromRtePos2().

◆ lsf()

Numeric lsf	(	VectorView	x,
		ConstMatrixView	A,
		ConstVectorView	y,
		bool	residual = `true`
	)

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>.

Parameters

[in]	x	As equation
[in]	A	As equation
[in]	y	As equation
[in]	residual	(optional) Returns the residual if true

Returns: Squared residual or 0

Definition at line 596 of file lin_alg.cc.

◆ lubacksub()

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.

Parameters

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 92 of file lin_alg.cc.

Referenced by matrix_exp(), solve(), test_lusolve1D(), test_lusolve4D(), and test_solve_linear_system().

◆ ludcmp()

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.)

Parameters

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 57 of file lin_alg.cc.

Referenced by matrix_exp(), solve(), test_lusolve1D(), test_lusolve4D(), and test_solve_linear_system().

◆ matrix_exp()

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.

Parameters

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 396 of file lin_alg.cc.

References ARTS_ASSERT, c, id_mat(), is_size(), joker, Constant::k, lubacksub(), ludcmp(), mult(), N, ConstMatrixView::ncols(), norm_inf(), pow(), and q.

Referenced by ext2trans(), test_matrix_exp1D(), test_matrix_exp3D(), test_matrix_exp4D(), and test_real_diagonalize().

◆ norm_inf()

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).

Parameters

A	Input: arbitrary matrix

Returns: Maximum absolute row sum norm

Definition at line 471 of file lin_alg.cc.

References abs, ConstMatrixView::ncols(), norm_inf(), and ConstMatrixView::nrows().

Referenced by matrix_exp(), and norm_inf().

◆ solve()

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.

Parameters

x	The solution vector x.
A	The matrix A defining the system.
b	The vector b.

Definition at line 139 of file lin_alg.cc.

References ARTS_ASSERT, b, lubacksub(), ludcmp(), ConstMatrixView::ncols(), ConstVectorView::nelem(), and ConstMatrixView::nrows().

Functions

Detailed Description

Function Documentation

◆ det()

◆ diagonalize() [1/2]

◆ diagonalize() [2/2]

◆ id_mat()

◆ inv() [1/2]

◆ inv() [2/2]

◆ linreg()

◆ lsf()

◆ lubacksub()

◆ ludcmp()

◆ matrix_exp()

◆ norm_inf()

◆ solve()