Library Functions Manual

NAME¶

gedmd - DMD driver, Dynamic Mode Decomposition

SYNOPSIS¶

Functions¶

subroutine CGEDMD (jobs, jobz, jobr, jobf, whtsvd, m, n, x, ldx, y, ldy, nrnk, tol, k, eigs, z, ldz, res, b, ldb, w, ldw, s, lds, zwork, lzwork, rwork, lrwork, iwork, liwork, info)
CGEDMD computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. subroutine CGEDMDQ (jobs, jobz, jobr, jobq, jobt, jobf, whtsvd, m, n, f, ldf, x, ldx, y, ldy, nrnk, tol, k, eigs, z, ldz, res, b, ldb, v, ldv, s, lds, zwork, lzwork, work, lwork, iwork, liwork, info)
CGEDMDQ computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. subroutine DGEDMD (jobs, jobz, jobr, jobf, whtsvd, m, n, x, ldx, y, ldy, nrnk, tol, k, reig, imeig, z, ldz, res, b, ldb, w, ldw, s, lds, work, lwork, iwork, liwork, info)
DGEDMD computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. subroutine DGEDMDQ (jobs, jobz, jobr, jobq, jobt, jobf, whtsvd, m, n, f, ldf, x, ldx, y, ldy, nrnk, tol, k, reig, imeig, z, ldz, res, b, ldb, v, ldv, s, lds, work, lwork, iwork, liwork, info)
DGEDMDQ computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. subroutine SGEDMD (jobs, jobz, jobr, jobf, whtsvd, m, n, x, ldx, y, ldy, nrnk, tol, k, reig, imeig, z, ldz, res, b, ldb, w, ldw, s, lds, work, lwork, iwork, liwork, info)
SGEDMD computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. subroutine SGEDMDQ (jobs, jobz, jobr, jobq, jobt, jobf, whtsvd, m, n, f, ldf, x, ldx, y, ldy, nrnk, tol, k, reig, imeig, z, ldz, res, b, ldb, v, ldv, s, lds, work, lwork, iwork, liwork, info)
SGEDMDQ computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. subroutine ZGEDMD (jobs, jobz, jobr, jobf, whtsvd, m, n, x, ldx, y, ldy, nrnk, tol, k, eigs, z, ldz, res, b, ldb, w, ldw, s, lds, zwork, lzwork, rwork, lrwork, iwork, liwork, info)
ZGEDMD computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices. subroutine ZGEDMDQ (jobs, jobz, jobr, jobq, jobt, jobf, whtsvd, m, n, f, ldf, x, ldx, y, ldy, nrnk, tol, k, eigs, z, ldz, res, b, ldb, v, ldv, s, lds, zwork, lzwork, work, lwork, iwork, liwork, info)
ZGEDMDQ computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.

Detailed Description¶

Function Documentation¶

subroutine CGEDMD (character, intent(in) jobs, character, intent(in) jobz, character, intent(in) jobr, character, intent(in) jobf, integer, intent(in) whtsvd, integer, intent(in) m, integer, intent(in) n, complex(kind=wp), dimension(ldx,), intent(inout) x, integer, intent(in) ldx, complex(kind=wp), dimension(ldy,), intent(inout) y, integer, intent(in) ldy, integer, intent(in) nrnk, real(kind=wp), intent(in) tol, integer, intent(out) k, complex(kind=wp), dimension(), intent(out) eigs, complex(kind=wp), dimension(ldz,), intent(out) z, integer, intent(in) ldz, real(kind=wp), dimension(), intent(out) res, complex(kind=wp), dimension(ldb,), intent(out) b, integer, intent(in) ldb, complex(kind=wp), dimension(ldw,), intent(out) w, integer, intent(in) ldw, complex(kind=wp), dimension(lds,), intent(out) s, integer, intent(in) lds, complex(kind=wp), dimension(), intent(out) zwork, integer, intent(in) lzwork, real(kind=wp), dimension(), intent(out) rwork, integer, intent(in) lrwork, integer, dimension(*), intent(out) iwork, integer, intent(in) liwork, integer, intent(out) info)¶

CGEDMD computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.

Purpose:

!>    CGEDMD computes the Dynamic Mode Decomposition (DMD) for
!>    a pair of data snapshot matrices. For the input matrices
!>    X and Y such that Y = A*X with an unaccessible matrix
!>    A, CGEDMD computes a certain number of Ritz pairs of A using
!>    the standard Rayleigh-Ritz extraction from a subspace of
!>    range(X) that is determined using the leading left singular
!>    vectors of X. Optionally, CGEDMD returns the residuals
!>    of the computed Ritz pairs, the information needed for
!>    a refinement of the Ritz vectors, or the eigenvectors of
!>    the Exact DMD.
!>    For further details see the references listed
!>    below. For more details of the implementation see [3].
!>

References:

!>    [1] P. Schmid: Dynamic mode decomposition of numerical
!>        and experimental data,
!>        Journal of Fluid Mechanics 656, 5-28, 2010.
!>    [2] Z. Drmac, I. Mezic, R. Mohr: Data driven modal
!>        decompositions: analysis and enhancements,
!>        SIAM J. on Sci. Comp. 40 (4), A2253-A2285, 2018.
!>    [3] Z. Drmac: A LAPACK implementation of the Dynamic
!>        Mode Decomposition I. Technical report. AIMDyn Inc.
!>        and LAPACK Working Note 298.
!>    [4] J. Tu, C. W. Rowley, D. M. Luchtenburg, S. L.
!>        Brunton, N. Kutz: On Dynamic Mode Decomposition:
!>        Theory and Applications, Journal of Computational
!>        Dynamics 1(2), 391 -421, 2014.
!>

Developed and supported by:

!>    Developed and coded by Zlatko Drmac, Faculty of Science,
!>    University of Zagreb;  drmac@math.hr
!>    In cooperation with
!>    AIMdyn Inc., Santa Barbara, CA.
!>    and supported by
!>    - DARPA SBIR project  Contract No: W31P4Q-21-C-0007
!>    - DARPA PAI project  Contract No: HR0011-18-9-0033
!>    - DARPA MoDyL project 
!>    Contract No: HR0011-16-C-0116
!>    Any opinions, findings and conclusions or recommendations
!>    expressed in this material are those of the author and
!>    do not necessarily reflect the views of the DARPA SBIR
!>    Program Office
!>

Distribution Statement A:

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

Parameters

JOBS

!>    JOBS (input) CHARACTER*1
!>    Determines whether the initial data snapshots are scaled
!>    by a diagonal matrix.
!>    'S' :: The data snapshots matrices X and Y are multiplied
!>           with a diagonal matrix D so that X*D has unit
!>           nonzero columns (in the Euclidean 2-norm)
!>    'C' :: The snapshots are scaled as with the 'S' option.
!>           If it is found that an i-th column of X is zero
!>           vector and the corresponding i-th column of Y is
!>           non-zero, then the i-th column of Y is set to
!>           zero and a warning flag is raised.
!>    'Y' :: The data snapshots matrices X and Y are multiplied
!>           by a diagonal matrix D so that Y*D has unit
!>           nonzero columns (in the Euclidean 2-norm)
!>    'N' :: No data scaling.
!>

JOBZ

!>    JOBZ (input) CHARACTER*1
!>    Determines whether the eigenvectors (Koopman modes) will
!>    be computed.
!>    'V' :: The eigenvectors (Koopman modes) will be computed
!>           and returned in the matrix Z.
!>           See the description of Z.
!>    'F' :: The eigenvectors (Koopman modes) will be returned
!>           in factored form as the product X(:,1:K)*W, where X
!>           contains a POD basis (leading left singular vectors
!>           of the data matrix X) and W contains the eigenvectors
!>           of the corresponding Rayleigh quotient.
!>           See the descriptions of K, X, W, Z.
!>    'N' :: The eigenvectors are not computed.
!>

JOBR

!>    JOBR (input) CHARACTER*1
!>    Determines whether to compute the residuals.
!>    'R' :: The residuals for the computed eigenpairs will be
!>           computed and stored in the array RES.
!>           See the description of RES.
!>           For this option to be legal, JOBZ must be 'V'.
!>    'N' :: The residuals are not computed.
!>

JOBF

!>    JOBF (input) CHARACTER*1
!>    Specifies whether to store information needed for post-
!>    processing (e.g. computing refined Ritz vectors)
!>    'R' :: The matrix needed for the refinement of the Ritz
!>           vectors is computed and stored in the array B.
!>           See the description of B.
!>    'E' :: The unscaled eigenvectors of the Exact DMD are
!>           computed and returned in the array B. See the
!>           description of B.
!>    'N' :: No eigenvector refinement data is computed.
!>

WHTSVD

!>    WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 }
!>    Allows for a selection of the SVD algorithm from the
!>    LAPACK library.
!>    1 :: CGESVD (the QR SVD algorithm)
!>    2 :: CGESDD (the Divide and Conquer algorithm; if enough
!>         workspace available, this is the fastest option)
!>    3 :: CGESVDQ (the preconditioned QR SVD  ; this and 4
!>         are the most accurate options)
!>    4 :: CGEJSV (the preconditioned Jacobi SVD; this and 3
!>         are the most accurate options)
!>    For the four methods above, a significant difference in
!>    the accuracy of small singular values is possible if
!>    the snapshots vary in norm so that X is severely
!>    ill-conditioned. If small (smaller than EPS*||X||)
!>    singular values are of interest and JOBS=='N',  then
!>    the options (3, 4) give the most accurate results, where
!>    the option 4 is slightly better and with stronger
!>    theoretical background.
!>    If JOBS=='S', i.e. the columns of X will be normalized,
!>    then all methods give nearly equally accurate results.
!>

