gedmd(3)

📦
version 3.12.1 · 9 versions
, in arch rolling

✨ NAME ✨

gedmd — DMD driver, Dynamic Mode Decomposition 💅

(also known as gedmd(3), straight out of your Library Functions Manual — hi bestie, welcome 🫶)

🌸 SYNOPSIS 🌸

The functions 🎀

Ok so there's a whole squad of them, one per precision/type combo. Same vibe, different letters. Here's the full lineup:

  • cgedmd — CGEDMD computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.
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)
  • cgedmdq — CGEDMDQ 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)
  • dgedmd — DGEDMD 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)
  • dgedmdq — DGEDMDQ 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)
  • sgedmd — SGEDMD 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)
  • sgedmdq — SGEDMDQ 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)
  • zgedmd — ZGEDMD 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)
  • zgedmdq — ZGEDMDQ 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)

💫 Detailed Description 💫

Not much prose lives here, the real content is all under the function docs below, so let's go 👇

🎀 Function Documentation 🎀

💜 subroutine cgedmd 💜

The whole signature, character for character, because I would never mess with your argument list:

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

(…and this is exactly where part 1 runs out of page, bestie — X keeps going in the next chunk 🫶)

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 — the little status receipt at the end, bestie. it tells you what went down 💌
    • -i < 0 — the i-th argument had an illegal value 🙈 (rude of it, honestly)
    • = 0 — successful return. slay 🎉
    • = 1 — void input, quick exit because M=0 or N=0. nothing to do, we go home 🏠
    • = 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 the warning flag is set with INFO=4. ⚠️

Author

Zlatko Drmac

Definition at line 496 of file cgedmd.f90.

🌀 subroutine cgedmdq — the QR-flavored one 🎀

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. yes bestie, it said what it said 💅

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: (yes, it shows up twice — not a glitch, we're just built different 🫶)

Approved for Public Release, Distribution Unlimited. Cleared by DARPA on September 29, 2022

💖 Parameters

  • 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. keepin' it raw 💀
  • 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 — 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 — 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 — 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 — 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 — 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 — input INTEGER, M >= 0. The state space dimension (the number of rows of F). 📏
  • 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 — okay bestie, this one's a big deal 💖 it's an 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. No cap, this is the main character of the whole routine. 👑
  • LDF — input INTEGER, LDF >= M. The leading dimension of the array F. Just give it the space it needs, bestie. ✨
  • 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. It's giving "glow-up era" for your snapshots. 🦋
  • LDX — input INTEGER, LDX >= N. The leading dimension of the array X. Yes, it's N, not M — stay woke. 💅
  • 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. Trailing snapshots need love too, okay? 🌸
  • LDY — input INTEGER, LDY >= N. The leading dimension of the array Y. Same energy as LDX, but for Y. 🎀
  • 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. Basically, you're the boss of how many singular values get to stay. 💪
  • TOL — input REAL(KIND=WP), 0 <= TOL < 1. The tolerance for truncating small singular values. See the description of NRNK. Small but mighty — don't sleep on this one. 🔥
  • 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. This is the actual count that comes out — no cap. 📊
  • 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. The Ritz values are literally the eigenvalues of your compressed operator — iconic. 👑
  • 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. Two modes, both slay — pick your fighter. 💅
  • LDZ — input INTEGER, LDZ >= M. The leading dimension of the array Z. Don't forget to size it right! 🎀
  • 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. Residuals are like the vibe check — how close did we actually get? 📏

  • 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. Note the semicolon in B(1:N,1;K) — that's in the original, bestie, don't come for me. 😭
  • LDB — input INTEGER, LDB >= MIN(M,N). The leading dimension of the array B. Size it like you mean it. ✨
  • 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. It's giving collaboration — Q, X, and V working together. 🫶
  • LDV — input INTEGER, LDV >= N-1. The leading dimension of the array V. Yep, N-1 — not N. 💅
  • 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. Heads up: S gets clobbered by CGEEV — don't get attached. 💀
  • LDS — input INTEGER, LDS >= N-1. The leading dimension of the array S. Match the N-1 vibe. 🎀
  • 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. Workspace query? ZWORK(1) and ZWORK(2) are your new besties. 🔍
  • 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)

