ARTS  2.3.1285(git:92a29ea9-dirty)
lin_alg.cc File Reference

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.

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...
 
void matrix_exp2 (MatrixView F, ConstMatrixView A)
 
void matrix_exp_4x4 (MatrixView arts_f, ConstMatrixView A, const Index &q)
 
void matrix_exp2_4x4 (MatrixView arts_f, ConstMatrixView A)
 
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. More...
 
Matrix4x4ViewMap MapToEigen4x4 (Tensor3View &A, const Index &i)
 
MatrixViewMap MapToEigen (Tensor3View &A, const Index &i)
 
void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__explicit (MatrixView F, ConstMatrixView A)
 
void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__eigen (MatrixView F, ConstMatrixView A)
 
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)
 
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)
 
void propmat4x4_to_transmat4x4 (MatrixView F, Tensor3View dF_upp, Tensor3View dF_low, ConstMatrixView A, ConstTensor3View dA_upp, ConstTensor3View dA_low, const Index &q)
 
void matrix_exp_dmatrix_exp (MatrixView F, Tensor3View dF, ConstMatrixView A, ConstTensor3View dA, const Index &q)
 General exponential of a Matrix with their derivatives. More...
 
void matrix_exp_dmatrix_exp (MatrixView F, MatrixView dF, ConstMatrixView A, ConstMatrixView dA, const Index &q)
 General exponential of a Matrix with their derivatives. 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...
 
Eigen::ComplexEigenSolver< Eigen::MatrixXcd > eig (const Eigen::Ref< Eigen::MatrixXcd > A)
 Return the Eigen decomposition of the eigen matrix. 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

◆ cayley_hamilton_fitted_method_4x4_propmat_to_transmat__eigen() [1/2]

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 |

| d -v -w a |

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