!>    M (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) 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 (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 (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 (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 (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 (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 (workspace/output) COMPLEX(KIND=WP) N-by-N array
!>    The array S(1:K,1:K) is used for the matrix Rayleigh
!>    quotient. This content is overwritten during
!>    the eigenvalue decomposition by CGEEV.
!>    See the description of K.
!>

LDS

!>    LDS (input) INTEGER, LDS >= N
!>    The leading dimension of the array S.
!>

ZWORK

!>    ZWORK (workspace/output) COMPLEX(KIND=WP) LZWORK-by-1 array
!>    ZWORK is used as complex workspace in the complex SVD, as
!>    specified by WHTSVD (1,2, 3 or 4) and for CGEEV for computing
!>    the eigenvalues of a Rayleigh quotient.
!>    If the call to CGEDMD is only workspace query, then
!>    ZWORK(1) contains the minimal complex workspace length and
!>    ZWORK(2) is the optimal complex workspace length.
!>    Hence, the length of work is at least 2.
!>    See the description of LZWORK.
!>

LZWORK

!>    LZWORK (input) INTEGER
!>    The minimal length of the workspace vector ZWORK.
!>    LZWORK is calculated as MAX(LZWORK_SVD, LZWORK_CGEEV),
!>    where LZWORK_CGEEV = MAX( 1, 2*N )  and the minimal
!>    LZWORK_SVD is calculated as follows
!>    If WHTSVD == 1 :: CGESVD ::
!>       LZWORK_SVD = MAX(1,2*MIN(M,N)+MAX(M,N))
!>    If WHTSVD == 2 :: CGESDD ::
!>       LZWORK_SVD = 2*MIN(M,N)*MIN(M,N)+2*MIN(M,N)+MAX(M,N)
!>    If WHTSVD == 3 :: CGESVDQ ::
!>       LZWORK_SVD = obtainable by a query
!>    If WHTSVD == 4 :: CGEJSV ::
!>       LZWORK_SVD = obtainable by a query
!>    If on entry LZWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace lengths and returns them in
!>    LZWORK(1) and LZWORK(2), respectively.
!>

RWORK

!>    RWORK (workspace/output) REAL(KIND=WP) LRWORK-by-1 array
!>    On exit, RWORK(1:N) contains the singular values of
!>    X (for JOBS=='N') or column scaled X (JOBS=='S', 'C').
!>    If WHTSVD==4, then RWORK(N+1) and RWORK(N+2) contain
!>    scaling factor RWORK(N+2)/RWORK(N+1) used to scale X
!>    and Y to avoid overflow in the SVD of X.
!>    This may be of interest if the scaling option is off
!>    and as many as possible smallest eigenvalues are
!>    desired to the highest feasible accuracy.
!>    If the call to CGEDMD is only workspace query, then
!>    RWORK(1) contains the minimal workspace length.
!>    See the description of LRWORK.
!>

LRWORK

!>    LRWORK (input) INTEGER
!>    The minimal length of the workspace vector RWORK.
!>    LRWORK is calculated as follows:
!>    LRWORK = MAX(1, N+LRWORK_SVD,N+LRWORK_CGEEV), where
!>    LRWORK_CGEEV = MAX(1,2*N) and RWORK_SVD is the real workspace
!>    for the SVD subroutine determined by the input parameter
!>    WHTSVD.
!>    If WHTSVD == 1 :: CGESVD ::
!>       LRWORK_SVD = 5*MIN(M,N)
!>    If WHTSVD == 2 :: CGESDD ::
!>       LRWORK_SVD =  MAX(5*MIN(M,N)*MIN(M,N)+7*MIN(M,N),
!>       2*MAX(M,N)*MIN(M,N)+2*MIN(M,N)*MIN(M,N)+MIN(M,N) ) )
!>    If WHTSVD == 3 :: CGESVDQ ::
!>       LRWORK_SVD = obtainable by a query
!>    If WHTSVD == 4 :: CGEJSV ::
!>       LRWORK_SVD = obtainable by a query
!>    If on entry LRWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    real workspace length and returns it in RWORK(1).
!>

IWORK

!>    IWORK (workspace/output) INTEGER LIWORK-by-1 array
!>    Workspace that is required only if WHTSVD equals
!>    2 , 3 or 4. (See the description of WHTSVD).
!>    If on entry LWORK =-1 or LIWORK=-1, then the
!>    minimal length of IWORK is computed and returned in
!>    IWORK(1). See the description of LIWORK.
!>

LIWORK

!>    LIWORK (input) INTEGER
!>    The minimal length of the workspace vector IWORK.
!>    If WHTSVD == 1, then only IWORK(1) is used; LIWORK >=1
!>    If WHTSVD == 2, then LIWORK >= MAX(1,8*MIN(M,N))
!>    If WHTSVD == 3, then LIWORK >= MAX(1,M+N-1)
!>    If WHTSVD == 4, then LIWORK >= MAX(3,M+3*N)
!>    If on entry LIWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace lengths for  ZWORK, RWORK and
!>    IWORK. See the descriptions of ZWORK, RWORK and IWORK.
!>

INFO

!>    INFO (output) INTEGER
!>    -i < 0 :: On entry, the i-th argument had an
!>              illegal value
!>       = 0 :: Successful return.
!>       = 1 :: Void input. Quick exit (M=0 or N=0).
!>       = 2 :: The SVD computation of X did not converge.
!>              Suggestion: Check the input data and/or
!>              repeat with different WHTSVD.
!>       = 3 :: The computation of the eigenvalues did not
!>              converge.
!>       = 4 :: If data scaling was requested on input and
!>              the procedure found inconsistency in the data
!>              such that for some column index i,
!>              X(:,i) = 0 but Y(:,i) /= 0, then Y(:,i) is set
!>              to zero if JOBS=='C'. The computation proceeds
!>              with original or modified data and warning
!>              flag is set with INFO=4.
!>

Author

Zlatko Drmac

Definition at line 496 of file cgedmd.f90.

subroutine CGEDMDQ (character, intent(in) jobs, character, intent(in) jobz, character, intent(in) jobr, character, intent(in) jobq, character, intent(in) jobt, character, intent(in) jobf, integer, intent(in) whtsvd, integer, intent(in) m, integer, intent(in) n, complex(kind=wp), dimension(ldf,), intent(inout) f, integer, intent(in) ldf, complex(kind=wp), dimension(ldx,), intent(out) x, integer, intent(in) ldx, complex(kind=wp), dimension(ldy,), intent(out) y, integer, intent(in) ldy, integer, intent(in) nrnk, real(kind=wp), intent(in) tol, integer, intent(out) k, complex(kind=wp), dimension(), intent(out) eigs, complex(kind=wp), dimension(ldz,), intent(out) z, integer, intent(in) ldz, real(kind=wp), dimension(), intent(out) res, complex(kind=wp), dimension(ldb,), intent(out) b, integer, intent(in) ldb, complex(kind=wp), dimension(ldv,), intent(out) v, integer, intent(in) ldv, complex(kind=wp), dimension(lds,), intent(out) s, integer, intent(in) lds, complex(kind=wp), dimension(), intent(out) zwork, integer, intent(in) lzwork, real(kind=wp), dimension(), intent(out) work, integer, intent(in) lwork, integer, dimension(), intent(out) iwork, integer, intent(in) liwork, integer, intent(out) info)¶

CGEDMDQ computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.

Purpose:

!>    CGEDMDQ computes the Dynamic Mode Decomposition (DMD) for
!>    a pair of data snapshot matrices, using a QR factorization
!>    based compression of the data. For the input matrices
!>    X and Y such that Y = A*X with an unaccessible matrix
!>    A, CGEDMDQ computes a certain number of Ritz pairs of A using
!>    the standard Rayleigh-Ritz extraction from a subspace of
!>    range(X) that is determined using the leading left singular 
!>    vectors of X. Optionally, CGEDMDQ returns the residuals 
!>    of the computed Ritz pairs, the information needed for
!>    a refinement of the Ritz vectors, or the eigenvectors of
!>    the Exact DMD.
!>    For further details see the references listed
!>    below. For more details of the implementation see [3].      
!>

References:

!>    [1] P. Schmid: Dynamic mode decomposition of numerical
!>        and experimental data,
!>        Journal of Fluid Mechanics 656, 5-28, 2010.
!>    [2] Z. Drmac, I. Mezic, R. Mohr: Data driven modal
!>        decompositions: analysis and enhancements,
!>        SIAM J. on Sci. Comp. 40 (4), A2253-A2285, 2018.
!>    [3] Z. Drmac: A LAPACK implementation of the Dynamic
!>        Mode Decomposition I. Technical report. AIMDyn Inc.
!>        and LAPACK Working Note 298.      
!>    [4] J. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. 
!>        Brunton, N. Kutz: On Dynamic Mode Decomposition:
!>        Theory and Applications, Journal of Computational
!>        Dynamics 1(2), 391 -421, 2014.
!>

Developed and supported by:

!>    Developed and coded by Zlatko Drmac, Faculty of Science,
!>    University of Zagreb;  drmac@math.hr
!>    In cooperation with
!>    AIMdyn Inc., Santa Barbara, CA.
!>    and supported by
!>    - DARPA SBIR project  Contract No: W31P4Q-21-C-0007
!>    - DARPA PAI project  Contract No: HR0011-18-9-0033
!>    - DARPA MoDyL project 
!>    Contract No: HR0011-16-C-0116
!>    Any opinions, findings and conclusions or recommendations 
!>    expressed in this material are those of the author and 
!>    do not necessarily reflect the views of the DARPA SBIR 
!>    Program Office.      
!>

Developed and supported by:

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

Parameters

JOBS

!>    JOBS (input) CHARACTER*1
!>    Determines whether the initial data snapshots are scaled
!>    by a diagonal matrix. The data snapshots are the columns
!>    of F. The leading N-1 columns of F are denoted X and the
!>    trailing N-1 columns are denoted Y. 
!>    'S' :: The data snapshots matrices X and Y are multiplied
!>           with a diagonal matrix D so that X*D has unit
!>           nonzero columns (in the Euclidean 2-norm)
!>    'C' :: The snapshots are scaled as with the 'S' option.
!>           If it is found that an i-th column of X is zero
!>           vector and the corresponding i-th column of Y is
!>           non-zero, then the i-th column of Y is set to
!>           zero and a warning flag is raised.
!>    'Y' :: The data snapshots matrices X and Y are multiplied
!>           by a diagonal matrix D so that Y*D has unit
!>           nonzero columns (in the Euclidean 2-norm)    
!>    'N' :: No data scaling.   
!>

JOBZ

!>    JOBZ (input) CHARACTER*1
!>    Determines whether the eigenvectors (Koopman modes) will
!>    be computed.
!>    'V' :: The eigenvectors (Koopman modes) will be computed
!>           and returned in the matrix Z.
!>           See the description of Z.
!>    'F' :: The eigenvectors (Koopman modes) will be returned
!>           in factored form as the product Z*V, where Z
!>           is orthonormal and V contains the eigenvectors
!>           of the corresponding Rayleigh quotient.
!>           See the descriptions of F, V, Z.
!>    'Q' :: The eigenvectors (Koopman modes) will be returned
!>           in factored form as the product Q*Z, where Z
!>           contains the eigenvectors of the compression of the
!>           underlying discretised operator onto the span of
!>           the data snapshots. See the descriptions of F, V, Z.   
!>           Q is from the inital QR facorization.    
!>    'N' :: The eigenvectors are not computed.  
!>

JOBR

!>    JOBR (input) CHARACTER*1 
!>    Determines whether to compute the residuals.
!>    'R' :: The residuals for the computed eigenpairs will
!>           be computed and stored in the array RES.
!>           See the description of RES.
!>           For this option to be legal, JOBZ must be 'V'.
!>    'N' :: The residuals are not computed.
!>

JOBQ

!>    JOBQ (input) CHARACTER*1 
!>    Specifies whether to explicitly compute and return the
!>    unitary matrix from the QR factorization.
!>    'Q' :: The matrix Q of the QR factorization of the data
!>           snapshot matrix is computed and stored in the
!>           array F. See the description of F.       
!>    'N' :: The matrix Q is not explicitly computed.
!>

JOBT

!>    JOBT (input) CHARACTER*1 
!>    Specifies whether to return the upper triangular factor
!>    from the QR factorization.
!>    'R' :: The matrix R of the QR factorization of the data 
!>           snapshot matrix F is returned in the array Y.
!>           See the description of Y and Further details.       
!>    'N' :: The matrix R is not returned. 
!>

JOBF

!>    JOBF (input) CHARACTER*1
!>    Specifies whether to store information needed for post-
!>    processing (e.g. computing refined Ritz vectors)
!>    'R' :: The matrix needed for the refinement of the Ritz
!>           vectors is computed and stored in the array B.
!>           See the description of B.
!>    'E' :: The unscaled eigenvectors of the Exact DMD are 
!>           computed and returned in the array B. See the
!>           description of B.
!>    'N' :: No eigenvector refinement data is computed.   
!>    To be useful on exit, this option needs JOBQ='Q'.    
!>

WHTSVD

!>    WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 }
!>    Allows for a selection of the SVD algorithm from the
!>    LAPACK library.
!>    1 :: CGESVD (the QR SVD algorithm)
!>    2 :: CGESDD (the Divide and Conquer algorithm; if enough
!>         workspace available, this is the fastest option)
!>    3 :: CGESVDQ (the preconditioned QR SVD  ; this and 4
!>         are the most accurate options)
!>    4 :: CGEJSV (the preconditioned Jacobi SVD; this and 3
!>         are the most accurate options)
!>    For the four methods above, a significant difference in
!>    the accuracy of small singular values is possible if
!>    the snapshots vary in norm so that X is severely
!>    ill-conditioned. If small (smaller than EPS*||X||)
!>    singular values are of interest and JOBS=='N',  then
!>    the options (3, 4) give the most accurate results, where
!>    the option 4 is slightly better and with stronger 
!>    theoretical background.
!>    If JOBS=='S', i.e. the columns of X will be normalized,
!>    then all methods give nearly equally accurate results.
!>

!>    M (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 (input/output) COMPLEX(KIND=WP) M-by-N array
!>    > On entry,
!>    the columns of F are the sequence of data snapshots 
!>    from a single trajectory, taken at equidistant discrete
!>    times. It is assumed that the column norms of F are 
!>    in the range of the normalized floating point numbers. 
!>    < On exit,
!>    If JOBQ == 'Q', the array F contains the orthogonal 
!>    matrix/factor of the QR factorization of the initial 
!>    data snapshots matrix F. See the description of JOBQ. 
!>    If JOBQ == 'N', the entries in F strictly below the main
!>    diagonal contain, column-wise, the information on the 
!>    Householder vectors, as returned by CGEQRF. The 
!>    remaining information to restore the orthogonal matrix
!>    of the initial QR factorization is stored in ZWORK(1:MIN(M,N)). 
!>    See the description of ZWORK.
!>

LDF

!>    LDF (input) INTEGER, LDF >= M 
!>    The leading dimension of the array F.
!>

!>    X (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 (workspace/output) COMPLEX(KIND=WP) MIN(M,N)-by-(N) array
!>    Y is used as workspace to hold representations of the
!>    trailing N-1 snapshots in the orthonormal basis computed
!>    in the QR factorization of F.
!>    On exit, 
!>    If JOBT == 'R', Y contains the MIN(M,N)-by-N upper
!>    triangular factor from the QR factorization of the data
!>    snapshot matrix F.
!>

LDY

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

NRNK

!>    NRNK (input) INTEGER
!>    Determines the mode how to compute the numerical rank,
!>    i.e. how to truncate small singular values of the input
!>    matrix X. On input, if
!>    NRNK = -1 :: i-th singular value sigma(i) is truncated
!>                 if sigma(i) <= TOL*sigma(1)
!>                 This option is recommended.
!>    NRNK = -2 :: i-th singular value sigma(i) is truncated
!>                 if sigma(i) <= TOL*sigma(i-1)
!>                 This option is included for R&D purposes.
!>                 It requires highly accurate SVD, which
!>                 may not be feasible.      
!>    The numerical rank can be enforced by using positive 
!>    value of NRNK as follows: 
!>    0 < NRNK <= N-1 :: at most NRNK largest singular values
!>    will be used. If the number of the computed nonzero
!>    singular values is less than NRNK, then only those
!>    nonzero values will be used and the actually used
!>    dimension is less than NRNK. The actual number of
!>    the nonzero singular values is returned in the variable
!>    K. See the description of K.
!>

TOL

!>    TOL (input) REAL(KIND=WP), 0 <= TOL < 1
!>    The tolerance for truncating small singular values.
!>    See the description of NRNK.  
!>

!>    K (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 (workspace/output) COMPLEX(KIND=WP)  M-by-(N-1) array
!>    If JOBZ =='V' then Z contains the Ritz vectors. Z(:,i)
!>    is an eigenvector of the i-th Ritz value; ||Z(:,i)||_2=1.
!>    If JOBZ == 'F', then the Z(:,i)'s are given implicitly as
!>    Z*V, where Z contains orthonormal matrix (the product of
!>    Q from the initial QR factorization and the SVD/POD_basis
!>    returned by CGEDMD in X) and the second factor (the 
!>    eigenvectors of the Rayleigh quotient) is in the array V, 
!>    as returned by CGEDMD. That is,  X(:,1:K)*V(:,i)
!>    is an eigenvector corresponding to EIGS(i). The columns 
!>    of V(1:K,1:K) are the computed eigenvectors of the 
!>    K-by-K Rayleigh quotient.  
!>    See the descriptions of EIGS, X and V.      
!>

LDZ

!>    LDZ (input) INTEGER , LDZ >= M
!>    The leading dimension of the array Z.
!>

RES

!>    RES (output) REAL(KIND=WP) (N-1)-by-1 array
!>    RES(1:K) contains the residuals for the K computed 
!>    Ritz pairs, 
!>    RES(i) = || A * Z(:,i) - EIGS(i)*Z(:,i))||_2.
!>    See the description of EIGS and Z.      
!>

!>    B (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 (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 (output) COMPLEX(KIND=WP) (N-1)-by-(N-1) array
!>    The array S(1:K,1:K) is used for the matrix Rayleigh
!>    quotient. This content is overwritten during
!>    the eigenvalue decomposition by CGEEV.
!>    See the description of K.
!>

LDS

!>    LDS (input) INTEGER, LDS >= N-1        
!>    The leading dimension of the array S.
!>

LZWORK

!>    ZWORK (workspace/output) COMPLEX(KIND=WP) LWORK-by-1 array
!>    On exit, 
!>    ZWORK(1:MIN(M,N)) contains the scalar factors of the 
!>    elementary reflectors as returned by CGEQRF of the 
!>    M-by-N input matrix F.   
!>    If the call to CGEDMDQ is only workspace query, then
!>    ZWORK(1) contains the minimal complex workspace length and
!>    ZWORK(2) is the optimal complex workspace length. 
!>    Hence, the length of work is at least 2.
!>    See the description of LZWORK.      
!>

LZWORK

!>    LZWORK (input) INTEGER
!>    The minimal length of the  workspace vector ZWORK.
!>    LZWORK is calculated as follows:
!>    Let MLWQR  = N (minimal workspace for CGEQRF[M,N])
!>        MLWDMD = minimal workspace for CGEDMD (see the
!>                 description of LWORK in CGEDMD)
!>        MLWMQR = N (minimal workspace for 
!>                   ZUNMQR['L','N',M,N,N])
!>        MLWGQR = N (minimal workspace for ZUNGQR[M,N,N])
!>        MINMN  = MIN(M,N)      
!>    Then
!>    LZWORK = MAX(2, MIN(M,N)+MLWQR, MINMN+MLWDMD)
!>    is further updated as follows:
!>       if   JOBZ == 'V' or JOBZ == 'F' THEN 
!>            LZWORK = MAX( LZWORK, MINMN+MLWMQR )
!>       if   JOBQ == 'Q' THEN
!>            LZWORK = MAX( ZLWORK, MINMN+MLWGQR)      
!>

WORK

!>    WORK (workspace/output) REAL(KIND=WP) LWORK-by-1 array
!>    On exit,
!>    WORK(1:N-1) contains the singular values of 
!>    the input submatrix F(1:M,1:N-1).
!>    If the call to CGEDMDQ is only workspace query, then
!>    WORK(1) contains the minimal workspace length and
!>    WORK(2) is the optimal workspace length. hence, the
!>    length of work is at least 2.
!>    See the description of LWORK.
!>

LWORK

!>    LWORK (input) INTEGER
!>    The minimal length of the  workspace vector WORK.
!>    LWORK is the same as in CGEDMD, because in CGEDMDQ
!>    only CGEDMD requires real workspace for snapshots
!>    of dimensions MIN(M,N)-by-(N-1).
!>    If on entry LWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace lengths for both WORK and
!>    IWORK. See the descriptions of WORK and IWORK.          
!>

IWORK

!>    IWORK (workspace/output) INTEGER LIWORK-by-1 array
!>    Workspace that is required only if WHTSVD equals
!>    2 , 3 or 4. (See the description of WHTSVD).
!>    If on entry LWORK =-1 or LIWORK=-1, then the
!>    minimal length of IWORK is computed and returned in
!>    IWORK(1). See the description of LIWORK.
!>

LIWORK

!>    LIWORK (input) INTEGER
!>    The minimal length of the workspace vector IWORK.
!>    If WHTSVD == 1, then only IWORK(1) is used; LIWORK >=1
!>    Let M1=MIN(M,N), N1=N-1. Then      
!>    If WHTSVD == 2, then LIWORK >= MAX(1,8*MIN(M,N))
!>    If WHTSVD == 3, then LIWORK >= MAX(1,M+N-1)
!>    If WHTSVD == 4, then LIWORK >= MAX(3,M+3*N)
!>    If on entry LIWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace lengths for both WORK and
!>    IWORK. See the descriptions of WORK and IWORK.
!>

INFO

!>    INFO (output) INTEGER
!>    -i < 0 :: On entry, the i-th argument had an
!>              illegal value
!>       = 0 :: Successful return.
!>       = 1 :: Void input. Quick exit (M=0 or N=0).
!>       = 2 :: The SVD computation of X did not converge.
!>              Suggestion: Check the input data and/or
!>              repeat with different WHTSVD.
!>       = 3 :: The computation of the eigenvalues did not
!>              converge.
!>       = 4 :: If data scaling was requested on input and
!>              the procedure found inconsistency in the data
!>              such that for some column index i,
!>              X(:,i) = 0 but Y(:,i) /= 0, then Y(:,i) is set
!>              to zero if JOBS=='C'. The computation proceeds
!>              with original or modified data and warning
!>              flag is set with INFO=4.  
!>

Author

Zlatko Drmac

Definition at line 552 of file cgedmdq.f90.

subroutine DGEDMD (character, intent(in) jobs, character, intent(in) jobz, character, intent(in) jobr, character, intent(in) jobf, integer, intent(in) whtsvd, integer, intent(in) m, integer, intent(in) n, real(kind=wp), dimension(ldx,), intent(inout) x, integer, intent(in) ldx, real(kind=wp), dimension(ldy,), intent(inout) y, integer, intent(in) ldy, integer, intent(in) nrnk, real(kind=wp), intent(in) tol, integer, intent(out) k, real(kind=wp), dimension(), intent(out) reig, real(kind=wp), dimension(), intent(out) imeig, real(kind=wp), dimension(ldz,), intent(out) z, integer, intent(in) ldz, real(kind=wp), dimension(), intent(out) res, real(kind=wp), dimension(ldb,), intent(out) b, integer, intent(in) ldb, real(kind=wp), dimension(ldw,), intent(out) w, integer, intent(in) ldw, real(kind=wp), dimension(lds,), intent(out) s, integer, intent(in) lds, real(kind=wp), dimension(), intent(out) work, integer, intent(in) lwork, integer, dimension(*), intent(out) iwork, integer, intent(in) liwork, integer, intent(out) info)¶

DGEDMD computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.

Purpose:

!>    DGEDMD computes the Dynamic Mode Decomposition (DMD) for
!>    a pair of data snapshot matrices. For the input matrices
!>    X and Y such that Y = A*X with an unaccessible matrix
!>    A, DGEDMD computes a certain number of Ritz pairs of A using
!>    the standard Rayleigh-Ritz extraction from a subspace of
!>    range(X) that is determined using the leading left singular
!>    vectors of X. Optionally, DGEDMD returns the residuals
!>    of the computed Ritz pairs, the information needed for
!>    a refinement of the Ritz vectors, or the eigenvectors of
!>    the Exact DMD.
!>    For further details see the references listed
!>    below. For more details of the implementation see [3].
!>

References:

!>    [1] P. Schmid: Dynamic mode decomposition of numerical
!>        and experimental data,
!>        Journal of Fluid Mechanics 656, 5-28, 2010.
!>    [2] Z. Drmac, I. Mezic, R. Mohr: Data driven modal
!>        decompositions: analysis and enhancements,
!>        SIAM J. on Sci. Comp. 40 (4), A2253-A2285, 2018.
!>    [3] Z. Drmac: A LAPACK implementation of the Dynamic
!>        Mode Decomposition I. Technical report. AIMDyn Inc.
!>        and LAPACK Working Note 298.
!>    [4] J. Tu, C. W. Rowley, D. M. Luchtenburg, S. L.
!>        Brunton, N. Kutz: On Dynamic Mode Decomposition:
!>        Theory and Applications, Journal of Computational
!>        Dynamics 1(2), 391 -421, 2014.
!>

Developed and supported by:

!>    Developed and coded by Zlatko Drmac, Faculty of Science,
!>    University of Zagreb;  drmac@math.hr
!>    In cooperation with
!>    AIMdyn Inc., Santa Barbara, CA.
!>    and supported by
!>    - DARPA SBIR project  Contract No: W31P4Q-21-C-0007
!>    - DARPA PAI project  Contract No: HR0011-18-9-0033
!>    - DARPA MoDyL project 
!>    Contract No: HR0011-16-C-0116
!>    Any opinions, findings and conclusions or recommendations
!>    expressed in this material are those of the author and
!>    do not necessarily reflect the views of the DARPA SBIR
!>    Program Office
!>

Distribution Statement A:

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

Parameters

JOBS

!>    JOBS (input) is CHARACTER*1
!>    Determines whether the initial data snapshots are scaled
!>    by a diagonal matrix.
!>    'S' :: The data snapshots matrices X and Y are multiplied
!>           with a diagonal matrix D so that X*D has unit
!>           nonzero columns (in the Euclidean 2-norm)
!>    'C' :: The snapshots are scaled as with the 'S' option.
!>           If it is found that an i-th column of X is zero
!>           vector and the corresponding i-th column of Y is
!>           non-zero, then the i-th column of Y is set to
!>           zero and a warning flag is raised.
!>    'Y' :: The data snapshots matrices X and Y are multiplied
!>           by a diagonal matrix D so that Y*D has unit
!>           nonzero columns (in the Euclidean 2-norm)
!>    'N' :: No data scaling.
!>

JOBZ

!>    JOBZ (input) CHARACTER*1
!>    Determines whether the eigenvectors (Koopman modes) will
!>    be computed.
!>    'V' :: The eigenvectors (Koopman modes) will be computed
!>           and returned in the matrix Z.
!>           See the description of Z.
!>    'F' :: The eigenvectors (Koopman modes) will be returned
!>           in factored form as the product X(:,1:K)*W, where X
!>           contains a POD basis (leading left singular vectors
!>           of the data matrix X) and W contains the eigenvectors
!>           of the corresponding Rayleigh quotient.
!>           See the descriptions of K, X, W, Z.
!>    'N' :: The eigenvectors are not computed.
!>

JOBR

!>    JOBR (input) CHARACTER*1
!>    Determines whether to compute the residuals.
!>    'R' :: The residuals for the computed eigenpairs will be
!>           computed and stored in the array RES.
!>           See the description of RES.
!>           For this option to be legal, JOBZ must be 'V'.
!>    'N' :: The residuals are not computed.
!>

JOBF

!>    JOBF (input) CHARACTER*1
!>    Specifies whether to store information needed for post-
!>    processing (e.g. computing refined Ritz vectors)
!>    'R' :: The matrix needed for the refinement of the Ritz
!>           vectors is computed and stored in the array B.
!>           See the description of B.
!>    'E' :: The unscaled eigenvectors of the Exact DMD are
!>           computed and returned in the array B. See the
!>           description of B.
!>    'N' :: No eigenvector refinement data is computed.
!>

WHTSVD

!>    WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 }
!>    Allows for a selection of the SVD algorithm from the
!>    LAPACK library.
!>    1 :: DGESVD (the QR SVD algorithm)
!>    2 :: DGESDD (the Divide and Conquer algorithm; if enough
!>         workspace available, this is the fastest option)
!>    3 :: DGESVDQ (the preconditioned QR SVD  ; this and 4
!>         are the most accurate options)
!>    4 :: DGEJSV (the preconditioned Jacobi SVD; this and 3
!>         are the most accurate options)
!>    For the four methods above, a significant difference in
!>    the accuracy of small singular values is possible if
!>    the snapshots vary in norm so that X is severely
!>    ill-conditioned. If small (smaller than EPS*||X||)
!>    singular values are of interest and JOBS=='N',  then
!>    the options (3, 4) give the most accurate results, where
!>    the option 4 is slightly better and with stronger
!>    theoretical background.
!>    If JOBS=='S', i.e. the columns of X will be normalized,
!>    then all methods give nearly equally accurate results.
!>

!>    M (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

!>    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

!>    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 (output) INTEGER,  0 <= K <= N
!>    The dimension of the POD basis for the data snapshot
!>    matrix X and the number of the computed Ritz pairs.
!>    The value of K is determined according to the rule set
!>    by the parameters NRNK and TOL.
!>    See the descriptions of NRNK and TOL.
!>

REIG

!>    REIG (output) REAL(KIND=WP) N-by-1 array
!>    The leading K (K<=N) entries of REIG contain
!>    the real parts of the computed eigenvalues
!>    REIG(1:K) + sqrt(-1)*IMEIG(1:K).
!>    See the descriptions of K, IMEIG, and Z.
!>

IMEIG

!>    IMEIG (output) REAL(KIND=WP) N-by-1 array
!>    The leading K (K<=N) entries of IMEIG contain
!>    the imaginary parts of the computed eigenvalues
!>    REIG(1:K) + sqrt(-1)*IMEIG(1:K).
!>    The eigenvalues are determined as follows:
!>    If IMEIG(i) == 0, then the corresponding eigenvalue is
!>    real, LAMBDA(i) = REIG(i).
!>    If IMEIG(i)>0, then the corresponding complex
!>    conjugate pair of eigenvalues reads
!>    LAMBDA(i)   = REIG(i) + sqrt(-1)*IMAG(i)
!>    LAMBDA(i+1) = REIG(i) - sqrt(-1)*IMAG(i)
!>    That is, complex conjugate pairs have consecutive
!>    indices (i,i+1), with the positive imaginary part
!>    listed first.
!>    See the descriptions of K, REIG, and Z.
!>

!>    Z (workspace/output) REAL(KIND=WP)  M-by-N array
!>    If JOBZ =='V' then
!>       Z contains real Ritz vectors as follows:
!>       If IMEIG(i)=0, then Z(:,i) is an eigenvector of
!>       the i-th Ritz value; ||Z(:,i)||_2=1.
!>       If IMEIG(i) > 0 (and IMEIG(i+1) < 0) then
!>       [Z(:,i) Z(:,i+1)] span an invariant subspace and
!>       the Ritz values extracted from this subspace are
!>       REIG(i) + sqrt(-1)*IMEIG(i) and
!>       REIG(i) - sqrt(-1)*IMEIG(i).
!>       The corresponding eigenvectors are
!>       Z(:,i) + sqrt(-1)*Z(:,i+1) and
!>       Z(:,i) - sqrt(-1)*Z(:,i+1), respectively.
!>       || Z(:,i:i+1)||_F = 1.
!>    If JOBZ == 'F', then the above descriptions hold for
!>    the columns of X(:,1:K)*W(1:K,1:K), where the columns
!>    of W(1:k,1:K) are the computed eigenvectors of the
!>    K-by-K Rayleigh quotient. The columns of W(1:K,1:K)
!>    are similarly structured: If IMEIG(i) == 0 then
!>    X(:,1:K)*W(:,i) is an eigenvector, and if IMEIG(i)>0
!>    then X(:,1:K)*W(:,i)+sqrt(-1)*X(:,1:K)*W(:,i+1) and
!>         X(:,1:K)*W(:,i)-sqrt(-1)*X(:,1:K)*W(:,i+1)
!>    are the eigenvectors of LAMBDA(i), LAMBDA(i+1).
!>    See the descriptions of REIG, IMEIG, X and W.
!>

LDZ

!>    LDZ (input) INTEGER , LDZ >= M
!>    The leading dimension of the array Z.
!>

RES

!>    RES (output) REAL(KIND=WP) N-by-1 array
!>    RES(1:K) contains the residuals for the K computed
!>    Ritz pairs.
!>    If LAMBDA(i) is real, then
!>       RES(i) = || A * Z(:,i) - LAMBDA(i)*Z(:,i))||_2.
!>    If [LAMBDA(i), LAMBDA(i+1)] is a complex conjugate pair
!>    then
!>    RES(i)=RES(i+1) = || A * Z(:,i:i+1) - Z(:,i:i+1) *B||_F
!>    where B = [ real(LAMBDA(i)) imag(LAMBDA(i)) ]
!>              [-imag(LAMBDA(i)) real(LAMBDA(i)) ].
!>    It holds that
!>    RES(i)   = || A*ZC(:,i)   - LAMBDA(i)  *ZC(:,i)   ||_2
!>    RES(i+1) = || A*ZC(:,i+1) - LAMBDA(i+1)*ZC(:,i+1) ||_2
!>    where ZC(:,i)   =  Z(:,i) + sqrt(-1)*Z(:,i+1)
!>          ZC(:,i+1) =  Z(:,i) - sqrt(-1)*Z(:,i+1)
!>    See the description of REIG, IMEIG and Z.
!>

!>    B (output) REAL(KIND=WP)  M-by-N array.
!>    IF JOBF =='R', B(1:M,1:K) contains A*U(:,1:K), and can
!>    be used for computing the refined vectors; see further
!>    details in the provided references.
!>    If JOBF == 'E', B(1:M,1;K) contains
!>    A*U(:,1:K)*W(1:K,1:K), which are the vectors from the
!>    Exact DMD, up to scaling by the inverse eigenvalues.
!>    If JOBF =='N', then B is not referenced.
!>    See the descriptions of X, W, K.
!>

LDB

!>    LDB (input) INTEGER, LDB >= M
!>    The leading dimension of the array B.
!>

!>    W (workspace/output) REAL(KIND=WP) N-by-N array
!>    On exit, W(1:K,1:K) contains the K computed
!>    eigenvectors of the matrix Rayleigh quotient (real and
!>    imaginary parts for each complex conjugate pair of the
!>    eigenvalues). The Ritz vectors (returned in Z) are the
!>    product of X (containing a POD basis for the input
!>    matrix X) and W. See the descriptions of K, S, X and Z.
!>    W is also used as a workspace to temporarily store the
!>    right singular vectors of X.
!>

LDW

!>    LDW (input) INTEGER, LDW >= N
!>    The leading dimension of the array W.
!>

!>    S (workspace/output) REAL(KIND=WP) N-by-N array
!>    The array S(1:K,1:K) is used for the matrix Rayleigh
!>    quotient. This content is overwritten during
!>    the eigenvalue decomposition by DGEEV.
!>    See the description of K.
!>

LDS

!>    LDS (input) INTEGER, LDS >= N
!>    The leading dimension of the array S.
!>

WORK

!>    WORK (workspace/output) REAL(KIND=WP) LWORK-by-1 array
!>    On exit, WORK(1:N) contains the singular values of
!>    X (for JOBS=='N') or column scaled X (JOBS=='S', 'C').
!>    If WHTSVD==4, then WORK(N+1) and WORK(N+2) contain
!>    scaling factor WORK(N+2)/WORK(N+1) used to scale X
!>    and Y to avoid overflow in the SVD of X.
!>    This may be of interest if the scaling option is off
!>    and as many as possible smallest eigenvalues are
!>    desired to the highest feasible accuracy.
!>    If the call to DGEDMD is only workspace query, then
!>    WORK(1) contains the minimal workspace length and
!>    WORK(2) is the optimal workspace length. Hence, the
!>    leng of work is at least 2.
!>    See the description of LWORK.
!>

LWORK

!>    LWORK (input) INTEGER
!>    The minimal length of the workspace vector WORK.
!>    LWORK is calculated as follows:
!>    If WHTSVD == 1 ::
!>       If JOBZ == 'V', then
!>       LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,4*N)).
!>       If JOBZ == 'N'  then
!>       LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,3*N)).
!>       Here LWORK_SVD = MAX(1,3*N+M,5*N) is the minimal
!>       workspace length of DGESVD.
!>    If WHTSVD == 2 ::
!>       If JOBZ == 'V', then
!>       LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,4*N))
!>       If JOBZ == 'N', then
!>       LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,3*N))
!>       Here LWORK_SVD = MAX(M, 5*N*N+4*N)+3*N*N is the
!>       minimal workspace length of DGESDD.
!>    If WHTSVD == 3 ::
!>       If JOBZ == 'V', then
!>       LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,4*N))
!>       If JOBZ == 'N', then
!>       LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,3*N))
!>       Here LWORK_SVD = N+M+MAX(3*N+1,
!>                       MAX(1,3*N+M,5*N),MAX(1,N))
!>       is the minimal workspace length of DGESVDQ.
!>    If WHTSVD == 4 ::
!>       If JOBZ == 'V', then
!>       LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,4*N))
!>       If JOBZ == 'N', then
!>       LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,3*N))
!>       Here LWORK_SVD = MAX(7,2*M+N,6*N+2*N*N) is the
!>       minimal workspace length of DGEJSV.
!>    The above expressions are not simplified in order to
!>    make the usage of WORK more transparent, and for
!>    easier checking. In any case, LWORK >= 2.
!>    If on entry LWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace lengths for both WORK and
!>    IWORK. See the descriptions of WORK and IWORK.
!>

