Interface to LAPACK library (should usually not directly be used but only indirectly via Modelica.Math.Matrices)
This package contains external Modelica functions as interface to the LAPACK library (http://www.netlib.org/lapack) that provides FORTRAN subroutines to solve linear algebra tasks. Usually, these functions are not directly called, but only via the much more convenient interface of Modelica.Math.Matrices. The documentation of the LAPACK functions is a copy of the original FORTRAN code. The details of LAPACK are described in:
See also http://en.wikipedia.org/wiki/Lapack.
This package contains a direct interface to the LAPACK subroutines
Extends from Modelica.Icons.Package (Icon for standard packages).
Name | Description |
---|---|
dgeev | Compute eigenvalues and (right) eigenvectors for real nonsymmetric matrix A |
dgeev_eigenValues | Compute eigenvalues for real nonsymmetric matrix A |
dgegv | Compute generalized eigenvalues for a (A,B) system |
dgelsx | Computes the minimum-norm solution to a real linear least squares problem with rank deficient A |
dgelsx_vec | Computes the minimum-norm solution to a real linear least squares problem with rank deficient A |
dgelsy | Computes the minimum-norm solution to a real linear least squares problem with rank deficient A |
dgelsy_vec | Computes the minimum-norm solution to a real linear least squares problem with rank deficient A |
dgels_vec | Solves overdetermined or underdetermined real linear equations A*x=b with a b vector |
dgesv | Solve real system of linear equations A*X=B with a B matrix |
dgesv_vec | Solve real system of linear equations A*x=b with a b vector |
dgglse_vec | Solve a linear equality constrained least squares problem |
dgtsv | Solve real system of linear equations A*X=B with B matrix and tridiagonal A |
dgtsv_vec | Solve real system of linear equations A*x=b with b vector and tridiagonal A |
dgbsv | Solve real system of linear equations A*X=B with a B matrix |
dgbsv_vec | Solve real system of linear equations A*x=b with a b vector |
dgesvd | Determine singular value decomposition |
dgesvd_sigma | Determine singular values |
dgetrf | Compute LU factorization of square or rectangular matrix A (A = P*L*U) |
dgetrs | Solves a system of linear equations with the LU decomposition from dgetrf(..) |
dgetrs_vec | Solves a system of linear equations with the LU decomposition from dgetrf(..) |
dgetri | Computes the inverse of a matrix using the LU factorization from dgetrf(..) |
dgeqpf | Compute QR factorization of square or rectangular matrix A with column pivoting (A(:,p) = Q*R) |
dorgqr | Generates a Real orthogonal matrix Q which is defined as the product of elementary reflectors as returned from dgeqpf |
dgees | Computes real Schur form T of real nonsymmetric matrix A, and, optionally, the matrix of Schur vectors Z as well as the eigenvalues |
dtrsen | Reorder the real Schur factorization of a real matrix |
dgesvx | Solve real system of linear equations op(A)*X=B, op(A) is A or A' according to the Boolean input transposed |
dtrsyl | Solve the real Sylvester matrix equation op(A)*X + X*op(B) = scale*C or op(A)*X - X*op(B) = scale*C |
dhseqr | Compute eigenvalues of a matrix H using lapack routine DHSEQR for Hessenberg form matrix |
dlange | Norm of a matrix |
dgecon | Estimates the reciprocal of the condition number of a general real matrix A |
dgehrd | reduces a real general matrix A to upper Hessenberg form H by an orthogonal similarity transformation: Q' * A * Q = H |
dgeqrf | computes a QR factorization without pivoting |
dgeevx | Compute the eigenvalues and the (real) left and right eigenvectors of matrix A, using lapack routine dgeevx |
dgesdd | Determine singular value decomposition |
dggev | Compute generalized eigenvalues, as well as the left and right eigenvectors for a (A,B) system |
dggevx | Compute generalized eigenvalues for a (A,B) system, using lapack routine dggevx |
dhgeqz | Compute generalized eigenvalues for a (A,B) system |
dormhr | overwrites the general real M-by-N matrix C with Q * C or C * Q or Q' * C or C * Q', where Q is an orthogonal matrix as returned by dgehrd |
dormqr | overwrites the general real M-by-N matrix C with Q * C or C * Q or Q' * C or C * Q', where Q is an orthogonal matrix of a QR factorization as returned by dgeqrf |
dtrevc | Compute the right and/or left eigenvectors of a real upper quasi-triangular matrix T |
dpotrf | Computes the Cholesky factorization of a real symmetric positive definite matrix A |
dtrsm | Solve one of the matrix equations op( A )*X = alpha*B, or X*op( A ) = alpha*B, where A is triangular matrix. BLAS routine |
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 |
Compute eigenvalues and (right) eigenvectors for real nonsymmetric matrix A
This function is not a full interface to the LAPACK function DGEEV, but calls it in such a way that only eigenvalues and right eigenvectors are computed. Lapack documentation Purpose ======= DGEEV computes for an N-by-N real nonsymmetric matrix A, the eigenvalues and, optionally, the left and/or right eigenvectors. The right eigenvector v(j) of A satisfies A * v(j) = lambda(j) * v(j) where lambda(j) is its eigenvalue. The left eigenvector u(j) of A satisfies u(j)**H * A = lambda(j) * u(j)**H where u(j)**H denotes the conjugate transpose of u(j). The computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real. Arguments ========= JOBVL (input) CHARACTER*1 = 'N': left eigenvectors of A are not computed; = 'V': left eigenvectors of A are computed. JOBVR (input) CHARACTER*1 = 'N': right eigenvectors of A are not computed; = 'V': right eigenvectors of A are computed. N (input) INTEGER The order of the matrix A. N >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the N-by-N matrix A. On exit, A has been overwritten. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,N). WR (output) DOUBLE PRECISION array, dimension (N) WI (output) DOUBLE PRECISION array, dimension (N) WR and WI contain the real and imaginary parts, respectively, of the computed eigenvalues. Complex conjugate pairs of eigenvalues appear consecutively with the eigenvalue having the positive imaginary part first. VL (output) DOUBLE PRECISION array, dimension (LDVL,N) If JOBVL = 'V', the left eigenvectors u(j) are stored one after another in the columns of VL, in the same order as their eigenvalues. If JOBVL = 'N', VL is not referenced. If the j-th eigenvalue is real, then u(j) = VL(:,j), the j-th column of VL. If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and u(j+1) = VL(:,j) - i*VL(:,j+1). LDVL (input) INTEGER The leading dimension of the array VL. LDVL >= 1; if JOBVL = 'V', LDVL >= N. VR (output) DOUBLE PRECISION array, dimension (LDVR,N) If JOBVR = 'V', the right eigenvectors v(j) are stored one after another in the columns of VR, in the same order as their eigenvalues. If JOBVR = 'N', VR is not referenced. If the j-th eigenvalue is real, then v(j) = VR(:,j), the j-th column of VR. If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and v(j+1) = VR(:,j) - i*VR(:,j+1). LDVR (input) INTEGER The leading dimension of the array VR. LDVR >= 1; if JOBVR = 'V', LDVR >= N. WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) INTEGER The dimension of the array WORK. LWORK >= max(1,3*N), and if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N. For good performance, LWORK must generally be larger. 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. > 0: if INFO = i, the QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed; elements i+1:N of WR and WI contain eigenvalues which have converged.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, size(A, 1)] |
Name | Description |
---|---|
eigenReal[size(A, 1)] | Real part of eigen values |
eigenImag[size(A, 1)] | Imaginary part of eigen values |
eigenVectors[size(A, 1), size(A, 1)] | Right eigen vectors |
info |
Compute eigenvalues for real nonsymmetric matrix A
Lapack documentation Purpose ======= DGEEV computes for an N-by-N real nonsymmetric matrix A, the eigenvalues and, optionally, the left and/or right eigenvectors. The right eigenvector v(j) of A satisfies A * v(j) = lambda(j) * v(j) where lambda(j) is its eigenvalue. The left eigenvector u(j) of A satisfies u(j)**H * A = lambda(j) * u(j)**H where u(j)**H denotes the conjugate transpose of u(j). The computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real. Arguments ========= JOBVL (input) CHARACTER*1 = 'N': left eigenvectors of A are not computed; = 'V': left eigenvectors of A are computed. JOBVR (input) CHARACTER*1 = 'N': right eigenvectors of A are not computed; = 'V': right eigenvectors of A are computed. N (input) INTEGER The order of the matrix A. N >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the N-by-N matrix A. On exit, A has been overwritten. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,N). WR (output) DOUBLE PRECISION array, dimension (N) WI (output) DOUBLE PRECISION array, dimension (N) WR and WI contain the real and imaginary parts, respectively, of the computed eigenvalues. Complex conjugate pairs of eigenvalues appear consecutively with the eigenvalue having the positive imaginary part first. VL (output) DOUBLE PRECISION array, dimension (LDVL,N) If JOBVL = 'V', the left eigenvectors u(j) are stored one after another in the columns of VL, in the same order as their eigenvalues. If JOBVL = 'N', VL is not referenced. If the j-th eigenvalue is real, then u(j) = VL(:,j), the j-th column of VL. If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and u(j+1) = VL(:,j) - i*VL(:,j+1). LDVL (input) INTEGER The leading dimension of the array VL. LDVL >= 1; if JOBVL = 'V', LDVL >= N. VR (output) DOUBLE PRECISION array, dimension (LDVR,N) If JOBVR = 'V', the right eigenvectors v(j) are stored one after another in the columns of VR, in the same order as their eigenvalues. If JOBVR = 'N', VR is not referenced. If the j-th eigenvalue is real, then v(j) = VR(:,j), the j-th column of VR. If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and v(j+1) = VR(:,j) - i*VR(:,j+1). LDVR (input) INTEGER The leading dimension of the array VR. LDVR >= 1; if JOBVR = 'V', LDVR >= N. WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) INTEGER The dimension of the array WORK. LWORK >= max(1,3*N), and if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N. For good performance, LWORK must generally be larger. 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. > 0: if INFO = i, the QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed; elements i+1:N of WR and WI contain eigenvalues which have converged.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, size(A, 1)] |
Name | Description |
---|---|
EigenReal[size(A, 1)] | |
EigenImag[size(A, 1)] | |
info |
Compute generalized eigenvalues for a (A,B) system
Lapack documentation Purpose ======= This routine is deprecated and has been replaced by routine DGGEV. DGEGV computes the eigenvalues and, optionally, the left and/or right eigenvectors of a real matrix pair (A,B). Given two square matrices A and B, the generalized nonsymmetric eigenvalue problem (GNEP) is to find the eigenvalues lambda and corresponding (non-zero) eigenvectors x such that A*x = lambda*B*x. An alternate form is to find the eigenvalues mu and corresponding eigenvectors y such that mu*A*y = B*y. These two forms are equivalent with mu = 1/lambda and x = y if neither lambda nor mu is zero. In order to deal with the case that lambda or mu is zero or small, two values alpha and beta are returned for each eigenvalue, such that lambda = alpha/beta and mu = beta/alpha. The vectors x and y in the above equations are right eigenvectors of the matrix pair (A,B). Vectors u and v satisfying u**H*A = lambda*u**H*B or mu*v**H*A = v**H*B are left eigenvectors of (A,B). Note: this routine performs "full balancing" on A and B -- see "Further Details", below. Arguments ========= JOBVL (input) CHARACTER*1 = 'N': do not compute the left generalized eigenvectors; = 'V': compute the left generalized eigenvectors (returned in VL). JOBVR (input) CHARACTER*1 = 'N': do not compute the right generalized eigenvectors; = 'V': compute the right generalized eigenvectors (returned in VR). N (input) INTEGER The order of the matrices A, B, VL, and VR. N >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA, N) On entry, the matrix A. If JOBVL = 'V' or JOBVR = 'V', then on exit A contains the real Schur form of A from the generalized Schur factorization of the pair (A,B) after balancing. If no eigenvectors were computed, then only the diagonal blocks from the Schur form will be correct. See DGGHRD and DHGEQZ for details. LDA (input) INTEGER The leading dimension of A. LDA >= max(1,N). B (input/output) DOUBLE PRECISION array, dimension (LDB, N) On entry, the matrix B. If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the upper triangular matrix obtained from B in the generalized Schur factorization of the pair (A,B) after balancing. If no eigenvectors were computed, then only those elements of B corresponding to the diagonal blocks from the Schur form of A will be correct. See DGGHRD and DHGEQZ for details. LDB (input) INTEGER The leading dimension of B. LDB >= max(1,N). ALPHAR (output) DOUBLE PRECISION array, dimension (N) The real parts of each scalar alpha defining an eigenvalue of GNEP. ALPHAI (output) DOUBLE PRECISION array, dimension (N) The imaginary parts of each scalar alpha defining an eigenvalue of GNEP. If ALPHAI(j) is zero, then the j-th eigenvalue is real; if positive, then the j-th and (j+1)-st eigenvalues are a complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j). BETA (output) DOUBLE PRECISION array, dimension (N) The scalars beta that define the eigenvalues of GNEP. Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and beta = BETA(j) represent the j-th eigenvalue of the matrix pair (A,B), in one of the forms lambda = alpha/beta or mu = beta/alpha. Since either lambda or mu may overflow, they should not, in general, be computed. VL (output) DOUBLE PRECISION array, dimension (LDVL,N) If JOBVL = 'V', the left eigenvectors u(j) are stored in the columns of VL, in the same order as their eigenvalues. If the j-th eigenvalue is real, then u(j) = VL(:,j). If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and u(j+1) = VL(:,j) - i*VL(:,j+1). Each eigenvector is scaled so that its largest component has abs(real part) + abs(imag. part) = 1, except for eigenvectors corresponding to an eigenvalue with alpha = beta = 0, which are set to zero. Not referenced if JOBVL = 'N'. LDVL (input) INTEGER The leading dimension of the matrix VL. LDVL >= 1, and if JOBVL = 'V', LDVL >= N. VR (output) DOUBLE PRECISION array, dimension (LDVR,N) If JOBVR = 'V', the right eigenvectors x(j) are stored in the columns of VR, in the same order as their eigenvalues. If the j-th eigenvalue is real, then x(j) = VR(:,j). If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then x(j) = VR(:,j) + i*VR(:,j+1) and x(j+1) = VR(:,j) - i*VR(:,j+1). Each eigenvector is scaled so that its largest component has abs(real part) + abs(imag. part) = 1, except for eigenvalues corresponding to an eigenvalue with alpha = beta = 0, which are set to zero. Not referenced if JOBVR = 'N'. LDVR (input) INTEGER The leading dimension of the matrix VR. LDVR >= 1, and if JOBVR = 'V', LDVR >= N. WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) INTEGER The dimension of the array WORK. LWORK >= max(1,8*N). For good performance, LWORK must generally be larger. To compute the optimal value of LWORK, call ILAENV to get blocksizes (for DGEQRF, DORMQR, and DORGQR.) Then compute: NB -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR; The optimal LWORK is: 2*N + MAX( 6*N, N*(NB+1) ). 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 failed. No eigenvectors have been calculated, but ALPHAR(j), ALPHAI(j), and BETA(j) should be correct for j=INFO+1,...,N. > N: errors that usually indicate LAPACK problems: =N+1: error return from DGGBAL =N+2: error return from DGEQRF =N+3: error return from DORMQR =N+4: error return from DORGQR =N+5: error return from DGGHRD =N+6: error return from DHGEQZ (other than failed iteration) =N+7: error return from DTGEVC =N+8: error return from DGGBAK (computing VL) =N+9: error return from DGGBAK (computing VR) =N+10: error return from DLASCL (various calls) Further Details =============== Balancing --------- This driver calls DGGBAL to both permute and scale rows and columns of A and B. The permutations PL and PR are chosen so that PL*A*PR and PL*B*R will be upper triangular except for the diagonal blocks A(i:j,i:j) and B(i:j,i:j), with i and j as close together as possible. The diagonal scaling matrices DL and DR are chosen so that the pair DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to one (except for the elements that start out zero.) After the eigenvalues and eigenvectors of the balanced matrices have been computed, DGGBAK transforms the eigenvectors back to what they would have been (in perfect arithmetic) if they had not been balanced. Contents of A and B on Exit -------- -- - --- - -- ---- If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or both), then on exit the arrays A and B will contain the real Schur form[*] of the "balanced" versions of A and B. If no eigenvectors are computed, then only the diagonal blocks will be correct. [*] See DHGEQZ, DGEGS, or read the book "Matrix Computations", by Golub & van Loan, pub. by Johns Hopkins U. Press.
Extends from Modelica.Icons.Function (Icon for functions).
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 |
Computes the minimum-norm solution to a real linear least squares problem with rank deficient A
Lapack documentation Purpose ======= This routine is deprecated and has been replaced by routine DGELSY. DGELSX computes the minimum-norm solution to a real linear least squares problem: minimize || A * X - B || using a complete orthogonal factorization of A. A is an M-by-N matrix which may be rank-deficient. Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix X. The routine first computes a QR factorization with column pivoting: A * P = Q * [ R11 R12 ] [ 0 R22 ] with R11 defined as the largest leading submatrix whose estimated condition number is less than 1/RCOND. The order of R11, RANK, is the effective rank of A. Then, R22 is considered to be negligible, and R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization: A * P = Q * [ T11 0 ] * Z [ 0 0 ] The minimum-norm solution is then X = P * Z' [ inv(T11)*Q1'*B ] [ 0 ] where Q1 consists of the first RANK columns of Q. Arguments ========= M (input) INTEGER The number of rows of the matrix A. M >= 0. N (input) INTEGER The number of columns of the matrix A. N >= 0. NRHS (input) INTEGER The number of right hand sides, i.e., the number of columns of matrices B and X. NRHS >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, A has been overwritten by details of its complete orthogonal factorization. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,M). B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) On entry, the M-by-NRHS right hand side matrix B. On exit, the N-by-NRHS solution matrix X. If m >= n and RANK = n, the residual sum-of-squares for the solution in the i-th column is given by the sum of squares of elements N+1:M in that column. LDB (input) INTEGER The leading dimension of the array B. LDB >= max(1,M,N). JPVT (input/output) INTEGER array, dimension (N) On entry, if JPVT(i) .ne. 0, the i-th column of A is an initial column, otherwise it is a free column. Before the QR factorization of A, all initial columns are permuted to the leading positions; only the remaining free columns are moved as a result of column pivoting during the factorization. On exit, if JPVT(i) = k, then the i-th column of A*P was the k-th column of A. RCOND (input) DOUBLE PRECISION RCOND is used to determine the effective rank of A, which is defined as the order of the largest leading triangular submatrix R11 in the QR factorization with pivoting of A, whose estimated condition number < 1/RCOND. RANK (output) INTEGER The effective rank of A, i.e., the order of the submatrix R11. This is the same as the order of the submatrix T11 in the complete orthogonal factorization of A. WORK (workspace) DOUBLE PRECISION array, dimension (max( min(M,N)+3*N, 2*min(M,N)+NRHS )), INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, :] | |
B[size(A, 1), :] | |
rcond | Reciprocal condition number to estimate rank |
Name | Description |
---|---|
X[max(size(A, 1), size(A, 2)), size(B, 2)] | Solution is in first size(A,2) rows |
info | |
rank | Effective rank of A |
Computes the minimum-norm solution to a real linear least squares problem with rank deficient A
Lapack documentation Purpose ======= This routine is deprecated and has been replaced by routine DGELSY. DGELSX computes the minimum-norm solution to a real linear least squares problem: minimize || A * X - B || using a complete orthogonal factorization of A. A is an M-by-N matrix which may be rank-deficient. Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix X. The routine first computes a QR factorization with column pivoting: A * P = Q * [ R11 R12 ] [ 0 R22 ] with R11 defined as the largest leading submatrix whose estimated condition number is less than 1/RCOND. The order of R11, RANK, is the effective rank of A. Then, R22 is considered to be negligible, and R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization: A * P = Q * [ T11 0 ] * Z [ 0 0 ] The minimum-norm solution is then X = P * Z' [ inv(T11)*Q1'*B ] [ 0 ] where Q1 consists of the first RANK columns of Q. Arguments ========= M (input) INTEGER The number of rows of the matrix A. M >= 0. N (input) INTEGER The number of columns of the matrix A. N >= 0. NRHS (input) INTEGER The number of right hand sides, i.e., the number of columns of matrices B and X. NRHS >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, A has been overwritten by details of its complete orthogonal factorization. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,M). B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) On entry, the M-by-NRHS right hand side matrix B. On exit, the N-by-NRHS solution matrix X. If m >= n and RANK = n, the residual sum-of-squares for the solution in the i-th column is given by the sum of squares of elements N+1:M in that column. LDB (input) INTEGER The leading dimension of the array B. LDB >= max(1,M,N). JPVT (input/output) INTEGER array, dimension (N) On entry, if JPVT(i) .ne. 0, the i-th column of A is an initial column, otherwise it is a free column. Before the QR factorization of A, all initial columns are permuted to the leading positions; only the remaining free columns are moved as a result of column pivoting during the factorization. On exit, if JPVT(i) = k, then the i-th column of A*P was the k-th column of A. RCOND (input) DOUBLE PRECISION RCOND is used to determine the effective rank of A, which is defined as the order of the largest leading triangular submatrix R11 in the QR factorization with pivoting of A, whose estimated condition number < 1/RCOND. RANK (output) INTEGER The effective rank of A, i.e., the order of the submatrix R11. This is the same as the order of the submatrix T11 in the complete orthogonal factorization of A. WORK (workspace) DOUBLE PRECISION array, dimension (max( min(M,N)+3*N, 2*min(M,N)+NRHS )), INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, :] | |
b[size(A, 1)] | |
rcond | Reciprocal condition number to estimate rank |
Name | Description |
---|---|
x[max(size(A, 1), size(A, 2))] | solution is in first size(A,2) rows |
info | |
rank | Effective rank of A |
Computes the minimum-norm solution to a real linear least squares problem with rank deficient A
Lapack documentation Purpose ======= DGELSY computes the minimum-norm solution to a real linear least squares problem: minimize || A * X - B || using a complete orthogonal factorization of A. A is an M-by-N matrix which may be rank-deficient. Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix X. The routine first computes a QR factorization with column pivoting: A * P = Q * [ R11 R12 ] [ 0 R22 ] with R11 defined as the largest leading submatrix whose estimated condition number is less than 1/RCOND. The order of R11, RANK, is the effective rank of A. Then, R22 is considered to be negligible, and R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization: A * P = Q * [ T11 0 ] * Z [ 0 0 ] The minimum-norm solution is then X = P * Z' [ inv(T11)*Q1'*B ] [ 0 ] where Q1 consists of the first RANK columns of Q. This routine is basically identical to the original xGELSX except three differences: o The call to the subroutine xGEQPF has been substituted by the call to the subroutine xGEQP3. This subroutine is a Blas-3 version of the QR factorization with column pivoting. o Matrix B (the right hand side) is updated with Blas-3. o The permutation of matrix B (the right hand side) is faster and more simple. Arguments ========= M (input) INTEGER The number of rows of the matrix A. M >= 0. N (input) INTEGER The number of columns of the matrix A. N >= 0. NRHS (input) INTEGER The number of right hand sides, i.e., the number of columns of matrices B and X. NRHS >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, A has been overwritten by details of its complete orthogonal factorization. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,M). B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) On entry, the M-by-NRHS right hand side matrix B. On exit, the N-by-NRHS solution matrix X. LDB (input) INTEGER The leading dimension of the array B. LDB >= max(1,M,N). JPVT (input/output) INTEGER array, dimension (N) On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted to the front of AP, otherwise column i is a free column. On exit, if JPVT(i) = k, then the i-th column of AP was the k-th column of A. RCOND (input) DOUBLE PRECISION RCOND is used to determine the effective rank of A, which is defined as the order of the largest leading triangular submatrix R11 in the QR factorization with pivoting of A, whose estimated condition number < 1/RCOND. RANK (output) INTEGER The effective rank of A, i.e., the order of the submatrix R11. This is the same as the order of the submatrix T11 in the complete orthogonal factorization of A. WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) INTEGER The dimension of the array WORK. The unblocked strategy requires that: LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ), where MN = min( M, N ). The block algorithm requires that: LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ), where NB is an upper bound on the blocksize returned by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR, and DORMRZ. 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.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, :] | |
B[size(A, 1), :] | |
rcond | Reciprocal condition number to estimate rank |
Name | Description |
---|---|
X[max(size(A, 1), size(A, 2)), size(B, 2)] | Solution is in first size(A,2) rows |
info | |
rank | Effective rank of A |
Computes the minimum-norm solution to a real linear least squares problem with rank deficient A
Lapack documentation Purpose ======= DGELSY computes the minimum-norm solution to a real linear least squares problem: minimize || A * X - B || using a complete orthogonal factorization of A. A is an M-by-N matrix which may be rank-deficient. Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix X. The routine first computes a QR factorization with column pivoting: A * P = Q * [ R11 R12 ] [ 0 R22 ] with R11 defined as the largest leading submatrix whose estimated condition number is less than 1/RCOND. The order of R11, RANK, is the effective rank of A. Then, R22 is considered to be negligible, and R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization: A * P = Q * [ T11 0 ] * Z [ 0 0 ] The minimum-norm solution is then X = P * Z' [ inv(T11)*Q1'*B ] [ 0 ] where Q1 consists of the first RANK columns of Q. This routine is basically identical to the original xGELSX except three differences: o The call to the subroutine xGEQPF has been substituted by the call to the subroutine xGEQP3. This subroutine is a Blas-3 version of the QR factorization with column pivoting. o Matrix B (the right hand side) is updated with Blas-3. o The permutation of matrix B (the right hand side) is faster and more simple. Arguments ========= M (input) INTEGER The number of rows of the matrix A. M >= 0. N (input) INTEGER The number of columns of the matrix A. N >= 0. NRHS (input) INTEGER The number of right hand sides, i.e., the number of columns of matrices B and X. NRHS >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, A has been overwritten by details of its complete orthogonal factorization. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,M). B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) On entry, the M-by-NRHS right hand side matrix B. On exit, the N-by-NRHS solution matrix X. LDB (input) INTEGER The leading dimension of the array B. LDB >= max(1,M,N). JPVT (input/output) INTEGER array, dimension (N) On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted to the front of AP, otherwise column i is a free column. On exit, if JPVT(i) = k, then the i-th column of AP was the k-th column of A. RCOND (input) DOUBLE PRECISION RCOND is used to determine the effective rank of A, which is defined as the order of the largest leading triangular submatrix R11 in the QR factorization with pivoting of A, whose estimated condition number < 1/RCOND. RANK (output) INTEGER The effective rank of A, i.e., the order of the submatrix R11. This is the same as the order of the submatrix T11 in the complete orthogonal factorization of A. WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) INTEGER The dimension of the array WORK. The unblocked strategy requires that: LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ), where MN = min( M, N ). The block algorithm requires that: LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ), where NB is an upper bound on the blocksize returned by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR, and DORMRZ. 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.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, :] | |
b[size(A, 1)] | |
rcond | Reciprocal condition number to estimate rank |
Name | Description |
---|---|
x[max(size(A, 1), size(A, 2))] | solution is in first size(A,2) rows |
info | |
rank | Effective rank of A |
Solves overdetermined or underdetermined real linear equations A*x=b with a b vector
Lapack documentation Purpose ======= DGELS solves overdetermined or underdetermined real linear systems involving an M-by-N matrix A, or its transpose, using a QR or LQ factorization of A. It is assumed that A has full rank. The following options are provided: 1. If TRANS = 'N' and m >= n: find the least squares solution of an overdetermined system, i.e., solve the least squares problem minimize || B - A*X ||. 2. If TRANS = 'N' and m < n: find the minimum norm solution of an underdetermined system A * X = B. 3. If TRANS = 'T' and m >= n: find the minimum norm solution of an undetermined system A**T * X = B. 4. If TRANS = 'T' and m < n: find the least squares solution of an overdetermined system, i.e., solve the least squares problem minimize || B - A**T * X ||. Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix X. Arguments ========= TRANS (input) CHARACTER*1 = 'N': the linear system involves A; = 'T': the linear system involves A**T. M (input) INTEGER The number of rows of the matrix A. M >= 0. N (input) INTEGER The number of columns of the matrix A. N >= 0. NRHS (input) INTEGER The number of right hand sides, i.e., the number of columns of the matrices B and X. NRHS >=0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, if M >= N, A is overwritten by details of its QR factorization as returned by DGEQRF; if M < N, A is overwritten by details of its LQ factorization as returned by DGELQF. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,M). B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) On entry, the matrix B of right hand side vectors, stored columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS if TRANS = 'T'. On exit, if INFO = 0, B is overwritten by the solution vectors, stored columnwise: if TRANS = 'N' and m >= n, rows 1 to n of B contain the least squares solution vectors; the residual sum of squares for the solution in each column is given by the sum of squares of elements N+1 to M in that column; if TRANS = 'N' and m < n, rows 1 to N of B contain the minimum norm solution vectors; if TRANS = 'T' and m >= n, rows 1 to M of B contain the minimum norm solution vectors; if TRANS = 'T' and m < n, rows 1 to M of B contain the least squares solution vectors; the residual sum of squares for the solution in each column is given by the sum of squares of elements M+1 to N in that column. LDB (input) INTEGER The leading dimension of the array B. LDB >= MAX(1,M,N). WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) INTEGER The dimension of the array WORK. LWORK >= max( 1, MN + max( MN, NRHS ) ). For optimal performance, LWORK >= max( 1, MN + max( MN, NRHS )*NB ). where MN = min(M,N) and NB is the optimum block size. 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 > 0: if INFO = i, the i-th diagonal element of the triangular factor of A is zero, so that A does not have full rank; the least squares solution could not be computed.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, :] | |
b[size(A, 1)] |
Name | Description |
---|---|
x[max(size(A, 1), size(A, 2))] | solution is in first size(A,2) rows |
info |
Solve real system of linear equations A*X=B with a B matrix
Lapack documentation Purpose ======= DGESV computes the solution to a real system of linear equations A * X = B, where A is an N-by-N matrix and X and B are N-by-NRHS matrices. The LU decomposition with partial pivoting and row interchanges is used to factor A as A = P * L * U, where P is a permutation matrix, L is unit lower triangular, and U is upper triangular. The factored form of A is then used to solve the system of equations A * X = B. Arguments ========= N (input) INTEGER The number of linear equations, i.e., the order of the matrix A. N >= 0. NRHS (input) INTEGER The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the N-by-N coefficient matrix A. On exit, the factors L and U from the factorization A = P*L*U; the unit diagonal elements of L are not stored. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,N). IPIV (output) INTEGER array, dimension (N) The pivot indices that define the permutation matrix P; row i of the matrix was interchanged with row IPIV(i). B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) On entry, the N-by-NRHS matrix of right hand side matrix B. On exit, if INFO = 0, the N-by-NRHS solution matrix X. LDB (input) INTEGER The leading dimension of the array B. LDB >= max(1,N). INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, U(i,i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, so the solution could not be computed.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, size(A, 1)] | |
B[size(A, 1), :] |
Name | Description |
---|---|
X[size(A, 1), size(B, 2)] | |
info |
Solve real system of linear equations A*x=b with a b vector
Same as function LAPACK.dgesv, but right hand side is a vector and not a matrix. For details of the arguments, see documentation of dgesv.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, size(A, 1)] | |
b[size(A, 1)] |
Name | Description |
---|---|
x[size(A, 1)] | |
info |
Solve a linear equality constrained least squares problem
Lapack documentation Purpose ======= DGGLSE solves the linear equality-constrained least squares (LSE) problem: minimize || c - A*x ||_2 subject to B*x = d where A is an M-by-N matrix, B is a P-by-N matrix, c is a given M-vector, and d is a given P-vector. It is assumed that P <= N <= M+P, and rank(B) = P and rank( (A) ) = N. ( (B) ) These conditions ensure that the LSE problem has a unique solution, which is obtained using a generalized RQ factorization of the matrices (B, A) given by B = (0 R)*Q, A = Z*T*Q. Arguments ========= M (input) INTEGER The number of rows of the matrix A. M >= 0. N (input) INTEGER The number of columns of the matrices A and B. N >= 0. P (input) INTEGER The number of rows of the matrix B. 0 <= P <= N <= M+P. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, the elements on and above the diagonal of the array contain the min(M,N)-by-N upper trapezoidal matrix T. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,M). B (input/output) DOUBLE PRECISION array, dimension (LDB,N) On entry, the P-by-N matrix B. On exit, the upper triangle of the subarray B(1:P,N-P+1:N) contains the P-by-P upper triangular matrix R. LDB (input) INTEGER The leading dimension of the array B. LDB >= max(1,P). C (input/output) DOUBLE PRECISION array, dimension (M) On entry, C contains the right hand side vector for the least squares part of the LSE problem. On exit, the residual sum of squares for the solution is given by the sum of squares of elements N-P+1 to M of vector C. D (input/output) DOUBLE PRECISION array, dimension (P) On entry, D contains the right hand side vector for the constrained equation. On exit, D is destroyed. X (output) DOUBLE PRECISION array, dimension (N) On exit, X is the solution of the LSE problem. WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) INTEGER The dimension of the array WORK. LWORK >= max(1,M+N+P). For optimum performance LWORK >= P+min(M,N)+max(M,N)*NB, where NB is an upper bound for the optimal blocksizes for DGEQRF, SGERQF, DORMQR and SORMRQ. 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: the upper triangular factor R associated with B in the generalized RQ factorization of the pair (B, A) is singular, so that rank(B) < P; the least squares solution could not be computed. = 2: the (N-P) by (N-P) part of the upper trapezoidal factor T associated with A in the generalized RQ factorization of the pair (B, A) is singular, so that rank( (A) ) < N; the least squares solution could not ( (B) ) be computed.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, :] | Minimize |A*x - c|^2 |
c[size(A, 1)] | |
B[:, size(A, 2)] | subject to B*x=d |
d[size(B, 1)] |
Name | Description |
---|---|
x[size(A, 2)] | solution vector |
info |
Solve real system of linear equations A*X=B with B matrix and tridiagonal A
Lapack documentation Purpose ======= DGTSV solves the equation A*X = B, where A is an n by n tridiagonal matrix, by Gaussian elimination with partial pivoting. Note that the equation A'*X = B may be solved by interchanging the order of the arguments DU and DL. Arguments ========= N (input) INTEGER The order of the matrix A. N >= 0. NRHS (input) INTEGER The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0. DL (input/output) DOUBLE PRECISION array, dimension (N-1) On entry, DL must contain the (n-1) sub-diagonal elements of A. On exit, DL is overwritten by the (n-2) elements of the second super-diagonal of the upper triangular matrix U from the LU factorization of A, in DL(1), ..., DL(n-2). D (input/output) DOUBLE PRECISION array, dimension (N) On entry, D must contain the diagonal elements of A. On exit, D is overwritten by the n diagonal elements of U. DU (input/output) DOUBLE PRECISION array, dimension (N-1) On entry, DU must contain the (n-1) super-diagonal elements of A. On exit, DU is overwritten by the (n-1) elements of the first super-diagonal of U. B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) On entry, the N by NRHS matrix of right hand side matrix B. On exit, if INFO = 0, the N by NRHS solution matrix X. LDB (input) INTEGER The leading dimension of the array B. LDB >= max(1,N). INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, U(i,i) is exactly zero, and the solution has not been computed. The factorization has not been completed unless i = N.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
superdiag[:] | |
diag[size(superdiag, 1) + 1] | |
subdiag[size(superdiag, 1)] | |
B[size(diag, 1), :] |
Name | Description |
---|---|
X[size(B, 1), size(B, 2)] | |
info |
Solve real system of linear equations A*x=b with b vector and tridiagonal A
Same as function LAPACK.dgtsv, but right hand side is a vector and not a matrix. For details of the arguments, see documentation of dgtsv.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
superdiag[:] | |
diag[size(superdiag, 1) + 1] | |
subdiag[size(superdiag, 1)] | |
b[size(diag, 1)] |
Name | Description |
---|---|
x[size(b, 1)] | |
info |
Solve real system of linear equations A*X=B with a B matrix
Lapack documentation Purpose ======= DGBSV computes the solution to a real system of linear equations A * X = B, where A is a band matrix of order N with KL subdiagonals and KU superdiagonals, and X and B are N-by-NRHS matrices. The LU decomposition with partial pivoting and row interchanges is used to factor A as A = L * U, where L is a product of permutation and unit lower triangular matrices with KL subdiagonals, and U is upper triangular with KL+KU superdiagonals. The factored form of A is then used to solve the system of equations A * X = B. Arguments ========= N (input) INTEGER The number of linear equations, i.e., the order of the matrix A. N >= 0. KL (input) INTEGER The number of subdiagonals within the band of A. KL >= 0. KU (input) INTEGER The number of superdiagonals within the band of A. KU >= 0. NRHS (input) INTEGER The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0. AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N) On entry, the matrix A in band storage, in rows KL+1 to 2*KL+KU+1; rows 1 to KL of the array need not be set. The j-th column of A is stored in the j-th column of the array AB as follows: AB(KL+KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+KL) On exit, details of the factorization: U is stored as an upper triangular band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and the multipliers used during the factorization are stored in rows KL+KU+2 to 2*KL+KU+1. See below for further details. LDAB (input) INTEGER The leading dimension of the array AB. LDAB >= 2*KL+KU+1. IPIV (output) INTEGER array, dimension (N) The pivot indices that define the permutation matrix P; row i of the matrix was interchanged with row IPIV(i). B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) On entry, the N-by-NRHS right hand side matrix B. On exit, if INFO = 0, the N-by-NRHS solution matrix X. LDB (input) INTEGER The leading dimension of the array B. LDB >= max(1,N). INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, U(i,i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, and the solution has not been computed. Further Details =============== The band storage scheme is illustrated by the following example, when M = N = 6, KL = 2, KU = 1: On entry: On exit: * * * + + + * * * u14 u25 u36 * * + + + + * * u13 u24 u35 u46 * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 * a31 a42 a53 a64 * * m31 m42 m53 m64 * * Array elements marked * are not used by the routine; elements marked + need not be set on entry, but are required by the routine to store elements of U because of fill-in resulting from the row interchanges.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
n | Number of equations |
kLower | Number of lower bands |
kUpper | Number of upper bands |
A[2*kLower + kUpper + 1, n] | |
B[n, :] |
Name | Description |
---|---|
X[n, size(B, 2)] | |
info |
Solve real system of linear equations A*x=b with a b vector
Same as function LAPACK.dgbsv, but right hand side is a vector and not a matrix. For details of the arguments, see documentation of dgbsv.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
n | Number of equations |
kLower | Number of lower bands |
kUpper | Number of upper bands |
A[2*kLower + kUpper + 1, n] | |
b[n] |
Name | Description |
---|---|
x[n] | |
info |
Determine singular value decomposition
Lapack documentation Purpose ======= DGESVD computes the singular value decomposition (SVD) of a real M-by-N matrix A, optionally computing the left and/or right singular vectors. The SVD is written A = U * SIGMA * transpose(V) where SIGMA is an M-by-N matrix which is zero except for its min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA are the singular values of A; they are real and non-negative, and are returned in descending order. The first min(m,n) columns of U and V are the left and right singular vectors of A. Note that the routine returns V**T, not V. Arguments ========= JOBU (input) CHARACTER*1 Specifies options for computing all or part of the matrix U: = 'A': all M columns of U are returned in array U: = 'S': the first min(m,n) columns of U (the left singular vectors) are returned in the array U; = 'O': the first min(m,n) columns of U (the left singular vectors) are overwritten on the array A; = 'N': no columns of U (no left singular vectors) are computed. JOBVT (input) CHARACTER*1 Specifies options for computing all or part of the matrix V**T: = 'A': all N rows of V**T are returned in the array VT; = 'S': the first min(m,n) rows of V**T (the right singular vectors) are returned in the array VT; = 'O': the first min(m,n) rows of V**T (the right singular vectors) are overwritten on the array A; = 'N': no rows of V**T (no right singular vectors) are computed. JOBVT and JOBU cannot both be 'O'. M (input) INTEGER The number of rows of the input matrix A. M >= 0. N (input) INTEGER The number of columns of the input matrix A. N >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, if JOBU = 'O', A is overwritten with the first min(m,n) columns of U (the left singular vectors, stored columnwise); if JOBVT = 'O', A is overwritten with the first min(m,n) rows of V**T (the right singular vectors, stored rowwise); if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A are destroyed. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,M). S (output) DOUBLE PRECISION array, dimension (min(M,N)) The singular values of A, sorted so that S(i) >= S(i+1). U (output) DOUBLE PRECISION array, dimension (LDU,UCOL) (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'. If JOBU = 'A', U contains the M-by-M orthogonal matrix U; if JOBU = 'S', U contains the first min(m,n) columns of U (the left singular vectors, stored columnwise); if JOBU = 'N' or 'O', U is not referenced. LDU (input) INTEGER The leading dimension of the array U. LDU >= 1; if JOBU = 'S' or 'A', LDU >= M. VT (output) DOUBLE PRECISION array, dimension (LDVT,N) If JOBVT = 'A', VT contains the N-by-N orthogonal matrix V**T; if JOBVT = 'S', VT contains the first min(m,n) rows of V**T (the right singular vectors, stored rowwise); if JOBVT = 'N' or 'O', VT is not referenced. LDVT (input) INTEGER The leading dimension of the array VT. LDVT >= 1; if JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N). WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK; if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged superdiagonal elements of an upper bidiagonal matrix B whose diagonal is in S (not necessarily sorted). B satisfies A = U * B * VT, so it has the same singular values as A, and singular vectors related by U and VT. LWORK (input) INTEGER The dimension of the array WORK. LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)). For good performance, LWORK should generally be larger. 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. > 0: if DBDSQR did not converge, INFO specifies how many superdiagonals of an intermediate bidiagonal form B did not converge to zero. See the description of WORK above for details.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, :] |
Name | Description |
---|---|
sigma[min(size(A, 1), size(A, 2))] | |
U[size(A, 1), size(A, 1)] | |
VT[size(A, 2), size(A, 2)] | |
info |
Determine singular values
Lapack documentation Purpose ======= DGESVD computes the singular value decomposition (SVD) of a real M-by-N matrix A, optionally computing the left and/or right singular vectors. The SVD is written A = U * SIGMA * transpose(V) where SIGMA is an M-by-N matrix which is zero except for its min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA are the singular values of A; they are real and non-negative, and are returned in descending order. The first min(m,n) columns of U and V are the left and right singular vectors of A. Note that the routine returns V**T, not V. Arguments ========= JOBU (input) CHARACTER*1 Specifies options for computing all or part of the matrix U: = 'A': all M columns of U are returned in array U: = 'S': the first min(m,n) columns of U (the left singular vectors) are returned in the array U; = 'O': the first min(m,n) columns of U (the left singular vectors) are overwritten on the array A; = 'N': no columns of U (no left singular vectors) are computed. JOBVT (input) CHARACTER*1 Specifies options for computing all or part of the matrix V**T: = 'A': all N rows of V**T are returned in the array VT; = 'S': the first min(m,n) rows of V**T (the right singular vectors) are returned in the array VT; = 'O': the first min(m,n) rows of V**T (the right singular vectors) are overwritten on the array A; = 'N': no rows of V**T (no right singular vectors) are computed. JOBVT and JOBU cannot both be 'O'. M (input) INTEGER The number of rows of the input matrix A. M >= 0. N (input) INTEGER The number of columns of the input matrix A. N >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, if JOBU = 'O', A is overwritten with the first min(m,n) columns of U (the left singular vectors, stored columnwise); if JOBVT = 'O', A is overwritten with the first min(m,n) rows of V**T (the right singular vectors, stored rowwise); if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A are destroyed. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,M). S (output) DOUBLE PRECISION array, dimension (min(M,N)) The singular values of A, sorted so that S(i) >= S(i+1). U (output) DOUBLE PRECISION array, dimension (LDU,UCOL) (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'. If JOBU = 'A', U contains the M-by-M orthogonal matrix U; if JOBU = 'S', U contains the first min(m,n) columns of U (the left singular vectors, stored columnwise); if JOBU = 'N' or 'O', U is not referenced. LDU (input) INTEGER The leading dimension of the array U. LDU >= 1; if JOBU = 'S' or 'A', LDU >= M. VT (output) DOUBLE PRECISION array, dimension (LDVT,N) If JOBVT = 'A', VT contains the N-by-N orthogonal matrix V**T; if JOBVT = 'S', VT contains the first min(m,n) rows of V**T (the right singular vectors, stored rowwise); if JOBVT = 'N' or 'O', VT is not referenced. LDVT (input) INTEGER The leading dimension of the array VT. LDVT >= 1; if JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N). WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK; if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged superdiagonal elements of an upper bidiagonal matrix B whose diagonal is in S (not necessarily sorted). B satisfies A = U * B * VT, so it has the same singular values as A, and singular vectors related by U and VT. LWORK (input) INTEGER The dimension of the array WORK. LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)). For good performance, LWORK should generally be larger. 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. > 0: if DBDSQR did not converge, INFO specifies how many superdiagonals of an intermediate bidiagonal form B did not converge to zero. See the description of WORK above for details.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, :] |
Name | Description |
---|---|
sigma[min(size(A, 1), size(A, 2))] | |
info |
Compute LU factorization of square or rectangular matrix A (A = P*L*U)
Lapack documentation Purpose ======= DGETRF computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges. The factorization has the form A = P * L * U where P is a permutation matrix, L is lower triangular with unit diagonal elements (lower trapezoidal if m > n), and U is upper triangular (upper trapezoidal if m < n). This is the right-looking Level 3 BLAS version of the algorithm. Arguments ========= M (input) INTEGER The number of rows of the matrix A. M >= 0. N (input) INTEGER The number of columns of the matrix A. N >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the M-by-N matrix to be factored. On exit, the factors L and U from the factorization A = P*L*U; the unit diagonal elements of L are not stored. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,M). IPIV (output) INTEGER array, dimension (min(M,N)) The pivot indices; for 1 <= i <= min(M,N), row i of the matrix was interchanged with row IPIV(i). INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 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.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, :] | Square or rectangular matrix |
Name | Description |
---|---|
LU[size(A, 1), size(A, 2)] | |
pivots[min(size(A, 1), size(A, 2))] | Pivot vector |
info | Information |
Solves a system of linear equations with the LU decomposition from dgetrf(..)
Lapack documentation Purpose ======= DGETRS solves a system of linear equations A * X = B or A' * X = B with a general N-by-N matrix A using the LU factorization computed by DGETRF. Arguments ========= TRANS (input) CHARACTER*1 Specifies the form of the system of equations: = 'N': A * X = B (No transpose) = 'T': A'* X = B (Transpose) = 'C': A'* X = B (Conjugate transpose = Transpose) N (input) INTEGER The order of the matrix A. N >= 0. NRHS (input) INTEGER The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0. A (input) DOUBLE PRECISION array, dimension (LDA,N) The factors L and U from the factorization A = P*L*U as computed by DGETRF. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,N). IPIV (input) INTEGER array, dimension (N) The pivot indices from DGETRF; for 1<=i<=N, row i of the matrix was interchanged with row IPIV(i). B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) On entry, the right hand side matrix B. On exit, the solution matrix X. LDB (input) INTEGER The leading dimension of the array B. LDB >= max(1,N). INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
LU[:, size(LU, 1)] | LU factorization of dgetrf of a square matrix |
pivots[size(LU, 1)] | Pivot vector of dgetrf |
B[size(LU, 1), :] | Right hand side matrix B |
Name | Description |
---|---|
X[size(B, 1), size(B, 2)] | Solution matrix X |
info |
Solves a system of linear equations with the LU decomposition from dgetrf(..)
Lapack documentation Purpose ======= DGETRS solves a system of linear equations A * X = B or A' * X = B with a general N-by-N matrix A using the LU factorization computed by DGETRF. Arguments ========= TRANS (input) CHARACTER*1 Specifies the form of the system of equations: = 'N': A * X = B (No transpose) = 'T': A'* X = B (Transpose) = 'C': A'* X = B (Conjugate transpose = Transpose) N (input) INTEGER The order of the matrix A. N >= 0. NRHS (input) INTEGER The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0. A (input) DOUBLE PRECISION array, dimension (LDA,N) The factors L and U from the factorization A = P*L*U as computed by DGETRF. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,N). IPIV (input) INTEGER array, dimension (N) The pivot indices from DGETRF; for 1<=i<=N, row i of the matrix was interchanged with row IPIV(i). B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) On entry, the right hand side matrix B. On exit, the solution matrix X. LDB (input) INTEGER The leading dimension of the array B. LDB >= max(1,N). INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
LU[:, size(LU, 1)] | LU factorization of dgetrf of a square matrix |
pivots[size(LU, 1)] | Pivot vector of dgetrf |
b[size(LU, 1)] | Right hand side vector b |
Name | Description |
---|---|
x[size(b, 1)] | |
info |
Computes the inverse of a matrix using the LU factorization from dgetrf(..)
Lapack documentation Purpose ======= DGETRI computes the inverse of a matrix using the LU factorization computed by DGETRF. This method inverts U and then computes inv(A) by solving the system inv(A)*L = inv(U) for inv(A). Arguments ========= N (input) INTEGER The order of the matrix A. N >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the factors L and U from the factorization A = P*L*U as computed by DGETRF. On exit, if INFO = 0, the inverse of the original matrix A. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,N). IPIV (input) INTEGER array, dimension (N) The pivot indices from DGETRF; for 1<=i<=N, row i of the matrix was interchanged with row IPIV(i). WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO=0, then WORK(1) returns the optimal LWORK. LWORK (input) INTEGER The dimension of the array WORK. LWORK >= max(1,N). For optimal performance LWORK >= N*NB, where NB is the optimal blocksize returned by ILAENV. 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 > 0: if INFO = i, U(i,i) is exactly zero; the matrix is singular and its inverse could not be computed.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
LU[:, size(LU, 1)] | LU factorization of dgetrf of a square matrix |
pivots[size(LU, 1)] | Pivot vector of dgetrf |
Name | Description |
---|---|
inv[size(LU, 1), size(LU, 2)] | Inverse of matrix P*L*U |
info |
Compute QR factorization of square or rectangular matrix A with column pivoting (A(:,p) = Q*R)
Lapack documentation Purpose ======= This routine is deprecated and has been replaced by routine DGEQP3. DGEQPF computes a QR factorization with column pivoting of a real M-by-N matrix A: A*P = Q*R. Arguments ========= M (input) INTEGER The number of rows of the matrix A. M >= 0. N (input) INTEGER The number of columns of the matrix A. N >= 0 A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, the upper triangle of the array contains the min(M,N)-by-N upper triangular matrix R; the elements below the diagonal, together with the array TAU, represent the orthogonal matrix Q as a product of min(m,n) elementary reflectors. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,M). JPVT (input/output) INTEGER array, dimension (N) On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted to the front of A*P (a leading column); if JPVT(i) = 0, the i-th column of A is a free column. On exit, if JPVT(i) = k, then the i-th column of A*P was the k-th column of A. TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) The scalar factors of the elementary reflectors. WORK (workspace) DOUBLE PRECISION array, dimension (3*N) 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 elementary reflectors Q = H(1) H(2) . . . H(n) Each H(i) has the form H = I - tau * v * v' where tau is a real scalar, and v is a real vector with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i). The matrix P is represented in jpvt as follows: If jpvt(j) = i then the jth column of P is the ith canonical unit vector.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, :] | Square or rectangular matrix |
Name | Description |
---|---|
QR[size(A, 1), size(A, 2)] | QR factorization in packed format |
tau[min(size(A, 1), size(A, 2))] | The scalar factors of the elementary reflectors of Q |
p[size(A, 2)] | Pivot vector |
info |
Generates a Real orthogonal matrix Q which is defined as the product of elementary reflectors as returned from dgeqpf
Lapack documentation Purpose ======= DORGQR generates an M-by-N real matrix Q with orthonormal columns, which is defined as the first N columns of a product of K elementary reflectors of order M Q = H(1) H(2) . . . H(k) as returned by DGEQRF. Arguments ========= M (input) INTEGER The number of rows of the matrix Q. M >= 0. N (input) INTEGER The number of columns of the matrix Q. M >= N >= 0. K (input) INTEGER The number of elementary reflectors whose product defines the matrix Q. N >= K >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the i-th column must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by DGEQRF in the first k columns of its array argument A. On exit, the M-by-N matrix Q. LDA (input) INTEGER The first dimension of the array A. LDA >= max(1,M). TAU (input) DOUBLE PRECISION array, dimension (K) TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by DGEQRF. WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,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). 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 has an illegal value
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
QR[:, :] | QR from dgeqpf |
tau[min(size(QR, 1), size(QR, 2))] | The scalar factors of the elementary reflectors of Q |
Name | Description |
---|---|
Q[size(QR, 1), size(QR, 2)] | Orthogonal matrix Q |
info |
Computes real Schur form T of real nonsymmetric matrix A, and, optionally, the matrix of Schur vectors Z as well as the eigenvalues
Lapack documentation Purpose ======= DGEES computes for an N-by-N real nonsymmetric matrix A, the eigenvalues, the real Schur form T, and, optionally, the matrix of Schur vectors Z. This gives the Schur factorization A = Z*T*(Z**T). Optionally, it also orders the eigenvalues on the diagonal of the real Schur form so that selected eigenvalues are at the top left. The leading columns of Z then form an orthonormal basis for the invariant subspace corresponding to the selected eigenvalues. A matrix is in real Schur form if it is upper quasi-triangular with 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in the form [ a b ] [ c a ] where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc). Arguments ========= JOBVS (input) CHARACTER*1 = 'N': Schur vectors are not computed; = 'V': Schur vectors are computed. SORT (input) CHARACTER*1 Specifies whether or not to order the eigenvalues on the diagonal of the Schur form. = 'N': Eigenvalues are not ordered; = 'S': Eigenvalues are ordered (see SELECT). SELECT (external procedure) LOGICAL FUNCTION of two DOUBLE PRECISION arguments SELECT must be declared EXTERNAL in the calling subroutine. If SORT = 'S', SELECT is used to select eigenvalues to sort to the top left of the Schur form. If SORT = 'N', SELECT is not referenced. An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if SELECT(WR(j),WI(j)) is true; i.e., if either one of a complex conjugate pair of eigenvalues is selected, then both complex eigenvalues are selected. Note that a selected complex eigenvalue may no longer satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since ordering may change the value of complex eigenvalues (especially if the eigenvalue is ill-conditioned); in this case INFO is set to N+2 (see INFO below). N (input) INTEGER The order of the matrix A. N >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the N-by-N matrix A. On exit, A has been overwritten by its real Schur form T. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,N). SDIM (output) INTEGER If SORT = 'N', SDIM = 0. If SORT = 'S', SDIM = number of eigenvalues (after sorting) for which SELECT is true. (Complex conjugate pairs for which SELECT is true for either eigenvalue count as 2.) WR (output) DOUBLE PRECISION array, dimension (N) WI (output) DOUBLE PRECISION array, dimension (N) WR and WI contain the real and imaginary parts, respectively, of the computed eigenvalues in the same order that they appear on the diagonal of the output Schur form T. Complex conjugate pairs of eigenvalues will appear consecutively with the eigenvalue having the positive imaginary part first. VS (output) DOUBLE PRECISION array, dimension (LDVS,N) If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur vectors. If JOBVS = 'N', VS is not referenced. LDVS (input) INTEGER The leading dimension of the array VS. LDVS >= 1; if JOBVS = 'V', LDVS >= N. WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) contains the optimal LWORK. LWORK (input) INTEGER The dimension of the array WORK. LWORK >= max(1,3*N). For good performance, LWORK must generally be larger. 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. BWORK (workspace) LOGICAL array, dimension (N) Not referenced if SORT = 'N'. INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if INFO = i, and i is <= N: the QR algorithm failed to compute all the eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI contain those eigenvalues which have converged; if JOBVS = 'V', VS contains the matrix which reduces A to its partially converged Schur form. = N+1: the eigenvalues could not be reordered because some eigenvalues were too close to separate (the problem is very ill-conditioned); = N+2: after reordering, roundoff changed values of some complex eigenvalues so that leading eigenvalues in the Schur form no longer satisfy SELECT=.TRUE. This could also be caused by underflow due to scaling.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, size(A, 1)] | Square matrix |
Name | Description |
---|---|
T[size(A, 1), size(A, 2)] | Real Schur form with A = Z*T*Z' |
Z[size(A, 1), size(A, 1)] | orthogonal matrix Z of Schur vectors |
eval_real[size(A, 1)] | real part of the eigenvectors of A |
eval_imag[size(A, 1)] | imaginary part of the eigenvectors of A |
info |
Reorder the real Schur factorization of a real matrix
Lapack documentation Purpose ======= DTRSEN reorders the real Schur factorization of a real matrix A = Q*T*Q**T, so that a selected cluster of eigenvalues appears in the leading diagonal blocks of the upper quasi-triangular matrix T, and the leading columns of Q form an orthonormal basis of the corresponding right invariant subspace. Optionally the routine computes the reciprocal condition numbers of the cluster of eigenvalues and/or the invariant subspace. T must be in Schur canonical form (as returned by DHSEQR), that is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block has its diagonal elements equal and its off-diagonal elements of opposite sign. Arguments ========= JOB (input) CHARACTER*1 Specifies whether condition numbers are required for the cluster of eigenvalues (S) or the invariant subspace (SEP): = 'N': none; = 'E': for eigenvalues only (S); = 'V': for invariant subspace only (SEP); = 'B': for both eigenvalues and invariant subspace (S and SEP). COMPQ (input) CHARACTER*1 = 'V': update the matrix Q of Schur vectors; = 'N': do not update Q. SELECT (input) LOGICAL array, dimension (N) SELECT specifies the eigenvalues in the selected cluster. To select a real eigenvalue w(j), SELECT(j) must be set to .TRUE.. To select a complex conjugate pair of eigenvalues w(j) and w(j+1), corresponding to a 2-by-2 diagonal block, either SELECT(j) or SELECT(j+1) or both must be set to .TRUE.; a complex conjugate pair of eigenvalues must be either both included in the cluster or both excluded. N (input) INTEGER The order of the matrix T. N >= 0. T (input/output) DOUBLE PRECISION array, dimension (LDT,N) On entry, the upper quasi-triangular matrix T, in Schur canonical form. On exit, T is overwritten by the reordered matrix T, again in Schur canonical form, with the selected eigenvalues in the leading diagonal blocks. LDT (input) INTEGER The leading dimension of the array T. LDT >= max(1,N). Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) On entry, if COMPQ = 'V', the matrix Q of Schur vectors. On exit, if COMPQ = 'V', Q has been postmultiplied by the orthogonal transformation matrix which reorders T; the leading M columns of Q form an orthonormal basis for the specified invariant subspace. If COMPQ = 'N', Q is not referenced. LDQ (input) INTEGER The leading dimension of the array Q. LDQ >= 1; and if COMPQ = 'V', LDQ >= N. WR (output) DOUBLE PRECISION array, dimension (N) WI (output) DOUBLE PRECISION array, dimension (N) The real and imaginary parts, respectively, of the reordered eigenvalues of T. The eigenvalues are stored in the same order as on the diagonal of T, with WR(i) = T(i,i) and, if T(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) > 0 and WI(i+1) = -WI(i). Note that if a complex eigenvalue is sufficiently ill-conditioned, then its value may differ significantly from its value before reordering. M (output) INTEGER The dimension of the specified invariant subspace. 0 < = M <= N. S (output) DOUBLE PRECISION If JOB = 'E' or 'B', S is a lower bound on the reciprocal condition number for the selected cluster of eigenvalues. S cannot underestimate the true reciprocal condition number by more than a factor of sqrt(N). If M = 0 or N, S = 1. If JOB = 'N' or 'V', S is not referenced. SEP (output) DOUBLE PRECISION If JOB = 'V' or 'B', SEP is the estimated reciprocal condition number of the specified invariant subspace. If M = 0 or N, SEP = norm(T). If JOB = 'N' or 'E', SEP is not referenced. WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) INTEGER The dimension of the array WORK. If JOB = 'N', LWORK >= max(1,N); if JOB = 'E', LWORK >= max(1,M*(N-M)); if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)). 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. IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK)) On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. LIWORK (input) INTEGER The dimension of the array IWORK. If JOB = 'N' or 'E', LIWORK >= 1; if JOB = 'V' or 'B', LIWORK >= max(1,M*(N-M)). If LIWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the IWORK array, returns this value as the first entry of the IWORK array, and no error message related to LIWORK is issued by XERBLA. INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value = 1: reordering of T failed because some eigenvalues are too close to separate (the problem is very ill-conditioned); T may have been partially reordered, and WR and WI contain the eigenvalues in the same order as in T; S and SEP (if requested) are set to zero. Further Details =============== DTRSEN first collects the selected eigenvalues by computing an orthogonal transformation Z to move them to the top left corner of T. In other words, the selected eigenvalues are the eigenvalues of T11 in: Z'*T*Z = ( T11 T12 ) n1 ( 0 T22 ) n2 n1 n2 where N = n1+n2 and Z' means the transpose of Z. The first n1 columns of Z span the specified invariant subspace of T. If T has been obtained from the real Schur factorization of a matrix A = Q*T*Q', then the reordered real Schur factorization of A is given by A = (Q*Z)*(Z'*T*Z)*(Q*Z)', and the first n1 columns of Q*Z span the corresponding invariant subspace of A. The reciprocal condition number of the average of the eigenvalues of T11 may be returned in S. S lies between 0 (very badly conditioned) and 1 (very well conditioned). It is computed as follows. First we compute R so that P = ( I R ) n1 ( 0 0 ) n2 n1 n2 is the projector on the invariant subspace associated with T11. R is the solution of the Sylvester equation: T11*R - R*T22 = T12. Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote the two-norm of M. Then S is computed as the lower bound (1 + F-norm(R)**2)**(-1/2) on the reciprocal of 2-norm(P), the true reciprocal condition number. S cannot underestimate 1 / 2-norm(P) by more than a factor of sqrt(N). An approximate error bound for the computed average of the eigenvalues of T11 is EPS * norm(T) / S where EPS is the machine precision. The reciprocal condition number of the right invariant subspace spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP. SEP is defined as the separation of T11 and T22: sep( T11, T22 ) = sigma-min( C ) where sigma-min(C) is the smallest singular value of the n1*n2-by-n1*n2 matrix C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) ) I(m) is an m by m identity matrix, and kprod denotes the Kronecker product. We estimate sigma-min(C) by the reciprocal of an estimate of the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C) cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2). When SEP is small, small changes in T can cause large changes in the invariant subspace. An approximate bound on the maximum angular error in the computed right invariant subspace is EPS * norm(T) / SEP
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
job | Specifies the usage of a condition number |
compq | Is "V" if Schur vector matrix is to be updated |
select[:] | Specifies the eigenvalues to reorder |
T[:, :] | Real Schur form to be reordered |
Q[:, size(T, 2)] | Matrix of the Schur vectors |
Name | Description |
---|---|
To[:, :] | Reordered Schur form |
Qo[:, :] | Reordered Schur vectors |
wr[size(T, 2)] | Reordered eigenvalues, real part |
wi[size(T, 2)] | Reordered eigenvalues, imaginary part |
m | Dimension of the invariant sub space spanned bei the selected eigenvalues |
s | Lower bound of the reciprocal condition number. Not referenced for job==V |
sep | Estimated reciprocal condition number of the specified invariant subspace |
info |
Solve real system of linear equations op(A)*X=B, op(A) is A or A' according to the Boolean input transposed
Lapack documentation Purpose ======= DGESVX uses the LU factorization to compute the solution to a real system of linear equations A * X = B, where A is an N-by-N matrix and X and B are N-by-NRHS matrices. Error bounds on the solution and a condition estimate are also provided. Description =========== The following steps are performed: 1. If FACT = 'E', real scaling factors are computed to equilibrate the system: TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B Whether or not the system will be equilibrated depends on the scaling of the matrix A, but if equilibration is used, A is overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N') or diag(C)*B (if TRANS = 'T' or 'C'). 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the matrix A (after equilibration if FACT = 'E') as A = P * L * U, where P is a permutation matrix, L is a unit lower triangular matrix, and U is upper triangular. 3. If some U(i,i)=0, so that U is exactly singular, then the routine returns with INFO = i. Otherwise, the factored form of A is used to estimate the condition number of the matrix A. If the reciprocal of the condition number is less than machine precision, INFO = N+1 is returned as a warning, but the routine still goes on to solve for X and compute error bounds as described below. 4. The system of equations is solved for X using the factored form of A. 5. Iterative refinement is applied to improve the computed solution matrix and calculate error bounds and backward error estimates for it. 6. If equilibration was used, the matrix X is premultiplied by diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so that it solves the original system before equilibration. Arguments ========= FACT (input) CHARACTER*1 Specifies whether or not the factored form of the matrix A is supplied on entry, and if not, whether the matrix A should be equilibrated before it is factored. = 'F': On entry, AF and IPIV contain the factored form of A. If EQUED is not 'N', the matrix A has been equilibrated with scaling factors given by R and C. A, AF, and IPIV are not modified. = 'N': The matrix A will be copied to AF and factored. = 'E': The matrix A will be equilibrated if necessary, then copied to AF and factored. TRANS (input) CHARACTER*1 Specifies the form of the system of equations: = 'N': A * X = B (No transpose) = 'T': A**T * X = B (Transpose) = 'C': A**H * X = B (Transpose) N (input) INTEGER The number of linear equations, i.e., the order of the matrix A. N >= 0. NRHS (input) INTEGER The number of right hand sides, i.e., the number of columns of the matrices B and X. NRHS >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the N-by-N matrix A. If FACT = 'F' and EQUED is not 'N', then A must have been equilibrated by the scaling factors in R and/or C. A is not modified if FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit. On exit, if EQUED .ne. 'N', A is scaled as follows: EQUED = 'R': A := diag(R) * A EQUED = 'C': A := A * diag(C) EQUED = 'B': A := diag(R) * A * diag(C). LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,N). AF (input or output) DOUBLE PRECISION array, dimension (LDAF,N) If FACT = 'F', then AF is an input argument and on entry contains the factors L and U from the factorization A = P*L*U as computed by DGETRF. If EQUED .ne. 'N', then AF is the factored form of the equilibrated matrix A. If FACT = 'N', then AF is an output argument and on exit returns the factors L and U from the factorization A = P*L*U of the original matrix A. If FACT = 'E', then AF is an output argument and on exit returns the factors L and U from the factorization A = P*L*U of the equilibrated matrix A (see the description of A for the form of the equilibrated matrix). LDAF (input) INTEGER The leading dimension of the array AF. LDAF >= max(1,N). IPIV (input or output) INTEGER array, dimension (N) If FACT = 'F', then IPIV is an input argument and on entry contains the pivot indices from the factorization A = P*L*U as computed by DGETRF; row i of the matrix was interchanged with row IPIV(i). If FACT = 'N', then IPIV is an output argument and on exit contains the pivot indices from the factorization A = P*L*U of the original matrix A. If FACT = 'E', then IPIV is an output argument and on exit contains the pivot indices from the factorization A = P*L*U of the equilibrated matrix A. EQUED (input or output) CHARACTER*1 Specifies the form of equilibration that was done. = 'N': No equilibration (always true if FACT = 'N'). = 'R': Row equilibration, i.e., A has been premultiplied by diag(R). = 'C': Column equilibration, i.e., A has been postmultiplied by diag(C). = 'B': Both row and column equilibration, i.e., A has been replaced by diag(R) * A * diag(C). EQUED is an input argument if FACT = 'F'; otherwise, it is an output argument. R (input or output) DOUBLE PRECISION array, dimension (N) The row scale factors for A. If EQUED = 'R' or 'B', A is multiplied on the left by diag(R); if EQUED = 'N' or 'C', R is not accessed. R is an input argument if FACT = 'F'; otherwise, R is an output argument. If FACT = 'F' and EQUED = 'R' or 'B', each element of R must be positive. C (input or output) DOUBLE PRECISION array, dimension (N) The column scale factors for A. If EQUED = 'C' or 'B', A is multiplied on the right by diag(C); if EQUED = 'N' or 'R', C is not accessed. C is an input argument if FACT = 'F'; otherwise, C is an output argument. If FACT = 'F' and EQUED = 'C' or 'B', each element of C must be positive. B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) On entry, the N-by-NRHS right hand side matrix B. On exit, if EQUED = 'N', B is not modified; if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by diag(R)*B; if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is overwritten by diag(C)*B. LDB (input) INTEGER The leading dimension of the array B. LDB >= max(1,N). X (output) DOUBLE PRECISION array, dimension (LDX,NRHS) If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to the original system of equations. Note that A and B are modified on exit if EQUED .ne. 'N', and the solution to the equilibrated system is inv(diag(C))*X if TRANS = 'N' and EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C' and EQUED = 'R' or 'B'. LDX (input) INTEGER The leading dimension of the array X. LDX >= max(1,N). RCOND (output) DOUBLE PRECISION The estimate of the reciprocal condition number of the matrix A after equilibration (if done). If RCOND is less than the machine precision (in particular, if RCOND = 0), the matrix is singular to working precision. This condition is indicated by a return code of INFO > 0. FERR (output) DOUBLE PRECISION array, dimension (NRHS) The estimated forward error bound for each solution vector X(j) (the j-th column of the solution matrix X). If XTRUE is the true solution corresponding to X(j), FERR(j) is an estimated upper bound for the magnitude of the largest element in (X(j) - XTRUE) divided by the magnitude of the largest element in X(j). The estimate is as reliable as the estimate for RCOND, and is almost always a slight overestimate of the true error. BERR (output) DOUBLE PRECISION array, dimension (NRHS) The componentwise relative backward error of each solution vector X(j) (i.e., the smallest relative change in any element of A or B that makes X(j) an exact solution). WORK (workspace/output) DOUBLE PRECISION array, dimension (4*N) On exit, WORK(1) contains the reciprocal pivot growth factor norm(A)/norm(U). The "max absolute element" norm is used. If WORK(1) is much less than 1, then the stability of the LU factorization of the (equilibrated) matrix A could be poor. This also means that the solution X, condition estimator RCOND, and forward error bound FERR could be unreliable. If factorization fails with 0<INFO<=N, then WORK(1) contains the reciprocal pivot growth factor for the leading INFO columns of A. IWORK (workspace) INTEGER array, dimension (N) INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, and i is <= N: U(i,i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, so the solution and error bounds could not be computed. RCOND = 0 is returned. = N+1: U is nonsingular, but RCOND is less than machine precision, meaning that the matrix is singular to working precision. Nevertheless, the solution and error bounds are computed because there are a number of situations where the computed solution can be more accurate than the value of RCOND would suggest.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, size(A, 1)] | Real square matrix A |
B[size(A, 1), :] | Real matrix B |
transposed | True if the equation to be solved is A'*X=B |
Name | Description |
---|---|
X[size(A, 1), size(B, 2)] | Solution matrix |
info | |
rcond | reciprocal condition number of the matrix A |
Solve the real Sylvester matrix equation op(A)*X + X*op(B) = scale*C or op(A)*X - X*op(B) = scale*C
Lapack documentation Purpose ======= DTRSYL solves the real Sylvester matrix equation: op(A)*X + X*op(B) = scale*C or op(A)*X - X*op(B) = scale*C, where op(A) = A or A**T, and A and B are both upper quasi- triangular. A is M-by-M and B is N-by-N; the right hand side C and the solution X are M-by-N; and scale is an output scale factor, set <= 1 to avoid overflow in X. A and B must be in Schur canonical form (as returned by DHSEQR), that is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block has its diagonal elements equal and its off-diagonal elements of opposite sign. Arguments ========= TRANA (input) CHARACTER*1 Specifies the option op(A): = 'N': op(A) = A (No transpose) = 'T': op(A) = A**T (Transpose) = 'C': op(A) = A**H (Conjugate transpose = Transpose) TRANB (input) CHARACTER*1 Specifies the option op(B): = 'N': op(B) = B (No transpose) = 'T': op(B) = B**T (Transpose) = 'C': op(B) = B**H (Conjugate transpose = Transpose) ISGN (input) INTEGER Specifies the sign in the equation: = +1: solve op(A)*X + X*op(B) = scale*C = -1: solve op(A)*X - X*op(B) = scale*C M (input) INTEGER The order of the matrix A, and the number of rows in the matrices X and C. M >= 0. N (input) INTEGER The order of the matrix B, and the number of columns in the matrices X and C. N >= 0. A (input) DOUBLE PRECISION array, dimension (LDA,M) The upper quasi-triangular matrix A, in Schur canonical form. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,M). B (input) DOUBLE PRECISION array, dimension (LDB,N) The upper quasi-triangular matrix B, in Schur canonical form. LDB (input) INTEGER The leading dimension of the array B. LDB >= max(1,N). C (input/output) DOUBLE PRECISION array, dimension (LDC,N) On entry, the M-by-N right hand side matrix C. On exit, C is overwritten by the solution matrix X. LDC (input) INTEGER The leading dimension of the array C. LDC >= max(1,M) SCALE (output) DOUBLE PRECISION The scale factor, scale, set <= 1 to avoid overflow in X. INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value = 1: A and B have common or very close eigenvalues; perturbed values were used to solve the equation (but the matrices A and B are unchanged).
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, :] | Upper quais-triangular matrix |
B[:, :] | Upper quais-triangular matrix |
C[if tranA then size(A, 1) else size(A, 2), if tranB then size(B, 1) else size(B, 2)] | Right side of the Sylvester equation |
tranA | True if op(A)=A' |
tranB | True if op(B)=B' |
isgn | Specifies the sign in the equation, +1 or -1 |
Name | Description |
---|---|
X[size(C, 1), size(C, 2)] | Solution of the Sylvester equation |
scale | Scale factor |
info |
Compute eigenvalues of a matrix H using lapack routine DHSEQR for Hessenberg form matrix
Lapack documentation Purpose ======= DHSEQR computes the eigenvalues of a Hessenberg matrix H and, optionally, the matrices T and Z from the Schur decomposition H = Z T Z**T, where T is an upper quasi-triangular matrix (the Schur form), and Z is the orthogonal matrix of Schur vectors. Optionally Z may be postmultiplied into an input orthogonal matrix Q so that this routine can give the Schur factorization of a matrix A which has been reduced to the Hessenberg form H by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T. Arguments ========= JOB (input) CHARACTER*1 = 'E': compute eigenvalues only; = 'S': compute eigenvalues and the Schur form T. COMPZ (input) CHARACTER*1 = 'N': no Schur vectors are computed; = 'I': Z is initialized to the unit matrix and the matrix Z of Schur vectors of H is returned; = 'V': Z must contain an orthogonal matrix Q on entry, and the product Q*Z is returned. N (input) INTEGER The order of the matrix H. N >= 0. ILO (input) INTEGER IHI (input) INTEGER It is assumed that H 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, and then passed to DGEHRD when the matrix output by DGEBAL is reduced to Hessenberg form. Otherwise ILO and IHI should be set to 1 and N respectively. If N>0, then 1<=ILO<=IHI<=N. If N = 0, then ILO = 1 and IHI = 0. H (input/output) DOUBLE PRECISION array, dimension (LDH,N) On entry, the upper Hessenberg matrix H. On exit, if INFO = 0 and JOB = 'S', then H contains the upper quasi-triangular matrix T from the Schur decomposition (the Schur form); 2-by-2 diagonal blocks (corresponding to complex conjugate pairs of eigenvalues) are returned in standard form, with H(i,i) = H(i+1,i+1) and H(i+1,i)*H(i,i+1)<0. If INFO = 0 and JOB = 'E', the contents of H are unspecified on exit. (The output value of H when INFO>0 is given under the description of INFO below.) Unlike earlier versions of DHSEQR, this subroutine may explicitly H(i,j) = 0 for i>j and j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N. LDH (input) INTEGER The leading dimension of the array H. LDH >= max(1,N). WR (output) DOUBLE PRECISION array, dimension (N) WI (output) DOUBLE PRECISION array, dimension (N) The real and imaginary parts, respectively, of the computed eigenvalues. If two eigenvalues are computed as a complex conjugate pair, they are stored in consecutive elements of WR and WI, say the i-th and (i+1)th, with WI(i) > 0 and WI(i+1) < 0. If JOB = 'S', the eigenvalues are stored in the same order as on the diagonal of the Schur form returned in H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i). Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N) If COMPZ = 'N', Z is not referenced. If COMPZ = 'I', on entry Z need not be set and on exit, if INFO = 0, Z contains the orthogonal matrix Z of the Schur vectors of H. If COMPZ = 'V', on entry Z must contain an N-by-N matrix Q, which is assumed to be equal to the unit matrix except for the submatrix Z(ILO:IHI,ILO:IHI). On exit, if INFO = 0, Z contains Q*Z. Normally Q is the orthogonal matrix generated by DORGHR after the call to DGEHRD which formed the Hessenberg matrix H. (The output value of Z when INFO>0 is given under the description of INFO below.) LDZ (input) INTEGER The leading dimension of the array Z. if COMPZ = 'I' or COMPZ = 'V', then LDZ>=MAX(1,N). Otherwise, LDZ>=1. WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns an estimate of the optimal value for LWORK. LWORK (input) INTEGER The dimension of the array WORK. LWORK >= max(1,N) is sufficient and delivers very good and sometimes optimal performance. However, LWORK as large as 11*N may be required for optimal performance. A workspace query is recommended to determine the optimal workspace size. If LWORK = -1, then DHSEQR does a workspace query. In this case, DHSEQR checks the input parameters and estimates the optimal workspace size for the given values of N, ILO and IHI. The estimate is returned in WORK(1). No error message related to LWORK is issued by XERBLA. Neither H nor Z are accessed. INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, DHSEQR failed to compute all of the eigenvalues. Elements 1:ilo-1 and i+1:n of WR and WI contain those eigenvalues which have been successfully computed. (Failures are rare.) If INFO > 0 and JOB = 'E', then on exit, the remaining unconverged eigenvalues are the eigen- values of the upper Hessenberg matrix rows and columns ILO through INFO of the final, output value of H. If INFO > 0 and JOB = 'S', then on exit (*) (initial value of H)*U = U*(final value of H) where U is an orthogonal matrix. The final value of H is upper Hessenberg and quasi-triangular in rows and columns INFO+1 through IHI. If INFO > 0 and COMPZ = 'V', then on exit (final value of Z) = (initial value of Z)*U where U is the orthogonal matrix in (*) (regard- less of the value of JOB.) If INFO > 0 and COMPZ = 'I', then on exit (final value of Z) = U where U is the orthogonal matrix in (*) (regard- less of the value of JOB.) If INFO > 0 and COMPZ = 'N', then Z is not accessed.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
H[:, size(H, 1)] | Matrix H with Hessenberg form |
eigenValuesOnly | True to compute the eigenvalues. False to compute the Schur form too |
compz | Specifies the computation of the Schur vectors |
Z[:, :] | Matrix Z |
Name | Description |
---|---|
alphaReal[size(H, 1)] | Real part of alpha (eigenvalue=(alphaReal+i*alphaImag)) |
alphaImag[size(H, 1)] | Imaginary part of alpha (eigenvalue=(alphaReal+i*alphaImag)) |
info | |
Ho[:, :] | Schur decomposition (if eigenValuesOnly==false, unspecified else) |
Zo[:, :] | |
work[3*max(1, size(H, 1))] |
Norm of a matrix
Lapack documentation Purpose ======= DLANGE returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real matrix A. Description =========== DLANGE returns the value DLANGE = ( max(abs(A(i,j))), NORM = 'M' or 'm' ( ( norm1(A), NORM = '1', 'O' or 'o' ( ( normI(A), NORM = 'I' or 'i' ( ( normF(A), NORM = 'F', 'f', 'E' or 'e' where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a consistent matrix norm. Arguments ========= NORM (input) CHARACTER*1 Specifies the value to be returned in DLANGE as described above. M (input) INTEGER The number of rows of the matrix A. M >= 0. When M = 0, DLANGE is set to zero. N (input) INTEGER The number of columns of the matrix A. N >= 0. When N = 0, DLANGE is set to zero. A (input) DOUBLE PRECISION array, dimension (LDA,N) The m by n matrix A. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(M,1). WORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)), where LWORK >= M when NORM = 'I'; otherwise, WORK is not referenced.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, :] | Real matrix A |
norm | specifies the norm, i.e., 1, I, F, M |
Name | Description |
---|---|
anorm | norm of A |
Estimates the reciprocal of the condition number of a general real matrix A
Lapack documentation Purpose ======= DGECON estimates the reciprocal of the condition number of a general real matrix A, in either the 1-norm or the infinity-norm, using the LU factorization computed by DGETRF. An estimate is obtained for norm(inv(A)), and the reciprocal of the condition number is computed as RCOND = 1 / ( norm(A) * norm(inv(A)) ). Arguments ========= NORM (input) CHARACTER*1 Specifies whether the 1-norm condition number or the infinity-norm condition number is required: = '1' or 'O': 1-norm; = 'I': Infinity-norm. N (input) INTEGER The order of the matrix A. N >= 0. A (input) DOUBLE PRECISION array, dimension (LDA,N) The factors L and U from the factorization A = P*L*U as computed by DGETRF. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,N). ANORM (input) DOUBLE PRECISION If NORM = '1' or 'O', the 1-norm of the original matrix A. If NORM = 'I', the infinity-norm of the original matrix A. RCOND (output) DOUBLE PRECISION The reciprocal of the condition number of the matrix A, computed as RCOND = 1/(norm(A) * norm(inv(A))). WORK (workspace) DOUBLE PRECISION array, dimension (4*N) IWORK (workspace) INTEGER array, dimension (N) INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
LU_of_A[:, :] | LU factorization of a real matrix A |
inf | Is true if infinity norm is used and false for 1-norm |
anorm | norm of A |
Name | Description |
---|---|
rcond | Reciprocal condition number of A |
info |
reduces a real general matrix A to upper Hessenberg form H by an orthogonal similarity transformation: Q' * A * Q = H
Lapack documentation 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).
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, size(A, 1)] | |
ilo | lowest index where the original matrix had been Hessenbergform |
ihi | highest index where the original matrix had been Hessenbergform |
Name | Description |
---|---|
Aout[size(A, 1), size(A, 2)] | contains the Hessenberg form in the upper triangle and the first subdiagonal and below the first subdiagonal it contains the elementary reflectors which represents (with array tau) as a product the orthogonal matrix Q |
tau[max(size(A, 1), 1) - 1] | scalar factors of the elementary reflectors |
info |
computes a QR factorization without pivoting
Lapack documentation Purpose ======= DGEQRF computes a QR factorization of a real M-by-N matrix A: A = Q * R. Arguments ========= M (input) INTEGER The number of rows of the matrix A. M >= 0. N (input) INTEGER The number of columns of the matrix A. N >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, the elements on and above the diagonal of the array contain the min(M,N)-by-N upper trapezoidal matrix R (R is upper triangular if m >= n); the elements below the diagonal, with the array TAU, represent the orthogonal matrix Q as a product of min(m,n) elementary reflectors (see Further Details). LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,M). TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) The scalar factors of the elementary reflectors (see Further Details). WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,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). 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 elementary reflectors Q = H(1) H(2) . . . H(k), where k = min(m,n). 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-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), and tau in TAU(i).
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, :] | Square or rectangular matrix |
Name | Description |
---|---|
Aout[size(A, 1), size(A, 2)] | the upper triangle of the array contains the upper trapezoidal matrix R; the elements below the diagonal, together with the array TAU, represent the orthogonal matrix Q as a product of elementary reflectors |
tau[min(size(A, 1), size(A, 2))] | scalar factors of the elementary reflectors |
info | |
work[3*max(1, size(A, 2))] |
Compute the eigenvalues and the (real) left and right eigenvectors of matrix A, using lapack routine dgeevx
Lapack documentation Purpose ======= DGEEVX computes for an N-by-N real nonsymmetric matrix A, the eigenvalues and, optionally, the left and/or right eigenvectors. Optionally also, it computes a balancing transformation to improve the conditioning of the eigenvalues and eigenvectors (ILO, IHI, SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues (RCONDE), and reciprocal condition numbers for the right eigenvectors (RCONDV). The right eigenvector v(j) of A satisfies A * v(j) = lambda(j) * v(j) where lambda(j) is its eigenvalue. The left eigenvector u(j) of A satisfies u(j)**H * A = lambda(j) * u(j)**H where u(j)**H denotes the conjugate transpose of u(j). The computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real. Balancing a matrix means permuting the rows and columns to make it more nearly upper triangular, and applying a diagonal similarity transformation D * A * D**(-1), where D is a diagonal matrix, to make its rows and columns closer in norm and the condition numbers of its eigenvalues and eigenvectors smaller. The computed reciprocal condition numbers correspond to the balanced matrix. Permuting rows and columns will not change the condition numbers (in exact arithmetic) but diagonal scaling will. For further explanation of balancing, see section 4.10.2 of the LAPACK Users' Guide. Arguments ========= BALANC (input) CHARACTER*1 Indicates how the input matrix should be diagonally scaled and/or permuted to improve the conditioning of its eigenvalues. = 'N': Do not diagonally scale or permute; = 'P': Perform permutations to make the matrix more nearly upper triangular. Do not diagonally scale; = 'S': Diagonally scale the matrix, i.e. replace A by D*A*D**(-1), where D is a diagonal matrix chosen to make the rows and columns of A more equal in norm. Do not permute; = 'B': Both diagonally scale and permute A. Computed reciprocal condition numbers will be for the matrix after balancing and/or permuting. Permuting does not change condition numbers (in exact arithmetic), but balancing does. JOBVL (input) CHARACTER*1 = 'N': left eigenvectors of A are not computed; = 'V': left eigenvectors of A are computed. If SENSE = 'E' or 'B', JOBVL must = 'V'. JOBVR (input) CHARACTER*1 = 'N': right eigenvectors of A are not computed; = 'V': right eigenvectors of A are computed. If SENSE = 'E' or 'B', JOBVR must = 'V'. SENSE (input) CHARACTER*1 Determines which reciprocal condition numbers are computed. = 'N': None are computed; = 'E': Computed for eigenvalues only; = 'V': Computed for right eigenvectors only; = 'B': Computed for eigenvalues and right eigenvectors. If SENSE = 'E' or 'B', both left and right eigenvectors must also be computed (JOBVL = 'V' and JOBVR = 'V'). N (input) INTEGER The order of the matrix A. N >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the N-by-N matrix A. On exit, A has been overwritten. If JOBVL = 'V' or JOBVR = 'V', A contains the real Schur form of the balanced version of the input matrix A. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,N). WR (output) DOUBLE PRECISION array, dimension (N) WI (output) DOUBLE PRECISION array, dimension (N) WR and WI contain the real and imaginary parts, respectively, of the computed eigenvalues. Complex conjugate pairs of eigenvalues will appear consecutively with the eigenvalue having the positive imaginary part first. VL (output) DOUBLE PRECISION array, dimension (LDVL,N) If JOBVL = 'V', the left eigenvectors u(j) are stored one after another in the columns of VL, in the same order as their eigenvalues. If JOBVL = 'N', VL is not referenced. If the j-th eigenvalue is real, then u(j) = VL(:,j), the j-th column of VL. If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and u(j+1) = VL(:,j) - i*VL(:,j+1). LDVL (input) INTEGER The leading dimension of the array VL. LDVL >= 1; if JOBVL = 'V', LDVL >= N. VR (output) DOUBLE PRECISION array, dimension (LDVR,N) If JOBVR = 'V', the right eigenvectors v(j) are stored one after another in the columns of VR, in the same order as their eigenvalues. If JOBVR = 'N', VR is not referenced. If the j-th eigenvalue is real, then v(j) = VR(:,j), the j-th column of VR. If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and v(j+1) = VR(:,j) - i*VR(:,j+1). LDVR (input) INTEGER The leading dimension of the array VR. LDVR >= 1, and if JOBVR = 'V', LDVR >= N. ILO (output) INTEGER IHI (output) INTEGER ILO and IHI are integer values determined when A was balanced. The balanced A(i,j) = 0 if I > J and J = 1,...,ILO-1 or I = IHI+1,...,N. SCALE (output) DOUBLE PRECISION array, dimension (N) Details of the permutations and scaling factors applied when balancing A. If P(j) is the index of the row and column interchanged with row and column j, and D(j) is the scaling factor applied to row and column j, then SCALE(J) = P(J), for J = 1,...,ILO-1 = D(J), for J = ILO,...,IHI = P(J) for J = IHI+1,...,N. The order in which the interchanges are made is N to IHI+1, then 1 to ILO-1. ABNRM (output) DOUBLE PRECISION The one-norm of the balanced matrix (the maximum of the sum of absolute values of elements of any column). RCONDE (output) DOUBLE PRECISION array, dimension (N) RCONDE(j) is the reciprocal condition number of the j-th eigenvalue. RCONDV (output) DOUBLE PRECISION array, dimension (N) RCONDV(j) is the reciprocal condition number of the j-th right eigenvector. WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) INTEGER The dimension of the array WORK. If SENSE = 'N' or 'E', LWORK >= max(1,2*N), and if JOBVL = 'V' or JOBVR = 'V', LWORK >= 3*N. If SENSE = 'V' or 'B', LWORK >= N*(N+6). For good performance, LWORK must generally be larger. 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. IWORK (workspace) INTEGER array, dimension (2*N-2) If SENSE = 'N' or 'E', not referenced. INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if INFO = i, the QR algorithm failed to compute all the eigenvalues, and no eigenvectors or condition numbers have been computed; elements 1:ILO-1 and i+1:N of WR and WI contain eigenvalues which have converged.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, size(A, 1)] |
Name | Description |
---|---|
alphaReal[size(A, 1)] | Real part of alpha (eigenvalue=(alphaReal+i*alphaImag)) |
alphaImag[size(A, 1)] | Imaginary part of alpha (eigenvalue=(alphaReal+i*alphaImag)) |
lEigenVectors[size(A, 1), size(A, 1)] | left eigenvectors of matrix A |
rEigenVectors[size(A, 1), size(A, 1)] | right eigenvectors of matrix A |
AS[size(A, 1), size(A, 2)] | AS iss the real Schur form of the balanced version of the input matrix A |
info |
Determine singular value decomposition
Lapack documentation Purpose ======= DGESDD computes the singular value decomposition (SVD) of a real M-by-N matrix A, optionally computing the left and right singular vectors. If singular vectors are desired, it uses a divide-and-conquer algorithm. The SVD is written A = U * SIGMA * transpose(V) where SIGMA is an M-by-N matrix which is zero except for its min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA are the singular values of A; they are real and non-negative, and are returned in descending order. The first min(m,n) columns of U and V are the left and right singular vectors of A. Note that the routine returns VT = V**T, not V. The divide and conquer algorithm makes very mild assumptions about floating point arithmetic. It will work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none. Arguments ========= JOBZ (input) CHARACTER*1 Specifies options for computing all or part of the matrix U: = 'A': all M columns of U and all N rows of V**T are returned in the arrays U and VT; = 'S': the first min(M,N) columns of U and the first min(M,N) rows of V**T are returned in the arrays U and VT; = 'O': If M >= N, the first N columns of U are overwritten on the array A and all rows of V**T are returned in the array VT; otherwise, all columns of U are returned in the array U and the first M rows of V**T are overwritten in the array A; = 'N': no columns of U or rows of V**T are computed. M (input) INTEGER The number of rows of the input matrix A. M >= 0. N (input) INTEGER The number of columns of the input matrix A. N >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, if JOBZ = 'O', A is overwritten with the first N columns of U (the left singular vectors, stored columnwise) if M >= N; A is overwritten with the first M rows of V**T (the right singular vectors, stored rowwise) otherwise. if JOBZ .ne. 'O', the contents of A are destroyed. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,M). S (output) DOUBLE PRECISION array, dimension (min(M,N)) The singular values of A, sorted so that S(i) >= S(i+1). U (output) DOUBLE PRECISION array, dimension (LDU,UCOL) UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N; UCOL = min(M,N) if JOBZ = 'S'. If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M orthogonal matrix U; if JOBZ = 'S', U contains the first min(M,N) columns of U (the left singular vectors, stored columnwise); if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced. LDU (input) INTEGER The leading dimension of the array U. LDU >= 1; if JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. VT (output) DOUBLE PRECISION array, dimension (LDVT,N) If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the N-by-N orthogonal matrix V**T; if JOBZ = 'S', VT contains the first min(M,N) rows of V**T (the right singular vectors, stored rowwise); if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced. LDVT (input) INTEGER The leading dimension of the array VT. LDVT >= 1; if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; if JOBZ = 'S', LDVT >= min(M,N). WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK; LWORK (input) INTEGER The dimension of the array WORK. LWORK >= 1. If JOBZ = 'N', LWORK >= 3*min(M,N) + max(max(M,N),7*min(M,N)). If JOBZ = 'O', LWORK >= 3*min(M,N) + max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)). If JOBZ = 'S' or 'A' LWORK >= 3*min(M,N) + max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)). For good performance, LWORK should generally be larger. If LWORK = -1 but other input arguments are legal, WORK(1) returns the optimal LWORK. IWORK (workspace) INTEGER array, dimension (8*min(M,N)) INFO (output) INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: DBDSDC did not converge, updating process failed. Further Details =============== Based on contributions by Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, :] |
Name | Description |
---|---|
sigma[min(size(A, 1), size(A, 2))] | |
U[size(A, 1), size(A, 1)] | |
VT[size(A, 2), size(A, 2)] | |
info |
Compute generalized eigenvalues, as well as the left and right eigenvectors for a (A,B) system
Lapack documentation Purpose ======= DGGEV computes for a pair of N-by-N real nonsymmetric matrices (A,B) the generalized eigenvalues, and optionally, the left and/or right generalized eigenvectors. A generalized eigenvalue for a pair of matrices (A,B) is a scalar lambda or a ratio alpha/beta = lambda, such that A - lambda*B is singular. It is usually represented as the pair (alpha,beta), as there is a reasonable interpretation for beta=0, and even for both being zero. The right eigenvector v(j) corresponding to the eigenvalue lambda(j) of (A,B) satisfies A * v(j) = lambda(j) * B * v(j). The left eigenvector u(j) corresponding to the eigenvalue lambda(j) of (A,B) satisfies u(j)**H * A = lambda(j) * u(j)**H * B . where u(j)**H is the conjugate-transpose of u(j). Arguments ========= JOBVL (input) CHARACTER*1 = 'N': do not compute the left generalized eigenvectors; = 'V': compute the left generalized eigenvectors. JOBVR (input) CHARACTER*1 = 'N': do not compute the right generalized eigenvectors; = 'V': compute the right generalized eigenvectors. N (input) INTEGER The order of the matrices A, B, VL, and VR. N >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA, N) On entry, the matrix A in the pair (A,B). On exit, A has been overwritten. LDA (input) INTEGER The leading dimension of A. LDA >= max(1,N). B (input/output) DOUBLE PRECISION array, dimension (LDB, N) On entry, the matrix B in the pair (A,B). On exit, B has been overwritten. LDB (input) INTEGER The leading dimension of B. LDB >= max(1,N). ALPHAR (output) DOUBLE PRECISION array, dimension (N) ALPHAI (output) DOUBLE PRECISION array, dimension (N) BETA (output) DOUBLE PRECISION array, dimension (N) On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will be the generalized eigenvalues. If ALPHAI(j) is zero, then the j-th eigenvalue is real; if positive, then the j-th and (j+1)-st eigenvalues are a complex conjugate pair, with ALPHAI(j+1) negative. Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j) may easily over- or underflow, and BETA(j) may even be zero. Thus, the user should avoid naively computing the ratio alpha/beta. However, ALPHAR and ALPHAI will be always less than and usually comparable with norm(A) in magnitude, and BETA always less than and usually comparable with norm(B). VL (output) DOUBLE PRECISION array, dimension (LDVL,N) If JOBVL = 'V', the left eigenvectors u(j) are stored one after another in the columns of VL, in the same order as their eigenvalues. If the j-th eigenvalue is real, then u(j) = VL(:,j), the j-th column of VL. If the j-th and (j+1)-th eigenvalues form a complex conjugate pair, then u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1). Each eigenvector is scaled so the largest component has abs(real part)+abs(imag. part)=1. Not referenced if JOBVL = 'N'. LDVL (input) INTEGER The leading dimension of the matrix VL. LDVL >= 1, and if JOBVL = 'V', LDVL >= N. VR (output) DOUBLE PRECISION array, dimension (LDVR,N) If JOBVR = 'V', the right eigenvectors v(j) are stored one after another in the columns of VR, in the same order as their eigenvalues. If the j-th eigenvalue is real, then v(j) = VR(:,j), the j-th column of VR. If the j-th and (j+1)-th eigenvalues form a complex conjugate pair, then v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1). Each eigenvector is scaled so the largest component has abs(real part)+abs(imag. part)=1. Not referenced if JOBVR = 'N'. LDVR (input) INTEGER The leading dimension of the matrix VR. LDVR >= 1, and if JOBVR = 'V', LDVR >= N. WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) INTEGER The dimension of the array WORK. LWORK >= max(1,8*N). For good performance, LWORK must generally be larger. 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 failed. No eigenvectors have been calculated, but ALPHAR(j), ALPHAI(j), and BETA(j) should be correct for j=INFO+1,...,N. > N: =N+1: other than QZ iteration failed in DHGEQZ. =N+2: error return from DTGEVC.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, size(A, 1)] | |
B[size(A, 1), size(A, 1)] | |
nA | The actual dimensions of matrices A and B (the computation is performed for A[1:nA,1:nA], B[1:nA,1:nA]) |
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 |
lEigenVectors[size(A, 1), size(A, 1)] | left eigenvectors of matrix A |
rEigenVectors[size(A, 1), size(A, 1)] | right eigenvectors of matrix A |
info |
Compute generalized eigenvalues for a (A,B) system, using lapack routine dggevx
Lapack documentation Purpose ======= DGGEVX computes for a pair of N-by-N real nonsymmetric matrices (A,B) the generalized eigenvalues, and optionally, the left and/or right generalized eigenvectors. Optionally also, it computes a balancing transformation to improve the conditioning of the eigenvalues and eigenvectors (ILO, IHI, LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for the eigenvalues (RCONDE), and reciprocal condition numbers for the right eigenvectors (RCONDV). A generalized eigenvalue for a pair of matrices (A,B) is a scalar lambda or a ratio alpha/beta = lambda, such that A - lambda*B is singular. It is usually represented as the pair (alpha,beta), as there is a reasonable interpretation for beta=0, and even for both being zero. The right eigenvector v(j) corresponding to the eigenvalue lambda(j) of (A,B) satisfies A * v(j) = lambda(j) * B * v(j) . The left eigenvector u(j) corresponding to the eigenvalue lambda(j) of (A,B) satisfies u(j)**H * A = lambda(j) * u(j)**H * B. where u(j)**H is the conjugate-transpose of u(j). Arguments ========= BALANC (input) CHARACTER*1 Specifies the balance option to be performed. = 'N': do not diagonally scale or permute; = 'P': permute only; = 'S': scale only; = 'B': both permute and scale. Computed reciprocal condition numbers will be for the matrices after permuting and/or balancing. Permuting does not change condition numbers (in exact arithmetic), but balancing does. JOBVL (input) CHARACTER*1 = 'N': do not compute the left generalized eigenvectors; = 'V': compute the left generalized eigenvectors. JOBVR (input) CHARACTER*1 = 'N': do not compute the right generalized eigenvectors; = 'V': compute the right generalized eigenvectors. SENSE (input) CHARACTER*1 Determines which reciprocal condition numbers are computed. = 'N': none are computed; = 'E': computed for eigenvalues only; = 'V': computed for eigenvectors only; = 'B': computed for eigenvalues and eigenvectors. N (input) INTEGER The order of the matrices A, B, VL, and VR. N >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA, N) On entry, the matrix A in the pair (A,B). On exit, A has been overwritten. If JOBVL='V' or JOBVR='V' or both, then A contains the first part of the real Schur form of the "balanced" versions of the input A and B. LDA (input) INTEGER The leading dimension of A. LDA >= max(1,N). B (input/output) DOUBLE PRECISION array, dimension (LDB, N) On entry, the matrix B in the pair (A,B). On exit, B has been overwritten. If JOBVL='V' or JOBVR='V' or both, then B contains the second part of the real Schur form of the "balanced" versions of the input A and B. LDB (input) INTEGER The leading dimension of B. LDB >= max(1,N). ALPHAR (output) DOUBLE PRECISION array, dimension (N) ALPHAI (output) DOUBLE PRECISION array, dimension (N) BETA (output) DOUBLE PRECISION array, dimension (N) On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will be the generalized eigenvalues. If ALPHAI(j) is zero, then the j-th eigenvalue is real; if positive, then the j-th and (j+1)-st eigenvalues are a complex conjugate pair, with ALPHAI(j+1) negative. Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j) may easily over- or underflow, and BETA(j) may even be zero. Thus, the user should avoid naively computing the ratio ALPHA/BETA. However, ALPHAR and ALPHAI will be always less than and usually comparable with norm(A) in magnitude, and BETA always less than and usually comparable with norm(B). VL (output) DOUBLE PRECISION array, dimension (LDVL,N) If JOBVL = 'V', the left eigenvectors u(j) are stored one after another in the columns of VL, in the same order as their eigenvalues. If the j-th eigenvalue is real, then u(j) = VL(:,j), the j-th column of VL. If the j-th and (j+1)-th eigenvalues form a complex conjugate pair, then u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1). Each eigenvector will be scaled so the largest component have abs(real part) + abs(imag. part) = 1. Not referenced if JOBVL = 'N'. LDVL (input) INTEGER The leading dimension of the matrix VL. LDVL >= 1, and if JOBVL = 'V', LDVL >= N. VR (output) DOUBLE PRECISION array, dimension (LDVR,N) If JOBVR = 'V', the right eigenvectors v(j) are stored one after another in the columns of VR, in the same order as their eigenvalues. If the j-th eigenvalue is real, then v(j) = VR(:,j), the j-th column of VR. If the j-th and (j+1)-th eigenvalues form a complex conjugate pair, then v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1). Each eigenvector will be scaled so the largest component have abs(real part) + abs(imag. part) = 1. Not referenced if JOBVR = 'N'. LDVR (input) INTEGER The leading dimension of the matrix VR. LDVR >= 1, and if JOBVR = 'V', LDVR >= N. ILO (output) INTEGER IHI (output) INTEGER ILO and IHI are integer values such that on exit A(i,j) = 0 and B(i,j) = 0 if i > j and j = 1,...,ILO-1 or i = IHI+1,...,N. If BALANC = 'N' or 'S', ILO = 1 and IHI = N. LSCALE (output) DOUBLE PRECISION array, dimension (N) Details of the permutations and scaling factors applied to the left side of A and B. If PL(j) is the index of the row interchanged with row j, and DL(j) is the scaling factor applied to row j, then LSCALE(j) = PL(j) for j = 1,...,ILO-1 = DL(j) for j = ILO,...,IHI = PL(j) for j = IHI+1,...,N. The order in which the interchanges are made is N to IHI+1, then 1 to ILO-1. RSCALE (output) DOUBLE PRECISION array, dimension (N) Details of the permutations and scaling factors applied to the right side of A and B. If PR(j) is the index of the column interchanged with column j, and DR(j) is the scaling factor applied to column j, then RSCALE(j) = PR(j) for j = 1,...,ILO-1 = DR(j) for j = ILO,...,IHI = PR(j) for j = IHI+1,...,N The order in which the interchanges are made is N to IHI+1, then 1 to ILO-1. ABNRM (output) DOUBLE PRECISION The one-norm of the balanced matrix A. BBNRM (output) DOUBLE PRECISION The one-norm of the balanced matrix B. RCONDE (output) DOUBLE PRECISION array, dimension (N) If SENSE = 'E' or 'B', the reciprocal condition numbers of the eigenvalues, stored in consecutive elements of the array. For a complex conjugate pair of eigenvalues two consecutive elements of RCONDE are set to the same value. Thus RCONDE(j), RCONDV(j), and the j-th columns of VL and VR all correspond to the j-th eigenpair. If SENSE = 'N or 'V', RCONDE is not referenced. RCONDV (output) DOUBLE PRECISION array, dimension (N) If SENSE = 'V' or 'B', the estimated reciprocal condition numbers of the eigenvectors, stored in consecutive elements of the array. For a complex eigenvector two consecutive elements of RCONDV are set to the same value. If the eigenvalues cannot be reordered to compute RCONDV(j), RCONDV(j) is set to 0; this can only occur when the true value would be very small anyway. If SENSE = 'N' or 'E', RCONDV is not referenced. WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) INTEGER The dimension of the array WORK. LWORK >= max(1,2*N). If BALANC = 'S' or 'B', or JOBVL = 'V', or JOBVR = 'V', LWORK >= max(1,6*N). If SENSE = 'E' or 'B', LWORK >= max(1,10*N). If SENSE = 'V' or 'B', LWORK >= 2*N*N+8*N+16. 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. IWORK (workspace) INTEGER array, dimension (N+6) If SENSE = 'E', IWORK is not referenced. BWORK (workspace) LOGICAL array, dimension (N) If SENSE = 'N', BWORK is not referenced. INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value. = 1,...,N: The QZ iteration failed. No eigenvectors have been calculated, but ALPHAR(j), ALPHAI(j), and BETA(j) should be correct for j=INFO+1,...,N. > N: =N+1: other than QZ iteration failed in DHGEQZ. =N+2: error return from DTGEVC. Further Details =============== Balancing a matrix pair (A,B) includes, first, permuting rows and columns to isolate eigenvalues, second, applying diagonal similarity transformation to the rows and columns to make the rows and columns as close in norm as possible. The computed reciprocal condition numbers correspond to the balanced matrix. Permuting rows and columns will not change the condition numbers (in exact arithmetic) but diagonal scaling will. For further explanation of balancing, see section 4.11.1.2 of LAPACK Users' Guide. An approximate error bound on the chordal distance between the i-th computed generalized eigenvalue w and the corresponding exact eigenvalue lambda is chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I) An approximate error bound for the angle between the i-th computed eigenvector VL(i) or VR(i) is given by EPS * norm(ABNRM, BBNRM) / DIF(i). For further explanation of the reciprocal condition numbers RCONDE and RCONDV, see section 4.11 of LAPACK User's Guide.
Extends from Modelica.Icons.Function (Icon for functions).
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 |
lEigenVectors[size(A, 1), size(A, 1)] | left eigenvectors of matrix A |
rEigenVectors[size(A, 1), size(A, 1)] | right eigenvectors of matrix A |
info |
Compute generalized eigenvalues for a (A,B) system
Lapack documentation Purpose ======= DHGEQZ computes the eigenvalues of a real matrix pair (H,T), where H is an upper Hessenberg matrix and T is upper triangular, using the double-shift QZ method. Matrix pairs of this type are produced by the reduction to generalized upper Hessenberg form of a real matrix pair (A,B): A = Q1*H*Z1**T, B = Q1*T*Z1**T, as computed by DGGHRD. If JOB='S', then the Hessenberg-triangular pair (H,T) is also reduced to generalized Schur form, H = Q*S*Z**T, T = Q*P*Z**T, where Q and Z are orthogonal matrices, P is an upper triangular matrix, and S is a quasi-triangular matrix with 1-by-1 and 2-by-2 diagonal blocks. The 1-by-1 blocks correspond to real eigenvalues of the matrix pair (H,T) and the 2-by-2 blocks correspond to complex conjugate pairs of eigenvalues. Additionally, the 2-by-2 upper triangular diagonal blocks of P corresponding to 2-by-2 blocks of S are reduced to positive diagonal form, i.e., if S(j+1,j) is non-zero, then P(j+1,j) = P(j,j+1) = 0, P(j,j) > 0, and P(j+1,j+1) > 0. Optionally, the orthogonal matrix Q from the generalized Schur factorization may be postmultiplied into an input matrix Q1, and the orthogonal matrix Z may be postmultiplied into an input matrix Z1. If Q1 and Z1 are the orthogonal matrices from DGGHRD that reduced the matrix pair (A,B) to generalized upper Hessenberg form, then the output matrices Q1*Q and Z1*Z are the orthogonal factors from the generalized Schur factorization of (A,B): A = (Q1*Q)*S*(Z1*Z)**T, B = (Q1*Q)*P*(Z1*Z)**T. To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently, of (A,B)) are computed as a pair of values (alpha,beta), where alpha is complex and beta real. If beta is nonzero, lambda = alpha / beta is an eigenvalue of the generalized nonsymmetric eigenvalue problem (GNEP) A*x = lambda*B*x and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the alternate form of the GNEP mu*A*y = B*y. Real eigenvalues can be read directly from the generalized Schur form: alpha = S(i,i), beta = P(i,i). 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 eigenvalues only; = 'S': Compute eigenvalues and the Schur form. COMPQ (input) CHARACTER*1 = 'N': Left Schur vectors (Q) are not computed; = 'I': Q is initialized to the unit matrix and the matrix Q of left Schur vectors of (H,T) is returned; = 'V': Q must contain an orthogonal matrix Q1 on entry and the product Q1*Q is returned. COMPZ (input) CHARACTER*1 = 'N': Right Schur vectors (Z) are not computed; = 'I': Z is initialized to the unit matrix and the matrix Z of right Schur vectors of (H,T) is returned; = 'V': Z must contain an orthogonal matrix Z1 on entry and the product Z1*Z is returned. N (input) INTEGER The order of the matrices H, T, Q, and Z. N >= 0. ILO (input) INTEGER IHI (input) INTEGER ILO and IHI mark the rows and columns of H which are in Hessenberg form. It is assumed that A is already upper triangular in rows and columns 1:ILO-1 and IHI+1:N. If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0. H (input/output) DOUBLE PRECISION array, dimension (LDH, N) On entry, the N-by-N upper Hessenberg matrix H. On exit, if JOB = 'S', H contains the upper quasi-triangular matrix S from the generalized Schur factorization; 2-by-2 diagonal blocks (corresponding to complex conjugate pairs of eigenvalues) are returned in standard form, with H(i,i) = H(i+1,i+1) and H(i+1,i)*H(i,i+1) < 0. If JOB = 'E', the diagonal blocks of H match those of S, but the rest of H is unspecified. LDH (input) INTEGER The leading dimension of the array H. LDH >= max( 1, N ). T (input/output) DOUBLE PRECISION array, dimension (LDT, N) On entry, the N-by-N upper triangular matrix T. On exit, if JOB = 'S', T contains the upper triangular matrix P from the generalized Schur factorization; 2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks of S are reduced to positive diagonal form, i.e., if H(j+1,j) is non-zero, then T(j+1,j) = T(j,j+1) = 0, T(j,j) > 0, and T(j+1,j+1) > 0. If JOB = 'E', the diagonal blocks of T match those of P, but the rest of T is unspecified. LDT (input) INTEGER The leading dimension of the array T. LDT >= max( 1, N ). ALPHAR (output) DOUBLE PRECISION array, dimension (N) The real parts of each scalar alpha defining an eigenvalue of GNEP. ALPHAI (output) DOUBLE PRECISION array, dimension (N) The imaginary parts of each scalar alpha defining an eigenvalue of GNEP. If ALPHAI(j) is zero, then the j-th eigenvalue is real; if positive, then the j-th and (j+1)-st eigenvalues are a complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j). BETA (output) DOUBLE PRECISION array, dimension (N) The scalars beta that define the eigenvalues of GNEP. Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and beta = BETA(j) represent the j-th eigenvalue of the matrix pair (A,B), in one of the forms lambda = alpha/beta or mu = beta/alpha. Since either lambda or mu may overflow, they should not, in general, be computed. Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N) On entry, if COMPZ = 'V', the orthogonal matrix Q1 used in the reduction of (A,B) to generalized Hessenberg form. On exit, if COMPZ = 'I', the orthogonal matrix of left Schur vectors of (H,T), and if COMPZ = 'V', the orthogonal matrix of left Schur vectors of (A,B). Not referenced if COMPZ = 'N'. 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) On entry, if COMPZ = 'V', the orthogonal matrix Z1 used in the reduction of (A,B) to generalized Hessenberg form. On exit, if COMPZ = 'I', the orthogonal matrix of right Schur vectors of (H,T), and if COMPZ = 'V', the orthogonal matrix of right Schur vectors of (A,B). Not referenced if COMPZ = 'N'. 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 (MAX(1,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. (H,T) 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. (H,T) is not in Schur form, but ALPHAR(i), ALPHAI(i), and BETA(i), i=INFO-N+1,...,N should be correct. 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.
Extends from Modelica.Icons.Function (Icon for functions).
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 |
overwrites the general real M-by-N matrix C with Q * C or C * Q or Q' * C or C * Q', where Q is an orthogonal matrix as returned by dgehrd
Lapack documentation Purpose ======= DORMHR overwrites the general real M-by-N matrix C with SIDE = 'L' SIDE = 'R' TRANS = 'N': Q * C C * Q TRANS = 'T': Q**T * C C * Q**T where Q is a real orthogonal matrix of order nq, with nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Q is defined as the product of IHI-ILO elementary reflectors, as returned by DGEHRD: Q = H(ilo) H(ilo+1) . . . H(ihi-1). Arguments ========= SIDE (input) CHARACTER*1 = 'L': apply Q or Q**T from the Left; = 'R': apply Q or Q**T from the Right. TRANS (input) CHARACTER*1 = 'N': No transpose, apply Q; = 'T': Transpose, apply Q**T. M (input) INTEGER The number of rows of the matrix C. M >= 0. N (input) INTEGER The number of columns of the matrix C. 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). If SIDE = 'L', then 1 <= ILO <= IHI <= M, if M > 0, and ILO = 1 and IHI = 0, if M = 0; if SIDE = 'R', then 1 <= ILO <= IHI <= N, if N > 0, and ILO = 1 and IHI = 0, if N = 0. A (input) DOUBLE PRECISION array, dimension (LDA,M) if SIDE = 'L' (LDA,N) if SIDE = 'R' The vectors which define the elementary reflectors, as returned by DGEHRD. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,M) if SIDE = 'L'; LDA >= max(1,N) if SIDE = 'R'. TAU (input) DOUBLE PRECISION array, dimension (M-1) if SIDE = 'L' (N-1) if SIDE = 'R' TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by DGEHRD. C (input/output) DOUBLE PRECISION array, dimension (LDC,N) On entry, the M-by-N matrix C. On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q. LDC (input) INTEGER The leading dimension of the array C. LDC >= max(1,M). WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) INTEGER The dimension of the array WORK. If SIDE = 'L', LWORK >= max(1,N); if SIDE = 'R', LWORK >= max(1,M). For optimum performance LWORK >= N*NB if SIDE = 'L', and LWORK >= M*NB if SIDE = 'R', 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
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
C[:, :] | |
A[:, :] | |
tau[if side == "L" then size(C, 2) - 1 else size(C, 1) - 1] | |
side | |
trans | |
ilo | lowest index where the original matrix had been Hessenbergform |
ihi | highest index where the original matrix had been Hessenbergform |
Name | Description |
---|---|
Cout[size(C, 1), size(C, 2)] | contains the Hessenberg form in the upper triangle and the first subdiagonal and below the first subdiagonal it contains the elementary reflectors which represents (with array tau) as a product the orthogonal matrix Q |
info |
overwrites the general real M-by-N matrix C with Q * C or C * Q or Q' * C or C * Q', where Q is an orthogonal matrix of a QR factorization as returned by dgeqrf
Lapack documentation Purpose ======= DORMQR overwrites the general real M-by-N matrix C with SIDE = 'L' SIDE = 'R' TRANS = 'N': Q * C C * Q TRANS = 'T': Q**T * C C * Q**T where Q is a real orthogonal matrix defined as the product of k elementary reflectors Q = H(1) H(2) . . . H(k) as returned by DGEQRF. Q is of order M if SIDE = 'L' and of order N if SIDE = 'R'. Arguments ========= SIDE (input) CHARACTER*1 = 'L': apply Q or Q**T from the Left; = 'R': apply Q or Q**T from the Right. TRANS (input) CHARACTER*1 = 'N': No transpose, apply Q; = 'T': Transpose, apply Q**T. M (input) INTEGER The number of rows of the matrix C. M >= 0. N (input) INTEGER The number of columns of the matrix C. N >= 0. K (input) INTEGER The number of elementary reflectors whose product defines the matrix Q. If SIDE = 'L', M >= K >= 0; if SIDE = 'R', N >= K >= 0. A (input) DOUBLE PRECISION array, dimension (LDA,K) The i-th column must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by DGEQRF in the first k columns of its array argument A. A is modified by the routine but restored on exit. LDA (input) INTEGER The leading dimension of the array A. If SIDE = 'L', LDA >= max(1,M); if SIDE = 'R', LDA >= max(1,N). TAU (input) DOUBLE PRECISION array, dimension (K) TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by DGEQRF. C (input/output) DOUBLE PRECISION array, dimension (LDC,N) On entry, the M-by-N matrix C. On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q. LDC (input) INTEGER The leading dimension of the array C. LDC >= max(1,M). WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) INTEGER The dimension of the array WORK. If SIDE = 'L', LWORK >= max(1,N); if SIDE = 'R', LWORK >= max(1,M). For optimum performance LWORK >= N*NB if SIDE = 'L', and LWORK >= M*NB if SIDE = 'R', 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
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
C[:, :] | |
A[:, :] | |
tau[:] | |
side | |
trans |
Name | Description |
---|---|
Cout[size(C, 1), size(C, 2)] | contains Q*C or Q**T*C or C*Q**T or C*Q |
info |
Compute the right and/or left eigenvectors of a real upper quasi-triangular matrix T
Lapack documentation Purpose ======= DTREVC computes some or all of the right and/or left eigenvectors of a real upper quasi-triangular matrix T. Matrices of this type are produced by the Schur factorization of a real general matrix: A = Q*T*Q**T, as computed by DHSEQR. The right eigenvector x and the left eigenvector y of T corresponding to an eigenvalue w are defined by: T*x = w*x, (y**H)*T = w*(y**H) where y**H denotes the conjugate transpose of y. The eigenvalues are not input to this routine, but are read directly from the diagonal blocks of T. This routine returns the matrices X and/or Y of right and left eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an input matrix. If Q is the orthogonal factor that reduces a matrix A to Schur form T, then Q*X and Q*Y are the matrices of right and left eigenvectors of A. Arguments ========= SIDE (input) CHARACTER*1 = 'R': compute right eigenvectors only; = 'L': compute left eigenvectors only; = 'B': compute both right and left eigenvectors. HOWMNY (input) CHARACTER*1 = 'A': compute all right and/or left eigenvectors; = 'B': compute all right and/or left eigenvectors, backtransformed by the matrices in VR and/or VL; = 'S': compute selected right and/or left eigenvectors, as indicated by the logical array SELECT. SELECT (input/output) LOGICAL array, dimension (N) If HOWMNY = 'S', SELECT specifies the eigenvectors to be computed. If w(j) is a real eigenvalue, the corresponding real eigenvector is computed if SELECT(j) is .TRUE.. If w(j) and w(j+1) are the real and imaginary parts of a complex eigenvalue, the corresponding complex eigenvector is computed if either SELECT(j) or SELECT(j+1) is .TRUE., and on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to .FALSE.. Not referenced if HOWMNY = 'A' or 'B'. N (input) INTEGER The order of the matrix T. N >= 0. T (input) DOUBLE PRECISION array, dimension (LDT,N) The upper quasi-triangular matrix T in Schur canonical form. LDT (input) INTEGER The leading dimension of the array T. LDT >= max(1,N). VL (input/output) DOUBLE PRECISION array, dimension (LDVL,MM) On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must contain an N-by-N matrix Q (usually the orthogonal matrix Q of Schur vectors returned by DHSEQR). On exit, if SIDE = 'L' or 'B', VL contains: if HOWMNY = 'A', the matrix Y of left eigenvectors of T; if HOWMNY = 'B', the matrix Q*Y; if HOWMNY = 'S', the left eigenvectors of T specified by SELECT, stored consecutively in the columns of VL, in the same order as their eigenvalues. A complex eigenvector corresponding to a complex eigenvalue is stored in two consecutive columns, the first holding the real part, and the second the imaginary part. Not referenced if SIDE = 'R'. LDVL (input) INTEGER The leading dimension of the array VL. LDVL >= 1, and if SIDE = 'L' or 'B', LDVL >= N. VR (input/output) DOUBLE PRECISION array, dimension (LDVR,MM) On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must contain an N-by-N matrix Q (usually the orthogonal matrix Q of Schur vectors returned by DHSEQR). On exit, if SIDE = 'R' or 'B', VR contains: if HOWMNY = 'A', the matrix X of right eigenvectors of T; if HOWMNY = 'B', the matrix Q*X; if HOWMNY = 'S', the right eigenvectors of T specified by SELECT, stored consecutively in the columns of VR, in the same order as their eigenvalues. A complex eigenvector corresponding to a complex eigenvalue is stored in two consecutive columns, the first holding the real part and the second the imaginary part. Not referenced if SIDE = 'L'. LDVR (input) INTEGER The leading dimension of the array VR. LDVR >= 1, and if SIDE = 'R' or 'B', LDVR >= N. MM (input) INTEGER The number of columns in the arrays VL and/or VR. MM >= M. M (output) INTEGER The number of columns in the arrays VL and/or VR actually used to store the eigenvectors. If HOWMNY = 'A' or 'B', M is set to N. Each selected real eigenvector occupies one column and each selected complex eigenvector occupies two columns. WORK (workspace) DOUBLE PRECISION array, dimension (3*N) INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value Further Details =============== The algorithm used in this program is basically backward (forward) substitution, with scaling to make the code robust against possible overflow. Each eigenvector is normalized so that the element of largest magnitude has magnitude 1; here the magnitude of a complex number (x,y) is taken to be |x| + |y|.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
T[:, size(T, 1)] | Upper quasi triangular matrix |
side | Specify which eigenvectors |
howmny | Specify how many eigenvectors |
Q[size(T, 1), size(T, 1)] | Orthogonal matrix Q of Schur vectors returned by DHSEQR |
Name | Description |
---|---|
lEigenVectors[size(T, 1), size(T, 1)] | left eigenvectors of matrix T |
rEigenVectors[size(T, 1), size(T, 1)] | right eigenvectors of matrix T |
info |
Computes the Cholesky factorization of a real symmetric positive definite matrix A
Lapack documentation Purpose ======= DPOTRF computes the Cholesky factorization of a real symmetric positive definite matrix A. The factorization has the form A = U**T * U, if UPLO = 'U', or A = L * L**T, if UPLO = 'L', where U is an upper triangular matrix and L is lower triangular. This is the block version of the algorithm, calling Level 3 BLAS. Arguments ========= UPLO (input) CHARACTER*1 = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A is stored. N (input) INTEGER The order of the matrix A. N >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the symmetric matrix A. If UPLO = 'U', the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. On exit, if INFO = 0, the factor U or L from the Cholesky factorization A = U**T*U or A = L*L**T. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,N). INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, the leading minor of order i is not positive definite, and the factorization could not be completed.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, size(A, 1)] | Real symmetric positive definite matrix A |
upper | True if the upper triangle of A is provided |
Name | Description |
---|---|
Acholesky[size(A, 1), size(A, 1)] | Cholesky factor |
info |
Solve one of the matrix equations op( A )*X = alpha*B, or X*op( A ) = alpha*B, where A is triangular matrix. BLAS routine
Lapack documentation Purpose ======= DTRSM solves one of the matrix equations op( A )*X = alpha*B, or X*op( A ) = alpha*B, where alpha is a scalar, X and B are m by n matrices, A is a unit, or non-unit, upper or lower triangular matrix and op( A ) is one of op( A ) = A or op( A ) = A'. The matrix X is overwritten on B. Arguments ========== SIDE - CHARACTER*1. On entry, SIDE specifies whether op( A ) appears on the left or right of X as follows: SIDE = 'L' or 'l' op( A )*X = alpha*B. SIDE = 'R' or 'r' X*op( A ) = alpha*B. Unchanged on exit. UPLO - CHARACTER*1. On entry, UPLO specifies whether the matrix A is an upper or lower triangular matrix as follows: UPLO = 'U' or 'u' A is an upper triangular matrix. UPLO = 'L' or 'l' A is a lower triangular matrix. Unchanged on exit. TRANSA - CHARACTER*1. On entry, TRANSA specifies the form of op( A ) to be used in the matrix multiplication as follows: TRANSA = 'N' or 'n' op( A ) = A. TRANSA = 'T' or 't' op( A ) = A'. TRANSA = 'C' or 'c' op( A ) = A'. Unchanged on exit. DIAG - CHARACTER*1. On entry, DIAG specifies whether or not A is unit triangular as follows: DIAG = 'U' or 'u' A is assumed to be unit triangular. DIAG = 'N' or 'n' A is not assumed to be unit triangular. Unchanged on exit. M - INTEGER. On entry, M specifies the number of rows of B. M must be at least zero. Unchanged on exit. N - INTEGER. On entry, N specifies the number of columns of B. N must be at least zero. Unchanged on exit. ALPHA - DOUBLE PRECISION. On entry, ALPHA specifies the scalar alpha. When alpha is zero then A is not referenced and B need not be set before entry. Unchanged on exit. A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. Before entry with UPLO = 'U' or 'u', the leading k by k upper triangular part of the array A must contain the upper triangular matrix and the strictly lower triangular part of A is not referenced. Before entry with UPLO = 'L' or 'l', the leading k by k lower triangular part of the array A must contain the lower triangular matrix and the strictly upper triangular part of A is not referenced. Note that when DIAG = 'U' or 'u', the diagonal elements of A are not referenced either, but are assumed to be unity. Unchanged on exit. LDA - INTEGER. On entry, LDA specifies the first dimension of A as declared in the calling (sub) program. When SIDE = 'L' or 'l' then LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' then LDA must be at least max( 1, n ). Unchanged on exit. B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). Before entry, the leading m by n part of the array B must contain the right-hand side matrix B, and on exit is overwritten by the solution matrix X. LDB - INTEGER. On entry, LDB specifies the first dimension of B as declared in the calling (sub) program. LDB must be at least max( 1, m ). Unchanged on exit. Level 3 Blas routine.
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, :] | Input matrix A |
B[:, :] | Input matrix B |
alpha | Factor alpha |
right | True if A is right multiplication |
upper | True if A is upper triangular |
trans | True if op(A) means transposed(A) |
unitTriangular | True if A is unit triangular, i.e., all diagonal elements of A are equal to 1 |
Name | Description |
---|---|
X[size(B, 1), size(B, 2)] | Matrix Bout=alpha*op( A )*B, or B := alpha*B*op( A ) |
Generates a real orthogonal matrix Q which is defined as the product of IHI-ILO elementary reflectors of order N, as returned by DGEHRD
Lapack documentation 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 (MAX(1,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
Extends from Modelica.Icons.Function (Icon for functions).
Name | Description |
---|---|
A[:, size(A, 1)] | Square matrix with 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 |
tau[max(0, size(A, 1) - 1)] | scalar factors of the elementary reflectors |
Name | Description |
---|---|
Aout[size(A, 1), size(A, 2)] | Orthogonal matrix as a result of elementary reflectors |
info |