Parameters
FOutput: The matrix exponential of A (Has to be initialized before calling the function.
dFOutput: The derivative of the matrix exponential of A (Has to be initialized before calling the function. Page dimension is for each derivative.
AInput: arbitrary square matrix.
dAInput: derivative of A. Page dimension is for each derivative.

Definition at line 963 of file lin_alg.cc.

References MapToEigen4x4(), and sqrt().

Referenced by ext2trans().

◆ cayley_hamilton_fitted_method_4x4_propmat_to_transmat__eigen() [2/2]

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 |

| d -v -w a |

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 |

| dd -dv -dw da |

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

Parameters
FOutput: The matrix exponential of A (Has to be initialized before calling the function.
dFOutput: The derivative of the matrix exponential of A (Has to be initialized before calling the function. Page dimension is for each derivative.
AInput: arbitrary square matrix.
dAInput: derivative of A. Page dimension is for each derivative.

Definition at line 1550 of file lin_alg.cc.

References MapToEigen4x4(), ConstTensor3View::npages(), and sqrt().

◆ cayley_hamilton_fitted_method_4x4_propmat_to_transmat__explicit() [1/2]

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 |

| d -v -w a |

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

Parameters
FOutput: The matrix exponential of A (Has to be initialized before calling the function.
dFOutput: The derivative of the matrix exponential of A (Has to be initialized before calling the function. Page dimension is for each derivative.
AInput: arbitrary square matrix.
dAInput: derivative of A. Page dimension is for each derivative.

Definition at line 789 of file lin_alg.cc.

References sqrt().

◆ cayley_hamilton_fitted_method_4x4_propmat_to_transmat__explicit() [2/2]

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 |

| d -v -w a |

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 |

| dd -dv -dw da |

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

Parameters
FOutput: The matrix exponential of A (Has to be initialized before calling the function.
dFOutput: The derivative of the matrix exponential of A (Has to be initialized before calling the function. Page dimension is for each derivative.
AInput: arbitrary square matrix.
dAInput: derivative of A. Page dimension is for each derivative.

Definition at line 1111 of file lin_alg.cc.

References ConstTensor3View::npages(), and sqrt().

◆ det()

Determinant of N by N matrix. Simple recursive method.

Parameters
AIn: Matrix of size NxN.
Author
Richard Larsson
Date
2012-08-03

Definition at line 2256 of file lin_alg.cc.

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

Referenced by det().

◆ diagonalize() [1/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]PThe right eigenvectors.
[out]WRThe real eigenvalues.
[out]WIThe imaginary eigenvalues.
[in]AThe matrix to diagonalize.

Definition at line 245 of file lin_alg.cc.

References lapack::dgeev_(), Range::get_extent(), i, ConstMatrixView::mcr, ConstVectorView::mdata, ConstMatrixView::mdata, n, ConstMatrixView::ncols(), ConstVectorView::nelem(), and ConstMatrixView::nrows().

Referenced by matrix_exp2(), test_complex_diagonalize(), and test_real_diagonalize().

◆ diagonalize() [2/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]PThe right eigenvectors.
[out]WThe eigenvalues.
[in]AThe matrix to diagonalize.

Definition at line 330 of file lin_alg.cc.

References ComplexVectorView::get_c_array(), VectorView::get_c_array(), ComplexMatrixView::get_c_array(), ConstComplexMatrixView::get_column_extent(), i, n, ConstComplexMatrixView::ncols(), ConstComplexVectorView::nelem(), ConstComplexMatrixView::nrows(), swap(), and lapack::zgeev_().

◆ eig()

Eigen::ComplexEigenSolver<Eigen::MatrixXcd> eig ( const Eigen::Ref< Eigen::MatrixXcd >  A)

Return the Eigen decomposition of the eigen matrix.

Parameters
[in]Aa matrix to eigen value decompose
Returns
Object with eigenvalues and eigenvectors computed

Definition at line 2373 of file lin_alg.cc.

◆ id_mat()

◆ inv() [1/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]AinvThe MatrixView to contain the inverse of A.
[in]AThe matrix to be inverted.

Definition at line 167 of file lin_alg.cc.

References lapack::dgetrf_(), lapack::dgetri_(), ConstMatrixView::mdata, n, ConstMatrixView::ncols(), and ConstMatrixView::nrows().

Referenced by benchmark_inv(), benchmark_oem_linear(), CovarianceMatrix::invert_correlation_block(), matrix_exp2(), oem_gauss_newton_n_form(), oem_linear_mform(), test_complex_diagonalize(), test_inv(), test_oem_gauss_newton(), test_oem_levenberg_marquardt(), and test_real_diagonalize().

◆ inv() [2/2]

◆ 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
pOut: Fitted coefficients.
xIn: x-value of data points
yIn: y-value of data points
Author
Patrick Eriksson
Date
2013-01-25

Definition at line 2302 of file lin_alg.cc.

References i, n, 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]xAs equation
[in]AAs equation
[in]yAs equation
[in]residual(optional) Returns the residual if true
Returns
Squared residual or 0

Definition at line 2350 of file lin_alg.cc.

References mult(), n, r, solve(), and transpose().

◆ 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
xOutput: Solution vector of the equation system.
LUInput: LU decomposition of the matrix (output of function ludcp).
bInput: Right-hand-side vector of equation system.
indxInput: Pivoting information (output of function ludcp).

Definition at line 91 of file lin_alg.cc.

References DEBUG_ONLY, lapack::dgetrs_(), Range::get_stride(), i, is_size(), ConstMatrixView::mcr, ConstMatrixView::mdata, ConstVectorView::mrange, n, and ConstMatrixView::nrows().

Referenced by scale_rows(), solve(), 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
LUOutput: returns L and U in one matrix
indxOutput: Vector that records the row permutation.
AInput: Matrix for which the LU decomposition is performed

Definition at line 56 of file lin_alg.cc.

References lapack::dgetrf_(), i, is_size(), ConstMatrixView::mdata, n, ConstMatrixView::nrows(), and transpose().

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

◆ MapToEigen()

MatrixViewMap MapToEigen ( Tensor3View A,
const Index i 
)

Definition at line 750 of file lin_alg.cc.

References joker.

◆ MapToEigen4x4()

Matrix4x4ViewMap MapToEigen4x4 ( Tensor3View A,
const Index i 
)

Definition at line 745 of file lin_alg.cc.

References joker.

Referenced by cayley_hamilton_fitted_method_4x4_propmat_to_transmat__eigen(), and matrix_exp2_4x4().

◆ 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
FOutput: The matrix exponential of A (Has to be initialized before calling the function.
AInput: arbitrary square matrix
qInput: Parameter for the accuracy of the computation

Definition at line 391 of file lin_alg.cc.

References is_size(), n, and ConstMatrixView::ncols().

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

◆ matrix_exp2()

void matrix_exp2 ( MatrixView  F,
ConstMatrixView  A 
)

Definition at line 458 of file lin_alg.cc.

References diagonalize(), inv(), is_size(), n, and ConstMatrixView::nrows().

Referenced by test_real_diagonalize().

◆ matrix_exp2_4x4()

void matrix_exp2_4x4 ( MatrixView  arts_f,
ConstMatrixView  A 
)

Definition at line 521 of file lin_alg.cc.

References F, is_size(), MapToEigen4x4(), and q.

◆ matrix_exp_4x4()

void matrix_exp_4x4 ( MatrixView  arts_f,
ConstMatrixView  A,
const Index q 
)

Definition at line 481 of file lin_alg.cc.

References is_size(), and norm_inf().

◆ matrix_exp_dmatrix_exp() [1/2]

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.

Parameters
FOutput: The matrix exponential of A (Has to be initialized before calling the function.
dFOutput: The derivative of the matrix exponential of A (Has to be initialized before calling the function. Page dimension is for each derivative.
AInput: arbitrary square matrix.
dAInput: derivative of A. Page dimension is for each derivative.
qInput: Parameter for the accuracy of the computation, Matlab default is 6.

Definition at line 1926 of file lin_alg.cc.

References is_size(), joker, n, ConstMatrixView::ncols(), and ConstTensor3View::npages().

Referenced by test_matrixexp().

◆ matrix_exp_dmatrix_exp() [2/2]

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.

Parameters
FOutput: The matrix exponential of A (Has to be initialized before calling the function.
dFOutput: The derivative of the matrix exponential of A (Has to be initialized before calling the function.
AInput: arbitrary square matrix
dAInput: derivative of arbitrary square matrix
qInput: Parameter for the accuracy of the computation

Definition at line 2087 of file lin_alg.cc.

References is_size(), n, and ConstMatrixView::ncols().

◆ 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
AInput: arbitrary matrix
Returns
Maximum absolute row sum norm

Definition at line 2223 of file lin_alg.cc.

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

Referenced by matrix_exp_4x4(), norm_inf(), and propmat4x4_to_transmat4x4().

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

References is_size(), joker, norm_inf(), and ConstTensor3View::npages().

◆ 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
xThe solution vector x.
AThe matrix A defining the system.
bThe vector b.

Definition at line 138 of file lin_alg.cc.

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

Referenced by lsf().

◆ special_matrix_exp_and_dmatrix_exp_dx_for_rt()

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.

Parameters
FOutput: The matrix exponential of A (Has to be initialized before calling the function.
dFOutput: The derivative of the matrix exponential of A (Has to be initialized before calling the function. Page dimension is for each derivative.
AInput: arbitrary square matrix.
dAInput: derivative of A. Page dimension is for each derivative.
qInput: Parameter for the accuracy of the computation, Matlab default is 6.

Definition at line 559 of file lin_alg.cc.

References is_size(), joker, n, ConstMatrixView::ncols(), and ConstTensor3View::npages().