 |
LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
|
|
LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
|
| subroutine dgedmdq | ( | character, intent(in) | jobs, |
| character, intent(in) | jobz, | ||
| character, intent(in) | jobr, | ||
| character, intent(in) | jobq, | ||
| character, intent(in) | jobt, | ||
| character, intent(in) | jobf, | ||
| integer, intent(in) | whtsvd, | ||
| integer, intent(in) | m, | ||
| integer, intent(in) | n, | ||
| real(kind=wp), dimension(ldf,*), intent(inout) | f, | ||
| integer, intent(in) | ldf, | ||
| real(kind=wp), dimension(ldx,*), intent(out) | x, | ||
| integer, intent(in) | ldx, | ||
| real(kind=wp), dimension(ldy,*), intent(out) | y, | ||
| integer, intent(in) | ldy, | ||
| integer, intent(in) | nrnk, | ||
| real(kind=wp), intent(in) | tol, | ||
| integer, intent(out) | k, | ||
| real(kind=wp), dimension(*), intent(out) | reig, | ||
| real(kind=wp), dimension(*), intent(out) | imeig, | ||
| real(kind=wp), dimension(ldz,*), intent(out) | z, | ||
| integer, intent(in) | ldz, | ||
| real(kind=wp), dimension(*), intent(out) | res, | ||
| real(kind=wp), dimension(ldb,*), intent(out) | b, | ||
| integer, intent(in) | ldb, | ||
| real(kind=wp), dimension(ldv,*), intent(out) | v, | ||
| integer, intent(in) | ldv, | ||
| real(kind=wp), dimension(lds,*), intent(out) | s, | ||
| integer, intent(in) | lds, | ||
| real(kind=wp), dimension(*), intent(out) | work, | ||
| integer, intent(in) | lwork, | ||
| integer, dimension(*), intent(out) | iwork, | ||
| integer, intent(in) | liwork, | ||
| integer, intent(out) | info | ||
| ) |
DGEDMDQ computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.
DGEDMDQ computes the Dynamic Mode Decomposition (DMD) for
a pair of data snapshot matrices, using a QR factorization
based compression of the data. For the input matrices
X and Y such that Y = A*X with an unaccessible matrix
A, DGEDMDQ computes a certain number of Ritz pairs of A using
the standard Rayleigh-Ritz extraction from a subspace of
range(X) that is determined using the leading left singular
vectors of X. Optionally, DGEDMDQ returns the residuals
of the computed Ritz pairs, the information needed for
a refinement of the Ritz vectors, or the eigenvectors of
the Exact DMD.
For further details see the references listed
below. For more details of the implementation see [3]. [1] P. Schmid: Dynamic mode decomposition of numerical
and experimental data,
Journal of Fluid Mechanics 656, 5-28, 2010.
[2] Z. Drmac, I. Mezic, R. Mohr: Data driven modal
decompositions: analysis and enhancements,
SIAM J. on Sci. Comp. 40 (4), A2253-A2285, 2018.
[3] Z. Drmac: A LAPACK implementation of the Dynamic
Mode Decomposition I. Technical report. AIMDyn Inc.
and LAPACK Working Note 298.
[4] J. Tu, C. W. Rowley, D. M. Luchtenburg, S. L.
Brunton, N. Kutz: On Dynamic Mode Decomposition:
Theory and Applications, Journal of Computational
Dynamics 1(2), 391 -421, 2014. Developed and coded by Zlatko Drmac, Faculty of Science,
University of Zagreb; drmac@math.hr
In cooperation with
AIMdyn Inc., Santa Barbara, CA.
and supported by
- DARPA SBIR project "Koopman Operator-Based Forecasting
for Nonstationary Processes from Near-Term, Limited
Observational Data" Contract No: W31P4Q-21-C-0007
- DARPA PAI project "Physics-Informed Machine Learning
Methodologies" Contract No: HR0011-18-9-0033
- DARPA MoDyL project "A Data-Driven, Operator-Theoretic
Framework for Space-Time Analysis of Process Dynamics"
Contract No: HR0011-16-C-0116
Any opinions, findings and conclusions or recommendations
expressed in this material are those of the author and
do not necessarily reflect the views of the DARPA SBIR
Program Office. Approved for Public Release, Distribution Unlimited.
Cleared by DARPA on September 29, 2022 | [in] | JOBS | JOBS (input) CHARACTER*1
Determines whether the initial data snapshots are scaled
by a diagonal matrix. The data snapshots are the columns
of F. The leading N-1 columns of F are denoted X and the
trailing N-1 columns are denoted Y.
'S' :: The data snapshots matrices X and Y are multiplied
with a diagonal matrix D so that X*D has unit
nonzero columns (in the Euclidean 2-norm)
'C' :: The snapshots are scaled as with the 'S' option.
If it is found that an i-th column of X is zero
vector and the corresponding i-th column of Y is
non-zero, then the i-th column of Y is set to
zero and a warning flag is raised.
'Y' :: The data snapshots matrices X and Y are multiplied
by a diagonal matrix D so that Y*D has unit
nonzero columns (in the Euclidean 2-norm)
'N' :: No data scaling. |
| [in] | JOBZ | JOBZ (input) CHARACTER*1
Determines whether the eigenvectors (Koopman modes) will
be computed.
'V' :: The eigenvectors (Koopman modes) will be computed
and returned in the matrix Z.
See the description of Z.
'F' :: The eigenvectors (Koopman modes) will be returned
in factored form as the product Z*V, where Z
is orthonormal and V contains the eigenvectors
of the corresponding Rayleigh quotient.
See the descriptions of F, V, Z.
'Q' :: The eigenvectors (Koopman modes) will be returned
in factored form as the product Q*Z, where Z
contains the eigenvectors of the compression of the
underlying discretized operator onto the span of
the data snapshots. See the descriptions of F, V, Z.
Q is from the initial QR factorization.
'N' :: The eigenvectors are not computed. |
| [in] | JOBR | JOBR (input) CHARACTER*1
Determines whether to compute the residuals.
'R' :: The residuals for the computed eigenpairs will
be computed and stored in the array RES.
See the description of RES.
For this option to be legal, JOBZ must be 'V'.
'N' :: The residuals are not computed. |
| [in] | JOBQ | JOBQ (input) CHARACTER*1
Specifies whether to explicitly compute and return the
orthogonal matrix from the QR factorization.
'Q' :: The matrix Q of the QR factorization of the data
snapshot matrix is computed and stored in the
array F. See the description of F.
'N' :: The matrix Q is not explicitly computed. |
| [in] | JOBT | JOBT (input) CHARACTER*1
Specifies whether to return the upper triangular factor
from the QR factorization.
'R' :: The matrix R of the QR factorization of the data
snapshot matrix F is returned in the array Y.
See the description of Y and Further details.
'N' :: The matrix R is not returned. |
| [in] | JOBF | JOBF (input) CHARACTER*1
Specifies whether to store information needed for post-
processing (e.g. computing refined Ritz vectors)
'R' :: The matrix needed for the refinement of the Ritz
vectors is computed and stored in the array B.
See the description of B.
'E' :: The unscaled eigenvectors of the Exact DMD are
computed and returned in the array B. See the
description of B.
'N' :: No eigenvector refinement data is computed.
To be useful on exit, this option needs JOBQ='Q'. |
| [in] | WHTSVD | WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 }
Allows for a selection of the SVD algorithm from the
LAPACK library.
1 :: DGESVD (the QR SVD algorithm)
2 :: DGESDD (the Divide and Conquer algorithm; if enough
workspace available, this is the fastest option)
3 :: DGESVDQ (the preconditioned QR SVD ; this and 4
are the most accurate options)
4 :: DGEJSV (the preconditioned Jacobi SVD; this and 3
are the most accurate options)
For the four methods above, a significant difference in
the accuracy of small singular values is possible if
the snapshots vary in norm so that X is severely
ill-conditioned. If small (smaller than EPS*||X||)
singular values are of interest and JOBS=='N', then
the options (3, 4) give the most accurate results, where
the option 4 is slightly better and with stronger
theoretical background.
If JOBS=='S', i.e. the columns of X will be normalized,
then all methods give nearly equally accurate results. |
| [in] | M | M (input) INTEGER, M >= 0
The state space dimension (the number of rows of F). |
| [in] | N | N (input) INTEGER, 0 <= N <= M
The number of data snapshots from a single trajectory,
taken at equidistant discrete times. This is the
number of columns of F. |
| [in,out] | F | F (input/output) REAL(KIND=WP) M-by-N array
> On entry,
the columns of F are the sequence of data snapshots
from a single trajectory, taken at equidistant discrete
times. It is assumed that the column norms of F are
in the range of the normalized floating point numbers.
< On exit,
If JOBQ == 'Q', the array F contains the orthogonal
matrix/factor of the QR factorization of the initial
data snapshots matrix F. See the description of JOBQ.
If JOBQ == 'N', the entries in F strictly below the main
diagonal contain, column-wise, the information on the
Householder vectors, as returned by DGEQRF. The
remaining information to restore the orthogonal matrix
of the initial QR factorization is stored in WORK(1:N).
See the description of WORK. |
| [in] | LDF | LDF (input) INTEGER, LDF >= M
The leading dimension of the array F. |
| [in,out] | X | X (workspace/output) REAL(KIND=WP) MIN(M,N)-by-(N-1) array
X is used as workspace to hold representations of the
leading N-1 snapshots in the orthonormal basis computed
in the QR factorization of F.
On exit, the leading K columns of X contain the leading
K left singular vectors of the above described content
of X. To lift them to the space of the left singular
vectors U(:,1:K)of the input data, pre-multiply with the
Q factor from the initial QR factorization.
See the descriptions of F, K, V and Z. |
| [in] | LDX | LDX (input) INTEGER, LDX >= N
The leading dimension of the array X. |
| [in,out] | Y | Y (workspace/output) REAL(KIND=WP) MIN(M,N)-by-(N-1) array
Y is used as workspace to hold representations of the
trailing N-1 snapshots in the orthonormal basis computed
in the QR factorization of F.
On exit,
If JOBT == 'R', Y contains the MIN(M,N)-by-N upper
triangular factor from the QR factorization of the data
snapshot matrix F. |
| [in] | LDY | LDY (input) INTEGER , LDY >= N
The leading dimension of the array Y. |
| [in] | NRNK | NRNK (input) INTEGER
Determines the mode how to compute the numerical rank,
i.e. how to truncate small singular values of the input
matrix X. On input, if
NRNK = -1 :: i-th singular value sigma(i) is truncated
if sigma(i) <= TOL*sigma(1)
This option is recommended.
NRNK = -2 :: i-th singular value sigma(i) is truncated
if sigma(i) <= TOL*sigma(i-1)
This option is included for R&D purposes.
It requires highly accurate SVD, which
may not be feasible.
The numerical rank can be enforced by using positive
value of NRNK as follows:
0 < NRNK <= N-1 :: at most NRNK largest singular values
will be used. If the number of the computed nonzero
singular values is less than NRNK, then only those
nonzero values will be used and the actually used
dimension is less than NRNK. The actual number of
the nonzero singular values is returned in the variable
K. See the description of K. |
| [in] | TOL | TOL (input) REAL(KIND=WP), 0 <= TOL < 1
The tolerance for truncating small singular values.
See the description of NRNK. |
| [out] | K | K (output) INTEGER, 0 <= K <= N
The dimension of the SVD/POD basis for the leading N-1
data snapshots (columns of F) and the number of the
computed Ritz pairs. The value of K is determined
according to the rule set by the parameters NRNK and
TOL. See the descriptions of NRNK and TOL. |
| [out] | REIG | REIG (output) REAL(KIND=WP) (N-1)-by-1 array
The leading K (K<=N) entries of REIG contain
the real parts of the computed eigenvalues
REIG(1:K) + sqrt(-1)*IMEIG(1:K).
See the descriptions of K, IMEIG, Z. |
| [out] | IMEIG | IMEIG (output) REAL(KIND=WP) (N-1)-by-1 array
The leading K (K<N) entries of REIG contain
the imaginary parts of the computed eigenvalues
REIG(1:K) + sqrt(-1)*IMEIG(1:K).
The eigenvalues are determined as follows:
If IMEIG(i) == 0, then the corresponding eigenvalue is
real, LAMBDA(i) = REIG(i).
If IMEIG(i)>0, then the corresponding complex
conjugate pair of eigenvalues reads
LAMBDA(i) = REIG(i) + sqrt(-1)*IMAG(i)
LAMBDA(i+1) = REIG(i) - sqrt(-1)*IMAG(i)
That is, complex conjugate pairs have consequtive
indices (i,i+1), with the positive imaginary part
listed first.
See the descriptions of K, REIG, Z. |
| [out] | Z | Z (workspace/output) REAL(KIND=WP) M-by-(N-1) array
If JOBZ =='V' then
Z contains real Ritz vectors as follows:
If IMEIG(i)=0, then Z(:,i) is an eigenvector of
the i-th Ritz value.
If IMEIG(i) > 0 (and IMEIG(i+1) < 0) then
[Z(:,i) Z(:,i+1)] span an invariant subspace and
the Ritz values extracted from this subspace are
REIG(i) + sqrt(-1)*IMEIG(i) and
REIG(i) - sqrt(-1)*IMEIG(i).
The corresponding eigenvectors are
Z(:,i) + sqrt(-1)*Z(:,i+1) and
Z(:,i) - sqrt(-1)*Z(:,i+1), respectively.
If JOBZ == 'F', then the above descriptions hold for
the columns of Z*V, where the columns of V are the
eigenvectors of the K-by-K Rayleigh quotient, and Z is
orthonormal. The columns of V are similarly structured:
If IMEIG(i) == 0 then Z*V(:,i) is an eigenvector, and if
IMEIG(i) > 0 then Z*V(:,i)+sqrt(-1)*Z*V(:,i+1) and
Z*V(:,i)-sqrt(-1)*Z*V(:,i+1)
are the eigenvectors of LAMBDA(i), LAMBDA(i+1).
See the descriptions of REIG, IMEIG, X and V. |
| [in] | LDZ | LDZ (input) INTEGER , LDZ >= M
The leading dimension of the array Z. |
| [out] | RES | RES (output) REAL(KIND=WP) (N-1)-by-1 array
RES(1:K) contains the residuals for the K computed
Ritz pairs.
If LAMBDA(i) is real, then
RES(i) = || A * Z(:,i) - LAMBDA(i)*Z(:,i))||_2.
If [LAMBDA(i), LAMBDA(i+1)] is a complex conjugate pair
then
RES(i)=RES(i+1) = || A * Z(:,i:i+1) - Z(:,i:i+1) *B||_F
where B = [ real(LAMBDA(i)) imag(LAMBDA(i)) ]
[-imag(LAMBDA(i)) real(LAMBDA(i)) ].
It holds that
RES(i) = || A*ZC(:,i) - LAMBDA(i) *ZC(:,i) ||_2
RES(i+1) = || A*ZC(:,i+1) - LAMBDA(i+1)*ZC(:,i+1) ||_2
where ZC(:,i) = Z(:,i) + sqrt(-1)*Z(:,i+1)
ZC(:,i+1) = Z(:,i) - sqrt(-1)*Z(:,i+1)
See the description of Z. |
| [out] | B | B (output) REAL(KIND=WP) MIN(M,N)-by-(N-1) array.
IF JOBF =='R', B(1:N,1:K) contains A*U(:,1:K), and can
be used for computing the refined vectors; see further
details in the provided references.
If JOBF == 'E', B(1:N,1;K) contains
A*U(:,1:K)*W(1:K,1:K), which are the vectors from the
Exact DMD, up to scaling by the inverse eigenvalues.
In both cases, the content of B can be lifted to the
original dimension of the input data by pre-multiplying
with the Q factor from the initial QR factorization.
Here A denotes a compression of the underlying operator.
See the descriptions of F and X.
If JOBF =='N', then B is not referenced. |
| [in] | LDB | LDB (input) INTEGER, LDB >= MIN(M,N)
The leading dimension of the array B. |
| [out] | V | V (workspace/output) REAL(KIND=WP) (N-1)-by-(N-1) array
On exit, V(1:K,1:K) contains the K eigenvectors of
the Rayleigh quotient. The eigenvectors of a complex
conjugate pair of eigenvalues are returned in real form
as explained in the description of Z. The Ritz vectors
(returned in Z) are the product of X and V; see
the descriptions of X and Z. |
| [in] | LDV | LDV (input) INTEGER, LDV >= N-1
The leading dimension of the array V. |
| [out] | S | S (output) REAL(KIND=WP) (N-1)-by-(N-1) array
The array S(1:K,1:K) is used for the matrix Rayleigh
quotient. This content is overwritten during
the eigenvalue decomposition by DGEEV.
See the description of K. |
| [in] | LDS | LDS (input) INTEGER, LDS >= N-1
The leading dimension of the array S. |
| [out] | WORK | WORK (workspace/output) REAL(KIND=WP) LWORK-by-1 array
On exit,
WORK(1:MIN(M,N)) contains the scalar factors of the
elementary reflectors as returned by DGEQRF of the
M-by-N input matrix F.
WORK(MIN(M,N)+1:MIN(M,N)+N-1) contains the singular values of
the input submatrix F(1:M,1:N-1).
If the call to DGEDMDQ is only workspace query, then
WORK(1) contains the minimal workspace length and
WORK(2) is the optimal workspace length. Hence, the
length of work is at least 2.
See the description of LWORK. |
| [in] | LWORK | LWORK (input) INTEGER
The minimal length of the workspace vector WORK.
LWORK is calculated as follows:
Let MLWQR = N (minimal workspace for DGEQRF[M,N])
MLWDMD = minimal workspace for DGEDMD (see the
description of LWORK in DGEDMD) for
snapshots of dimensions MIN(M,N)-by-(N-1)
MLWMQR = N (minimal workspace for
DORMQR['L','N',M,N,N])
MLWGQR = N (minimal workspace for DORGQR[M,N,N])
Then
LWORK = MAX(N+MLWQR, N+MLWDMD)
is updated as follows:
if JOBZ == 'V' or JOBZ == 'F' THEN
LWORK = MAX( LWORK, MIN(M,N)+N-1+MLWMQR )
if JOBQ == 'Q' THEN
LWORK = MAX( LWORK, MIN(M,N)+N-1+MLWGQR)
If on entry LWORK = -1, then a workspace query is
assumed and the procedure only computes the minimal
and the optimal workspace lengths for both WORK and
IWORK. See the descriptions of WORK and IWORK. |
| [out] | IWORK | IWORK (workspace/output) INTEGER LIWORK-by-1 array
Workspace that is required only if WHTSVD equals
2 , 3 or 4. (See the description of WHTSVD).
If on entry LWORK =-1 or LIWORK=-1, then the
minimal length of IWORK is computed and returned in
IWORK(1). See the description of LIWORK. |
| [in] | LIWORK | LIWORK (input) INTEGER
The minimal length of the workspace vector IWORK.
If WHTSVD == 1, then only IWORK(1) is used; LIWORK >=1
Let M1=MIN(M,N), N1=N-1. Then
If WHTSVD == 2, then LIWORK >= MAX(1,8*MIN(M1,N1))
If WHTSVD == 3, then LIWORK >= MAX(1,M1+N1-1)
If WHTSVD == 4, then LIWORK >= MAX(3,M1+3*N1)
If on entry LIWORK = -1, then a workspace query is
assumed and the procedure only computes the minimal
and the optimal workspace lengths for both WORK and
IWORK. See the descriptions of WORK and IWORK. |
| [out] | INFO | INFO (output) INTEGER
-i < 0 :: On entry, the i-th argument had an
illegal value
= 0 :: Successful return.
= 1 :: Void input. Quick exit (M=0 or N=0).
= 2 :: The SVD computation of X did not converge.
Suggestion: Check the input data and/or
repeat with different WHTSVD.
= 3 :: The computation of the eigenvalues did not
converge.
= 4 :: If data scaling was requested on input and
the procedure found inconsistency in the data
such that for some column index i,
X(:,i) = 0 but Y(:,i) /= 0, then Y(:,i) is set
to zero if JOBS=='C'. The computation proceeds
with original or modified data and warning
flag is set with INFO=4. |
Definition at line 571 of file dgedmdq.f90.
1.9.7