IWORK

!>    IWORK (workspace/output) INTEGER LIWORK-by-1 array
!>    Workspace that is required only if WHTSVD equals
!>    2 , 3 or 4. (See the description of WHTSVD).
!>    If on entry LWORK =-1 or LIWORK=-1, then the
!>    minimal length of IWORK is computed and returned in
!>    IWORK(1). See the description of LIWORK.
!>

LIWORK

!>    LIWORK (input) INTEGER
!>    The minimal length of the workspace vector IWORK.
!>    If WHTSVD == 1, then only IWORK(1) is used; LIWORK >=1
!>    If WHTSVD == 2, then LIWORK >= MAX(1,8*MIN(M,N))
!>    If WHTSVD == 3, then LIWORK >= MAX(1,M+N-1)
!>    If WHTSVD == 4, then LIWORK >= MAX(3,M+3*N)
!>    If on entry LIWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace lengths for both WORK and
!>    IWORK. See the descriptions of WORK and IWORK.
!>

INFO

!>    INFO (output) INTEGER
!>    -i < 0 :: On entry, the i-th argument had an
!>              illegal value
!>       = 0 :: Successful return.
!>       = 1 :: Void input. Quick exit (M=0 or N=0).
!>       = 2 :: The SVD computation of X did not converge.
!>              Suggestion: Check the input data and/or
!>              repeat with different WHTSVD.
!>       = 3 :: The computation of the eigenvalues did not
!>              converge.
!>       = 4 :: If data scaling was requested on input and
!>              the procedure found inconsistency in the data
!>              such that for some column index i,
!>              X(:,i) = 0 but Y(:,i) /= 0, then Y(:,i) is set
!>              to zero if JOBS=='C'. The computation proceeds
!>              with original or modified data and warning
!>              flag is set with INFO=4.
!>

Author

Zlatko Drmac

Definition at line 531 of file dgedmd.f90.

subroutine DGEDMDQ (character, intent(in) jobs, character, intent(in) jobz, character, intent(in) jobr, character, intent(in) jobq, character, intent(in) jobt, character, intent(in) jobf, integer, intent(in) whtsvd, integer, intent(in) m, integer, intent(in) n, real(kind=wp), dimension(ldf,), intent(inout) f, integer, intent(in) ldf, real(kind=wp), dimension(ldx,), intent(out) x, integer, intent(in) ldx, real(kind=wp), dimension(ldy,), intent(out) y, integer, intent(in) ldy, integer, intent(in) nrnk, real(kind=wp), intent(in) tol, integer, intent(out) k, real(kind=wp), dimension(), intent(out) reig, real(kind=wp), dimension(), intent(out) imeig, real(kind=wp), dimension(ldz,), intent(out) z, integer, intent(in) ldz, real(kind=wp), dimension(), intent(out) res, real(kind=wp), dimension(ldb,), intent(out) b, integer, intent(in) ldb, real(kind=wp), dimension(ldv,), intent(out) v, integer, intent(in) ldv, real(kind=wp), dimension(lds,), intent(out) s, integer, intent(in) lds, real(kind=wp), dimension(), intent(out) work, integer, intent(in) lwork, integer, dimension(), intent(out) iwork, integer, intent(in) liwork, integer, intent(out) info)¶

