Name | Description |
---|---|
![]() | |
![]() | |
![]() | Solution of continuous-time algebraic Riccati equations |
![]() | Solution of discrete-time algebraic Riccati equations |
![]() | Determinant of a matrix (computed by LU decomposition) |
![]() | Compute eigenvalues and eigenvectors for a real, nonsymmetric matrix |
![]() | Solve a linear equality constrained least squares problem |
![]() | flip the columns of a matrix in left/right direction |
![]() | flip the columns of a matrix in up/down direction |
![]() | Return the Frobenius norm of a matrix |
![]() | Read matrix from a matlab file |
![]() | Compute invariant zeros of linear state space system with a genralized system matrix [A, B, C, D] which is of upper Hessenberg form |
![]() | Compute an upper Hessenberg matrix by repeatedly applicated householder similarity transformation |
![]() | reflect each of the vectors ai of matrix A=[a1, a2, ..., an] on a plane with orthogonal vector u |
![]() | Calculate the similarity transformation SAS of matrix A with householder matrix S = I - 2u*u' |
![]() | Solve overdetermined or underdetermined real system of linear equations A*x=b in a least squares sense (A may be rank deficient) |
![]() | LU decomposition of square or rectangular matrix |
![]() | Solve real system of linear equations P*L*U*x=b with a b vector and an LU decomposition (from LU(..)) |
![]() | Solve real system of linear equations P*L*U*X=B with a B vector and an LU decomposition (from LU(..)) |
![]() | Solution of continuous-time Lyapunov equation X*A + A'*X = C |
![]() | Returns the norm of a matrix |
![]() | generates a real orthogonal matrix Q which is defined as the product of IHI-ILO elementary reflectors of order N, as returned by DGEHRD |
![]() | print matrix |
![]() | QR decomposition of a square matrix without column pivoting (A = Q*R) |
![]() | Computes the real Schur form (RSF) of a square matrix |
![]() | Solve real system of linear equations A*x=b with a b vector (Gaussian elemination with partial pivoting) |
![]() | Solve real system of linear equations A*X=B with a B matrix (Gaussian elemination with partial pivoting) |
![]() | transform a real general matrix A to upper Hessenberg form H by an orthogonal similarity transformation: Q' * A * Q = H |
![]() | tarce(A) is the sum of the diagonal elements of A |
withQ + A'*X + X*A - X*G*X = 0
using the Schur vector approach proposed by Laub [1].-1 G = B*R *B'
It is assumed that Q is symmetric and positve semidefinite and R is symmetric, nonsingular and positive definite, (A,B) is stabilizable and (A,Q) is detectable.
The assumptions are not checked in this function
The assumptions guarantee that Hamiltonian matrix
has no pure imaginary eigenvalue and can be put to an ordered real Schur formH = [A, -G; -Q, -A']
with orthogonal similarity transformation U. S is ordered in such a way, that S11 contains the n stable eigenvalues of the closed loop system with system matrixU'*H*U = S = [S11, S12; 0, S22]
If U is partitioned to-1 A - B*R *B'*X
with dimenstions according to S, the solution X can be calculated byU = [U11, U12; U21, U22]
The algorithm uses LAPACK routines dgehrd (to compute the upper Hessenberg matrix of H), dorghr (to calculate the orthogonal matrix from the elementary reflectors as returned from dgehrd), dhseqr (to put transformed H to Schur form and to calculate the eigenvalues of the closed loop system) and dtrsen (to compute the ordered real Schur form and matrix U).X*U11 = U21.
[1] Laub, A.J. A Schur Method for Solving Algebraic Riccati equations. IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.
Name | Description |
---|---|
A[:, size(A, 1)] | |
B[size(A, 1), :] | |
R[size(B, 2), size(B, 2)] | |
Q[size(A, 1), size(A, 1)] | |
refine |
Name | Description |
---|---|
X[size(A, 1), size(A, 2)] | stabilizing solution of CARE |
ev[size(A, 1)] | eigenvalues of the closed loop system |
using the Schur vector approach proposed by Laub [1].-1 X = A'*X*A - A'*X*B*(R + B'*X*B) *B'*X*A + Q
It is assumed that Q is symmetric and positve semidefinite and R is symmetric, nonsingular and positive definite, (A,B) is stabilizable and (A,Q) is detectable.
The assumptions are not checked in this function
The assumptions guarantee that the Hamiltonian matrix
with-1 -1 -1 -1 H = [A, -A *G; Q*A, A' + Q*A *G ]
has no eigenvalue on the unit circle and can be put to an ordered real Schur form-1 G = B*R *B'
with orthogonal similarity transformation U. S is ordered in such a way, that S11 contains the n stable eigenvalues of the closed loop system with system matrixU'*H*U = S = [S11, S12; 0, S22]
If U is partitioned to-1 A - B*(R + B'*X*B) *B'*X*A
according to S, the solution X can be calculated byU = [U11, U12; U21, U22]
The algorithm uses LAPACK routines dgehrd (to compute the upper Hessenberg matrix of H), dorghr (to calculate the orthogonal matrix from the elementary reflectors as returned from dgehrd), dhseqr (to put transformed H to Schur form and to calculate the eigenvalues of the closed loop system) and dtrsen (to compute the ordered real Schur form and matrix U).X*U11 = U21.
[1] Laub, A.J. A Schur Method for Solving Algebraic Riccati equations. IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.
Name | Description |
---|---|
A[:, size(A, 1)] | |
B[size(A, 1), :] | |
R[size(B, 2), size(B, 2)] | |
Q[size(A, 1), size(A, 1)] |
Name | Description |
---|---|
S[size(A, 1), size(A, 2)] | orthogonal matrix of the Schur vectors associated to ordered rsf |
ev[size(A, 1)] | eigenvalues of the closed loop system |
Matrices.det(A);
This function call returns the determinant of matrix A computed by a LU decomposition. Usally, this function should never be used, because there are nearly always better numerical algorithms as by computing the determinant. E.g., use function Matrices.rank to compute the rank of a matrix.
Extends from Modelica.Icons.Function (Icon for a function).
Name | Description |
---|---|
A[:, size(A, 1)] |
Name | Description |
---|---|
result | Determinant of matrix A |
eigenvalues = Matrices.eigenValues(A); (eigenvalues, eigenvectors) = Matrices.eigenValues(A);
This function call returns the eigenvalues and optionally the (right) eigenvectors of a square matrix A. The first column of "eigenvalues" contains the real and the second column contains the imaginary part of the eigenvalues. If the i-th eigenvalue has no imaginary part, then eigenvectors[:,i] is the corresponding real eigenvector. If the i-th eigenvalue has an imaginary part, then eigenvalues[i+1,:] is the conjugate complex eigenvalue and eigenvectors[:,i] is the real and eigenvectors[:,i+1] is the imaginary part of the eigenvector of the i-th eigenvalue. With function Matrices.eigenValueMatrix, a real block diagonal matrix is constructed from the eigenvalues such that
A = eigenvectors * eigenValueMatrix(eigenvalues) * inv(eigenvectors)
provided the eigenvector matrix "eigenvectors" can be inverted (an inversion is possible, if all eigenvalues are different and no eigenvalue is zero).
Real A[3,3] = [1,2,3; 3,4,5; 2,1,4]; Real eval; algorithm eval := Matrices.eigenValues(A); // eval = [-0.618, 0; // 8.0 , 0; // 1.618, 0];
i.e., matrix A has the 3 real eigenvalues -0.618, 8, 1.618.
Extends from Modelica.Icons.Function (Icon for a function).
Name | Description |
---|---|
A[:, size(A, 1)] | Matrix |
Name | Description |
---|---|
eigenvalues[size(A, 1), 2] | Eigenvalues of matrix A (Re: first column, Im: second column) |
leftEigenvectors[size(A, 1), size(A, 2)] | Real-valued eigenvector matrix |
rightEigenvectors[size(A, 1), size(A, 2)] | Real-valued eigenvector matrix |
x = Matrices.equalityLeastSquares(A,a,B,b);
This function returns the solution x of the linear equality-constrained least squares problem:
min|A*x - a|^2 over x, subject to B*x = b
It is required that the dimensions of A and B fulfill the following relationship:
size(B,1) ≤ size(A,2) ≤ size(A,1) + size(B,1)
Note, the solution is computed with the LAPACK function "dgglse" using the generalized RQ factorization under the assumptions that B has full row rank (= size(B,1)) and the matrix [A;B] has full column rank (= size(A,2)). In this case, the problem has a unique solution.
Extends from Modelica.Icons.Function (Icon for a function).
Name | Description |
---|---|
A[:, :] | Minimize |A*x - a|^2 |
a[size(A, 1)] | |
B[:, size(A, 2)] | subject to B*x=b |
b[size(B, 1)] |
Name | Description |
---|---|
x[size(A, 2)] | solution vector |
Name | Description |
---|---|
A[:, :] | Matrix to be fliped |
Name | Description |
---|---|
Aflip[size(A, 1), size(A, 2)] | fliped matrix |
Name | Description |
---|---|
A[:, :] | Matrix to be fliped |
Name | Description |
---|---|
Aflip[size(A, 1), size(A, 2)] | fliped matrix |
Name | Description |
---|---|
A[:, :] | Input matrix |
Name | Description |
---|---|
result | frobenius norm of matrix A |
Name | Description |
---|---|
fileName | |
matrixName | Name of the matrix |
Name | Description |
---|---|
A[n, m] |
This function is an interface to LAPACK routine DHGEQZ to calculate invariant zeros of systems with generalized system matrices of upper Hessenberg form. DHGEQZ is described below: Purpose ========================================================== DHGEQZ implements a single-/double-shift version of the QZ method for finding the generalized eigenvalues w(j)=(ALPHAR(j) + i*ALPHAI(j))/BETAR(j) of the equation det( A - w(i) B ) = 0 In addition, the pair A,B may be reduced to generalized Schur form: B is upper triangular, and A is block upper triangular, where the diagonal blocks are either 1-by-1 or 2-by-2, the 2-by-2 blocks having complex generalized eigenvalues (see the description of the argument JOB.) If JOB='S', then the pair (A,B) is simultaneously reduced to Schur form by applying one orthogonal tranformation (usually called Q) on the left and another (usually called Z) on the right. The 2-by-2 upper-triangular diagonal blocks of B corresponding to 2-by-2 blocks of A will be reduced to positive diagonal matrices. (I.e., if A(j+1,j) is non-zero, then B(j+1,j)=B(j,j+1)=0 and B(j,j) and B(j+1,j+1) will be positive.) If JOB='E', then at each iteration, the same transformations are computed, but they are only applied to those parts of A and B which are needed to compute ALPHAR, ALPHAI, and BETAR. If JOB='S' and COMPQ and COMPZ are 'V' or 'I', then the orthogonal transformations used to reduce (A,B) are accumulated into the arrays Q and Z s.t.: Q(in) A(in) Z(in)* = Q(out) A(out) Z(out)* Q(in) B(in) Z(in)* = Q(out) B(out) Z(out)* Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973), pp. 241--256. Arguments ========= JOB (input) CHARACTER*1 = 'E': compute only ALPHAR, ALPHAI, and BETA. A and B will not necessarily be put into generalized Schur form. = 'S': put A and B into generalized Schur form, as well as computing ALPHAR, ALPHAI, and BETA. COMPQ (input) CHARACTER*1 = 'N': do not modify Q. = 'V': multiply the array Q on the right by the transpose of the orthogonal tranformation that is applied to the left side of A and B to reduce them to Schur form. = 'I': like COMPQ='V', except that Q will be initialized to the identity first. COMPZ (input) CHARACTER*1 = 'N': do not modify Z. = 'V': multiply the array Z on the right by the orthogonal tranformation that is applied to the right side of A and B to reduce them to Schur form. = 'I': like COMPZ='V', except that Z will be initialized to the identity first. N (input) INTEGER The order of the matrices A, B, Q, and Z. N >= 0. ILO (input) INTEGER IHI (input) INTEGER It is assumed that A is already upper triangular in rows and columns 1:ILO-1 and IHI+1:N. 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. A (input/output) DOUBLE PRECISION array, dimension (LDA, N) On entry, the N-by-N upper Hessenberg matrix A. Elements below the subdiagonal must be zero. If JOB='S', then on exit A and B will have been simultaneously reduced to generalized Schur form. If JOB='E', then on exit A will have been destroyed. The diagonal blocks will be correct, but the off-diagonal portion will be meaningless. LDA (input) INTEGER The leading dimension of the array A. LDA >= max( 1, N ). B (input/output) DOUBLE PRECISION array, dimension (LDB, N) On entry, the N-by-N upper triangular matrix B. Elements below the diagonal must be zero. 2-by-2 blocks in B corresponding to 2-by-2 blocks in A will be reduced to positive diagonal form. (I.e., if A(j+1,j) is non-zero, then B(j+1,j)=B(j,j+1)=0 and B(j,j) and B(j+1,j+1) will be positive.) If JOB='S', then on exit A and B will have been simultaneously reduced to Schur form. If JOB='E', then on exit B will have been destroyed. Elements corresponding to diagonal blocks of A will be correct, but the off-diagonal portion will be meaningless. LDB (input) INTEGER The leading dimension of the array B. LDB >= max( 1, N ). ALPHAR (output) DOUBLE PRECISION array, dimension (N) ALPHAR(1:N) will be set to real parts of the diagonal elements of A that would result from reducing A and B to Schur form and then further reducing them both to triangular form using unitary transformations s.t. the diagonal of B was non-negative real. Thus, if A(j,j) is in a 1-by-1 block (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=A(j,j). Note that the (real or complex) values (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the generalized eigenvalues of the matrix pencil A - wB. ALPHAI (output) DOUBLE PRECISION array, dimension (N) ALPHAI(1:N) will be set to imaginary parts of the diagonal elements of A that would result from reducing A and B to Schur form and then further reducing them both to triangular form using unitary transformations s.t. the diagonal of B was non-negative real. Thus, if A(j,j) is in a 1-by-1 block (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=0. Note that the (real or complex) values (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the generalized eigenvalues of the matrix pencil A - wB. BETA (output) DOUBLE PRECISION array, dimension (N) BETA(1:N) will be set to the (real) diagonal elements of B that would result from reducing A and B to Schur form and then further reducing them both to triangular form using unitary transformations s.t. the diagonal of B was non-negative real. Thus, if A(j,j) is in a 1-by-1 block (i.e., A(j+1,j)=A(j,j+1)=0), then BETA(j)=B(j,j). Note that the (real or complex) values (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the generalized eigenvalues of the matrix pencil A - wB. (Note that BETA(1:N) will always be non-negative, and no BETAI is necessary.) Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N) If COMPQ='N', then Q will not be referenced. If COMPQ='V' or 'I', then the transpose of the orthogonal transformations which are applied to A and B on the left will be applied to the array Q on the right. LDQ (input) INTEGER The leading dimension of the array Q. LDQ >= 1. If COMPQ='V' or 'I', then LDQ >= N. Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N) If COMPZ='N', then Z will not be referenced. If COMPZ='V' or 'I', then the orthogonal transformations which are applied to A and B on the right will be applied to the array Z on the right. LDZ (input) INTEGER The leading dimension of the array Z. LDZ >= 1. If COMPZ='V' or 'I', then LDZ >= N. WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) On exit, if INFO >= 0, WORK(1) returns the optimal LWORK. LWORK (input) INTEGER The dimension of the array WORK. LWORK >= max(1,N). If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA. INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value = 1,...,N: the QZ iteration did not converge. (A,B) is not in Schur form, but ALPHAR(i), ALPHAI(i), and BETA(i), i=INFO+1,...,N should be correct. = N+1,...,2*N: the shift calculation failed. (A,B) is not in Schur form, but ALPHAR(i), ALPHAI(i), and BETA(i), i=INFO-N+1,...,N should be correct. > 2*N: various "impossible" errors. Further Details =============== Iteration counters: JITER -- counts iterations. IITER -- counts iterations run since ILAST was last changed. This is therefore reset only when a 1-by-1 or 2-by-2 block deflates off the bottom. =====================================================================
Name | Description |
---|---|
A[:, size(A, 1)] | |
B[size(A, 1), size(A, 1)] |
Name | Description |
---|---|
alphaReal[size(A, 1)] | Real part of alpha (eigenvalue=(alphaReal+i*alphaImag)/beta) |
alphaImag[size(A, 1)] | Imaginary part of alpha |
beta[size(A, 1)] | Denominator of eigenvalue |
info |
Ai+1 = (I-2*u_i*u_i')*Ai*(I-2*u_i*u_i')with Householder vector u_i.
The elementary transformations can be subsumed under
A -> Q*A*Qand Q*A*Q is Hessenberg matrix.
Modelica_LinearSystems2.Math.Matrices.hess uses LAPACK routine dgehrd. In contrast to this function Modelica_LinearSystems2.Math.Matrices.Internal.hessenberg does not use any LAPACK routine.
Name | Description |
---|---|
H[size(H, 1), size(H, 2)] |
Name | Description |
---|---|
Ht[size(H, 1), size(H, 2)] |
Name | Description |
---|---|
A[:, :] | |
u[size(A, 1)] | householder vector |
Name | Description |
---|---|
RA[size(A, 1), size(A, 2)] | reflexion of A |
Name | Description |
---|---|
A[:, size(A, 1)] | |
u[size(A, 1)] | householder vector |
Name | Description |
---|---|
SAS[size(A, 1), size(A, 1)] |
x = Matrices.leastSquares(A,b);
A linear system of equations A*x = b has no solutions or infinitely many solutions if A is not square. Function "leastSquares" returns a solution in a least squarse sense:
size(A,1) > size(A,2): returns x such that |A*x - b|^2 is a minimum size(A,1) = size(A,2): returns x such that A*x = b size(A,1) < size(A,2): returns x such that |x|^2 is a minimum for all vectors x that fulfill A*x = b
Note, the solution is computed with the LAPACK function "dgelsx", i.e., QR or LQ factorization of A with column pivoting. If A does not have full rank, the solution is not unique and from the infinitely many solutions the one is selected that minimizes both |x|^2 and |A*x - b|^2.
Extends from Modelica.Icons.Function (Icon for a function).
Name | Description |
---|---|
A[:, :] | Matrix A |
b[size(A, 1)] | Vector b |
Name | Description |
---|---|
x[size(A, 2)] | Vector x such that min|A*x-b|^2 if size(A,1) >= size(A,2) or min|x|^2 and A*x=b, if size(A,1) < size(A,2) |
(LU, pivots) = Matrices.LU(A); (LU, pivots, info) = Matrices.LU(A);
This function call returns the LU decomposition of a "Real[m,n]" matrix A, i.e.,
P*L*U = A
where P is a permutation matrix (implicitely
defined by vector pivots
),
L is a lower triangular matrix with unit
diagonal elements (lower trapezoidal if m > n), and
U is an upper triangular matrix (upper trapezoidal if m < n).
Matrices L and U are stored in the returned
matrix LU
(the diagonal of L is not stored).
With the companion function
Matrices.LU_solve,
this decomposition can be used to solve
linear systems (P*L*U)*x = b with different right
hand side vectors b. If a linear system of equations with
just one right hand side vector b shall be solved, it is
more convenient to just use the function
Matrices.solve.
The optional third (Integer) output argument has the following meaning:
info = 0: | successful exit|
info > 0: | if info = i, U[i,i] is exactly zero. The factorization
has been completed, but the factor U is exactly singular, and division by zero will occur if it is used to solve a system of equations. |
The LU factorization is computed with the LAPACK function "dgetrf", i.e., by Gaussian elemination using partial pivoting with row interchanges. Vector "pivots" are the pivot indices, i.e., for 1 ≤ i ≤ min(m,n), row i of matrix A was interchanged with row pivots[i].
Real A[3,3] = [1,2,3; 3,4,5; 2,1,4]; Real b1[3] = {10,22,12}; Real b2[3] = { 7,13,10}; Real LU[3,3]; Integer pivots[3]; Real x1[3]; Real x2[3]; algorithm (LU, pivots) := Matrices.LU(A); x1 := Matrices.LU_solve(LU, pivots, b1); // x1 = {3,2,1} x2 := Matrices.LU_solve(LU, pivots, b2); // x2 = {1,0,2}
Extends from Modelica.Icons.Function (Icon for a function).
Name | Description |
---|---|
A[:, :] | Square or rectangular matrix |
Name | Description |
---|---|
LU[size(A, 1), size(A, 2)] | L,U factors (used with LU_solve(..)) |
pivots[min(size(A, 1), size(A, 2))] | pivot indices (used with LU_solve(..)) |
info | Information |
Matrices.LU_solve(LU, pivots, b);
This function call returns the solution x of the linear systems of equations
P*L*U*x = b;
where P is a permutation matrix (implicitely
defined by vector pivots
),
L is a lower triangular matrix with unit
diagonal elements (lower trapezoidal if m > n), and
U is an upper triangular matrix (upper trapezoidal if m < n).
The matrices of this decomposition are computed with function
Matrices.LU that
returns arguments LU
and pivots
used as input arguments of Matrices.LU_solve
.
With Matrices.LU
and Matrices.LU_solve
it is possible to efficiently solve linear systems
with different right hand side vectors. If a linear system of equations with
just one right hand side vector shall be solved, it is
more convenient to just use the function
Matrices.solve.
If a unique solution x does not exist (since the LU decomposition is singular), an exception is raised.
The LU factorization is computed with the LAPACK function "dgetrf", i.e., by Gaussian elemination using partial pivoting with row interchanges. Vector "pivots" are the pivot indices, i.e., for 1 ≤ i ≤ min(m,n), row i of matrix A was interchanged with row pivots[i].
Real A[3,3] = [1,2,3; 3,4,5; 2,1,4]; Real b1[3] = {10,22,12}; Real b2[3] = { 7,13,10}; Real LU[3,3]; Integer pivots[3]; Real x1[3]; Real x2[3]; algorithm (LU, pivots) := Matrices.LU(A); x1 := Matrices.LU_solve(LU, pivots, b1); // x1 = {3,2,1} x2 := Matrices.LU_solve(LU, pivots, b2); // x2 = {1,0,2}
Extends from Modelica.Icons.Function (Icon for a function).
Name | Description |
---|---|
LU[:, size(LU, 1)] | L,U factors of Matrices.LU(..) for a square matrix |
pivots[size(LU, 1)] | Pivots indices of Matrices.LU(..) |
b[size(LU, 1)] | Right hand side vector of P*L*U*x=b |
Name | Description |
---|---|
x[size(b, 1)] | Solution vector such that P*L*U*x = b |
Matrices.LU_solve(LU, pivots, B);
This function call returns the solution x of the linear systems of equations
P*L*U*X = B;
where P is a permutation matrix (implicitely
defined by vector pivots
),
L is a lower triangular matrix with unit
diagonal elements (lower trapezoidal if m > n), and
U is an upper triangular matrix (upper trapezoidal if m < n).
The matrices of this decomposition are computed with function
Matrices.LU that
returns arguments LU
and pivots
used as input arguments of Matrices.LU_solve
.
With Matrices.LU
and Matrices.LU_solve
it is possible to efficiently solve linear systems
with different right hand side matrices. If a linear system of equations with
just one right hand side matrix shall be solved, it is
more convenient to just use the function
Matrices.solve2.
If a unique solution X does not exist (since the LU decomposition is singular), an exception is raised.
The LU factorization is computed with the LAPACK function "dgetrf", i.e., by Gaussian elemination using partial pivoting with row interchanges. Vector "pivots" are the pivot indices, i.e., for 1 ≤ i ≤ min(m,n), row i of matrix A was interchanged with row pivots[i].
Real A[3,3] = [1,2,3; 3,4,5; 2,1,4]; Real B1[3] = [10, 20; 22, 44; 12, 24]; Real B2[3] = [ 7, 14; 13, 26; 10, 20]; Real LU[3,3]; Integer pivots[3]; Real X1[3,2]; Real X2[3,2]; algorithm (LU, pivots) := Matrices.LU(A); X1 := Matrices.LU_solve2(LU, pivots, B1); /* X1 = [3, 6; 2, 4; 1, 2] */ X2 := Matrices.LU_solve2(LU, pivots, B2); /* X2 = [1, 2; 0, 0; 2, 4] */
Name | Description |
---|---|
LU[:, size(LU, 1)] | L,U factors of Matrices.LU(..) for a square matrix |
pivots[size(LU, 1)] | Pivots indices of Matrices.LU(..) |
B[size(LU, 1), :] | Right hand side matrix of P*L*U*X=B |
Name | Description |
---|---|
X[size(B, 1), size(B, 2)] | Solution matrix such that P*L*U*X = B |
using the Schur method for Lyapunov equations proposed by Bartels and Stewart [1].XA + A'*X = C.
[1] Bartels, R.H. and Stewart G.W. Algorithm 432: Solution of the matrix equation AX + XB = C. Comm. ACM., Vol. 15, pp. 820-826, 1972.
Name | Description |
---|---|
A[:, size(A, 1)] | |
C[size(A, 1), size(A, 2)] | |
eps |
Name | Description |
---|---|
X[size(A, 1), size(A, 2)] | solution of the Riccati equation |
Matrices.norm(A); Matrices.norm(A, p=2);
The function call "Matrices.norm(A)
" returns the
2-norm of matrix A, i.e., the largest singular value of A.
The function call "Matrices.norm(A, p)
" returns the
p-norm of matrix A. The only allowed values for p are
Note, for any matrices A1, A2 the following inequality holds:
Matrices.norm(A1+A2,p) ≤ Matrices.norm(A1,p) + Matrices.norm(A2,p)
Note, for any matrix A and vector v the following inequality holds:
Vectors.norm(A*v,p) ≤ Matrices.norm(A,p)*Vectors.norm(A,p)
Extends from Modelica.Icons.Function (Icon for a function).
Name | Description |
---|---|
A[:, :] | Input matrix |
p | Type of p-norm (only allowed: 1, 2 or Modelica.Constants.inf) |
Name | Description |
---|---|
result | p-norm of matrix A |
Purpose /* ======= */ /* DORGHR generates a real orthogonal matrix Q which is defined as the */ /* product of IHI-ILO elementary reflectors of order N, as returned by */ /* DGEHRD: */ /* Q = H(ilo) H(ilo+1) . . . H(ihi-1). */ /* Arguments */ /* ========= */ /* N (input) INTEGER */ /* The order of the matrix Q. N >= 0. */ /* ILO (input) INTEGER */ /* IHI (input) INTEGER */ /* ILO and IHI must have the same values as in the previous call */ /* of DGEHRD. Q is equal to the unit matrix except in the */ /* submatrix Q(ilo+1:ihi,ilo+1:ihi). */ /* 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. */ /* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */ /* On entry, the vectors which define the elementary reflectors, */ /* as returned by DGEHRD. */ /* On exit, the N-by-N orthogonal matrix Q. */ /* LDA (input) INTEGER */ /* The leading dimension of the array A. LDA >= max(1,N). */ /* TAU (input) DOUBLE PRECISION array, dimension (N-1) */ /* TAU(i) must contain the scalar factor of the elementary */ /* reflector H(i), as returned by DGEHRD. */ /* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) */ /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */ /* LWORK (input) INTEGER */ /* The dimension of the array WORK. LWORK >= IHI-ILO. */ /* For optimum performance LWORK >= (IHI-ILO)*NB, where NB is */ /* the optimal blocksize. */ /* If LWORK = -1, then a workspace query is assumed; the routine */ /* only calculates the optimal size of the WORK array, returns */ /* this value as the first entry of the WORK array, and no error */ /* message related to LWORK is issued by XERBLA. */ /* INFO (output) INTEGER */ /* = 0: successful exit */ /* < 0: if INFO = -i, the i-th argument had an illegal value */ /* ===================================================================== */
Name | Description |
---|---|
A[:, size(A, 1)] | |
tau[size(A, 1) - 1] | scalar factors of the elementary reflectors |
ilo | lowest index where the original matrix had been Hessenbergform - ilo must have the same value as in the previous call of DGEHRD |
ihi | highest index where the original matrix had been Hessenbergform - ihi must have the same value as in the previous call of DGEHRD |
Name | Description |
---|---|
Q[size(A, 1), size(A, 2)] | Orthogonal matrix as a result of elementary reflectors |
info |
Name | Description |
---|---|
M[:, :] | |
significantDigits | Number of significant digits that are shown |
name | Independent variable name used for printing |
Name | Description |
---|---|
s |
Name | Description |
---|---|
A[:, :] | Rectangular matrix with size(A,1) >= size(A,2) |
Name | Description |
---|---|
Q[size(A, 1), size(A, 2)] | Rectangular matrix with orthonormal columns such that Q*R=A[:,p] |
R[min(size(A, 1), size(A, 2)), size(A, 2)] | Square upper triangular matrix |
tau[min(size(A, 1), size(A, 2))] | |
Q2[:, :] |
Name | Description |
---|---|
A[:, size(A, 1)] |
Name | Description |
---|---|
S[size(A, 1), size(A, 2)] | |
QZ[size(A, 1), size(A, 2)] | |
alphaReal[size(A, 1)] | Real part of eigenvalue=alphaReal+i*alphaImag |
alphaImag[size(A, 1)] | Imaginary part of eigenvalue=(alphaReal+i*alphaImag |
Matrices.solve(A,b);
This function call returns the solution x of the linear system of equations
A*x = b
If a unique solution x does not exist (since A is singular), an exception is raised.
Note, the solution is computed with the LAPACK function "dgesv", i.e., by Gaussian elemination with partial pivoting.
Real A[3,3] = [1,2,3; 3,4,5; 2,1,4]; Real b[3] = {10,22,12}; Real x[3]; algorithm x := Matrices.solve(A,b); // x = {3,2,1}
Extends from Modelica.Icons.Function (Icon for a function).
Name | Description |
---|---|
A[:, size(A, 1)] | Matrix A of A*x = b |
b[size(A, 1)] | Vector b of A*x = b |
Name | Description |
---|---|
x[size(b, 1)] | Vector x such that A*x = b |
Matrices.solve2(A,b);
This function call returns the solution X of the linear system of equations
A*X = B
If a unique solution X does not exist (since A is singular), an exception is raised.
Note, the solution is computed with the LAPACK function "dgesv", i.e., by Gaussian elemination with partial pivoting.
Real A[3,3] = [1,2,3; 3,4,5; 2,1,4]; Real B[3,2] = [10, 20; 22, 44; 12, 24]; Real X[3,2]; algorithm (LU, pivots) := Matrices.LU(A); X := Matrices.solve2(A, B1); /* X = [3, 6; 2, 4; 1, 2] */
Extends from Modelica.Icons.Function (Icon for a function).
Name | Description |
---|---|
A[:, size(A, 1)] | Matrix A of A*X = B |
B[size(A, 1), :] | Matrix B of A*X = B |
Name | Description |
---|---|
X[size(B, 1), size(B, 2)] | Matrix X such that A*X = B |
Inputs are
Outputs are
===================================================================== -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University June 30, 1999 Purpose ======= DGEHRD reduces a real general matrix A to upper Hessenberg form H by an orthogonal similarity transformation: Q' * A * Q = H . Arguments ========= N (input) INTEGER The order of the matrix A. N >= 0. ILO (input) INTEGER IHI (input) INTEGER It is assumed that A is already upper triangular in rows and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally set by a previous call to DGEBAL; otherwise they should be set to 1 and N respectively. See Further Details. 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the N-by-N general matrix to be reduced. On exit, the upper triangle and the first subdiagonal of A are overwritten with the upper Hessenberg matrix H, and the elements below the first subdiagonal, with the array TAU, represent the orthogonal matrix Q as a product of elementary reflectors. See Further Details. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,N). TAU (output) DOUBLE PRECISION array, dimension (N-1) The scalar factors of the elementary reflectors (see Further Details). Elements 1:ILO-1 and IHI:N-1 of TAU are set to zero. WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) INTEGER The length of the array WORK. LWORK >= max(1,N). For optimum performance LWORK >= N*NB, where NB is the optimal blocksize. If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA. INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value. Further Details =============== The matrix Q is represented as a product of (ihi-ilo) elementary reflectors Q = H(ilo) H(ilo+1) . . . H(ihi-1). Each H(i) has the form H(i) = I - tau * v * v' where tau is a real scalar, and v is a real vector with v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on exit in A(i+2:ihi,i), and tau in TAU(i). The contents of A are illustrated by the following example, with n = 7, ilo = 2 and ihi = 6: on entry, on exit, ( a a a a a a a ) ( a a h h h h a ) ( a a a a a a ) ( a h h h h a ) ( a a a a a a ) ( h h h h h h ) ( a a a a a a ) ( v2 h h h h h ) ( a a a a a a ) ( v2 v3 h h h h ) ( a a a a a a ) ( v2 v3 v4 h h h ) ( a ) ( a ) where a denotes an element of the original matrix A, h denotes a modified element of the upper Hessenberg matrix H, and vi denotes an element of the vector defining H(i). =====================================================================
Name | Description |
---|---|
A[:, size(A, 1)] | |
ilo | lowest index where the original matrix had been Hessenbergform |
ihi |
Name | Description |
---|---|
H[size(A, 1), size(A, 2)] | Upper Hessenberg form |
V[size(A, 1), size(A, 2)] | V=[v1,v2,..vn-1,0] with vi are vectors which define the elementary reflectors |
tau[size(A, 1) - 1] | |
info |
Name | Description |
---|---|
A[:, size(A, 1)] |
Name | Description |
---|---|
result |