Note: the original says ZLWORK in that last line — I'm not fixing it, I'm just the messenger. 😭

  • 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. Yes, "hence" is lowercase in the original — keep it cute. 🎀
  • 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. Set it to -1 to ask "how much space do I need?" — it's the polite thing to do. 💅
  • 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. Integer workspace, but make it fashion. ✨
  • LIWORK — (the description continues in the next part, but LIWORK is the length of IWORK — stay tuned, bestie) 🎀
  • LIWORK (input) INTEGER — bestie, this is 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 (output) INTEGER — the receipts, no cap:
    • -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

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. [1] P. Schmid: Dynamic mode decomposition of numerical and experimental data, Journal of Fluid Mechanics 656, 5-28, 2010.
  2. [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. [3] Z. Drmac: A LAPACK implementation of the Dynamic Mode Decomposition I. Technical report. AIMDyn Inc. and LAPACK Working Note 298.
  4. [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 (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 (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 (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 (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 (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 (input) INTEGER, M>= 0 — the state space dimension (the row dimension of X, Y).
  • N (input) INTEGER, 0 <= N <= M — the number of data snapshot pairs (the number of columns of X and Y).
  • 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 (input) INTEGER, LDX >= M — the leading dimension of the array X.
  • 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 (input) INTEGER , LDY >= M — the leading dimension of the array Y.
  • NRNK
  • NRNK — (input) INTEGER — this one sets the mode for computing the numerical rank, aka how we truncate the tiny singular values of the input matrix X. On input, if:
    • NRNK = -1 :: the i-th singular value sigma(i) is truncated if sigma(i) <= TOL*sigma(1). This option is recommended. 👑
    • NRNK = -2 :: the 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. iykyk, it's giving "lab experiment" energy 🧪

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 — (input) REAL(KIND=WP), 0 <= TOL < 1 — the tolerance for truncating small singular values. See the description of NRNK. 🎀
  • 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 — (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 — (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 — (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 — (input) INTEGER, LDZ >= M — the leading dimension of the array Z. 💅
  • 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 — (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 — (input) INTEGER, LDB >= M — the leading dimension of the array B.
  • 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 — (input) INTEGER, LDW >= N — the leading dimension of the array W.
  • 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 — (input) INTEGER, LDS >= N — the leading dimension of the array S.
  • 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. (source spelling, not mine 💀)

See the description of 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 — (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 — (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 — (output) INTEGER

okay so picking right back up with the INFO deets, because bestie this one's a whole mood board of possible outcomes 💅

  • -i < 0 — on entry, the i-th argument had an illegal value. classic. somebody passed something cursed.
  • = 0 — successful return. we love to see it 🎀
  • = 1 — void input. quick exit (M=0 or N=0). literally nothing to do, so it dipped immediately.
  • = 2 — the SVD computation of X did not converge 😭 Suggestion: check the input data and/or repeat with a different WHTSVD.
  • = 3 — the computation of the eigenvalues did not converge. sad.
  • = 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 the warning flag is set with INFO=4. so it just shrugs and keeps going 🫶

Author

Zlatko Drmac

Definition at line 530 of file dgedmd.f90.

okay now we get to the queen of the page, the one with the extra Q energy — DGEDMDQ 👑

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. that's the whole vibe, but like, let her explain.

✨ Purpose ✨

DGEDMDQ computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices, using a QR factorization based compression of the data.

so — 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. it does it all, no cap 💖

For further details see the references listed below. For more details of the implementation see [3].

📚 References 📚

  1. [1] P. Schmid: Dynamic mode decomposition of numerical and experimental data, Journal of Fluid Mechanics 656, 5-28, 2010.
  2. [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. [3] Z. Drmac: A LAPACK implementation of the Dynamic Mode Decomposition I. Technical report. AIMDyn Inc. and LAPACK Working Note 298.
  4. [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.

the reading list ate, honestly 🔥

🌸 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. (legal said what legal said 🫶)

🦋 Distribution Statement A 🦋

Approved for Public Release, Distribution Unlimited. Cleared by DARPA on September 29, 2022

iconic. cleared. periodt. 👑

✨ 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. (yes, that's how the slicing works, don't @ me)

  • '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. (the cleanup crew option 🧹)
  • '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. the main character decision, fr.

  • '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. (respectfully, we're skipping that part)

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'. ⚠️ so like, don't try it with anything else.
  • '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'. 👀 keep that in mind bestie.

WHTSVD

WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 } — allows for a selection of the SVD algorithm from the LAPACK library. pick your fighter 💅

  • 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. so like, if you scaled, it's basically a tie 🎀

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 — okay so this is the big one, and she has so much to say…

  • 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’s assumed that the column norms of F are in the range of the normalized floating point numbers. ✨
    • On exit:
      • If JOBQ == 'Q', 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 — input INTEGER, LDF >= M. The leading dimension of the array F. 🎀
  • 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 — input INTEGER, LDX >= N. The leading dimension of the array X. 🌸
  • 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 — input INTEGER, LDY >= N. The leading dimension of the array Y. 🫶
  • 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 a 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 — input REAL(KIND=WP), 0 <= TOL < 1. The tolerance for truncating small singular values. See the description of NRNK. 🌈
  • 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 — 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 — output REAL(KIND=WP) (N-1)-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, 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 — input INTEGER, LDZ >= M. The leading dimension of the array Z. 🍓
  • 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 — 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 — input INTEGER, LDB >= MIN(M,N). The leading dimension of the array B. 🌸
  • 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 — input INTEGER, LDV >= N-1. The leading dimension of the array V. 🫶
  • 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 — input INTEGER, LDS >= N-1. The leading dimension of the array S. 🌈
  • 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 —

okay so LWORK is (input) INTEGER — literally the minimal length of the workspace vector WORK, and the way it gets calculated is a whole little recipe, so stay with me bestie 💅

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)

and if LWORK = -1 on entry, it's a workspace query era: the procedure only computes the minimal and the optimal workspace lengths for both WORK and IWORK. See the descriptions of WORK and IWORK ✨ iykyk

IWORK 🫶

IWORK is (workspace/output) INTEGER LIWORK-by-1 array — not me with a whole array 😭. It's 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 is (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 the tea is:

  • 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 is (output) INTEGER — the vibes report, and no cap you need to read it:

  • -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 🌈

new function, new main character energy — and here's the whole signature, verbatim, because we do not mess with function signatures in this house 💅

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. periodt 💅

Purpose 💖

okay so 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 📚

the receipts, bestie:

  • [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 (input) CHARACTER1* — Determines whether the initial data snapshots are scaled by a diagonal matrix. vibe check:

  • '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 (input) CHARACTER1* — Determines whether the eigenvectors (Koopman modes) will be computed. Koopman mode era:

  • '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 (input) CHARACTER1* — 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 (input) CHARACTER1* — 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 (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 (input) INTEGER, M>= 0 — The state space dimension (the row dimension of X, Y).

N

N

N (input) INTEGER, 0 <= N <= M — the number of data snapshot pairs, i.e. the number of columns in X and Y. That's it, that's the whole job bestie 💅

X

X (input/output) REAL(KIND=WP) M-by-N array.

On entry, X holds the data snapshot matrix X. The vibe check here: it's assumed 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, aka the leading K left singular vectors of the input data matrix X, U(:,1:K). And all N columns of X contain all left singular vectors of the input matrix X. She ate, fr fr 👑

See the descriptions of K, Z and W.

LDX

LDX (input) INTEGER, LDX >= M — the leading dimension of the array X. Respect the constraint or it's giving undefined behaviour energy 😭

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. Same deal as LDX, no shortcuts 🎀

NRNK

NRNK (input) INTEGER — this one decides the mode for computing the numerical rank, meaning how we truncate the small singular values of the input matrix X. On input:

  • NRNK = -1 :: the i-th singular value sigma(i) is truncated if sigma(i) <= TOL*sigma(1). This option is recommended. 👑 (it's the main character of rank detection)
  • NRNK = -2 :: the 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. 🙈 (iconic but highkey risky)

And you can enforce the numerical rank by using a positive value of NRNK, like this:

  • 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. No cap on the lore, 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. They're literally holding hands, iykyk 🫶

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. Yes bestie, another one, keep up 💀

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. The leading-dimension saga continues, and honestly? Obsessed 🎀

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. She's doing double duty, main character energy fr 💅

LDW

LDW (input) INTEGER, LDW >= N — the leading dimension of the array W. (Note the N this time, not M — yes, they switch it up on you 🙈)

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. So don't get attached to whatever you put in there, it's gone 😭

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. (Niche? Yes. Iconic? Also yes 🌈)

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

LWORK (input) INTEGER — the minimal length of the workspace vector WORK, bestie. Here's exactly how it's calculated:

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.

These expressions are not simplified — and that's on purpose, no cap. It makes the usage of WORK more transparent and easier to check. 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

   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. Periodt. 💅

🎯 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 💫

  • JOBT (input) CHARACTER*1 — tells the routine whether you want the upper triangular factor back, bestie. Say less.
    • 'R' — the matrix R of the QR factorization of the data snapshot matrix F is returned in the array Y. Peep the description of Y and Further details for the receipts.
    • 'N' — the matrix R is not returned. We're gatekeeping it. 🫣

JOBF 🎀

  • JOBF (input) CHARACTER*1 — controls whether we store the info you need for post-processing (like, say, computing refined Ritz vectors).
    • 'R' — the matrix needed for the refinement of the Ritz vectors gets 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. We simply ✨ do not ✨.
Heads up: to be useful on exit, this option needs JOBQ='Q'. No JOBQ='Q', no vibe. 💅

WHTSVD 🌈

  • WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 } — lets you pick which SVD algorithm from the LAPACK library gets the aux.
    • 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)

Okay so for the four methods above, there can be a significant difference in the accuracy of small singular values if the snapshots vary in norm so that X is severely ill-conditioned. So if small (smaller than EPS*||X||) singular values are of interest and JOBS=='N', then options (3, 4) give the most accurate results, where option 4 is slightly better and comes with stronger theoretical background. Science slay. 🔬

And if JOBS=='S', i.e. the columns of X will be normalized, then all methods give nearly equally accurate results. Everyone wins, we love to see it. 💖

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. Basically it's 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. Noted, no cap.

< 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. Respect the leading dimension, it's literally the main character.

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. There it is. 🎁

LDY 📏

  • LDY (input) INTEGER, LDY >= N — the leading dimension of the array Y.

NRNK 👑

  • NRNK (input) INTEGER — determines the mode for how to compute the numerical rank, i.e. how to truncate small singular values of the input matrix X. On input, if:
    • NRNK = -1 — the i-th singular value sigma(i) is truncated if sigma(i) <= TOL*sigma(1). This option is recommended. ⭐ Use this one, bestie.
    • NRNK = -2 — the 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. 🔬 (so like, careful.)

The numerical rank can be enforced by using a 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<N) entries of REIG contain the imaginary parts of the computed eigenvalues
REIG(1:K) + sqrt(-1)*IMEIG(1:K).

The eigenvalues are determined as follows — okay so listen:

  • 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. Positivity first, always. 💅

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. 😭 (the formulas are giving "trust the process")

B

  • B — (output) REAL(KIND=WP) MIN(M,N)-by-(N-1) array. Bestie, here’s the lowkey important part: if JOBF =='R', then B(1:N,1:K) contains A*U(:,1:K), and it can be used for computing the refined vectors — see the references for the full lore. If JOBF == 'E', then 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 — (input) INTEGER, LDB >= MIN(M,N). The leading dimension of the array B.
  • 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 — (input) INTEGER, LDV >= N-1. The leading dimension of the array V.
  • 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 — (input) INTEGER, LDS >= N-1. The leading dimension of the array S.
  • 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 — (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 — (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 — (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 — (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

So 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. [1] P. Schmid: Dynamic mode decomposition of numerical and experimental data, Journal of Fluid Mechanics 656, 5-28, 2010.
  2. [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. [3] Z. Drmac: A LAPACK implementation of the Dynamic Mode Decomposition I. Technical report. AIMDyn Inc. and LAPACK Working Note 298.
  4. [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 — (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

okay so we're still deep in the argument list and it's giving "the lore continues" 💅 — no skipping, bestie, every flag is coming with us.

  • 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 — 🫶 this is literally where the excerpt stops, bestie; no type or description is given here, and I'm not gonna make one up. The next part picks it right up.
  • RWORK (REAL(KIND=WP) LRWORK-by-1 array, workspace/output) — okay so on exit, RWORK(1:N) is holding the singular values of X (that's when JOBS=='N') or of the column scaled X (when JOBS=='S' or 'C'). And if WHTSVD==4, then RWORK(N+1) and RWORK(N+2) contain the scaling factor RWORK(N+2)/RWORK(N+1) that got used to scale X and Y so the SVD of X doesn't overflow 💅
    • ngl that scaling factor can be interesting to you if the scaling option is off and you want as many as possible smallest eigenvalues to the highest feasible accuracy 👀
    • and if the call to ZGEDMD is only a workspace query, then RWORK(1) contains the minimal workspace length. See the description of LRWORK ✨
  • LRWORK (INTEGER, input) — the minimal length of the workspace vector RWORK, bestie. It's calculated like this:
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. So:

  • 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

And 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). slay 🌸

  • IWORK (INTEGER LIWORK-by-1 array, workspace/output) — workspace that's only required if WHTSVD equals 2, 3 or 4 (peep the description of WHTSVD for that one 🫶). And 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 (INTEGER, input) — the minimal length of the workspace vector IWORK. The tea ☕:
    • 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)

And 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 (INTEGER, output) — the vibe check report at the end 💀 it's giving:
    • -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

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. periodt ✨

🎯 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

JOBT (input) CHARACTER*1 — tells us whether we're handing back the upper triangular factor from the QR factorization, bestie. The deets:

  • 'R' — the matrix R of the QR factorization of the data snapshot matrix F gets returned in the array Y. Peep the description of Y and Further details.
  • 'N' — nope, R is not returned. Not today.

JOBF

JOBF (input) CHARACTER*1 — controls whether we stash the info needed for post-processing (like, say, computing refined Ritz vectors 👑). Options:

  • '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. Zero. Nada.
Heads up 🫶: to actually be useful on exit, this option needs JOBQ='Q'.

WHTSVD

WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 } — this is your SVD algorithm picker from the LAPACK library ✨:

  1. 1 :: ZGESVD (the QR SVD algorithm)
  2. 2 :: ZGESDD (the Divide and Conquer algorithm; if enough workspace available, this is the fastest option)
  3. 3 :: ZGESVDQ (the preconditioned QR SVD ; this and 4 are the most accurate options)
  4. 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. Iconic.

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. Basically 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. Don't skip leg day.

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 — decides the mode for computing the numerical rank, i.e. how we truncate small singular values of the input matrix X. On input:

  • NRNK = -1 :: the i-th singular value sigma(i) is truncated if sigma(i) <= TOL*sigma(1). This option is recommended. 👑
  • NRNK = -2 :: the 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 a 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. Be nice to your tiny sigmas 🌸.

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. 😭 (residuals never lie)

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 — okay so this is the minimal length of the ZWORK workspace vector, bestie. It gets calculated like this:

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 (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 a 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 (input) INTEGER — the minimal length of the WORK workspace vector. 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 (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 (input) INTEGER — the minimal length of the IWORK workspace vector. 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 (output) INTEGER — the exit-status tea, spilled:
    • -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.

That’s the LAPACK Version 3.12.0 gedmd(3) manual page, bestie — source/version footer included, date not included, and we’re out 💖✨💅

📦 where this page comes from