DGEDMDQ computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.

Purpose:

!>     DGEDMDQ computes the Dynamic Mode Decomposition (DMD) for
!>     a pair of data snapshot matrices, using a QR factorization
!>     based compression of the data. For the input matrices
!>     X and Y such that Y = A*X with an unaccessible matrix
!>     A, DGEDMDQ computes a certain number of Ritz pairs of A using
!>     the standard Rayleigh-Ritz extraction from a subspace of
!>     range(X) that is determined using the leading left singular 
!>     vectors of X. Optionally, DGEDMDQ returns the residuals 
!>     of the computed Ritz pairs, the information needed for
!>     a refinement of the Ritz vectors, or the eigenvectors of
!>     the Exact DMD.
!>     For further details see the references listed
!>     below. For more details of the implementation see [3].      
!>

References:

!>    [1] P. Schmid: Dynamic mode decomposition of numerical
!>        and experimental data,
!>        Journal of Fluid Mechanics 656, 5-28, 2010.
!>    [2] Z. Drmac, I. Mezic, R. Mohr: Data driven modal
!>        decompositions: analysis and enhancements,
!>        SIAM J. on Sci. Comp. 40 (4), A2253-A2285, 2018.
!>    [3] Z. Drmac: A LAPACK implementation of the Dynamic
!>        Mode Decomposition I. Technical report. AIMDyn Inc.
!>        and LAPACK Working Note 298.      
!>    [4] J. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. 
!>        Brunton, N. Kutz: On Dynamic Mode Decomposition:
!>        Theory and Applications, Journal of Computational
!>        Dynamics 1(2), 391 -421, 2014.
!>

Developed and supported by:

!>    Developed and coded by Zlatko Drmac, Faculty of Science,
!>    University of Zagreb;  drmac@math.hr
!>    In cooperation with
!>    AIMdyn Inc., Santa Barbara, CA.
!>    and supported by
!>    - DARPA SBIR project  Contract No: W31P4Q-21-C-0007
!>    - DARPA PAI project  Contract No: HR0011-18-9-0033
!>    - DARPA MoDyL project 
!>    Contract No: HR0011-16-C-0116
!>    Any opinions, findings and conclusions or recommendations 
!>    expressed in this material are those of the author and 
!>    do not necessarily reflect the views of the DARPA SBIR 
!>    Program Office.      
!>

Distribution Statement A:

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

Parameters

JOBS

!>    JOBS (input) CHARACTER*1
!>    Determines whether the initial data snapshots are scaled
!>    by a diagonal matrix. The data snapshots are the columns
!>    of F. The leading N-1 columns of F are denoted X and the
!>    trailing N-1 columns are denoted Y. 
!>    'S' :: The data snapshots matrices X and Y are multiplied
!>           with a diagonal matrix D so that X*D has unit
!>           nonzero columns (in the Euclidean 2-norm)
!>    'C' :: The snapshots are scaled as with the 'S' option.
!>           If it is found that an i-th column of X is zero
!>           vector and the corresponding i-th column of Y is
!>           non-zero, then the i-th column of Y is set to
!>           zero and a warning flag is raised.
!>    'Y' :: The data snapshots matrices X and Y are multiplied
!>           by a diagonal matrix D so that Y*D has unit
!>           nonzero columns (in the Euclidean 2-norm)    
!>    'N' :: No data scaling.   
!>

JOBZ

!>    JOBZ (input) CHARACTER*1
!>    Determines whether the eigenvectors (Koopman modes) will
!>    be computed.
!>    'V' :: The eigenvectors (Koopman modes) will be computed
!>           and returned in the matrix Z.
!>           See the description of Z.
!>    'F' :: The eigenvectors (Koopman modes) will be returned
!>           in factored form as the product Z*V, where Z
!>           is orthonormal and V contains the eigenvectors
!>           of the corresponding Rayleigh quotient.
!>           See the descriptions of F, V, Z.
!>    'Q' :: The eigenvectors (Koopman modes) will be returned
!>           in factored form as the product Q*Z, where Z
!>           contains the eigenvectors of the compression of the
!>           underlying discretized operator onto the span of
!>           the data snapshots. See the descriptions of F, V, Z.  
!>           Q is from the initial QR factorization.      
!>    'N' :: The eigenvectors are not computed.  
!>

JOBR

!>    JOBR (input) CHARACTER*1 
!>    Determines whether to compute the residuals.
!>    'R' :: The residuals for the computed eigenpairs will
!>           be computed and stored in the array RES.
!>           See the description of RES.
!>           For this option to be legal, JOBZ must be 'V'.
!>    'N' :: The residuals are not computed.
!>

JOBQ

!>    JOBQ (input) CHARACTER*1 
!>    Specifies whether to explicitly compute and return the
!>    orthogonal matrix from the QR factorization.
!>    'Q' :: The matrix Q of the QR factorization of the data
!>           snapshot matrix is computed and stored in the
!>           array F. See the description of F.       
!>    'N' :: The matrix Q is not explicitly computed.
!>

JOBT

!>    JOBT (input) CHARACTER*1 
!>    Specifies whether to return the upper triangular factor
!>    from the QR factorization.
!>    'R' :: The matrix R of the QR factorization of the data 
!>           snapshot matrix F is returned in the array Y.
!>           See the description of Y and Further details.       
!>    'N' :: The matrix R is not returned.    
!>

JOBF

!>    JOBF (input) CHARACTER*1
!>    Specifies whether to store information needed for post-
!>    processing (e.g. computing refined Ritz vectors)
!>    'R' :: The matrix needed for the refinement of the Ritz
!>           vectors is computed and stored in the array B.
!>           See the description of B.
!>    'E' :: The unscaled eigenvectors of the Exact DMD are 
!>           computed and returned in the array B. See the
!>           description of B.
!>    'N' :: No eigenvector refinement data is computed.
!>    To be useful on exit, this option needs JOBQ='Q'.
!>

WHTSVD

!>    WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 }
!>    Allows for a selection of the SVD algorithm from the
!>    LAPACK library.
!>    1 :: DGESVD (the QR SVD algorithm)
!>    2 :: DGESDD (the Divide and Conquer algorithm; if enough
!>         workspace available, this is the fastest option)
!>    3 :: DGESVDQ (the preconditioned QR SVD  ; this and 4
!>         are the most accurate options)
!>    4 :: DGEJSV (the preconditioned Jacobi SVD; this and 3
!>         are the most accurate options)
!>    For the four methods above, a significant difference in
!>    the accuracy of small singular values is possible if
!>    the snapshots vary in norm so that X is severely
!>    ill-conditioned. If small (smaller than EPS*||X||)
!>    singular values are of interest and JOBS=='N',  then
!>    the options (3, 4) give the most accurate results, where
!>    the option 4 is slightly better and with stronger 
!>    theoretical background.
!>    If JOBS=='S', i.e. the columns of X will be normalized,
!>    then all methods give nearly equally accurate results.
!>

