--- name: "gedmd" section: "3" platform: "linux" package: "lapack-doc" version: "3.12.1-2" release: "Arch Linux 2026-11-01" platforms: ["linux"] url: "https://new.carta.tech/gedmd" --- # gedmd(3) DMD driver, Dynamic Mode Decomposition ## NAME gedmd - DMD driver, Dynamic Mode Decomposition ## SYNOPSIS ### Functions subroutine **cgedmd** (jobs, jobz, jobr, jobf, whtsvd, m, n, x, ldx, y, ldy, nrnk, tol, k, eigs, z, ldz, res, b, ldb, w, ldw, s, lds, zwork, lzwork, rwork, lrwork, iwork, liwork, info)\ **CGEDMD** computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. subroutine **cgedmdq** (jobs, jobz, jobr, jobq, jobt, jobf, whtsvd, m, n, f, ldf, x, ldx, y, ldy, nrnk, tol, k, eigs, z, ldz, res, b, ldb, v, ldv, s, lds, zwork, lzwork, work, lwork, iwork, liwork, info)\ **CGEDMDQ** computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. subroutine **dgedmd** (jobs, jobz, jobr, jobf, whtsvd, m, n, x, ldx, y, ldy, nrnk, tol, k, reig, imeig, z, ldz, res, b, ldb, w, ldw, s, lds, work, lwork, iwork, liwork, info)\ **DGEDMD** computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. subroutine **dgedmdq** (jobs, jobz, jobr, jobq, jobt, jobf, whtsvd, m, n, f, ldf, x, ldx, y, ldy, nrnk, tol, k, reig, imeig, z, ldz, res, b, ldb, v, ldv, s, lds, work, lwork, iwork, liwork, info)\ **DGEDMDQ** computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. subroutine **sgedmd** (jobs, jobz, jobr, jobf, whtsvd, m, n, x, ldx, y, ldy, nrnk, tol, k, reig, imeig, z, ldz, res, b, ldb, w, ldw, s, lds, work, lwork, iwork, liwork, info)\ **SGEDMD** computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. subroutine **sgedmdq** (jobs, jobz, jobr, jobq, jobt, jobf, whtsvd, m, n, f, ldf, x, ldx, y, ldy, nrnk, tol, k, reig, imeig, z, ldz, res, b, ldb, v, ldv, s, lds, work, lwork, iwork, liwork, info)\ **SGEDMDQ** computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. subroutine **zgedmd** (jobs, jobz, jobr, jobf, whtsvd, m, n, x, ldx, y, ldy, nrnk, tol, k, eigs, z, ldz, res, b, ldb, w, ldw, s, lds, zwork, lzwork, rwork, lrwork, iwork, liwork, info)\ **ZGEDMD** computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. subroutine **zgedmdq** (jobs, jobz, jobr, jobq, jobt, jobf, whtsvd, m, n, f, ldf, x, ldx, y, ldy, nrnk, tol, k, eigs, z, ldz, res, b, ldb, v, ldv, s, lds, zwork, lzwork, work, lwork, iwork, liwork, info)\ **ZGEDMDQ** computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. ## Detailed Description ## Function Documentation ### subroutine cgedmd (character, intent(in) jobs, character, intent(in) jobz, character, intent(in) jobr, character, intent(in) jobf, integer, intent(in) whtsvd, integer, intent(in) m, integer, intent(in) n, complex(kind=wp), dimension(ldx,\*), intent(inout) x, integer, intent(in) ldx, complex(kind=wp), dimension(ldy,\*), intent(inout) y, integer, intent(in) ldy, integer, intent(in) nrnk, real(kind=wp), intent(in) tol, integer, intent(out) k, complex(kind=wp), dimension(\*), intent(out) eigs, complex(kind=wp), dimension(ldz,\*), intent(out) z, integer, intent(in) ldz, real(kind=wp), dimension(\*), intent(out) res, complex(kind=wp), dimension(ldb,\*), intent(out) b, integer, intent(in) ldb, complex(kind=wp), dimension(ldw,\*), intent(out) w, integer, intent(in) ldw, complex(kind=wp), dimension(lds,\*), intent(out) s, integer, intent(in) lds, complex(kind=wp), dimension(\*), intent(out) zwork, integer, intent(in) lzwork, real(kind=wp), dimension(\*), intent(out) rwork, integer, intent(in) lrwork, integer, dimension(\*), intent(out) iwork, integer, intent(in) liwork, integer, intent(out) info) **CGEDMD** computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. **Purpose:** ``` !> CGEDMD computes the Dynamic Mode Decomposition (DMD) for !> a pair of data snapshot matrices. For the input matrices !> X and Y such that Y = A*X with an unaccessible matrix !> A, CGEDMD 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, CGEDMD 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]. !> ``` **References:** ``` !> [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 supported by:** ``` !> 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 Contract No: W31P4Q-21-C-0007 !> - DARPA PAI project Contract No: HR0011-18-9-0033 !> - DARPA MoDyL project !> 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 !> ``` **Distribution Statement A:** ``` !> Approved for Public Release, Distribution Unlimited. !> Cleared by DARPA on September 29, 2022 !> ``` **Parameters** *JOBS* ``` !> JOBS (input) CHARACTER*1 !> Determines whether the initial data snapshots are scaled !> by a diagonal matrix. !> '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. !> ``` *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 X(:,1:K)*W, where X !> contains a POD basis (leading left singular vectors !> of the data matrix X) and W contains the eigenvectors !> of the corresponding Rayleigh quotient. !> See the descriptions of K, X, W, Z. !> 'N' :: The eigenvectors are not computed. !> ``` *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. !> ``` *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. !> ``` *WHTSVD* ``` !> WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 } !> Allows for a selection of the SVD algorithm from the !> LAPACK library. !> 1 :: CGESVD (the QR SVD algorithm) !> 2 :: CGESDD (the Divide and Conquer algorithm; if enough !> workspace available, this is the fastest option) !> 3 :: CGESVDQ (the preconditioned QR SVD ; this and 4 !> are the most accurate options) !> 4 :: CGEJSV (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. !> ``` *M* ``` !> M (input) INTEGER, M>= 0 !> The state space dimension (the row dimension of X, Y). !> ``` *N* ``` !> N (input) INTEGER, 0 <= N <= M !> The number of data snapshot pairs !> (the number of columns of X and Y). !> ``` *X* ``` !> X (input/output) COMPLEX(KIND=WP) M-by-N array !> > On entry, X contains the data snapshot matrix X. It is !> assumed that the column norms of X are in the range of !> the normalized floating point numbers. !> < On exit, the leading K columns of X contain a POD basis, !> i.e. the leading K left singular vectors of the input !> data matrix X, U(:,1:K). All N columns of X contain all !> left singular vectors of the input matrix X. !> See the descriptions of K, Z and W. !> ``` *LDX* ``` !> LDX (input) INTEGER, LDX >= M !> The leading dimension of the array X. !> ``` *Y* ``` !> Y (input/workspace/output) COMPLEX(KIND=WP) M-by-N array !> > On entry, Y contains the data snapshot matrix Y !> < On exit, !> If JOBR == 'R', the leading K columns of Y contain !> the residual vectors for the computed Ritz pairs. !> See the description of RES. !> If JOBR == 'N', Y contains the original input data, !> scaled according to the value of JOBS. !> ``` *LDY* ``` !> LDY (input) INTEGER , LDY >= M !> The leading dimension of the array Y. !> ``` *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 :: 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 descriptions of TOL and K. !> ``` *TOL* ``` !> TOL (input) REAL(KIND=WP), 0 <= TOL < 1 !> The tolerance for truncating small singular values. !> See the description of NRNK. !> ``` *K* ``` !> K (output) INTEGER, 0 <= K <= N !> The dimension of the POD basis for the data snapshot !> matrix X 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. !> ``` *EIGS* ``` !> EIGS (output) COMPLEX(KIND=WP) N-by-1 array !> The leading K (K<=N) entries of EIGS contain !> the computed eigenvalues (Ritz values). !> See the descriptions of K, and Z. !> ``` *Z* ``` !> Z (workspace/output) COMPLEX(KIND=WP) M-by-N array !> If JOBZ =='V' then Z contains the Ritz vectors. Z(:,i) !> is an eigenvector of the i-th Ritz value; ||Z(:,i)||_2=1. !> If JOBZ == 'F', then the Z(:,i)'s are given implicitly as !> the columns of X(:,1:K)*W(1:K,1:K), i.e. X(:,1:K)*W(:,i) !> is an eigenvector corresponding to EIGS(i). The columns !> of W(1:k,1:K) are the computed eigenvectors of the !> K-by-K Rayleigh quotient. !> See the descriptions of EIGS, X and W. !> ``` *LDZ* ``` !> LDZ (input) INTEGER , LDZ >= M !> The leading dimension of the array Z. !> ``` *RES* ``` !> RES (output) REAL(KIND=WP) N-by-1 array !> RES(1:K) contains the residuals for the K computed !> Ritz pairs, !> RES(i) = || A * Z(:,i) - EIGS(i)*Z(:,i))||_2. !> See the description of EIGS and Z. !> ``` *B* ``` !> B (output) COMPLEX(KIND=WP) M-by-N array. !> IF JOBF =='R', B(1:M,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:M,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. !> If JOBF =='N', then B is not referenced. !> See the descriptions of X, W, K. !> ``` *LDB* ``` !> LDB (input) INTEGER, LDB >= M !> The leading dimension of the array B. !> ``` *W* ``` !> W (workspace/output) COMPLEX(KIND=WP) N-by-N array !> On exit, W(1:K,1:K) contains the K computed !> eigenvectors of the matrix Rayleigh quotient. !> The Ritz vectors (returned in Z) are the !> product of X (containing a POD basis for the input !> matrix X) and W. See the descriptions of K, S, X and Z. !> W is also used as a workspace to temporarily store the !> right singular vectors of X. !> ``` *LDW* ``` !> LDW (input) INTEGER, LDW >= N !> The leading dimension of the array W. !> ``` *S* ``` !> S (workspace/output) COMPLEX(KIND=WP) N-by-N array !> The array S(1:K,1:K) is used for the matrix Rayleigh !> quotient. This content is overwritten during !> the eigenvalue decomposition by CGEEV. !> See the description of K. !> ``` *LDS* ``` !> LDS (input) INTEGER, LDS >= N !> The leading dimension of the array S. !> ``` *ZWORK* ``` !> ZWORK (workspace/output) COMPLEX(KIND=WP) LZWORK-by-1 array !> ZWORK is used as complex workspace in the complex SVD, as !> specified by WHTSVD (1,2, 3 or 4) and for CGEEV for computing !> the eigenvalues of a Rayleigh quotient. !> If the call to CGEDMD is only workspace query, then !> ZWORK(1) contains the minimal complex workspace length and !> ZWORK(2) is the optimal complex workspace length. !> Hence, the length of work is at least 2. !> See the description of LZWORK. !> ``` *LZWORK* ``` !> LZWORK (input) INTEGER !> The minimal length of the workspace vector ZWORK. !> LZWORK is calculated as MAX(LZWORK_SVD, LZWORK_CGEEV), !> where LZWORK_CGEEV = MAX( 1, 2*N ) and the minimal !> LZWORK_SVD is calculated as follows !> If WHTSVD == 1 :: CGESVD :: !> LZWORK_SVD = MAX(1,2*MIN(M,N)+MAX(M,N)) !> If WHTSVD == 2 :: CGESDD :: !> LZWORK_SVD = 2*MIN(M,N)*MIN(M,N)+2*MIN(M,N)+MAX(M,N) !> If WHTSVD == 3 :: CGESVDQ :: !> LZWORK_SVD = obtainable by a query !> If WHTSVD == 4 :: CGEJSV :: !> LZWORK_SVD = obtainable by a query !> If on entry LZWORK = -1, then a workspace query is !> assumed and the procedure only computes the minimal !> and the optimal workspace lengths and returns them in !> LZWORK(1) and LZWORK(2), respectively. !> ``` *RWORK* ``` !> RWORK (workspace/output) REAL(KIND=WP) LRWORK-by-1 array !> On exit, RWORK(1:N) contains the singular values of !> X (for JOBS=='N') or column scaled X (JOBS=='S', 'C'). !> If WHTSVD==4, then RWORK(N+1) and RWORK(N+2) contain !> scaling factor RWORK(N+2)/RWORK(N+1) used to scale X !> and Y to avoid overflow in the SVD of X. !> This may be of interest if the scaling option is off !> and as many as possible smallest eigenvalues are !> desired to the highest feasible accuracy. !> If the call to CGEDMD is only workspace query, then !> RWORK(1) contains the minimal workspace length. !> See the description of LRWORK. !> ``` *LRWORK* ``` !> LRWORK (input) INTEGER !> The minimal length of the workspace vector RWORK. !> LRWORK is calculated as follows: !> LRWORK = MAX(1, N+LRWORK_SVD,N+LRWORK_CGEEV), where !> LRWORK_CGEEV = MAX(1,2*N) and RWORK_SVD is the real workspace !> for the SVD subroutine determined by the input parameter !> WHTSVD. !> If WHTSVD == 1 :: CGESVD :: !> LRWORK_SVD = 5*MIN(M,N) !> If WHTSVD == 2 :: CGESDD :: !> LRWORK_SVD = MAX(5*MIN(M,N)*MIN(M,N)+7*MIN(M,N), !> 2*MAX(M,N)*MIN(M,N)+2*MIN(M,N)*MIN(M,N)+MIN(M,N) ) ) !> If WHTSVD == 3 :: CGESVDQ :: !> LRWORK_SVD = obtainable by a query !> If WHTSVD == 4 :: CGEJSV :: !> LRWORK_SVD = obtainable by a query !> If on entry LRWORK = -1, then a workspace query is !> assumed and the procedure only computes the minimal !> real workspace length and returns it in RWORK(1). !> ``` *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. !> ``` *LIWORK* ``` !> LIWORK (input) INTEGER !> The minimal length of the workspace vector IWORK. !> If WHTSVD == 1, then only IWORK(1) is used; LIWORK >=1 !> If WHTSVD == 2, then LIWORK >= MAX(1,8*MIN(M,N)) !> If WHTSVD == 3, then LIWORK >= MAX(1,M+N-1) !> If WHTSVD == 4, then LIWORK >= MAX(3,M+3*N) !> If on entry LIWORK = -1, then a workspace query is !> assumed and the procedure only computes the minimal !> and the optimal workspace lengths for ZWORK, RWORK and !> IWORK. See the descriptions of ZWORK, RWORK and IWORK. !> ``` *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. !> ``` **Author** Zlatko Drmac Definition at line **496** of file **cgedmd.f90**. ### subroutine cgedmdq (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, complex(kind=wp), dimension(ldf,\*), intent(inout) f, integer, intent(in) ldf, complex(kind=wp), dimension(ldx,\*), intent(out) x, integer, intent(in) ldx, complex(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, complex(kind=wp), dimension(\*), intent(out) eigs, complex(kind=wp), dimension(ldz,\*), intent(out) z, integer, intent(in) ldz, real(kind=wp), dimension(\*), intent(out) res, complex(kind=wp), dimension(ldb,\*), intent(out) b, integer, intent(in) ldb, complex(kind=wp), dimension(ldv,\*), intent(out) v, integer, intent(in) ldv, complex(kind=wp), dimension(lds,\*), intent(out) s, integer, intent(in) lds, complex(kind=wp), dimension(\*), intent(out) zwork, integer, intent(in) lzwork, real(kind=wp), dimension(\*), intent(out) work, integer, intent(in) lwork, integer, dimension(\*), intent(out) iwork, integer, intent(in) liwork, integer, intent(out) info) **CGEDMDQ** computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. **Purpose:** ``` !> CGEDMDQ 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, CGEDMDQ 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, CGEDMDQ 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]. !> ``` **References:** ``` !> [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 supported by:** ``` !> 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 Contract No: W31P4Q-21-C-0007 !> - DARPA PAI project Contract No: HR0011-18-9-0033 !> - DARPA MoDyL project !> 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. !> ``` **Developed and supported by:** ``` !> Approved for Public Release, Distribution Unlimited. !> Cleared by DARPA on September 29, 2022 !> ``` **Parameters** *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. !> ``` *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 discretised operator onto the span of !> the data snapshots. See the descriptions of F, V, Z. !> Q is from the inital QR facorization. !> 'N' :: The eigenvectors are not computed. !> ``` *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. !> ``` *JOBQ* ``` !> JOBQ (input) CHARACTER*1 !> Specifies whether to explicitly compute and return the !> unitary 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. !> ``` *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. !> ``` *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'. !> ``` *WHTSVD* ``` !> WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 } !> Allows for a selection of the SVD algorithm from the !> LAPACK library. !> 1 :: CGESVD (the QR SVD algorithm) !> 2 :: CGESDD (the Divide and Conquer algorithm; if enough !> workspace available, this is the fastest option) !> 3 :: CGESVDQ (the preconditioned QR SVD ; this and 4 !> are the most accurate options) !> 4 :: CGEJSV (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. !> ``` *M* ``` !> M (input) INTEGER, M >= 0 !> The state space dimension (the number of rows of F). !> ``` *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. !> ``` *F* ``` !> F (input/output) COMPLEX(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 CGEQRF. The !> remaining information to restore the orthogonal matrix !> of the initial QR factorization is stored in ZWORK(1:MIN(M,N)). !> See the description of ZWORK. !> ``` *LDF* ``` !> LDF (input) INTEGER, LDF >= M !> The leading dimension of the array F. !> ``` *X* ``` !> X (workspace/output) COMPLEX(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. !> ``` *LDX* ``` !> LDX (input) INTEGER, LDX >= N !> The leading dimension of the array X. !> ``` *Y* ``` !> Y (workspace/output) COMPLEX(KIND=WP) MIN(M,N)-by-(N) 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. !> ``` *LDY* ``` !> LDY (input) INTEGER , LDY >= N !> The leading dimension of the array Y. !> ``` *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. !> ``` *TOL* ``` !> TOL (input) REAL(KIND=WP), 0 <= TOL < 1 !> The tolerance for truncating small singular values. !> See the description of NRNK. !> ``` *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. !> ``` *EIGS* ``` !> EIGS (output) COMPLEX(KIND=WP) (N-1)-by-1 array !> The leading K (K<=N-1) entries of EIGS contain !> the computed eigenvalues (Ritz values). !> See the descriptions of K, and Z. !> ``` *Z* ``` !> Z (workspace/output) COMPLEX(KIND=WP) M-by-(N-1) array !> If JOBZ =='V' then Z contains the Ritz vectors. Z(:,i) !> is an eigenvector of the i-th Ritz value; ||Z(:,i)||_2=1. !> If JOBZ == 'F', then the Z(:,i)'s are given implicitly as !> Z*V, where Z contains orthonormal matrix (the product of !> Q from the initial QR factorization and the SVD/POD_basis !> returned by CGEDMD in X) and the second factor (the !> eigenvectors of the Rayleigh quotient) is in the array V, !> as returned by CGEDMD. That is, X(:,1:K)*V(:,i) !> is an eigenvector corresponding to EIGS(i). The columns !> of V(1:K,1:K) are the computed eigenvectors of the !> K-by-K Rayleigh quotient. !> See the descriptions of EIGS, X and V. !> ``` *LDZ* ``` !> LDZ (input) INTEGER , LDZ >= M !> The leading dimension of the array Z. !> ``` *RES* ``` !> RES (output) REAL(KIND=WP) (N-1)-by-1 array !> RES(1:K) contains the residuals for the K computed !> Ritz pairs, !> RES(i) = || A * Z(:,i) - EIGS(i)*Z(:,i))||_2. !> See the description of EIGS and Z. !> ``` *B* ``` !> B (output) COMPLEX(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. !> ``` *LDB* ``` !> LDB (input) INTEGER, LDB >= MIN(M,N) !> The leading dimension of the array B. !> ``` *V* ``` !> V (workspace/output) COMPLEX(KIND=WP) (N-1)-by-(N-1) array !> On exit, V(1:K,1:K) V contains the K eigenvectors of !> the Rayleigh quotient. The Ritz vectors !> (returned in Z) are the product of Q from the initial QR !> factorization (see the description of F) X (see the !> description of X) and V. !> ``` *LDV* ``` !> LDV (input) INTEGER, LDV >= N-1 !> The leading dimension of the array V. !> ``` *S* ``` !> S (output) COMPLEX(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 CGEEV. !> See the description of K. !> ``` *LDS* ``` !> LDS (input) INTEGER, LDS >= N-1 !> The leading dimension of the array S. !> ``` *ZWORK* ``` !> ZWORK (workspace/output) COMPLEX(KIND=WP) LWORK-by-1 array !> On exit, !> ZWORK(1:MIN(M,N)) contains the scalar factors of the !> elementary reflectors as returned by CGEQRF of the !> M-by-N input matrix F. !> If the call to CGEDMDQ is only workspace query, then !> ZWORK(1) contains the minimal complex workspace length and !> ZWORK(2) is the optimal complex workspace length. !> Hence, the length of work is at least 2. !> See the description of LZWORK. !> ``` *LZWORK* ``` !> LZWORK (input) INTEGER !> The minimal length of the workspace vector ZWORK. !> LZWORK is calculated as follows: !> Let MLWQR = N (minimal workspace for CGEQRF[M,N]) !> MLWDMD = minimal workspace for CGEDMD (see the !> description of LWORK in CGEDMD) !> MLWMQR = N (minimal workspace for !> ZUNMQR['L','N',M,N,N]) !> MLWGQR = N (minimal workspace for ZUNGQR[M,N,N]) !> MINMN = MIN(M,N) !> Then !> LZWORK = MAX(2, MIN(M,N)+MLWQR, MINMN+MLWDMD) !> is further updated as follows: !> if JOBZ == 'V' or JOBZ == 'F' THEN !> LZWORK = MAX( LZWORK, MINMN+MLWMQR ) !> if JOBQ == 'Q' THEN !> LZWORK = MAX( ZLWORK, MINMN+MLWGQR) !> ``` *WORK* ``` !> WORK (workspace/output) REAL(KIND=WP) LWORK-by-1 array !> On exit, !> WORK(1:N-1) contains the singular values of !> the input submatrix F(1:M,1:N-1). !> If the call to CGEDMDQ 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. !> ``` *LWORK* ``` !> LWORK (input) INTEGER !> The minimal length of the workspace vector WORK. !> LWORK is the same as in CGEDMD, because in CGEDMDQ !> only CGEDMD requires real workspace for snapshots !> of dimensions MIN(M,N)-by-(N-1). !> 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. !> ``` *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. !> ``` *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(M,N)) !> If WHTSVD == 3, then LIWORK >= MAX(1,M+N-1) !> If WHTSVD == 4, then LIWORK >= MAX(3,M+3*N) !> 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. !> ``` *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. !> ``` **Author** Zlatko Drmac Definition at line **552** of file **cgedmdq.f90**. ### subroutine dgedmd (character, intent(in) jobs, character, intent(in) jobz, character, intent(in) jobr, character, intent(in) jobf, integer, intent(in) whtsvd, integer, intent(in) m, integer, intent(in) n, real(kind=wp), dimension(ldx,\*), intent(inout) x, integer, intent(in) ldx, real(kind=wp), dimension(ldy,\*), intent(inout) 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(ldw,\*), intent(out) w, integer, intent(in) ldw, 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) **DGEDMD** computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. **Purpose:** ``` !> DGEDMD computes the Dynamic Mode Decomposition (DMD) for !> a pair of data snapshot matrices. For the input matrices !> X and Y such that Y = A*X with an unaccessible matrix !> A, DGEDMD 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, DGEDMD 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]. !> ``` **References:** ``` !> [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 supported by:** ``` !> 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 Contract No: W31P4Q-21-C-0007 !> - DARPA PAI project Contract No: HR0011-18-9-0033 !> - DARPA MoDyL project !> 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 !> ``` **Distribution Statement A:** ``` !> Approved for Public Release, Distribution Unlimited. !> Cleared by DARPA on September 29, 2022 !> ``` **Parameters** *JOBS* ``` !> JOBS (input) is CHARACTER*1 !> Determines whether the initial data snapshots are scaled !> by a diagonal matrix. !> '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. !> ``` *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 X(:,1:K)*W, where X !> contains a POD basis (leading left singular vectors !> of the data matrix X) and W contains the eigenvectors !> of the corresponding Rayleigh quotient. !> See the descriptions of K, X, W, Z. !> 'N' :: The eigenvectors are not computed. !> ``` *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. !> ``` *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. !> ``` *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. !> ``` *M* ``` !> M (input) INTEGER, M>= 0 !> The state space dimension (the row dimension of X, Y). !> ``` *N* ``` !> N (input) INTEGER, 0 <= N <= M !> The number of data snapshot pairs !> (the number of columns of X and Y). !> ``` *X* ``` !> X (input/output) REAL(KIND=WP) M-by-N array !> > On entry, X contains the data snapshot matrix X. It is !> assumed that the column norms of X are in the range of !> the normalized floating point numbers. !> < On exit, the leading K columns of X contain a POD basis, !> i.e. the leading K left singular vectors of the input !> data matrix X, U(:,1:K). All N columns of X contain all !> left singular vectors of the input matrix X. !> See the descriptions of K, Z and W. !> ``` *LDX* ``` !> LDX (input) INTEGER, LDX >= M !> The leading dimension of the array X. !> ``` *Y* ``` !> Y (input/workspace/output) REAL(KIND=WP) M-by-N array !> > On entry, Y contains the data snapshot matrix Y !> < On exit, !> If JOBR == 'R', the leading K columns of Y contain !> the residual vectors for the computed Ritz pairs. !> See the description of RES. !> If JOBR == 'N', Y contains the original input data, !> scaled according to the value of JOBS. !> ``` *LDY* ``` !> LDY (input) INTEGER , LDY >= M !> The leading dimension of the array Y. !> ``` *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 :: 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 descriptions of TOL and K. !> ``` *TOL* ``` !> TOL (input) REAL(KIND=WP), 0 <= TOL < 1 !> The tolerance for truncating small singular values. !> See the description of NRNK. !> ``` *K* ``` !> K (output) INTEGER, 0 <= K <= N !> The dimension of the POD basis for the data snapshot !> matrix X 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. !> ``` *REIG* ``` !> REIG (output) REAL(KIND=WP) N-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, and Z. !> ``` *IMEIG* ``` !> IMEIG (output) REAL(KIND=WP) N-by-1 array !> The leading K (K<=N) entries of IMEIG 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 consecutive !> indices (i,i+1), with the positive imaginary part !> listed first. !> See the descriptions of K, REIG, and Z. !> ``` *Z* ``` !> Z (workspace/output) REAL(KIND=WP) M-by-N 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; ||Z(:,i)||_2=1. !> 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. !> || Z(:,i:i+1)||_F = 1. !> If JOBZ == 'F', then the above descriptions hold for !> the columns of X(:,1:K)*W(1:K,1:K), where the columns !> of W(1:k,1:K) are the computed eigenvectors of the !> K-by-K Rayleigh quotient. The columns of W(1:K,1:K) !> are similarly structured: If IMEIG(i) == 0 then !> X(:,1:K)*W(:,i) is an eigenvector, and if IMEIG(i)>0 !> then X(:,1:K)*W(:,i)+sqrt(-1)*X(:,1:K)*W(:,i+1) and !> X(:,1:K)*W(:,i)-sqrt(-1)*X(:,1:K)*W(:,i+1) !> are the eigenvectors of LAMBDA(i), LAMBDA(i+1). !> See the descriptions of REIG, IMEIG, X and W. !> ``` *LDZ* ``` !> LDZ (input) INTEGER , LDZ >= M !> The leading dimension of the array Z. !> ``` *RES* ``` !> RES (output) REAL(KIND=WP) N-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 REIG, IMEIG and Z. !> ``` *B* ``` !> B (output) REAL(KIND=WP) M-by-N array. !> IF JOBF =='R', B(1:M,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:M,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. !> If JOBF =='N', then B is not referenced. !> See the descriptions of X, W, K. !> ``` *LDB* ``` !> LDB (input) INTEGER, LDB >= M !> The leading dimension of the array B. !> ``` *W* ``` !> W (workspace/output) REAL(KIND=WP) N-by-N array !> On exit, W(1:K,1:K) contains the K computed !> eigenvectors of the matrix Rayleigh quotient (real and !> imaginary parts for each complex conjugate pair of the !> eigenvalues). The Ritz vectors (returned in Z) are the !> product of X (containing a POD basis for the input !> matrix X) and W. See the descriptions of K, S, X and Z. !> W is also used as a workspace to temporarily store the !> right singular vectors of X. !> ``` *LDW* ``` !> LDW (input) INTEGER, LDW >= N !> The leading dimension of the array W. !> ``` *S* ``` !> S (workspace/output) REAL(KIND=WP) N-by-N 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. !> ``` *LDS* ``` !> LDS (input) INTEGER, LDS >= N !> The leading dimension of the array S. !> ``` *WORK* ``` !> WORK (workspace/output) REAL(KIND=WP) LWORK-by-1 array !> On exit, WORK(1:N) contains the singular values of !> X (for JOBS=='N') or column scaled X (JOBS=='S', 'C'). !> If WHTSVD==4, then WORK(N+1) and WORK(N+2) contain !> scaling factor WORK(N+2)/WORK(N+1) used to scale X !> and Y to avoid overflow in the SVD of X. !> This may be of interest if the scaling option is off !> and as many as possible smallest eigenvalues are !> desired to the highest feasible accuracy. !> If the call to DGEDMD is only workspace query, then !> WORK(1) contains the minimal workspace length and !> WORK(2) is the optimal workspace length. Hence, the !> leng of work is at least 2. !> See the description of LWORK. !> ``` *LWORK* ``` !> LWORK (input) INTEGER !> The minimal length of the workspace vector WORK. !> LWORK is calculated as follows: !> If WHTSVD == 1 :: !> If JOBZ == 'V', then !> LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,4*N)). !> If JOBZ == 'N' then !> LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,3*N)). !> Here LWORK_SVD = MAX(1,3*N+M,5*N) is the minimal !> workspace length of DGESVD. !> If WHTSVD == 2 :: !> If JOBZ == 'V', then !> LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,4*N)) !> If JOBZ == 'N', then !> LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,3*N)) !> Here LWORK_SVD = MAX(M, 5*N*N+4*N)+3*N*N is the !> minimal workspace length of DGESDD. !> If WHTSVD == 3 :: !> If JOBZ == 'V', then !> LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,4*N)) !> If JOBZ == 'N', then !> LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,3*N)) !> Here LWORK_SVD = N+M+MAX(3*N+1, !> MAX(1,3*N+M,5*N),MAX(1,N)) !> is the minimal workspace length of DGESVDQ. !> If WHTSVD == 4 :: !> If JOBZ == 'V', then !> LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,4*N)) !> If JOBZ == 'N', then !> LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,3*N)) !> Here LWORK_SVD = MAX(7,2*M+N,6*N+2*N*N) is the !> minimal workspace length of DGEJSV. !> The above expressions are not simplified in order to !> make the usage of WORK more transparent, and for !> easier checking. In any case, LWORK >= 2. !> 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. !> ``` *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. !> ``` *LIWORK* ``` !> LIWORK (input) INTEGER !> The minimal length of the workspace vector IWORK. !> If WHTSVD == 1, then only IWORK(1) is used; LIWORK >=1 !> If WHTSVD == 2, then LIWORK >= MAX(1,8*MIN(M,N)) !> If WHTSVD == 3, then LIWORK >= MAX(1,M+N-1) !> If WHTSVD == 4, then LIWORK >= MAX(3,M+3*N) !> 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. !> ``` *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. !> ``` **Author** Zlatko Drmac Definition at line **530** of file **dgedmd.f90**. ### 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. **Purpose:** ``` !> 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]. !> ``` **References:** ``` !> [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 supported by:** ``` !> 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 Contract No: W31P4Q-21-C-0007 !> - DARPA PAI project Contract No: HR0011-18-9-0033 !> - DARPA MoDyL project !> 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. !> ``` **Distribution Statement A:** ``` !> Approved for Public Release, Distribution Unlimited. !> Cleared by DARPA on September 29, 2022 !> ``` **Parameters** *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. !> ``` *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. !> ``` *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. !> ``` *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. !> ``` *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. !> ``` *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'. !> ``` *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. !> ``` *M* ``` !> M (input) INTEGER, M >= 0 !> The state space dimension (the number of rows of F). !> ``` *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. !> ``` *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. !> ``` *LDF* ``` !> LDF (input) INTEGER, LDF >= M !> The leading dimension of the array F. !> ``` *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. !> ``` *LDX* ``` !> LDX (input) INTEGER, LDX >= N !> The leading dimension of the array X. !> ``` *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. !> ``` *LDY* ``` !> LDY (input) INTEGER , LDY >= N !> The leading dimension of the array Y. !> ``` *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. !> ``` *TOL* ``` !> TOL (input) REAL(KIND=WP), 0 <= TOL < 1 !> The tolerance for truncating small singular values. !> See the description of NRNK. !> ``` *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. !> ``` *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. !> ``` *IMEIG* ``` !> IMEIG (output) REAL(KIND=WP) (N-1)-by-1 array !> The leading K (K 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. !> ``` *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. !> ``` *LDZ* ``` !> LDZ (input) INTEGER , LDZ >= M !> The leading dimension of the array Z. !> ``` *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. !> ``` *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. !> ``` *LDB* ``` !> LDB (input) INTEGER, LDB >= MIN(M,N) !> The leading dimension of the array B. !> ``` *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. !> ``` *LDV* ``` !> LDV (input) INTEGER, LDV >= N-1 !> The leading dimension of the array V. !> ``` *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. !> ``` *LDS* ``` !> LDS (input) INTEGER, LDS >= N-1 !> The leading dimension of the array S. !> ``` *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. !> ``` *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. !> ``` *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. !> ``` *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. !> ``` *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. !> ``` **Author** Zlatko Drmac Definition at line **571** of file **dgedmdq.f90**. ### subroutine sgedmd (character, intent(in) jobs, character, intent(in) jobz, character, intent(in) jobr, character, intent(in) jobf, integer, intent(in) whtsvd, integer, intent(in) m, integer, intent(in) n, real(kind=wp), dimension(ldx,\*), intent(inout) x, integer, intent(in) ldx, real(kind=wp), dimension(ldy,\*), intent(inout) 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(ldw,\*), intent(out) w, integer, intent(in) ldw, 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) **SGEDMD** computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. **Purpose:** ``` !> SGEDMD computes the Dynamic Mode Decomposition (DMD) for !> a pair of data snapshot matrices. For the input matrices !> X and Y such that Y = A*X with an unaccessible matrix !> A, SGEDMD 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, SGEDMD 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]. !> ``` **References:** ``` !> [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 supported by:** ``` !> 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 Contract No: W31P4Q-21-C-0007 !> - DARPA PAI project Contract No: HR0011-18-9-0033 !> - DARPA MoDyL project !> 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 !> ``` **Distribution Statement A:** ``` !> Distribution Statement A: !> Approved for Public Release, Distribution Unlimited. !> Cleared by DARPA on September 29, 2022 !> ``` **Parameters** *JOBS* ``` !> JOBS (input) CHARACTER*1 !> Determines whether the initial data snapshots are scaled !> by a diagonal matrix. !> '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. !> ``` *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 X(:,1:K)*W, where X !> contains a POD basis (leading left singular vectors !> of the data matrix X) and W contains the eigenvectors !> of the corresponding Rayleigh quotient. !> See the descriptions of K, X, W, Z. !> 'N' :: The eigenvectors are not computed. !> ``` *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. !> ``` *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. !> ``` *WHTSVD* ``` !> WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 } !> Allows for a selection of the SVD algorithm from the !> LAPACK library. !> 1 :: SGESVD (the QR SVD algorithm) !> 2 :: SGESDD (the Divide and Conquer algorithm; if enough !> workspace available, this is the fastest option) !> 3 :: SGESVDQ (the preconditioned QR SVD ; this and 4 !> are the most accurate options) !> 4 :: SGEJSV (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. !> ``` *M* ``` !> M (input) INTEGER, M>= 0 !> The state space dimension (the row dimension of X, Y). !> ``` *N* ``` !> N (input) INTEGER, 0 <= N <= M !> The number of data snapshot pairs !> (the number of columns of X and Y). !> ``` *X* ``` !> X (input/output) REAL(KIND=WP) M-by-N array !> > On entry, X contains the data snapshot matrix X. It is !> assumed that the column norms of X are in the range of !> the normalized floating point numbers. !> < On exit, the leading K columns of X contain a POD basis, !> i.e. the leading K left singular vectors of the input !> data matrix X, U(:,1:K). All N columns of X contain all !> left singular vectors of the input matrix X. !> See the descriptions of K, Z and W. !> ``` *LDX* ``` !> LDX (input) INTEGER, LDX >= M !> The leading dimension of the array X. !> ``` *Y* ``` !> Y (input/workspace/output) REAL(KIND=WP) M-by-N array !> > On entry, Y contains the data snapshot matrix Y !> < On exit, !> If JOBR == 'R', the leading K columns of Y contain !> the residual vectors for the computed Ritz pairs. !> See the description of RES. !> If JOBR == 'N', Y contains the original input data, !> scaled according to the value of JOBS. !> ``` *LDY* ``` !> LDY (input) INTEGER , LDY >= M !> The leading dimension of the array Y. !> ``` *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 :: 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 descriptions of TOL and K. !> ``` *TOL* ``` !> TOL (input) REAL(KIND=WP), 0 <= TOL < 1 !> The tolerance for truncating small singular values. !> See the description of NRNK. !> ``` *K* ``` !> K (output) INTEGER, 0 <= K <= N !> The dimension of the POD basis for the data snapshot !> matrix X 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. !> ``` *REIG* ``` !> REIG (output) REAL(KIND=WP) N-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, and Z. !> ``` *IMEIG* ``` !> IMEIG (output) REAL(KIND=WP) N-by-1 array !> The leading K (K<=N) entries of IMEIG 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 consecutive !> indices (i,i+1), with the positive imaginary part !> listed first. !> See the descriptions of K, REIG, and Z. !> ``` *Z* ``` !> Z (workspace/output) REAL(KIND=WP) M-by-N 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; ||Z(:,i)||_2=1. !> 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. !> || Z(:,i:i+1)||_F = 1. !> If JOBZ == 'F', then the above descriptions hold for !> the columns of X(:,1:K)*W(1:K,1:K), where the columns !> of W(1:k,1:K) are the computed eigenvectors of the !> K-by-K Rayleigh quotient. The columns of W(1:K,1:K) !> are similarly structured: If IMEIG(i) == 0 then !> X(:,1:K)*W(:,i) is an eigenvector, and if IMEIG(i)>0 !> then X(:,1:K)*W(:,i)+sqrt(-1)*X(:,1:K)*W(:,i+1) and !> X(:,1:K)*W(:,i)-sqrt(-1)*X(:,1:K)*W(:,i+1) !> are the eigenvectors of LAMBDA(i), LAMBDA(i+1). !> See the descriptions of REIG, IMEIG, X and W. !> ``` *LDZ* ``` !> LDZ (input) INTEGER , LDZ >= M !> The leading dimension of the array Z. !> ``` *RES* ``` !> RES (output) REAL(KIND=WP) N-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 REIG, IMEIG and Z. !> ``` *B* ``` !> B (output) REAL(KIND=WP) M-by-N array. !> IF JOBF =='R', B(1:M,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:M,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. !> If JOBF =='N', then B is not referenced. !> See the descriptions of X, W, K. !> ``` *LDB* ``` !> LDB (input) INTEGER, LDB >= M !> The leading dimension of the array B. !> ``` *W* ``` !> W (workspace/output) REAL(KIND=WP) N-by-N array !> On exit, W(1:K,1:K) contains the K computed !> eigenvectors of the matrix Rayleigh quotient (real and !> imaginary parts for each complex conjugate pair of the !> eigenvalues). The Ritz vectors (returned in Z) are the !> product of X (containing a POD basis for the input !> matrix X) and W. See the descriptions of K, S, X and Z. !> W is also used as a workspace to temporarily store the !> left singular vectors of X. !> ``` *LDW* ``` !> LDW (input) INTEGER, LDW >= N !> The leading dimension of the array W. !> ``` *S* ``` !> S (workspace/output) REAL(KIND=WP) N-by-N array !> The array S(1:K,1:K) is used for the matrix Rayleigh !> quotient. This content is overwritten during !> the eigenvalue decomposition by SGEEV. !> See the description of K. !> ``` *LDS* ``` !> LDS (input) INTEGER, LDS >= N !> The leading dimension of the array S. !> ``` *WORK* ``` !> WORK (workspace/output) REAL(KIND=WP) LWORK-by-1 array !> On exit, WORK(1:N) contains the singular values of !> X (for JOBS=='N') or column scaled X (JOBS=='S', 'C'). !> If WHTSVD==4, then WORK(N+1) and WORK(N+2) contain !> scaling factor WORK(N+2)/WORK(N+1) used to scale X !> and Y to avoid overflow in the SVD of X. !> This may be of interest if the scaling option is off !> and as many as possible smallest eigenvalues are !> desired to the highest feasible accuracy. !> If the call to SGEDMD 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. !> ``` *LWORK* ``` !> LWORK (input) INTEGER !> The minimal length of the workspace vector WORK. !> LWORK is calculated as follows: !> If WHTSVD == 1 :: !> If JOBZ == 'V', then !> LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,4*N)). !> If JOBZ == 'N' then !> LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,3*N)). !> Here LWORK_SVD = MAX(1,3*N+M,5*N) is the minimal !> workspace length of SGESVD. !> If WHTSVD == 2 :: !> If JOBZ == 'V', then !> LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,4*N)) !> If JOBZ == 'N', then !> LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,3*N)) !> Here LWORK_SVD = MAX(M, 5*N*N+4*N)+3*N*N is the !> minimal workspace length of SGESDD. !> If WHTSVD == 3 :: !> If JOBZ == 'V', then !> LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,4*N)) !> If JOBZ == 'N', then !> LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,3*N)) !> Here LWORK_SVD = N+M+MAX(3*N+1, !> MAX(1,3*N+M,5*N),MAX(1,N)) !> is the minimal workspace length of SGESVDQ. !> If WHTSVD == 4 :: !> If JOBZ == 'V', then !> LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,4*N)) !> If JOBZ == 'N', then !> LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,3*N)) !> Here LWORK_SVD = MAX(7,2*M+N,6*N+2*N*N) is the !> minimal workspace length of SGEJSV. !> The above expressions are not simplified in order to !> make the usage of WORK more transparent, and for !> easier checking. In any case, LWORK >= 2. !> 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. !> ``` *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. !> ``` *LIWORK* ``` !> LIWORK (input) INTEGER !> The minimal length of the workspace vector IWORK. !> If WHTSVD == 1, then only IWORK(1) is used; LIWORK >=1 !> If WHTSVD == 2, then LIWORK >= MAX(1,8*MIN(M,N)) !> If WHTSVD == 3, then LIWORK >= MAX(1,M+N-1) !> If WHTSVD == 4, then LIWORK >= MAX(3,M+3*N) !> 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. !> ``` *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. !> ``` **Author** Zlatko Drmac Definition at line **530** of file **sgedmd.f90**. ### subroutine sgedmdq (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) **SGEDMDQ** computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. **Purpose:** ``` !> SGEDMDQ 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, SGEDMDQ 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, SGEDMDQ 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]. !> ``` **References:** ``` !> [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 supported by:** ``` !> 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 Contract No: W31P4Q-21-C-0007 !> - DARPA PAI project Contract No: HR0011-18-9-0033 !> - DARPA MoDyL project !> 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. !> ``` **Distribution Statement A:** ``` !> Approved for Public Release, Distribution Unlimited. !> Cleared by DARPA on September 29, 2022 !> ``` **Parameters** *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. !> ``` *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. !> ``` *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. !> ``` *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. !> ``` *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. !> ``` *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'. !> ``` *WHTSVD* ``` !> WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 } !> Allows for a selection of the SVD algorithm from the !> LAPACK library. !> 1 :: SGESVD (the QR SVD algorithm) !> 2 :: SGESDD (the Divide and Conquer algorithm; if enough !> workspace available, this is the fastest option) !> 3 :: SGESVDQ (the preconditioned QR SVD ; this and 4 !> are the most accurate options) !> 4 :: SGEJSV (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. !> ``` *M* ``` !> M (input) INTEGER, M >= 0 !> The state space dimension (the number of rows of F) !> ``` *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. !> ``` *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 SGEQRF. The !> remaining information to restore the orthogonal matrix !> of the initial QR factorization is stored in WORK(1:N). !> See the description of WORK. !> ``` *LDF* ``` !> LDF (input) INTEGER, LDF >= M !> The leading dimension of the array F. !> ``` *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. !> ``` *LDX* ``` !> LDX (input) INTEGER, LDX >= N !> The leading dimension of the array X !> ``` *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. !> ``` *LDY* ``` !> LDY (input) INTEGER , LDY >= N !> The leading dimension of the array Y !> ``` *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. !> ``` *TOL* ``` !> TOL (input) REAL(KIND=WP), 0 <= TOL < 1 !> The tolerance for truncating small singular values. !> See the description of NRNK. !> ``` *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. !> ``` *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. !> ``` *IMEIG* ``` !> IMEIG (output) REAL(KIND=WP) (N-1)-by-1 array !> The leading K (K 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 consecutive !> indices (i,i+1), with the positive imaginary part !> listed first. !> See the descriptions of K, REIG, Z. !> ``` *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. !> ``` *LDZ* ``` !> LDZ (input) INTEGER , LDZ >= M !> The leading dimension of the array Z. !> ``` *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. !> ``` *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. !> ``` *LDB* ``` !> LDB (input) INTEGER, LDB >= MIN(M,N) !> The leading dimension of the array B. !> ``` *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. !> ``` *LDV* ``` !> LDV (input) INTEGER, LDV >= N-1 !> The leading dimension of the array V. !> ``` *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 SGEEV. !> See the description of K. !> ``` *LDS* ``` !> LDS (input) INTEGER, LDS >= N-1 !> The leading dimension of the array S. !> ``` *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 SGEQRF 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 SGEDMDQ 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. !> ``` *LWORK* ``` !> LWORK (input) INTEGER !> The minimal length of the workspace vector WORK. !> LWORK is calculated as follows: !> Let MLWQR = N (minimal workspace for SGEQRF[M,N]) !> MLWDMD = minimal workspace for SGEDMD (see the !> description of LWORK in SGEDMD) for !> snapshots of dimensions MIN(M,N)-by-(N-1) !> MLWMQR = N (minimal workspace for !> SORMQR['L','N',M,N,N]) !> MLWGQR = N (minimal workspace for SORGQR[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. !> ``` *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. !> ``` *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 worskpace 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. !> ``` *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. !> ``` **Author** Zlatko Drmac Definition at line **571** of file **sgedmdq.f90**. ### subroutine zgedmd (character, intent(in) jobs, character, intent(in) jobz, character, intent(in) jobr, character, intent(in) jobf, integer, intent(in) whtsvd, integer, intent(in) m, integer, intent(in) n, complex(kind=wp), dimension(ldx,\*), intent(inout) x, integer, intent(in) ldx, complex(kind=wp), dimension(ldy,\*), intent(inout) y, integer, intent(in) ldy, integer, intent(in) nrnk, real(kind=wp), intent(in) tol, integer, intent(out) k, complex(kind=wp), dimension(\*), intent(out) eigs, complex(kind=wp), dimension(ldz,\*), intent(out) z, integer, intent(in) ldz, real(kind=wp), dimension(\*), intent(out) res, complex(kind=wp), dimension(ldb,\*), intent(out) b, integer, intent(in) ldb, complex(kind=wp), dimension(ldw,\*), intent(out) w, integer, intent(in) ldw, complex(kind=wp), dimension(lds,\*), intent(out) s, integer, intent(in) lds, complex(kind=wp), dimension(\*), intent(out) zwork, integer, intent(in) lzwork, real(kind=wp), dimension(\*), intent(out) rwork, integer, intent(in) lrwork, integer, dimension(\*), intent(out) iwork, integer, intent(in) liwork, integer, intent(out) info) **ZGEDMD** computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. **Purpose:** ``` !> ZGEDMD computes the Dynamic Mode Decomposition (DMD) for !> a pair of data snapshot matrices. For the input matrices !> X and Y such that Y = A*X with an unaccessible matrix !> A, ZGEDMD 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, ZGEDMD 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]. !> ``` **References:** ``` !> [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 supported by:** ``` !> 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 Contract No: W31P4Q-21-C-0007 !> - DARPA PAI project Contract No: HR0011-18-9-0033 !> - DARPA MoDyL project !> 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 !> ``` **Distribution Statement A:** ``` !> Approved for Public Release, Distribution Unlimited. !> Cleared by DARPA on September 29, 2022 !> ``` **Parameters** *JOBS* ``` !> JOBS (input) CHARACTER*1 !> Determines whether the initial data snapshots are scaled !> by a diagonal matrix. !> '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. !> ``` *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 X(:,1:K)*W, where X !> contains a POD basis (leading left singular vectors !> of the data matrix X) and W contains the eigenvectors !> of the corresponding Rayleigh quotient. !> See the descriptions of K, X, W, Z. !> 'N' :: The eigenvectors are not computed. !> ``` *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. !> ``` *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. !> ``` *WHTSVD* ``` !> WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 } !> Allows for a selection of the SVD algorithm from the !> LAPACK library. !> 1 :: ZGESVD (the QR SVD algorithm) !> 2 :: ZGESDD (the Divide and Conquer algorithm; if enough !> workspace available, this is the fastest option) !> 3 :: ZGESVDQ (the preconditioned QR SVD ; this and 4 !> are the most accurate options) !> 4 :: ZGEJSV (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. !> ``` *M* ``` !> M (input) INTEGER, M>= 0 !> The state space dimension (the row dimension of X, Y). !> ``` *N* ``` !> N (input) INTEGER, 0 <= N <= M !> The number of data snapshot pairs !> (the number of columns of X and Y). !> ``` *X* ``` !> X (input/output) COMPLEX(KIND=WP) M-by-N array !> > On entry, X contains the data snapshot matrix X. It is !> assumed that the column norms of X are in the range of !> the normalized floating point numbers. !> < On exit, the leading K columns of X contain a POD basis, !> i.e. the leading K left singular vectors of the input !> data matrix X, U(:,1:K). All N columns of X contain all !> left singular vectors of the input matrix X. !> See the descriptions of K, Z and W. !> ``` *LDX* ``` !> LDX (input) INTEGER, LDX >= M !> The leading dimension of the array X. !> ``` *Y* ``` !> Y (input/workspace/output) COMPLEX(KIND=WP) M-by-N array !> > On entry, Y contains the data snapshot matrix Y !> < On exit, !> If JOBR == 'R', the leading K columns of Y contain !> the residual vectors for the computed Ritz pairs. !> See the description of RES. !> If JOBR == 'N', Y contains the original input data, !> scaled according to the value of JOBS. !> ``` *LDY* ``` !> LDY (input) INTEGER , LDY >= M !> The leading dimension of the array Y. !> ``` *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 :: 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 descriptions of TOL and K. !> ``` *TOL* ``` !> TOL (input) REAL(KIND=WP), 0 <= TOL < 1 !> The tolerance for truncating small singular values. !> See the description of NRNK. !> ``` *K* ``` !> K (output) INTEGER, 0 <= K <= N !> The dimension of the POD basis for the data snapshot !> matrix X 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. !> ``` *EIGS* ``` !> EIGS (output) COMPLEX(KIND=WP) N-by-1 array !> The leading K (K<=N) entries of EIGS contain !> the computed eigenvalues (Ritz values). !> See the descriptions of K, and Z. !> ``` *Z* ``` !> Z (workspace/output) COMPLEX(KIND=WP) M-by-N array !> If JOBZ =='V' then Z contains the Ritz vectors. Z(:,i) !> is an eigenvector of the i-th Ritz value; ||Z(:,i)||_2=1. !> If JOBZ == 'F', then the Z(:,i)'s are given implicitly as !> the columns of X(:,1:K)*W(1:K,1:K), i.e. X(:,1:K)*W(:,i) !> is an eigenvector corresponding to EIGS(i). The columns !> of W(1:k,1:K) are the computed eigenvectors of the !> K-by-K Rayleigh quotient. !> See the descriptions of EIGS, X and W. !> ``` *LDZ* ``` !> LDZ (input) INTEGER , LDZ >= M !> The leading dimension of the array Z. !> ``` *RES* ``` !> RES (output) REAL(KIND=WP) N-by-1 array !> RES(1:K) contains the residuals for the K computed !> Ritz pairs, !> RES(i) = || A * Z(:,i) - EIGS(i)*Z(:,i))||_2. !> See the description of EIGS and Z. !> ``` *B* ``` !> B (output) COMPLEX(KIND=WP) M-by-N array. !> IF JOBF =='R', B(1:M,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:M,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. !> If JOBF =='N', then B is not referenced. !> See the descriptions of X, W, K. !> ``` *LDB* ``` !> LDB (input) INTEGER, LDB >= M !> The leading dimension of the array B. !> ``` *W* ``` !> W (workspace/output) COMPLEX(KIND=WP) N-by-N array !> On exit, W(1:K,1:K) contains the K computed !> eigenvectors of the matrix Rayleigh quotient. !> The Ritz vectors (returned in Z) are the !> product of X (containing a POD basis for the input !> matrix X) and W. See the descriptions of K, S, X and Z. !> W is also used as a workspace to temporarily store the !> right singular vectors of X. !> ``` *LDW* ``` !> LDW (input) INTEGER, LDW >= N !> The leading dimension of the array W. !> ``` *S* ``` !> S (workspace/output) COMPLEX(KIND=WP) N-by-N array !> The array S(1:K,1:K) is used for the matrix Rayleigh !> quotient. This content is overwritten during !> the eigenvalue decomposition by ZGEEV. !> See the description of K. !> ``` *LDS* ``` !> LDS (input) INTEGER, LDS >= N !> The leading dimension of the array S. !> ``` *ZWORK* ``` !> ZWORK (workspace/output) COMPLEX(KIND=WP) LZWORK-by-1 array !> ZWORK is used as complex workspace in the complex SVD, as !> specified by WHTSVD (1,2, 3 or 4) and for ZGEEV for computing !> the eigenvalues of a Rayleigh quotient. !> If the call to ZGEDMD is only workspace query, then !> ZWORK(1) contains the minimal complex workspace length and !> ZWORK(2) is the optimal complex workspace length. !> Hence, the length of work is at least 2. !> See the description of LZWORK. !> ``` *LZWORK* ``` !> LZWORK (input) INTEGER !> The minimal length of the workspace vector ZWORK. !> LZWORK is calculated as MAX(LZWORK_SVD, LZWORK_ZGEEV), !> where LZWORK_ZGEEV = MAX( 1, 2*N ) and the minimal !> LZWORK_SVD is calculated as follows !> If WHTSVD == 1 :: ZGESVD :: !> LZWORK_SVD = MAX(1,2*MIN(M,N)+MAX(M,N)) !> If WHTSVD == 2 :: ZGESDD :: !> LZWORK_SVD = 2*MIN(M,N)*MIN(M,N)+2*MIN(M,N)+MAX(M,N) !> If WHTSVD == 3 :: ZGESVDQ :: !> LZWORK_SVD = obtainable by a query !> If WHTSVD == 4 :: ZGEJSV :: !> LZWORK_SVD = obtainable by a query !> If on entry LZWORK = -1, then a workspace query is !> assumed and the procedure only computes the minimal !> and the optimal workspace lengths and returns them in !> LZWORK(1) and LZWORK(2), respectively. !> ``` *RWORK* ``` !> RWORK (workspace/output) REAL(KIND=WP) LRWORK-by-1 array !> On exit, RWORK(1:N) contains the singular values of !> X (for JOBS=='N') or column scaled X (JOBS=='S', 'C'). !> If WHTSVD==4, then RWORK(N+1) and RWORK(N+2) contain !> scaling factor RWORK(N+2)/RWORK(N+1) used to scale X !> and Y to avoid overflow in the SVD of X. !> This may be of interest if the scaling option is off !> and as many as possible smallest eigenvalues are !> desired to the highest feasible accuracy. !> If the call to ZGEDMD is only workspace query, then !> RWORK(1) contains the minimal workspace length. !> See the description of LRWORK. !> ``` *LRWORK* ``` !> LRWORK (input) INTEGER !> The minimal length of the workspace vector RWORK. !> LRWORK is calculated as follows: !> LRWORK = MAX(1, N+LRWORK_SVD,N+LRWORK_ZGEEV), where !> LRWORK_ZGEEV = MAX(1,2*N) and RWORK_SVD is the real workspace !> for the SVD subroutine determined by the input parameter !> WHTSVD. !> If WHTSVD == 1 :: ZGESVD :: !> LRWORK_SVD = 5*MIN(M,N) !> If WHTSVD == 2 :: ZGESDD :: !> LRWORK_SVD = MAX(5*MIN(M,N)*MIN(M,N)+7*MIN(M,N), !> 2*MAX(M,N)*MIN(M,N)+2*MIN(M,N)*MIN(M,N)+MIN(M,N) ) ) !> If WHTSVD == 3 :: ZGESVDQ :: !> LRWORK_SVD = obtainable by a query !> If WHTSVD == 4 :: ZGEJSV :: !> LRWORK_SVD = obtainable by a query !> If on entry LRWORK = -1, then a workspace query is !> assumed and the procedure only computes the minimal !> real workspace length and returns it in RWORK(1). !> ``` *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. !> ``` *LIWORK* ``` !> LIWORK (input) INTEGER !> The minimal length of the workspace vector IWORK. !> If WHTSVD == 1, then only IWORK(1) is used; LIWORK >=1 !> If WHTSVD == 2, then LIWORK >= MAX(1,8*MIN(M,N)) !> If WHTSVD == 3, then LIWORK >= MAX(1,M+N-1) !> If WHTSVD == 4, then LIWORK >= MAX(3,M+3*N) !> If on entry LIWORK = -1, then a workspace query is !> assumed and the procedure only computes the minimal !> and the optimal workspace lengths for ZWORK, RWORK and !> IWORK. See the descriptions of ZWORK, RWORK and IWORK. !> ``` *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. !> ``` **Author** Zlatko Drmac Definition at line **496** of file **zgedmd.f90**. ### subroutine zgedmdq (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, complex(kind=wp), dimension(ldf,\*), intent(inout) f, integer, intent(in) ldf, complex(kind=wp), dimension(ldx,\*), intent(out) x, integer, intent(in) ldx, complex(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, complex(kind=wp), dimension(\*), intent(out) eigs, complex(kind=wp), dimension(ldz,\*), intent(out) z, integer, intent(in) ldz, real(kind=wp), dimension(\*), intent(out) res, complex(kind=wp), dimension(ldb,\*), intent(out) b, integer, intent(in) ldb, complex(kind=wp), dimension(ldv,\*), intent(out) v, integer, intent(in) ldv, complex(kind=wp), dimension(lds,\*), intent(out) s, integer, intent(in) lds, complex(kind=wp), dimension(\*), intent(out) zwork, integer, intent(in) lzwork, real(kind=wp), dimension(\*), intent(out) work, integer, intent(in) lwork, integer, dimension(\*), intent(out) iwork, integer, intent(in) liwork, integer, intent(out) info) **ZGEDMDQ** computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. **Purpose:** ``` !> ZGEDMDQ 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, ZGEDMDQ 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, ZGEDMDQ 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]. !> ``` **References:** ``` !> [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 supported by:** ``` !> 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 Contract No: W31P4Q-21-C-0007 !> - DARPA PAI project Contract No: HR0011-18-9-0033 !> - DARPA MoDyL project !> 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. !> ``` **Developed and supported by:** ``` !> Distribution Statement A: !> Approved for Public Release, Distribution Unlimited. !> Cleared by DARPA on September 29, 2022 !> ``` **Parameters** *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. !> ``` *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. !> ``` *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. !> ``` *JOBQ* ``` !> JOBQ (input) CHARACTER*1 !> Specifies whether to explicitly compute and return the !> unitary 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. !> ``` *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. !> ``` *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'. !> ``` *WHTSVD* ``` !> WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 } !> Allows for a selection of the SVD algorithm from the !> LAPACK library. !> 1 :: ZGESVD (the QR SVD algorithm) !> 2 :: ZGESDD (the Divide and Conquer algorithm; if enough !> workspace available, this is the fastest option) !> 3 :: ZGESVDQ (the preconditioned QR SVD ; this and 4 !> are the most accurate options) !> 4 :: ZGEJSV (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. !> ``` *M* ``` !> M (input) INTEGER, M >= 0 !> The state space dimension (the number of rows of F). !> ``` *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. !> ``` *F* ``` !> F (input/output) COMPLEX(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 ZGEQRF. The !> remaining information to restore the orthogonal matrix !> of the initial QR factorization is stored in ZWORK(1:MIN(M,N)). !> See the description of ZWORK. !> ``` *LDF* ``` !> LDF (input) INTEGER, LDF >= M !> The leading dimension of the array F. !> ``` *X* ``` !> X (workspace/output) COMPLEX(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. !> ``` *LDX* ``` !> LDX (input) INTEGER, LDX >= N !> The leading dimension of the array X. !> ``` *Y* ``` !> Y (workspace/output) COMPLEX(KIND=WP) MIN(M,N)-by-(N) 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. !> ``` *LDY* ``` !> LDY (input) INTEGER , LDY >= N !> The leading dimension of the array Y. !> ``` *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. !> ``` *TOL* ``` !> TOL (input) REAL(KIND=WP), 0 <= TOL < 1 !> The tolerance for truncating small singular values. !> See the description of NRNK. !> ``` *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. !> ``` *EIGS* ``` !> EIGS (output) COMPLEX(KIND=WP) (N-1)-by-1 array !> The leading K (K<=N-1) entries of EIGS contain !> the computed eigenvalues (Ritz values). !> See the descriptions of K, and Z. !> ``` *Z* ``` !> Z (workspace/output) COMPLEX(KIND=WP) M-by-(N-1) array !> If JOBZ =='V' then Z contains the Ritz vectors. Z(:,i) !> is an eigenvector of the i-th Ritz value; ||Z(:,i)||_2=1. !> If JOBZ == 'F', then the Z(:,i)'s are given implicitly as !> Z*V, where Z contains orthonormal matrix (the product of !> Q from the initial QR factorization and the SVD/POD_basis !> returned by ZGEDMD in X) and the second factor (the !> eigenvectors of the Rayleigh quotient) is in the array V, !> as returned by ZGEDMD. That is, X(:,1:K)*V(:,i) !> is an eigenvector corresponding to EIGS(i). The columns !> of V(1:K,1:K) are the computed eigenvectors of the !> K-by-K Rayleigh quotient. !> See the descriptions of EIGS, X and V. !> ``` *LDZ* ``` !> LDZ (input) INTEGER , LDZ >= M !> The leading dimension of the array Z. !> ``` *RES* ``` !> RES (output) REAL(KIND=WP) (N-1)-by-1 array !> RES(1:K) contains the residuals for the K computed !> Ritz pairs, !> RES(i) = || A * Z(:,i) - EIGS(i)*Z(:,i))||_2. !> See the description of EIGS and Z. !> ``` *B* ``` !> B (output) COMPLEX(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. !> ``` *LDB* ``` !> LDB (input) INTEGER, LDB >= MIN(M,N) !> The leading dimension of the array B. !> ``` *V* ``` !> V (workspace/output) COMPLEX(KIND=WP) (N-1)-by-(N-1) array !> On exit, V(1:K,1:K) V contains the K eigenvectors of !> the Rayleigh quotient. The Ritz vectors !> (returned in Z) are the product of Q from the initial QR !> factorization (see the description of F) X (see the !> description of X) and V. !> ``` *LDV* ``` !> LDV (input) INTEGER, LDV >= N-1 !> The leading dimension of the array V. !> ``` *S* ``` !> S (output) COMPLEX(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 ZGEEV. !> See the description of K. !> ``` *LDS* ``` !> LDS (input) INTEGER, LDS >= N-1 !> The leading dimension of the array S. !> ``` *ZWORK* ``` !> ZWORK (workspace/output) COMPLEX(KIND=WP) LWORK-by-1 array !> On exit, !> ZWORK(1:MIN(M,N)) contains the scalar factors of the !> elementary reflectors as returned by ZGEQRF of the !> M-by-N input matrix F. !> If the call to ZGEDMDQ is only workspace query, then !> ZWORK(1) contains the minimal complex workspace length and !> ZWORK(2) is the optimal complex workspace length. !> Hence, the length of work is at least 2. !> See the description of LZWORK. !> ``` *LZWORK* ``` !> LZWORK (input) INTEGER !> The minimal length of the workspace vector ZWORK. !> LZWORK is calculated as follows: !> Let MLWQR = N (minimal workspace for ZGEQRF[M,N]) !> MLWDMD = minimal workspace for ZGEDMD (see the !> description of LWORK in ZGEDMD) !> MLWMQR = N (minimal workspace for !> ZUNMQR['L','N',M,N,N]) !> MLWGQR = N (minimal workspace for ZUNGQR[M,N,N]) !> MINMN = MIN(M,N) !> Then !> LZWORK = MAX(2, MIN(M,N)+MLWQR, MINMN+MLWDMD) !> is further updated as follows: !> if JOBZ == 'V' or JOBZ == 'F' THEN !> LZWORK = MAX(LZWORK, MINMN+MLWMQR) !> if JOBQ == 'Q' THEN !> LZWORK = MAX(ZLWORK, MINMN+MLWGQR) !> ``` *WORK* ``` !> WORK (workspace/output) REAL(KIND=WP) LWORK-by-1 array !> On exit, !> WORK(1:N-1) contains the singular values of !> the input submatrix F(1:M,1:N-1). !> If the call to ZGEDMDQ 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. !> ``` *LWORK* ``` !> LWORK (input) INTEGER !> The minimal length of the workspace vector WORK. !> LWORK is the same as in ZGEDMD, because in ZGEDMDQ !> only ZGEDMD requires real workspace for snapshots !> of dimensions MIN(M,N)-by-(N-1). !> If on entry LWORK = -1, then a workspace query is !> assumed and the procedure only computes the minimal !> and the optimal workspace length for WORK. !> ``` *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. !> ``` *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. !> ``` *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. !> ``` **Author** Zlatko Drmac Definition at line **551** of file **zgedmdq.f90**. ## Author Generated automatically by Doxygen for LAPACK from the source code.