!>    M (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 (input/output) REAL(KIND=WP) M-by-N array
!>    > On entry,
!>    the columns of F are the sequence of data snapshots 
!>    from a single trajectory, taken at equidistant discrete
!>    times. It is assumed that the column norms of F are 
!>    in the range of the normalized floating point numbers. 
!>    < On exit,
!>    If JOBQ == 'Q', the array F contains the orthogonal 
!>    matrix/factor of the QR factorization of the initial 
!>    data snapshots matrix F. See the description of JOBQ. 
!>    If JOBQ == 'N', the entries in F strictly below the main
!>    diagonal contain, column-wise, the information on the 
!>    Householder vectors, as returned by DGEQRF. The 
!>    remaining information to restore the orthogonal matrix
!>    of the initial QR factorization is stored in WORK(1:N). 
!>    See the description of WORK.
!>

LDF

!>    LDF (input) INTEGER, LDF >= M 
!>    The leading dimension of the array F.
!>

!>    X (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 (workspace/output) REAL(KIND=WP) MIN(M,N)-by-(N-1) array
!>    Y is used as workspace to hold representations of the
!>    trailing N-1 snapshots in the orthonormal basis computed
!>    in the QR factorization of F.
!>    On exit, 
!>    If JOBT == 'R', Y contains the MIN(M,N)-by-N upper
!>    triangular factor from the QR factorization of the data
!>    snapshot matrix F.
!>

LDY

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

NRNK

!>    NRNK (input) INTEGER
!>    Determines the mode how to compute the numerical rank,
!>    i.e. how to truncate small singular values of the input
!>    matrix X. On input, if
!>    NRNK = -1 :: i-th singular value sigma(i) is truncated
!>                 if sigma(i) <= TOL*sigma(1)
!>                 This option is recommended.  
!>    NRNK = -2 :: i-th singular value sigma(i) is truncated
!>                 if sigma(i) <= TOL*sigma(i-1)
!>                 This option is included for R&D purposes.
!>                 It requires highly accurate SVD, which
!>                 may not be feasible.       
!>    The numerical rank can be enforced by using positive 
!>    value of NRNK as follows: 
!>    0 < NRNK <= N-1 :: at most NRNK largest singular values
!>    will be used. If the number of the computed nonzero
!>    singular values is less than NRNK, then only those
!>    nonzero values will be used and the actually used
!>    dimension is less than NRNK. The actual number of
!>    the nonzero singular values is returned in the variable
!>    K. See the description of K.
!>

TOL

!>     TOL (input) REAL(KIND=WP), 0 <= TOL < 1
!>     The tolerance for truncating small singular values.
!>     See the description of NRNK.  
!>

!>     K (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:
!>    If IMEIG(i) == 0, then the corresponding eigenvalue is
!>    real, LAMBDA(i) = REIG(i).
!>    If IMEIG(i)>0, then the corresponding complex
!>    conjugate pair of eigenvalues reads
!>    LAMBDA(i)   = REIG(i) + sqrt(-1)*IMAG(i)
!>    LAMBDA(i+1) = REIG(i) - sqrt(-1)*IMAG(i)
!>    That is, complex conjugate pairs have consequtive
!>    indices (i,i+1), with the positive imaginary part
!>    listed first.
!>    See the descriptions of K, REIG, Z.     
!>

!>    Z (workspace/output) REAL(KIND=WP)  M-by-(N-1) array
!>    If JOBZ =='V' then
!>       Z contains real Ritz vectors as follows:
!>       If IMEIG(i)=0, then Z(:,i) is an eigenvector of
!>       the i-th Ritz value.
!>       If IMEIG(i) > 0 (and IMEIG(i+1) < 0) then
!>       [Z(:,i) Z(:,i+1)] span an invariant subspace and
!>       the Ritz values extracted from this subspace are
!>       REIG(i) + sqrt(-1)*IMEIG(i) and
!>       REIG(i) - sqrt(-1)*IMEIG(i).
!>       The corresponding eigenvectors are
!>       Z(:,i) + sqrt(-1)*Z(:,i+1) and
!>       Z(:,i) - sqrt(-1)*Z(:,i+1), respectively.
!>    If JOBZ == 'F', then the above descriptions hold for
!>    the columns of Z*V, where the columns of V are the
!>    eigenvectors of the K-by-K Rayleigh quotient, and Z is
!>    orthonormal. The columns of V are similarly structured:
!>    If IMEIG(i) == 0 then Z*V(:,i) is an eigenvector, and if 
!>    IMEIG(i) > 0 then Z*V(:,i)+sqrt(-1)*Z*V(:,i+1) and
!>                      Z*V(:,i)-sqrt(-1)*Z*V(:,i+1)
!>    are the eigenvectors of LAMBDA(i), LAMBDA(i+1).
!>    See the descriptions of REIG, IMEIG, X and V.
!>

LDZ

!>    LDZ (input) INTEGER , LDZ >= M
!>    The leading dimension of the array Z.
!>

RES

!>    RES (output) REAL(KIND=WP) (N-1)-by-1 array
!>    RES(1:K) contains the residuals for the K computed 
!>    Ritz pairs.       
!>    If LAMBDA(i) is real, then
!>       RES(i) = || A * Z(:,i) - LAMBDA(i)*Z(:,i))||_2.
!>    If [LAMBDA(i), LAMBDA(i+1)] is a complex conjugate pair
!>    then
!>    RES(i)=RES(i+1) = || A * Z(:,i:i+1) - Z(:,i:i+1) *B||_F
!>    where B = [ real(LAMBDA(i)) imag(LAMBDA(i)) ]
!>              [-imag(LAMBDA(i)) real(LAMBDA(i)) ].
!>    It holds that
!>    RES(i)   = || A*ZC(:,i)   - LAMBDA(i)  *ZC(:,i)   ||_2
!>    RES(i+1) = || A*ZC(:,i+1) - LAMBDA(i+1)*ZC(:,i+1) ||_2
!>    where ZC(:,i)   =  Z(:,i) + sqrt(-1)*Z(:,i+1)
!>          ZC(:,i+1) =  Z(:,i) - sqrt(-1)*Z(:,i+1)
!>    See the description of Z.
!>

!>    B (output) REAL(KIND=WP)  MIN(M,N)-by-(N-1) array.
!>    IF JOBF =='R', B(1:N,1:K) contains A*U(:,1:K), and can
!>    be used for computing the refined vectors; see further 
!>    details in the provided references. 
!>    If JOBF == 'E', B(1:N,1;K) contains 
!>    A*U(:,1:K)*W(1:K,1:K), which are the vectors from the
!>    Exact DMD, up to scaling by the inverse eigenvalues.   
!>    In both cases, the content of B can be lifted to the 
!>    original dimension of the input data by pre-multiplying
!>    with the Q factor from the initial QR factorization.
!>    Here A denotes a compression of the underlying operator.
!>    See the descriptions of F and X.
!>    If JOBF =='N', then B is not referenced.
!>

LDB

!>    LDB (input) INTEGER, LDB >= MIN(M,N)
!>    The leading dimension of the array B.
!>

!>    V (workspace/output) REAL(KIND=WP) (N-1)-by-(N-1) array
!>    On exit, V(1:K,1:K) contains the K eigenvectors of
!>    the Rayleigh quotient. The eigenvectors of a complex
!>    conjugate pair of eigenvalues are returned in real form
!>    as explained in the description of Z. The Ritz vectors
!>    (returned in Z) are the product of X and V; see
!>    the descriptions of X and Z.
!>

LDV

!>    LDV (input) INTEGER, LDV >= N-1
!>    The leading dimension of the array V.
!>

!>    S (output) REAL(KIND=WP) (N-1)-by-(N-1) array
!>    The array S(1:K,1:K) is used for the matrix Rayleigh
!>    quotient. This content is overwritten during
!>    the eigenvalue decomposition by DGEEV.
!>    See the description of K.
!>

LDS

!>    LDS (input) INTEGER, LDS >= N-1        
!>    The leading dimension of the array S.
!>

WORK

!>    WORK (workspace/output) REAL(KIND=WP) LWORK-by-1 array
!>    On exit, 
!>    WORK(1:MIN(M,N)) contains the scalar factors of the 
!>    elementary reflectors as returned by DGEQRF of the 
!>    M-by-N input matrix F.
!>    WORK(MIN(M,N)+1:MIN(M,N)+N-1) contains the singular values of 
!>    the input submatrix F(1:M,1:N-1).
!>    If the call to DGEDMDQ is only workspace query, then
!>    WORK(1) contains the minimal workspace length and
!>    WORK(2) is the optimal workspace length. Hence, the
!>    length of work is at least 2.
!>    See the description of LWORK.
!>

LWORK

!>    LWORK (input) INTEGER
!>    The minimal length of the  workspace vector WORK.
!>    LWORK is calculated as follows:
!>    Let MLWQR  = N (minimal workspace for DGEQRF[M,N])
!>        MLWDMD = minimal workspace for DGEDMD (see the
!>                 description of LWORK in DGEDMD) for 
!>                 snapshots of dimensions MIN(M,N)-by-(N-1)
!>        MLWMQR = N (minimal workspace for 
!>                   DORMQR['L','N',M,N,N])
!>        MLWGQR = N (minimal workspace for DORGQR[M,N,N])
!>    Then
!>    LWORK = MAX(N+MLWQR, N+MLWDMD)
!>    is updated as follows:
!>       if   JOBZ == 'V' or JOBZ == 'F' THEN 
!>            LWORK = MAX( LWORK, MIN(M,N)+N-1+MLWMQR )
!>       if   JOBQ == 'Q' THEN
!>            LWORK = MAX( LWORK, MIN(M,N)+N-1+MLWGQR)
!>    If on entry LWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace lengths for both WORK and
!>    IWORK. See the descriptions of WORK and IWORK.          
!>

IWORK

!>    IWORK (workspace/output) INTEGER LIWORK-by-1 array
!>    Workspace that is required only if WHTSVD equals
!>    2 , 3 or 4. (See the description of WHTSVD).
!>    If on entry LWORK =-1 or LIWORK=-1, then the
!>    minimal length of IWORK is computed and returned in
!>    IWORK(1). See the description of LIWORK.
!>

LIWORK

!>    LIWORK (input) INTEGER
!>    The minimal length of the workspace vector IWORK.
!>    If WHTSVD == 1, then only IWORK(1) is used; LIWORK >=1
!>    Let M1=MIN(M,N), N1=N-1. Then    
!>    If WHTSVD == 2, then LIWORK >= MAX(1,8*MIN(M1,N1))
!>    If WHTSVD == 3, then LIWORK >= MAX(1,M1+N1-1)
!>    If WHTSVD == 4, then LIWORK >= MAX(3,M1+3*N1)
!>    If on entry LIWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace lengths for both WORK and
!>    IWORK. See the descriptions of WORK and IWORK.
!>

INFO

!>    INFO (output) INTEGER
!>    -i < 0 :: On entry, the i-th argument had an
!>              illegal value
!>       = 0 :: Successful return.
!>       = 1 :: Void input. Quick exit (M=0 or N=0).
!>       = 2 :: The SVD computation of X did not converge.
!>              Suggestion: Check the input data and/or
!>              repeat with different WHTSVD.
!>       = 3 :: The computation of the eigenvalues did not
!>              converge.
!>       = 4 :: If data scaling was requested on input and
!>              the procedure found inconsistency in the data
!>              such that for some column index i,
!>              X(:,i) = 0 but Y(:,i) /= 0, then Y(:,i) is set
!>              to zero if JOBS=='C'. The computation proceeds
!>              with original or modified data and warning
!>              flag is set with INFO=4.  
!>

Author

Zlatko Drmac

Definition at line 571 of file dgedmdq.f90.

subroutine SGEDMD (character, intent(in) jobs, character, intent(in) jobz, character, intent(in) jobr, character, intent(in) jobf, integer, intent(in) whtsvd, integer, intent(in) m, integer, intent(in) n, real(kind=wp), dimension(ldx,), intent(inout) x, integer, intent(in) ldx, real(kind=wp), dimension(ldy,), intent(inout) y, integer, intent(in) ldy, integer, intent(in) nrnk, real(kind=wp), intent(in) tol, integer, intent(out) k, real(kind=wp), dimension(), intent(out) reig, real(kind=wp), dimension(), intent(out) imeig, real(kind=wp), dimension(ldz,), intent(out) z, integer, intent(in) ldz, real(kind=wp), dimension(), intent(out) res, real(kind=wp), dimension(ldb,), intent(out) b, integer, intent(in) ldb, real(kind=wp), dimension(ldw,), intent(out) w, integer, intent(in) ldw, real(kind=wp), dimension(lds,), intent(out) s, integer, intent(in) lds, real(kind=wp), dimension(), intent(out) work, integer, intent(in) lwork, integer, dimension(*), intent(out) iwork, integer, intent(in) liwork, integer, intent(out) info)¶

SGEDMD computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.

Purpose:

!>    SGEDMD computes the Dynamic Mode Decomposition (DMD) for
!>    a pair of data snapshot matrices. For the input matrices
!>    X and Y such that Y = A*X with an unaccessible matrix
!>    A, SGEDMD computes a certain number of Ritz pairs of A using
!>    the standard Rayleigh-Ritz extraction from a subspace of
!>    range(X) that is determined using the leading left singular
!>    vectors of X. Optionally, SGEDMD returns the residuals
!>    of the computed Ritz pairs, the information needed for
!>    a refinement of the Ritz vectors, or the eigenvectors of
!>    the Exact DMD.
!>    For further details see the references listed
!>    below. For more details of the implementation see [3].
!>

References:

!>    [1] P. Schmid: Dynamic mode decomposition of numerical
!>        and experimental data,
!>        Journal of Fluid Mechanics 656, 5-28, 2010.
!>    [2] Z. Drmac, I. Mezic, R. Mohr: Data driven modal
!>        decompositions: analysis and enhancements,
!>        SIAM J. on Sci. Comp. 40 (4), A2253-A2285, 2018.
!>    [3] Z. Drmac: A LAPACK implementation of the Dynamic
!>        Mode Decomposition I. Technical report. AIMDyn Inc.
!>        and LAPACK Working Note 298.
!>    [4] J. Tu, C. W. Rowley, D. M. Luchtenburg, S. L.
!>        Brunton, N. Kutz: On Dynamic Mode Decomposition:
!>        Theory and Applications, Journal of Computational
!>        Dynamics 1(2), 391 -421, 2014.
!>

Developed and supported by:

!>    Developed and coded by Zlatko Drmac, Faculty of Science,
!>    University of Zagreb;  drmac@math.hr
!>    In cooperation with
!>    AIMdyn Inc., Santa Barbara, CA.
!>    and supported by
!>    - DARPA SBIR project  Contract No: W31P4Q-21-C-0007
!>    - DARPA PAI project  Contract No: HR0011-18-9-0033
!>    - DARPA MoDyL project 
!>    Contract No: HR0011-16-C-0116
!>    Any opinions, findings and conclusions or recommendations
!>    expressed in this material are those of the author and
!>    do not necessarily reflect the views of the DARPA SBIR
!>    Program Office
!>

Distribution Statement A:

!>    Distribution Statement A:
!>    Approved for Public Release, Distribution Unlimited.
!>    Cleared by DARPA on September 29, 2022
!>

Parameters

JOBS

!>    JOBS (input) CHARACTER*1
!>    Determines whether the initial data snapshots are scaled
!>    by a diagonal matrix.
!>    'S' :: The data snapshots matrices X and Y are multiplied
!>           with a diagonal matrix D so that X*D has unit
!>           nonzero columns (in the Euclidean 2-norm)
!>    'C' :: The snapshots are scaled as with the 'S' option.
!>           If it is found that an i-th column of X is zero
!>           vector and the corresponding i-th column of Y is
!>           non-zero, then the i-th column of Y is set to
!>           zero and a warning flag is raised.
!>    'Y' :: The data snapshots matrices X and Y are multiplied
!>           by a diagonal matrix D so that Y*D has unit
!>           nonzero columns (in the Euclidean 2-norm)
!>    'N' :: No data scaling.
!>

JOBZ

!>    JOBZ (input) CHARACTER*1
!>    Determines whether the eigenvectors (Koopman modes) will
!>    be computed.
!>    'V' :: The eigenvectors (Koopman modes) will be computed
!>           and returned in the matrix Z.
!>           See the description of Z.
!>    'F' :: The eigenvectors (Koopman modes) will be returned
!>           in factored form as the product X(:,1:K)*W, where X
!>           contains a POD basis (leading left singular vectors
!>           of the data matrix X) and W contains the eigenvectors
!>           of the corresponding Rayleigh quotient.
!>           See the descriptions of K, X, W, Z.
!>    'N' :: The eigenvectors are not computed.
!>

JOBR

!>    JOBR (input) CHARACTER*1
!>    Determines whether to compute the residuals.
!>    'R' :: The residuals for the computed eigenpairs will be
!>           computed and stored in the array RES.
!>           See the description of RES.
!>           For this option to be legal, JOBZ must be 'V'.
!>    'N' :: The residuals are not computed.
!>

JOBF

!>    JOBF (input) CHARACTER*1
!>    Specifies whether to store information needed for post-
!>    processing (e.g. computing refined Ritz vectors)
!>    'R' :: The matrix needed for the refinement of the Ritz
!>           vectors is computed and stored in the array B.
!>           See the description of B.
!>    'E' :: The unscaled eigenvectors of the Exact DMD are
!>           computed and returned in the array B. See the
!>           description of B.
!>    'N' :: No eigenvector refinement data is computed.
!>

WHTSVD

!>    WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 }
!>    Allows for a selection of the SVD algorithm from the
!>    LAPACK library.
!>    1 :: SGESVD (the QR SVD algorithm)
!>    2 :: SGESDD (the Divide and Conquer algorithm; if enough
!>         workspace available, this is the fastest option)
!>    3 :: SGESVDQ (the preconditioned QR SVD  ; this and 4
!>         are the most accurate options)
!>    4 :: SGEJSV (the preconditioned Jacobi SVD; this and 3
!>         are the most accurate options)
!>    For the four methods above, a significant difference in
!>    the accuracy of small singular values is possible if
!>    the snapshots vary in norm so that X is severely
!>    ill-conditioned. If small (smaller than EPS*||X||)
!>    singular values are of interest and JOBS=='N',  then
!>    the options (3, 4) give the most accurate results, where
!>    the option 4 is slightly better and with stronger
!>    theoretical background.
!>    If JOBS=='S', i.e. the columns of X will be normalized,
!>    then all methods give nearly equally accurate results.
!>

!>    M (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

!>    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

!>    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 (output) INTEGER,  0 <= K <= N
!>    The dimension of the POD basis for the data snapshot
!>    matrix X and the number of the computed Ritz pairs.
!>    The value of K is determined according to the rule set
!>    by the parameters NRNK and TOL.
!>    See the descriptions of NRNK and TOL.
!>

REIG

!>    REIG (output) REAL(KIND=WP) N-by-1 array
!>    The leading K (K<=N) entries of REIG contain
!>    the real parts of the computed eigenvalues
!>    REIG(1:K) + sqrt(-1)*IMEIG(1:K).
!>    See the descriptions of K, IMEIG, and Z.
!>

IMEIG

!>    IMEIG (output) REAL(KIND=WP) N-by-1 array
!>    The leading K (K<=N) entries of IMEIG contain
!>    the imaginary parts of the computed eigenvalues
!>    REIG(1:K) + sqrt(-1)*IMEIG(1:K).
!>    The eigenvalues are determined as follows:
!>    If IMEIG(i) == 0, then the corresponding eigenvalue is
!>    real, LAMBDA(i) = REIG(i).
!>    If IMEIG(i)>0, then the corresponding complex
!>    conjugate pair of eigenvalues reads
!>    LAMBDA(i)   = REIG(i) + sqrt(-1)*IMAG(i)
!>    LAMBDA(i+1) = REIG(i) - sqrt(-1)*IMAG(i)
!>    That is, complex conjugate pairs have consecutive
!>    indices (i,i+1), with the positive imaginary part
!>    listed first.
!>    See the descriptions of K, REIG, and Z.
!>

!>    Z (workspace/output) REAL(KIND=WP)  M-by-N array
!>    If JOBZ =='V' then
!>       Z contains real Ritz vectors as follows:
!>       If IMEIG(i)=0, then Z(:,i) is an eigenvector of
!>       the i-th Ritz value; ||Z(:,i)||_2=1.
!>       If IMEIG(i) > 0 (and IMEIG(i+1) < 0) then
!>       [Z(:,i) Z(:,i+1)] span an invariant subspace and
!>       the Ritz values extracted from this subspace are
!>       REIG(i) + sqrt(-1)*IMEIG(i) and
!>       REIG(i) - sqrt(-1)*IMEIG(i).
!>       The corresponding eigenvectors are
!>       Z(:,i) + sqrt(-1)*Z(:,i+1) and
!>       Z(:,i) - sqrt(-1)*Z(:,i+1), respectively.
!>       || Z(:,i:i+1)||_F = 1.
!>    If JOBZ == 'F', then the above descriptions hold for
!>    the columns of X(:,1:K)*W(1:K,1:K), where the columns
!>    of W(1:k,1:K) are the computed eigenvectors of the
!>    K-by-K Rayleigh quotient. The columns of W(1:K,1:K)
!>    are similarly structured: If IMEIG(i) == 0 then
!>    X(:,1:K)*W(:,i) is an eigenvector, and if IMEIG(i)>0
!>    then X(:,1:K)*W(:,i)+sqrt(-1)*X(:,1:K)*W(:,i+1) and
!>         X(:,1:K)*W(:,i)-sqrt(-1)*X(:,1:K)*W(:,i+1)
!>    are the eigenvectors of LAMBDA(i), LAMBDA(i+1).
!>    See the descriptions of REIG, IMEIG, X and W.
!>

LDZ

!>    LDZ (input) INTEGER , LDZ >= M
!>    The leading dimension of the array Z.
!>

RES

!>    RES (output) REAL(KIND=WP) N-by-1 array
!>    RES(1:K) contains the residuals for the K computed
!>    Ritz pairs.
!>    If LAMBDA(i) is real, then
!>       RES(i) = || A * Z(:,i) - LAMBDA(i)*Z(:,i))||_2.
!>    If [LAMBDA(i), LAMBDA(i+1)] is a complex conjugate pair
!>    then
!>    RES(i)=RES(i+1) = || A * Z(:,i:i+1) - Z(:,i:i+1) *B||_F
!>    where B = [ real(LAMBDA(i)) imag(LAMBDA(i)) ]
!>              [-imag(LAMBDA(i)) real(LAMBDA(i)) ].
!>    It holds that
!>    RES(i)   = || A*ZC(:,i)   - LAMBDA(i)  *ZC(:,i)   ||_2
!>    RES(i+1) = || A*ZC(:,i+1) - LAMBDA(i+1)*ZC(:,i+1) ||_2
!>    where ZC(:,i)   =  Z(:,i) + sqrt(-1)*Z(:,i+1)
!>          ZC(:,i+1) =  Z(:,i) - sqrt(-1)*Z(:,i+1)
!>    See the description of REIG, IMEIG and Z.
!>

!>    B (output) REAL(KIND=WP)  M-by-N array.
!>    IF JOBF =='R', B(1:M,1:K) contains A*U(:,1:K), and can
!>    be used for computing the refined vectors; see further
!>    details in the provided references.
!>    If JOBF == 'E', B(1:M,1;K) contains
!>    A*U(:,1:K)*W(1:K,1:K), which are the vectors from the
!>    Exact DMD, up to scaling by the inverse eigenvalues.
!>    If JOBF =='N', then B is not referenced.
!>    See the descriptions of X, W, K.
!>

LDB

!>    LDB (input) INTEGER, LDB >= M
!>    The leading dimension of the array B.
!>

!>    W (workspace/output) REAL(KIND=WP) N-by-N array
!>    On exit, W(1:K,1:K) contains the K computed
!>    eigenvectors of the matrix Rayleigh quotient (real and
!>    imaginary parts for each complex conjugate pair of the
!>    eigenvalues). The Ritz vectors (returned in Z) are the
!>    product of X (containing a POD basis for the input
!>    matrix X) and W. See the descriptions of K, S, X and Z.
!>    W is also used as a workspace to temporarily store the
!>    left singular vectors of X.
!>

LDW

!>    LDW (input) INTEGER, LDW >= N
!>    The leading dimension of the array W.
!>

!>    S (workspace/output) REAL(KIND=WP) N-by-N array
!>    The array S(1:K,1:K) is used for the matrix Rayleigh
!>    quotient. This content is overwritten during
!>    the eigenvalue decomposition by SGEEV.
!>    See the description of K.
!>

LDS

!>    LDS (input) INTEGER, LDS >= N
!>    The leading dimension of the array S.
!>

WORK

!>    WORK (workspace/output) REAL(KIND=WP) LWORK-by-1 array
!>    On exit, WORK(1:N) contains the singular values of
!>    X (for JOBS=='N') or column scaled X (JOBS=='S', 'C').
!>    If WHTSVD==4, then WORK(N+1) and WORK(N+2) contain
!>    scaling factor WORK(N+2)/WORK(N+1) used to scale X
!>    and Y to avoid overflow in the SVD of X.
!>    This may be of interest if the scaling option is off
!>    and as many as possible smallest eigenvalues are
!>    desired to the highest feasible accuracy.
!>    If the call to SGEDMD is only workspace query, then
!>    WORK(1) contains the minimal workspace length and
!>    WORK(2) is the optimal workspace length. Hence, the
!>    length of work is at least 2.
!>    See the description of LWORK.
!>

LWORK

!>    LWORK (input) INTEGER
!>    The minimal length of the workspace vector WORK.
!>    LWORK is calculated as follows:
!>    If WHTSVD == 1 ::
!>       If JOBZ == 'V', then
!>       LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,4*N)).
!>       If JOBZ == 'N'  then
!>       LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,3*N)).
!>       Here LWORK_SVD = MAX(1,3*N+M,5*N) is the minimal
!>       workspace length of SGESVD.
!>    If WHTSVD == 2 ::
!>       If JOBZ == 'V', then
!>       LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,4*N))
!>       If JOBZ == 'N', then
!>       LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,3*N))
!>       Here LWORK_SVD = MAX(M, 5*N*N+4*N)+3*N*N is the
!>       minimal workspace length of SGESDD.
!>    If WHTSVD == 3 ::
!>       If JOBZ == 'V', then
!>       LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,4*N))
!>       If JOBZ == 'N', then
!>       LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,3*N))
!>       Here LWORK_SVD = N+M+MAX(3*N+1,
!>                       MAX(1,3*N+M,5*N),MAX(1,N))
!>       is the minimal workspace length of SGESVDQ.
!>    If WHTSVD == 4 ::
!>       If JOBZ == 'V', then
!>       LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,4*N))
!>       If JOBZ == 'N', then
!>       LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,3*N))
!>       Here LWORK_SVD = MAX(7,2*M+N,6*N+2*N*N) is the
!>       minimal workspace length of SGEJSV.
!>    The above expressions are not simplified in order to
!>    make the usage of WORK more transparent, and for
!>    easier checking. In any case, LWORK >= 2.
!>    If on entry LWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace lengths for both WORK and
!>    IWORK. See the descriptions of WORK and IWORK.
!>

IWORK

!>    IWORK (workspace/output) INTEGER LIWORK-by-1 array
!>    Workspace that is required only if WHTSVD equals
!>    2 , 3 or 4. (See the description of WHTSVD).
!>    If on entry LWORK =-1 or LIWORK=-1, then the
!>    minimal length of IWORK is computed and returned in
!>    IWORK(1). See the description of LIWORK.
!>

LIWORK

!>    LIWORK (input) INTEGER
!>    The minimal length of the workspace vector IWORK.
!>    If WHTSVD == 1, then only IWORK(1) is used; LIWORK >=1
!>    If WHTSVD == 2, then LIWORK >= MAX(1,8*MIN(M,N))
!>    If WHTSVD == 3, then LIWORK >= MAX(1,M+N-1)
!>    If WHTSVD == 4, then LIWORK >= MAX(3,M+3*N)
!>    If on entry LIWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace lengths for both WORK and
!>    IWORK. See the descriptions of WORK and IWORK.
!>

INFO

!>    INFO (output) INTEGER
!>    -i < 0 :: On entry, the i-th argument had an
!>              illegal value
!>       = 0 :: Successful return.
!>       = 1 :: Void input. Quick exit (M=0 or N=0).
!>       = 2 :: The SVD computation of X did not converge.
!>              Suggestion: Check the input data and/or
!>              repeat with different WHTSVD.
!>       = 3 :: The computation of the eigenvalues did not
!>              converge.
!>       = 4 :: If data scaling was requested on input and
!>              the procedure found inconsistency in the data
!>              such that for some column index i,
!>              X(:,i) = 0 but Y(:,i) /= 0, then Y(:,i) is set
!>              to zero if JOBS=='C'. The computation proceeds
!>              with original or modified data and warning
!>              flag is set with INFO=4.
!>

Author

Zlatko Drmac

Definition at line 530 of file sgedmd.f90.

subroutine SGEDMDQ (character, intent(in) jobs, character, intent(in) jobz, character, intent(in) jobr, character, intent(in) jobq, character, intent(in) jobt, character, intent(in) jobf, integer, intent(in) whtsvd, integer, intent(in) m, integer, intent(in) n, real(kind=wp), dimension(ldf,), intent(inout) f, integer, intent(in) ldf, real(kind=wp), dimension(ldx,), intent(out) x, integer, intent(in) ldx, real(kind=wp), dimension(ldy,), intent(out) y, integer, intent(in) ldy, integer, intent(in) nrnk, real(kind=wp), intent(in) tol, integer, intent(out) k, real(kind=wp), dimension(), intent(out) reig, real(kind=wp), dimension(), intent(out) imeig, real(kind=wp), dimension(ldz,), intent(out) z, integer, intent(in) ldz, real(kind=wp), dimension(), intent(out) res, real(kind=wp), dimension(ldb,), intent(out) b, integer, intent(in) ldb, real(kind=wp), dimension(ldv,), intent(out) v, integer, intent(in) ldv, real(kind=wp), dimension(lds,), intent(out) s, integer, intent(in) lds, real(kind=wp), dimension(), intent(out) work, integer, intent(in) lwork, integer, dimension(), intent(out) iwork, integer, intent(in) liwork, integer, intent(out) info)¶

SGEDMDQ computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.

Purpose:

!>    SGEDMDQ computes the Dynamic Mode Decomposition (DMD) for
!>    a pair of data snapshot matrices, using a QR factorization
!>    based compression of the data. For the input matrices
!>    X and Y such that Y = A*X with an unaccessible matrix
!>    A, SGEDMDQ computes a certain number of Ritz pairs of A using
!>    the standard Rayleigh-Ritz extraction from a subspace of
!>    range(X) that is determined using the leading left singular 
!>    vectors of X. Optionally, SGEDMDQ returns the residuals 
!>    of the computed Ritz pairs, the information needed for
!>    a refinement of the Ritz vectors, or the eigenvectors of
!>    the Exact DMD.
!>    For further details see the references listed
!>    below. For more details of the implementation see [3].      
!>

References:

!>    [1] P. Schmid: Dynamic mode decomposition of numerical
!>        and experimental data,
!>        Journal of Fluid Mechanics 656, 5-28, 2010.
!>    [2] Z. Drmac, I. Mezic, R. Mohr: Data driven modal
!>        decompositions: analysis and enhancements,
!>        SIAM J. on Sci. Comp. 40 (4), A2253-A2285, 2018.
!>    [3] Z. Drmac: A LAPACK implementation of the Dynamic
!>        Mode Decomposition I. Technical report. AIMDyn Inc.
!>        and LAPACK Working Note 298.      
!>    [4] J. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. 
!>        Brunton, N. Kutz: On Dynamic Mode Decomposition:
!>        Theory and Applications, Journal of Computational
!>        Dynamics 1(2), 391 -421, 2014.
!>

Developed and supported by:

!>    Developed and coded by Zlatko Drmac, Faculty of Science,
!>    University of Zagreb;  drmac@math.hr
!>    In cooperation with
!>    AIMdyn Inc., Santa Barbara, CA.
!>    and supported by
!>    - DARPA SBIR project  Contract No: W31P4Q-21-C-0007
!>    - DARPA PAI project  Contract No: HR0011-18-9-0033
!>    - DARPA MoDyL project 
!>    Contract No: HR0011-16-C-0116
!>    Any opinions, findings and conclusions or recommendations 
!>    expressed in this material are those of the author and 
!>    do not necessarily reflect the views of the DARPA SBIR 
!>    Program Office.      
!>

Distribution Statement A:

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

Parameters

JOBS

!>    JOBS (input) CHARACTER*1
!>    Determines whether the initial data snapshots are scaled
!>    by a diagonal matrix. The data snapshots are the columns
!>    of F. The leading N-1 columns of F are denoted X and the
!>    trailing N-1 columns are denoted Y. 
!>    'S' :: The data snapshots matrices X and Y are multiplied
!>           with a diagonal matrix D so that X*D has unit
!>           nonzero columns (in the Euclidean 2-norm)
!>    'C' :: The snapshots are scaled as with the 'S' option.
!>           If it is found that an i-th column of X is zero
!>           vector and the corresponding i-th column of Y is
!>           non-zero, then the i-th column of Y is set to
!>           zero and a warning flag is raised.
!>    'Y' :: The data snapshots matrices X and Y are multiplied
!>           by a diagonal matrix D so that Y*D has unit
!>           nonzero columns (in the Euclidean 2-norm)    
!>    'N' :: No data scaling.   
!>

JOBZ

!>    JOBZ (input) CHARACTER*1
!>    Determines whether the eigenvectors (Koopman modes) will
!>    be computed.
!>    'V' :: The eigenvectors (Koopman modes) will be computed
!>           and returned in the matrix Z.
!>           See the description of Z.
!>    'F' :: The eigenvectors (Koopman modes) will be returned
!>           in factored form as the product Z*V, where Z
!>           is orthonormal and V contains the eigenvectors
!>           of the corresponding Rayleigh quotient.
!>           See the descriptions of F, V, Z.
!>    'Q' :: The eigenvectors (Koopman modes) will be returned
!>           in factored form as the product Q*Z, where Z
!>           contains the eigenvectors of the compression of the
!>           underlying discretized operator onto the span of
!>           the data snapshots. See the descriptions of F, V, Z. 
!>           Q is from the initial QR factorization.  
!>    'N' :: The eigenvectors are not computed.  
!>

JOBR

!>    JOBR (input) CHARACTER*1 
!>    Determines whether to compute the residuals.
!>    'R' :: The residuals for the computed eigenpairs will
!>           be computed and stored in the array RES.
!>           See the description of RES.
!>           For this option to be legal, JOBZ must be 'V'.
!>    'N' :: The residuals are not computed.
!>

JOBQ

!>    JOBQ (input) CHARACTER*1 
!>    Specifies whether to explicitly compute and return the
!>    orthogonal matrix from the QR factorization.
!>    'Q' :: The matrix Q of the QR factorization of the data
!>           snapshot matrix is computed and stored in the
!>           array F. See the description of F.       
!>    'N' :: The matrix Q is not explicitly computed.
!>

JOBT

!>    JOBT (input) CHARACTER*1 
!>    Specifies whether to return the upper triangular factor
!>    from the QR factorization.
!>    'R' :: The matrix R of the QR factorization of the data 
!>           snapshot matrix F is returned in the array Y.
!>           See the description of Y and Further details.       
!>    'N' :: The matrix R is not returned.    
!>

JOBF

!>    JOBF (input) CHARACTER*1
!>    Specifies whether to store information needed for post-
!>    processing (e.g. computing refined Ritz vectors)
!>    'R' :: The matrix needed for the refinement of the Ritz
!>           vectors is computed and stored in the array B.
!>           See the description of B.
!>    'E' :: The unscaled eigenvectors of the Exact DMD are 
!>           computed and returned in the array B. See the
!>           description of B.
!>    'N' :: No eigenvector refinement data is computed.   
!>    To be useful on exit, this option needs JOBQ='Q'.      
!>

WHTSVD

!>    WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 }
!>    Allows for a selection of the SVD algorithm from the
!>    LAPACK library.
!>    1 :: SGESVD (the QR SVD algorithm)
!>    2 :: SGESDD (the Divide and Conquer algorithm; if enough
!>         workspace available, this is the fastest option)
!>    3 :: SGESVDQ (the preconditioned QR SVD  ; this and 4
!>         are the most accurate options)
!>    4 :: SGEJSV (the preconditioned Jacobi SVD; this and 3
!>         are the most accurate options)
!>    For the four methods above, a significant difference in
!>    the accuracy of small singular values is possible if
!>    the snapshots vary in norm so that X is severely
!>    ill-conditioned. If small (smaller than EPS*||X||)
!>    singular values are of interest and JOBS=='N',  then
!>    the options (3, 4) give the most accurate results, where
!>    the option 4 is slightly better and with stronger 
!>    theoretical background.
!>    If JOBS=='S', i.e. the columns of X will be normalized,
!>    then all methods give nearly equally accurate results.
!>

!>    M (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 (input/output) REAL(KIND=WP) M-by-N array
!>    > On entry,
!>    the columns of F are the sequence of data snapshots 
!>    from a single trajectory, taken at equidistant discrete
!>    times. It is assumed that the column norms of F are 
!>    in the range of the normalized floating point numbers. 
!>    < On exit,
!>    If JOBQ == 'Q', the array F contains the orthogonal 
!>    matrix/factor of the QR factorization of the initial 
!>    data snapshots matrix F. See the description of JOBQ. 
!>    If JOBQ == 'N', the entries in F strictly below the main
!>    diagonal contain, column-wise, the information on the 
!>    Householder vectors, as returned by SGEQRF. The 
!>    remaining information to restore the orthogonal matrix
!>    of the initial QR factorization is stored in WORK(1:N). 
!>    See the description of WORK.
!>

LDF

!>    LDF (input) INTEGER, LDF >= M 
!>    The leading dimension of the array F.
!>

!>    X (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 (workspace/output) REAL(KIND=WP) MIN(M,N)-by-(N-1) array
!>    Y is used as workspace to hold representations of the
!>    trailing N-1 snapshots in the orthonormal basis computed
!>    in the QR factorization of F.
!>    On exit, 
!>    If JOBT == 'R', Y contains the MIN(M,N)-by-N upper
!>    triangular factor from the QR factorization of the data
!>    snapshot matrix F.
!>

LDY

!>    LDY (input) INTEGER , LDY >= N
!>    The leading dimension of the array Y   
!>

NRNK

!>    NRNK (input) INTEGER
!>    Determines the mode how to compute the numerical rank,
!>    i.e. how to truncate small singular values of the input
!>    matrix X. On input, if
!>    NRNK = -1 :: i-th singular value sigma(i) is truncated
!>                 if sigma(i) <= TOL*sigma(1)
!>                 This option is recommended.   
!>    NRNK = -2 :: i-th singular value sigma(i) is truncated
!>                 if sigma(i) <= TOL*sigma(i-1)
!>                 This option is included for R&D purposes.
!>                 It requires highly accurate SVD, which
!>                 may not be feasible.     
!>    The numerical rank can be enforced by using positive 
!>    value of NRNK as follows: 
!>    0 < NRNK <= N-1 :: at most NRNK largest singular values
!>    will be used. If the number of the computed nonzero
!>    singular values is less than NRNK, then only those
!>    nonzero values will be used and the actually used
!>    dimension is less than NRNK. The actual number of
!>    the nonzero singular values is returned in the variable
!>    K. See the description of K.
!>

TOL

!>    TOL (input) REAL(KIND=WP), 0 <= TOL < 1
!>    The tolerance for truncating small singular values.
!>    See the description of NRNK.  
!>

!>    K (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:
!>    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

!>    LDZ (input) INTEGER , LDZ >= M
!>    The leading dimension of the array Z.
!>

RES

!>    RES (output) REAL(KIND=WP) (N-1)-by-1 array
!>    RES(1:K) contains the residuals for the K computed 
!>    Ritz pairs.       
!>    If LAMBDA(i) is real, then
!>       RES(i) = || A * Z(:,i) - LAMBDA(i)*Z(:,i))||_2.
!>    If [LAMBDA(i), LAMBDA(i+1)] is a complex conjugate pair
!>    then
!>    RES(i)=RES(i+1) = || A * Z(:,i:i+1) - Z(:,i:i+1) *B||_F
!>    where B = [ real(LAMBDA(i)) imag(LAMBDA(i)) ]
!>              [-imag(LAMBDA(i)) real(LAMBDA(i)) ].
!>    It holds that
!>    RES(i)   = || A*ZC(:,i)   - LAMBDA(i)  *ZC(:,i)   ||_2
!>    RES(i+1) = || A*ZC(:,i+1) - LAMBDA(i+1)*ZC(:,i+1) ||_2
!>    where ZC(:,i)   =  Z(:,i) + sqrt(-1)*Z(:,i+1)
!>          ZC(:,i+1) =  Z(:,i) - sqrt(-1)*Z(:,i+1)
!>    See the description of Z.
!>

!>    B (output) REAL(KIND=WP)  MIN(M,N)-by-(N-1) array.
!>    IF JOBF =='R', B(1:N,1:K) contains A*U(:,1:K), and can
!>    be used for computing the refined vectors; see further 
!>    details in the provided references. 
!>    If JOBF == 'E', B(1:N,1;K) contains 
!>    A*U(:,1:K)*W(1:K,1:K), which are the vectors from the
!>    Exact DMD, up to scaling by the inverse eigenvalues.   
!>    In both cases, the content of B can be lifted to the 
!>    original dimension of the input data by pre-multiplying
!>    with the Q factor from the initial QR factorization.     
!>    Here A denotes a compression of the underlying operator.      
!>    See the descriptions of F and X.
!>    If JOBF =='N', then B is not referenced.
!>

LDB

!>    LDB (input) INTEGER, LDB >= MIN(M,N)
!>    The leading dimension of the array B.
!>

!>    V (workspace/output) REAL(KIND=WP) (N-1)-by-(N-1) array
!>    On exit, V(1:K,1:K) contains the K eigenvectors of
!>    the Rayleigh quotient. The eigenvectors of a complex
!>    conjugate pair of eigenvalues are returned in real form
!>    as explained in the description of Z. The Ritz vectors
!>    (returned in Z) are the product of X and V; see
!>    the descriptions of X and Z.
!>

LDV

!>    LDV (input) INTEGER, LDV >= N-1
!>    The leading dimension of the array V.
!>

!>    S (output) REAL(KIND=WP) (N-1)-by-(N-1) array
!>    The array S(1:K,1:K) is used for the matrix Rayleigh
!>    quotient. This content is overwritten during
!>    the eigenvalue decomposition by SGEEV.
!>    See the description of K.
!>

LDS

!>    LDS (input) INTEGER, LDS >= N-1        
!>    The leading dimension of the array S.
!>

WORK

!>    WORK (workspace/output) REAL(KIND=WP) LWORK-by-1 array
!>    On exit, 
!>    WORK(1:MIN(M,N)) contains the scalar factors of the 
!>    elementary reflectors as returned by SGEQRF of the 
!>    M-by-N input matrix F.
!>    WORK(MIN(M,N)+1:MIN(M,N)+N-1) contains the singular values of 
!>    the input submatrix F(1:M,1:N-1).
!>    If the call to SGEDMDQ is only workspace query, then
!>    WORK(1) contains the minimal workspace length and
!>    WORK(2) is the optimal workspace length. Hence, the
!>    length of work is at least 2.
!>    See the description of LWORK.
!>

LWORK

!>    LWORK (input) INTEGER
!>    The minimal length of the  workspace vector WORK.
!>    LWORK is calculated as follows:
!>    Let MLWQR  = N (minimal workspace for SGEQRF[M,N])
!>        MLWDMD = minimal workspace for SGEDMD (see the
!>                 description of LWORK in SGEDMD) for 
!>                 snapshots of dimensions MIN(M,N)-by-(N-1)
!>        MLWMQR = N (minimal workspace for 
!>                   SORMQR['L','N',M,N,N])
!>        MLWGQR = N (minimal workspace for SORGQR[M,N,N])
!>    Then
!>    LWORK = MAX(N+MLWQR, N+MLWDMD)
!>    is updated as follows:
!>       if   JOBZ == 'V' or JOBZ == 'F' THEN 
!>            LWORK = MAX( LWORK,MIN(M,N)+N-1 +MLWMQR )
!>       if   JOBQ == 'Q' THEN
!>            LWORK = MAX( LWORK,MIN(M,N)+N-1+MLWGQR)
!>    If on entry LWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace lengths for both WORK and
!>    IWORK. See the descriptions of WORK and IWORK.          
!>

IWORK

!>    IWORK (workspace/output) INTEGER LIWORK-by-1 array
!>    Workspace that is required only if WHTSVD equals
!>    2 , 3 or 4. (See the description of WHTSVD).
!>    If on entry LWORK =-1 or LIWORK=-1, then the
!>    minimal length of IWORK is computed and returned in
!>    IWORK(1). See the description of LIWORK.
!>

LIWORK

!>    LIWORK (input) INTEGER
!>    The minimal length of the workspace vector IWORK.
!>    If WHTSVD == 1, then only IWORK(1) is used; LIWORK >=1
!>    Let M1=MIN(M,N), N1=N-1. Then
!>    If WHTSVD == 2, then LIWORK >= MAX(1,8*MIN(M1,N1))
!>    If WHTSVD == 3, then LIWORK >= MAX(1,M1+N1-1)
!>    If WHTSVD == 4, then LIWORK >= MAX(3,M1+3*N1)
!>    If on entry LIWORK = -1, then a worskpace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace lengths for both WORK and
!>    IWORK. See the descriptions of WORK and IWORK.
!>

INFO

!>    INFO (output) INTEGER
!>    -i < 0 :: On entry, the i-th argument had an
!>              illegal value
!>       = 0 :: Successful return.
!>       = 1 :: Void input. Quick exit (M=0 or N=0).
!>       = 2 :: The SVD computation of X did not converge.
!>              Suggestion: Check the input data and/or
!>              repeat with different WHTSVD.
!>       = 3 :: The computation of the eigenvalues did not
!>              converge.
!>       = 4 :: If data scaling was requested on input and
!>              the procedure found inconsistency in the data
!>              such that for some column index i,
!>              X(:,i) = 0 but Y(:,i) /= 0, then Y(:,i) is set
!>              to zero if JOBS=='C'. The computation proceeds
!>              with original or modified data and warning
!>              flag is set with INFO=4.  
!>

Author

Zlatko Drmac

Definition at line 571 of file sgedmdq.f90.

subroutine ZGEDMD (character, intent(in) jobs, character, intent(in) jobz, character, intent(in) jobr, character, intent(in) jobf, integer, intent(in) whtsvd, integer, intent(in) m, integer, intent(in) n, complex(kind=wp), dimension(ldx,), intent(inout) x, integer, intent(in) ldx, complex(kind=wp), dimension(ldy,), intent(inout) y, integer, intent(in) ldy, integer, intent(in) nrnk, real(kind=wp), intent(in) tol, integer, intent(out) k, complex(kind=wp), dimension(), intent(out) eigs, complex(kind=wp), dimension(ldz,), intent(out) z, integer, intent(in) ldz, real(kind=wp), dimension(), intent(out) res, complex(kind=wp), dimension(ldb,), intent(out) b, integer, intent(in) ldb, complex(kind=wp), dimension(ldw,), intent(out) w, integer, intent(in) ldw, complex(kind=wp), dimension(lds,), intent(out) s, integer, intent(in) lds, complex(kind=wp), dimension(), intent(out) zwork, integer, intent(in) lzwork, real(kind=wp), dimension(), intent(out) rwork, integer, intent(in) lrwork, integer, dimension(*), intent(out) iwork, integer, intent(in) liwork, integer, intent(out) info)¶

ZGEDMD computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.

Purpose:

!>    ZGEDMD computes the Dynamic Mode Decomposition (DMD) for
!>    a pair of data snapshot matrices. For the input matrices
!>    X and Y such that Y = A*X with an unaccessible matrix
!>    A, ZGEDMD computes a certain number of Ritz pairs of A using
!>    the standard Rayleigh-Ritz extraction from a subspace of
!>    range(X) that is determined using the leading left singular
!>    vectors of X. Optionally, ZGEDMD returns the residuals
!>    of the computed Ritz pairs, the information needed for
!>    a refinement of the Ritz vectors, or the eigenvectors of
!>    the Exact DMD.
!>    For further details see the references listed
!>    below. For more details of the implementation see [3].
!>

References:

!>    [1] P. Schmid: Dynamic mode decomposition of numerical
!>        and experimental data,
!>        Journal of Fluid Mechanics 656, 5-28, 2010.
!>    [2] Z. Drmac, I. Mezic, R. Mohr: Data driven modal
!>        decompositions: analysis and enhancements,
!>        SIAM J. on Sci. Comp. 40 (4), A2253-A2285, 2018.
!>    [3] Z. Drmac: A LAPACK implementation of the Dynamic
!>        Mode Decomposition I. Technical report. AIMDyn Inc.
!>        and LAPACK Working Note 298.
!>    [4] J. Tu, C. W. Rowley, D. M. Luchtenburg, S. L.
!>        Brunton, N. Kutz: On Dynamic Mode Decomposition:
!>        Theory and Applications, Journal of Computational
!>        Dynamics 1(2), 391 -421, 2014.
!>

Developed and supported by:

!>    Developed and coded by Zlatko Drmac, Faculty of Science,
!>    University of Zagreb;  drmac@math.hr
!>    In cooperation with
!>    AIMdyn Inc., Santa Barbara, CA.
!>    and supported by
!>    - DARPA SBIR project  Contract No: W31P4Q-21-C-0007
!>    - DARPA PAI project  Contract No: HR0011-18-9-0033
!>    - DARPA MoDyL project 
!>    Contract No: HR0011-16-C-0116
!>    Any opinions, findings and conclusions or recommendations
!>    expressed in this material are those of the author and
!>    do not necessarily reflect the views of the DARPA SBIR
!>    Program Office
!>

Distribution Statement A:

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

Parameters

JOBS

!>    JOBS (input) CHARACTER*1
!>    Determines whether the initial data snapshots are scaled
!>    by a diagonal matrix.
!>    'S' :: The data snapshots matrices X and Y are multiplied
!>           with a diagonal matrix D so that X*D has unit
!>           nonzero columns (in the Euclidean 2-norm)
!>    'C' :: The snapshots are scaled as with the 'S' option.
!>           If it is found that an i-th column of X is zero
!>           vector and the corresponding i-th column of Y is
!>           non-zero, then the i-th column of Y is set to
!>           zero and a warning flag is raised.
!>    'Y' :: The data snapshots matrices X and Y are multiplied
!>           by a diagonal matrix D so that Y*D has unit
!>           nonzero columns (in the Euclidean 2-norm)
!>    'N' :: No data scaling.
!>

JOBZ

!>    JOBZ (input) CHARACTER*1
!>    Determines whether the eigenvectors (Koopman modes) will
!>    be computed.
!>    'V' :: The eigenvectors (Koopman modes) will be computed
!>           and returned in the matrix Z.
!>           See the description of Z.
!>    'F' :: The eigenvectors (Koopman modes) will be returned
!>           in factored form as the product X(:,1:K)*W, where X
!>           contains a POD basis (leading left singular vectors
!>           of the data matrix X) and W contains the eigenvectors
!>           of the corresponding Rayleigh quotient.
!>           See the descriptions of K, X, W, Z.
!>    'N' :: The eigenvectors are not computed.
!>

JOBR

!>    JOBR (input) CHARACTER*1
!>    Determines whether to compute the residuals.
!>    'R' :: The residuals for the computed eigenpairs will be
!>           computed and stored in the array RES.
!>           See the description of RES.
!>           For this option to be legal, JOBZ must be 'V'.
!>    'N' :: The residuals are not computed.
!>

JOBF

!>    JOBF (input) CHARACTER*1
!>    Specifies whether to store information needed for post-
!>    processing (e.g. computing refined Ritz vectors)
!>    'R' :: The matrix needed for the refinement of the Ritz
!>           vectors is computed and stored in the array B.
!>           See the description of B.
!>    'E' :: The unscaled eigenvectors of the Exact DMD are
!>           computed and returned in the array B. See the
!>           description of B.
!>    'N' :: No eigenvector refinement data is computed.
!>

WHTSVD

!>    WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 }
!>    Allows for a selection of the SVD algorithm from the
!>    LAPACK library.
!>    1 :: ZGESVD (the QR SVD algorithm)
!>    2 :: ZGESDD (the Divide and Conquer algorithm; if enough
!>         workspace available, this is the fastest option)
!>    3 :: ZGESVDQ (the preconditioned QR SVD  ; this and 4
!>         are the most accurate options)
!>    4 :: ZGEJSV (the preconditioned Jacobi SVD; this and 3
!>         are the most accurate options)
!>    For the four methods above, a significant difference in
!>    the accuracy of small singular values is possible if
!>    the snapshots vary in norm so that X is severely
!>    ill-conditioned. If small (smaller than EPS*||X||)
!>    singular values are of interest and JOBS=='N',  then
!>    the options (3, 4) give the most accurate results, where
!>    the option 4 is slightly better and with stronger
!>    theoretical background.
!>    If JOBS=='S', i.e. the columns of X will be normalized,
!>    then all methods give nearly equally accurate results.
!>

!>    M (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).
!>

LDX

!>    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 (input) INTEGER, LDX >= M
!>    The leading dimension of the array X.
!>

!>    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 (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 (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 (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 (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 (workspace/output) COMPLEX(KIND=WP) N-by-N array
!>    The array S(1:K,1:K) is used for the matrix Rayleigh
!>    quotient. This content is overwritten during
!>    the eigenvalue decomposition by ZGEEV.
!>    See the description of K.
!>

LDS

!>    LDS (input) INTEGER, LDS >= N
!>    The leading dimension of the array S.
!>

ZWORK

!>    ZWORK (workspace/output) COMPLEX(KIND=WP) LZWORK-by-1 array
!>    ZWORK is used as complex workspace in the complex SVD, as
!>    specified by WHTSVD (1,2, 3 or 4) and for ZGEEV for computing
!>    the eigenvalues of a Rayleigh quotient.
!>    If the call to ZGEDMD is only workspace query, then
!>    ZWORK(1) contains the minimal complex workspace length and
!>    ZWORK(2) is the optimal complex workspace length.
!>    Hence, the length of work is at least 2.
!>    See the description of LZWORK.
!>

LZWORK

!>    LZWORK (input) INTEGER
!>    The minimal length of the workspace vector ZWORK.
!>    LZWORK is calculated as MAX(LZWORK_SVD, LZWORK_ZGEEV),
!>    where LZWORK_ZGEEV = MAX( 1, 2*N )  and the minimal
!>    LZWORK_SVD is calculated as follows
!>    If WHTSVD == 1 :: ZGESVD ::
!>       LZWORK_SVD = MAX(1,2*MIN(M,N)+MAX(M,N))
!>    If WHTSVD == 2 :: ZGESDD ::
!>       LZWORK_SVD = 2*MIN(M,N)*MIN(M,N)+2*MIN(M,N)+MAX(M,N)
!>    If WHTSVD == 3 :: ZGESVDQ ::
!>       LZWORK_SVD = obtainable by a query
!>    If WHTSVD == 4 :: ZGEJSV ::
!>       LZWORK_SVD = obtainable by a query
!>    If on entry LZWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace lengths and returns them in
!>    LZWORK(1) and LZWORK(2), respectively.
!>

RWORK

!>    RWORK (workspace/output) REAL(KIND=WP) LRWORK-by-1 array
!>    On exit, RWORK(1:N) contains the singular values of
!>    X (for JOBS=='N') or column scaled X (JOBS=='S', 'C').
!>    If WHTSVD==4, then RWORK(N+1) and RWORK(N+2) contain
!>    scaling factor RWORK(N+2)/RWORK(N+1) used to scale X
!>    and Y to avoid overflow in the SVD of X.
!>    This may be of interest if the scaling option is off
!>    and as many as possible smallest eigenvalues are
!>    desired to the highest feasible accuracy.
!>    If the call to ZGEDMD is only workspace query, then
!>    RWORK(1) contains the minimal workspace length.
!>    See the description of LRWORK.
!>

LRWORK

!>    LRWORK (input) INTEGER
!>    The minimal length of the workspace vector RWORK.
!>    LRWORK is calculated as follows:
!>    LRWORK = MAX(1, N+LRWORK_SVD,N+LRWORK_ZGEEV), where
!>    LRWORK_ZGEEV = MAX(1,2*N) and RWORK_SVD is the real workspace
!>    for the SVD subroutine determined by the input parameter
!>    WHTSVD.
!>    If WHTSVD == 1 :: ZGESVD ::
!>       LRWORK_SVD = 5*MIN(M,N)
!>    If WHTSVD == 2 :: ZGESDD ::
!>       LRWORK_SVD =  MAX(5*MIN(M,N)*MIN(M,N)+7*MIN(M,N),
!>       2*MAX(M,N)*MIN(M,N)+2*MIN(M,N)*MIN(M,N)+MIN(M,N) ) )
!>    If WHTSVD == 3 :: ZGESVDQ ::
!>       LRWORK_SVD = obtainable by a query
!>    If WHTSVD == 4 :: ZGEJSV ::
!>       LRWORK_SVD = obtainable by a query
!>    If on entry LRWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    real workspace length and returns it in RWORK(1).
!>

IWORK

!>    IWORK (workspace/output) INTEGER LIWORK-by-1 array
!>    Workspace that is required only if WHTSVD equals
!>    2 , 3 or 4. (See the description of WHTSVD).
!>    If on entry LWORK =-1 or LIWORK=-1, then the
!>    minimal length of IWORK is computed and returned in
!>    IWORK(1). See the description of LIWORK.
!>

LIWORK

!>    LIWORK (input) INTEGER
!>    The minimal length of the workspace vector IWORK.
!>    If WHTSVD == 1, then only IWORK(1) is used; LIWORK >=1
!>    If WHTSVD == 2, then LIWORK >= MAX(1,8*MIN(M,N))
!>    If WHTSVD == 3, then LIWORK >= MAX(1,M+N-1)
!>    If WHTSVD == 4, then LIWORK >= MAX(3,M+3*N)
!>    If on entry LIWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace lengths for  ZWORK, RWORK and
!>    IWORK. See the descriptions of ZWORK, RWORK and IWORK.
!>

INFO

!>    INFO (output) INTEGER
!>    -i < 0 :: On entry, the i-th argument had an
!>              illegal value
!>       = 0 :: Successful return.
!>       = 1 :: Void input. Quick exit (M=0 or N=0).
!>       = 2 :: The SVD computation of X did not converge.
!>              Suggestion: Check the input data and/or
!>              repeat with different WHTSVD.
!>       = 3 :: The computation of the eigenvalues did not
!>              converge.
!>       = 4 :: If data scaling was requested on input and
!>              the procedure found inconsistency in the data
!>              such that for some column index i,
!>              X(:,i) = 0 but Y(:,i) /= 0, then Y(:,i) is set
!>              to zero if JOBS=='C'. The computation proceeds
!>              with original or modified data and warning
!>              flag is set with INFO=4.
!>

Author

Zlatko Drmac

Definition at line 493 of file zgedmd.f90.

subroutine ZGEDMDQ (character, intent(in) jobs, character, intent(in) jobz, character, intent(in) jobr, character, intent(in) jobq, character, intent(in) jobt, character, intent(in) jobf, integer, intent(in) whtsvd, integer, intent(in) m, integer, intent(in) n, complex(kind=wp), dimension(ldf,), intent(inout) f, integer, intent(in) ldf, complex(kind=wp), dimension(ldx,), intent(out) x, integer, intent(in) ldx, complex(kind=wp), dimension(ldy,), intent(out) y, integer, intent(in) ldy, integer, intent(in) nrnk, real(kind=wp), intent(in) tol, integer, intent(out) k, complex(kind=wp), dimension(), intent(out) eigs, complex(kind=wp), dimension(ldz,), intent(out) z, integer, intent(in) ldz, real(kind=wp), dimension(), intent(out) res, complex(kind=wp), dimension(ldb,), intent(out) b, integer, intent(in) ldb, complex(kind=wp), dimension(ldv,), intent(out) v, integer, intent(in) ldv, complex(kind=wp), dimension(lds,), intent(out) s, integer, intent(in) lds, complex(kind=wp), dimension(), intent(out) zwork, integer, intent(in) lzwork, real(kind=wp), dimension(), intent(out) work, integer, intent(in) lwork, integer, dimension(), intent(out) iwork, integer, intent(in) liwork, integer, intent(out) info)¶

ZGEDMDQ computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.

Purpose:

!>    ZGEDMDQ computes the Dynamic Mode Decomposition (DMD) for
!>    a pair of data snapshot matrices, using a QR factorization
!>    based compression of the data. For the input matrices
!>    X and Y such that Y = A*X with an unaccessible matrix
!>    A, ZGEDMDQ computes a certain number of Ritz pairs of A using
!>    the standard Rayleigh-Ritz extraction from a subspace of
!>    range(X) that is determined using the leading left singular 
!>    vectors of X. Optionally, ZGEDMDQ returns the residuals 
!>    of the computed Ritz pairs, the information needed for
!>    a refinement of the Ritz vectors, or the eigenvectors of
!>    the Exact DMD.
!>    For further details see the references listed
!>    below. For more details of the implementation see [3].      
!>

References:

!>    [1] P. Schmid: Dynamic mode decomposition of numerical
!>        and experimental data,
!>        Journal of Fluid Mechanics 656, 5-28, 2010.
!>    [2] Z. Drmac, I. Mezic, R. Mohr: Data driven modal
!>        decompositions: analysis and enhancements,
!>        SIAM J. on Sci. Comp. 40 (4), A2253-A2285, 2018.
!>    [3] Z. Drmac: A LAPACK implementation of the Dynamic
!>        Mode Decomposition I. Technical report. AIMDyn Inc.
!>        and LAPACK Working Note 298.      
!>    [4] J. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. 
!>        Brunton, N. Kutz: On Dynamic Mode Decomposition:
!>        Theory and Applications, Journal of Computational
!>        Dynamics 1(2), 391 -421, 2014.
!>

Developed and supported by:

!>    Developed and coded by Zlatko Drmac, Faculty of Science,
!>    University of Zagreb;  drmac@math.hr
!>    In cooperation with
!>    AIMdyn Inc., Santa Barbara, CA.
!>    and supported by
!>    - DARPA SBIR project  Contract No: W31P4Q-21-C-0007
!>    - DARPA PAI project  Contract No: HR0011-18-9-0033
!>    - DARPA MoDyL project 
!>    Contract No: HR0011-16-C-0116
!>    Any opinions, findings and conclusions or recommendations 
!>    expressed in this material are those of the author and 
!>    do not necessarily reflect the views of the DARPA SBIR 
!>    Program Office.      
!>

Developed and supported by:

!>    Distribution Statement A: 
!>    Approved for Public Release, Distribution Unlimited.
!>    Cleared by DARPA on September 29, 2022
!>

Parameters

JOBS

!>    JOBS (input) CHARACTER*1
!>    Determines whether the initial data snapshots are scaled
!>    by a diagonal matrix. The data snapshots are the columns
!>    of F. The leading N-1 columns of F are denoted X and the
!>    trailing N-1 columns are denoted Y. 
!>    'S' :: The data snapshots matrices X and Y are multiplied
!>           with a diagonal matrix D so that X*D has unit
!>           nonzero columns (in the Euclidean 2-norm)
!>    'C' :: The snapshots are scaled as with the 'S' option.
!>           If it is found that an i-th column of X is zero
!>           vector and the corresponding i-th column of Y is
!>           non-zero, then the i-th column of Y is set to
!>           zero and a warning flag is raised.
!>    'Y' :: The data snapshots matrices X and Y are multiplied
!>           by a diagonal matrix D so that Y*D has unit
!>           nonzero columns (in the Euclidean 2-norm)    
!>    'N' :: No data scaling.   
!>

JOBZ

!>    JOBZ (input) CHARACTER*1
!>    Determines whether the eigenvectors (Koopman modes) will
!>    be computed.
!>    'V' :: The eigenvectors (Koopman modes) will be computed
!>           and returned in the matrix Z.
!>           See the description of Z.
!>    'F' :: The eigenvectors (Koopman modes) will be returned
!>           in factored form as the product Z*V, where Z
!>           is orthonormal and V contains the eigenvectors
!>           of the corresponding Rayleigh quotient.
!>           See the descriptions of F, V, Z.
!>    'Q' :: The eigenvectors (Koopman modes) will be returned
!>           in factored form as the product Q*Z, where Z
!>           contains the eigenvectors of the compression of the
!>           underlying discretized operator onto the span of
!>           the data snapshots. See the descriptions of F, V, Z.
!>           Q is from the initial QR factorization.  
!>    'N' :: The eigenvectors are not computed.  
!>

JOBR

!>    JOBR (input) CHARACTER*1 
!>    Determines whether to compute the residuals.
!>    'R' :: The residuals for the computed eigenpairs will
!>           be computed and stored in the array RES.
!>           See the description of RES.
!>           For this option to be legal, JOBZ must be 'V'.
!>    'N' :: The residuals are not computed.
!>

JOBQ

!>    JOBQ (input) CHARACTER*1 
!>    Specifies whether to explicitly compute and return the
!>    unitary matrix from the QR factorization.
!>    'Q' :: The matrix Q of the QR factorization of the data
!>           snapshot matrix is computed and stored in the
!>           array F. See the description of F.       
!>    'N' :: The matrix Q is not explicitly computed.
!>

JOBT

!>    JOBT (input) CHARACTER*1 
!>    Specifies whether to return the upper triangular factor
!>    from the QR factorization.
!>    'R' :: The matrix R of the QR factorization of the data 
!>           snapshot matrix F is returned in the array Y.
!>           See the description of Y and Further details.       
!>    'N' :: The matrix R is not returned. 
!>

JOBF

!>    JOBF (input) CHARACTER*1
!>    Specifies whether to store information needed for post-
!>    processing (e.g. computing refined Ritz vectors)
!>    'R' :: The matrix needed for the refinement of the Ritz
!>           vectors is computed and stored in the array B.
!>           See the description of B.
!>    'E' :: The unscaled eigenvectors of the Exact DMD are 
!>           computed and returned in the array B. See the
!>           description of B.
!>    'N' :: No eigenvector refinement data is computed.   
!>    To be useful on exit, this option needs JOBQ='Q'.    
!>

WHTSVD WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 } Allows for a selection of the SVD algorithm from the LAPACK library. 1 :: ZGESVD (the QR SVD algorithm) 2 :: ZGESDD (the Divide and Conquer algorithm; if enough workspace available, this is the fastest option) 3 :: ZGESVDQ (the preconditioned QR SVD ; this and 4 are the most accurate options) 4 :: ZGEJSV (the preconditioned Jacobi SVD; this and 3 are the most accurate options) For the four methods above, a significant difference in the accuracy of small singular values is possible if the snapshots vary in norm so that X is severely ill-conditioned. If small (smaller than EPS*||X||) singular values are of interest and JOBS=='N', then the options (3, 4) give the most accurate results, where the option 4 is slightly better and with stronger theoretical background. If JOBS=='S', i.e. the columns of X will be normalized, then all methods give nearly equally accurate results.
M

!>    M (input) INTEGER, M >= 0 
!>    The state space dimension (the number of rows of F).
!>

!>    N (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 (input/output) COMPLEX(KIND=WP) M-by-N array
!>    > On entry,
!>    the columns of F are the sequence of data snapshots 
!>    from a single trajectory, taken at equidistant discrete
!>    times. It is assumed that the column norms of F are 
!>    in the range of the normalized floating point numbers. 
!>    < On exit,
!>    If JOBQ == 'Q', the array F contains the orthogonal 
!>    matrix/factor of the QR factorization of the initial 
!>    data snapshots matrix F. See the description of JOBQ. 
!>    If JOBQ == 'N', the entries in F strictly below the main
!>    diagonal contain, column-wise, the information on the 
!>    Householder vectors, as returned by ZGEQRF. The 
!>    remaining information to restore the orthogonal matrix
!>    of the initial QR factorization is stored in ZWORK(1:MIN(M,N)). 
!>    See the description of ZWORK.
!>

LDF

!>    LDF (input) INTEGER, LDF >= M 
!>    The leading dimension of the array F.
!>

!>    X (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 (workspace/output) COMPLEX(KIND=WP) MIN(M,N)-by-(N) array
!>    Y is used as workspace to hold representations of the
!>    trailing N-1 snapshots in the orthonormal basis computed
!>    in the QR factorization of F.
!>    On exit, 
!>    If JOBT == 'R', Y contains the MIN(M,N)-by-N upper
!>    triangular factor from the QR factorization of the data
!>    snapshot matrix F.
!>

LDY

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

NRNK

!>    NRNK (input) INTEGER
!>    Determines the mode how to compute the numerical rank,
!>    i.e. how to truncate small singular values of the input
!>    matrix X. On input, if
!>    NRNK = -1 :: i-th singular value sigma(i) is truncated
!>                 if sigma(i) <= TOL*sigma(1)
!>                 This option is recommended.  
!>    NRNK = -2 :: i-th singular value sigma(i) is truncated
!>                 if sigma(i) <= TOL*sigma(i-1)
!>                 This option is included for R&D purposes.
!>                 It requires highly accurate SVD, which
!>                 may not be feasible.      
!>    The numerical rank can be enforced by using positive 
!>    value of NRNK as follows: 
!>    0 < NRNK <= N-1 :: at most NRNK largest singular values
!>    will be used. If the number of the computed nonzero
!>    singular values is less than NRNK, then only those
!>    nonzero values will be used and the actually used
!>    dimension is less than NRNK. The actual number of
!>    the nonzero singular values is returned in the variable
!>    K. See the description of K.
!>

TOL

!>    TOL (input) REAL(KIND=WP), 0 <= TOL < 1
!>    The tolerance for truncating small singular values.
!>    See the description of NRNK.  
!>

!>    K (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 (workspace/output) COMPLEX(KIND=WP)  M-by-(N-1) array
!>    If JOBZ =='V' then Z contains the Ritz vectors. Z(:,i)
!>    is an eigenvector of the i-th Ritz value; ||Z(:,i)||_2=1.
!>    If JOBZ == 'F', then the Z(:,i)'s are given implicitly as
!>    Z*V, where Z contains orthonormal matrix (the product of
!>    Q from the initial QR factorization and the SVD/POD_basis
!>    returned by ZGEDMD in X) and the second factor (the 
!>    eigenvectors of the Rayleigh quotient) is in the array V, 
!>    as returned by ZGEDMD. That is,  X(:,1:K)*V(:,i)
!>    is an eigenvector corresponding to EIGS(i). The columns 
!>    of V(1:K,1:K) are the computed eigenvectors of the 
!>    K-by-K Rayleigh quotient.  
!>    See the descriptions of EIGS, X and V.      
!>

LDZ

!>    LDZ (input) INTEGER , LDZ >= M
!>    The leading dimension of the array Z.
!>

RES

!>    RES (output) REAL(KIND=WP) (N-1)-by-1 array
!>    RES(1:K) contains the residuals for the K computed 
!>    Ritz pairs, 
!>    RES(i) = || A * Z(:,i) - EIGS(i)*Z(:,i))||_2.
!>    See the description of EIGS and Z.      
!>

!>    B (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 (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 (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.
!>

LZWORK

!>    ZWORK (workspace/output) COMPLEX(KIND=WP) LWORK-by-1 array
!>    On exit, 
!>    ZWORK(1:MIN(M,N)) contains the scalar factors of the 
!>    elementary reflectors as returned by ZGEQRF of the 
!>    M-by-N input matrix F.   
!>    If the call to ZGEDMDQ is only workspace query, then
!>    ZWORK(1) contains the minimal complex workspace length and
!>    ZWORK(2) is the optimal complex workspace length. 
!>    Hence, the length of work is at least 2.
!>    See the description of LZWORK.      
!>

LZWORK

!>    LZWORK (input) INTEGER
!>    The minimal length of the  workspace vector ZWORK.
!>    LZWORK is calculated as follows:
!>    Let MLWQR  = N (minimal workspace for ZGEQRF[M,N])
!>        MLWDMD = minimal workspace for ZGEDMD (see the
!>                 description of LWORK in ZGEDMD)
!>        MLWMQR = N (minimal workspace for 
!>                   ZUNMQR['L','N',M,N,N])
!>        MLWGQR = N (minimal workspace for ZUNGQR[M,N,N])
!>        MINMN  = MIN(M,N)      
!>    Then
!>    LZWORK = MAX(2, MIN(M,N)+MLWQR, MINMN+MLWDMD)
!>    is further updated as follows:
!>       if   JOBZ == 'V' or JOBZ == 'F' THEN 
!>            LZWORK = MAX(LZWORK, MINMN+MLWMQR)
!>       if   JOBQ == 'Q' THEN
!>            LZWORK = MAX(ZLWORK, MINMN+MLWGQR)      
!>

WORK

!>    WORK (workspace/output) REAL(KIND=WP) LWORK-by-1 array
!>    On exit,
!>    WORK(1:N-1) contains the singular values of 
!>    the input submatrix F(1:M,1:N-1).
!>    If the call to ZGEDMDQ is only workspace query, then
!>    WORK(1) contains the minimal workspace length and
!>    WORK(2) is the optimal workspace length. hence, the
!>    length of work is at least 2.
!>    See the description of LWORK.
!>

LWORK

!>    LWORK (input) INTEGER
!>    The minimal length of the  workspace vector WORK.
!>    LWORK is the same as in ZGEDMD, because in ZGEDMDQ
!>    only ZGEDMD requires real workspace for snapshots
!>    of dimensions MIN(M,N)-by-(N-1). 
!>    If on entry LWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace length for WORK.          
!>

IWORK

!>    IWORK (workspace/output) INTEGER LIWORK-by-1 array
!>    Workspace that is required only if WHTSVD equals
!>    2 , 3 or 4. (See the description of WHTSVD).
!>    If on entry LWORK =-1 or LIWORK=-1, then the
!>    minimal length of IWORK is computed and returned in
!>    IWORK(1). See the description of LIWORK.
!>

LIWORK

!>    LIWORK (input) INTEGER
!>    The minimal length of the workspace vector IWORK.
!>    If WHTSVD == 1, then only IWORK(1) is used; LIWORK >=1
!>    Let M1=MIN(M,N), N1=N-1. Then
!>    If WHTSVD == 2, then LIWORK >= MAX(1,8*MIN(M1,N1))
!>    If WHTSVD == 3, then LIWORK >= MAX(1,M1+N1-1)
!>    If WHTSVD == 4, then LIWORK >= MAX(3,M1+3*N1)
!>    If on entry LIWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace lengths for both WORK and
!>    IWORK. See the descriptions of WORK and IWORK.
!>

INFO

!>    INFO (output) INTEGER
!>    -i < 0 :: On entry, the i-th argument had an
!>              illegal value
!>       = 0 :: Successful return.
!>       = 1 :: Void input. Quick exit (M=0 or N=0).
!>       = 2 :: The SVD computation of X did not converge.
!>              Suggestion: Check the input data and/or
!>              repeat with different WHTSVD.
!>       = 3 :: The computation of the eigenvalues did not
!>              converge.
!>       = 4 :: If data scaling was requested on input and
!>              the procedure found inconsistency in the data
!>              such that for some column index i,
!>              X(:,i) = 0 but Y(:,i) /= 0, then Y(:,i) is set
!>              to zero if JOBS=='C'. The computation proceeds
!>              with original or modified data and warning
!>              flag is set with INFO=4.  
!>

Author

Zlatko Drmac

Definition at line 550 of file zgedmdq.f90.

Author¶

Generated automatically by Doxygen for LAPACK from the source code.

Version 3.12.0

LAPACK

Source file:	gedmd.3.en.gz (from lapack-man 3.12.0-2.1)
Source last updated:	2024-07-04T11:54:09Z
Converted to HTML:	2024-09-04T01:09:27Z