LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dgedmdq()

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

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

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
[in]JOBS
!>    JOBS (input) CHARACTER*1
!>    Determines whether the initial data snapshots are scaled
!>    by a diagonal matrix. The data snapshots are the columns
!>    of F. The leading N-1 columns of F are denoted X and the
!>    trailing N-1 columns are denoted Y.
!>    'S' :: The data snapshots matrices X and Y are multiplied
!>           with a diagonal matrix D so that X*D has unit
!>           nonzero columns (in the Euclidean 2-norm)
!>    'C' :: The snapshots are scaled as with the 'S' option.
!>           If it is found that an i-th column of X is zero
!>           vector and the corresponding i-th column of Y is
!>           non-zero, then the i-th column of Y is set to
!>           zero and a warning flag is raised.
!>    'Y' :: The data snapshots matrices X and Y are multiplied
!>           by a diagonal matrix D so that Y*D has unit
!>           nonzero columns (in the Euclidean 2-norm)
!>    'N' :: No data scaling.
!>
[in]JOBZ
!>    JOBZ (input) CHARACTER*1
!>    Determines whether the eigenvectors (Koopman modes) will
!>    be computed.
!>    'V' :: The eigenvectors (Koopman modes) will be computed
!>           and returned in the matrix Z.
!>           See the description of Z.
!>    'F' :: The eigenvectors (Koopman modes) will be returned
!>           in factored form as the product Z*V, where Z
!>           is orthonormal and V contains the eigenvectors
!>           of the corresponding Rayleigh quotient.
!>           See the descriptions of F, V, Z.
!>    'Q' :: The eigenvectors (Koopman modes) will be returned
!>           in factored form as the product Q*Z, where Z
!>           contains the eigenvectors of the compression of the
!>           underlying discretized operator onto the span of
!>           the data snapshots. See the descriptions of F, V, Z.
!>           Q is from the initial QR factorization.
!>    'N' :: The eigenvectors are not computed.
!>
[in]JOBR
!>    JOBR (input) CHARACTER*1
!>    Determines whether to compute the residuals.
!>    'R' :: The residuals for the computed eigenpairs will
!>           be computed and stored in the array RES.
!>           See the description of RES.
!>           For this option to be legal, JOBZ must be 'V'.
!>    'N' :: The residuals are not computed.
!>
[in]JOBQ
!>    JOBQ (input) CHARACTER*1
!>    Specifies whether to explicitly compute and return the
!>    orthogonal matrix from the QR factorization.
!>    'Q' :: The matrix Q of the QR factorization of the data
!>           snapshot matrix is computed and stored in the
!>           array F. See the description of F.
!>    'N' :: The matrix Q is not explicitly computed.
!>
[in]JOBT
!>    JOBT (input) CHARACTER*1
!>    Specifies whether to return the upper triangular factor
!>    from the QR factorization.
!>    'R' :: The matrix R of the QR factorization of the data
!>           snapshot matrix F is returned in the array Y.
!>           See the description of Y and Further details.
!>    'N' :: The matrix R is not returned.
!>
[in]JOBF
!>    JOBF (input) CHARACTER*1
!>    Specifies whether to store information needed for post-
!>    processing (e.g. computing refined Ritz vectors)
!>    'R' :: The matrix needed for the refinement of the Ritz
!>           vectors is computed and stored in the array B.
!>           See the description of B.
!>    'E' :: The unscaled eigenvectors of the Exact DMD are
!>           computed and returned in the array B. See the
!>           description of B.
!>    'N' :: No eigenvector refinement data is computed.
!>    To be useful on exit, this option needs JOBQ='Q'.
!>
[in]WHTSVD
!>    WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 }
!>    Allows for a selection of the SVD algorithm from the
!>    LAPACK library.
!>    1 :: DGESVD (the QR SVD algorithm)
!>    2 :: DGESDD (the Divide and Conquer algorithm; if enough
!>         workspace available, this is the fastest option)
!>    3 :: DGESVDQ (the preconditioned QR SVD  ; this and 4
!>         are the most accurate options)
!>    4 :: DGEJSV (the preconditioned Jacobi SVD; this and 3
!>         are the most accurate options)
!>    For the four methods above, a significant difference in
!>    the accuracy of small singular values is possible if
!>    the snapshots vary in norm so that X is severely
!>    ill-conditioned. If small (smaller than EPS*||X||)
!>    singular values are of interest and JOBS=='N',  then
!>    the options (3, 4) give the most accurate results, where
!>    the option 4 is slightly better and with stronger
!>    theoretical background.
!>    If JOBS=='S', i.e. the columns of X will be normalized,
!>    then all methods give nearly equally accurate results.
!>
[in]M
!>    M (input) INTEGER, M >= 0
!>    The state space dimension (the number of rows of F).
!>
[in]N
!>    N (input) INTEGER, 0 <= N <= M
!>    The number of data snapshots from a single trajectory,
!>    taken at equidistant discrete times. This is the
!>    number of columns of F.
!>
[in,out]F
!>    F (input/output) REAL(KIND=WP) M-by-N array
!>    > On entry,
!>    the columns of F are the sequence of data snapshots
!>    from a single trajectory, taken at equidistant discrete
!>    times. It is assumed that the column norms of F are
!>    in the range of the normalized floating point numbers.
!>    < On exit,
!>    If JOBQ == 'Q', the array F contains the orthogonal
!>    matrix/factor of the QR factorization of the initial
!>    data snapshots matrix F. See the description of JOBQ.
!>    If JOBQ == 'N', the entries in F strictly below the main
!>    diagonal contain, column-wise, the information on the
!>    Householder vectors, as returned by DGEQRF. The
!>    remaining information to restore the orthogonal matrix
!>    of the initial QR factorization is stored in WORK(1:N).
!>    See the description of WORK.
!>
[in]LDF
!>    LDF (input) INTEGER, LDF >= M
!>    The leading dimension of the array F.
!>
[in,out]X
!>    X (workspace/output) REAL(KIND=WP) MIN(M,N)-by-(N-1) array
!>    X is used as workspace to hold representations of the
!>    leading N-1 snapshots in the orthonormal basis computed
!>    in the QR factorization of F.
!>    On exit, the leading K columns of X contain the leading
!>    K left singular vectors of the above described content
!>    of X. To lift them to the space of the left singular
!>    vectors U(:,1:K)of the input data, pre-multiply with the
!>    Q factor from the initial QR factorization.
!>    See the descriptions of F, K, V  and Z.
!>
[in]LDX
!>    LDX (input) INTEGER, LDX >= N
!>    The leading dimension of the array X.
!>
[in,out]Y
!>    Y (workspace/output) REAL(KIND=WP) MIN(M,N)-by-(N-1) array
!>    Y is used as workspace to hold representations of the
!>    trailing N-1 snapshots in the orthonormal basis computed
!>    in the QR factorization of F.
!>    On exit,
!>    If JOBT == 'R', Y contains the MIN(M,N)-by-N upper
!>    triangular factor from the QR factorization of the data
!>    snapshot matrix F.
!>
[in]LDY
!>    LDY (input) INTEGER , LDY >= N
!>    The leading dimension of the array Y.
!>
[in]NRNK
!>    NRNK (input) INTEGER
!>    Determines the mode how to compute the numerical rank,
!>    i.e. how to truncate small singular values of the input
!>    matrix X. On input, if
!>    NRNK = -1 :: i-th singular value sigma(i) is truncated
!>                 if sigma(i) <= TOL*sigma(1)
!>                 This option is recommended.
!>    NRNK = -2 :: i-th singular value sigma(i) is truncated
!>                 if sigma(i) <= TOL*sigma(i-1)
!>                 This option is included for R&D purposes.
!>                 It requires highly accurate SVD, which
!>                 may not be feasible.
!>    The numerical rank can be enforced by using positive
!>    value of NRNK as follows:
!>    0 < NRNK <= N-1 :: at most NRNK largest singular values
!>    will be used. If the number of the computed nonzero
!>    singular values is less than NRNK, then only those
!>    nonzero values will be used and the actually used
!>    dimension is less than NRNK. The actual number of
!>    the nonzero singular values is returned in the variable
!>    K. See the description of K.
!>
[in]TOL
!>     TOL (input) REAL(KIND=WP), 0 <= TOL < 1
!>     The tolerance for truncating small singular values.
!>     See the description of NRNK.
!>
[out]K
!>     K (output) INTEGER,  0 <= K <= N
!>     The dimension of the SVD/POD basis for the leading N-1
!>     data snapshots (columns of F) and the number of the
!>     computed Ritz pairs. The value of K is determined
!>     according to the rule set by the parameters NRNK and
!>     TOL. See the descriptions of NRNK and TOL.
!>
[out]REIG
!>    REIG (output) REAL(KIND=WP) (N-1)-by-1 array
!>    The leading K (K<=N) entries of REIG contain
!>    the real parts of the computed eigenvalues
!>    REIG(1:K) + sqrt(-1)*IMEIG(1:K).
!>    See the descriptions of K, IMEIG, Z.
!>
[out]IMEIG
!>    IMEIG (output) REAL(KIND=WP) (N-1)-by-1 array
!>    The leading K (K<N) entries of REIG contain
!>    the imaginary parts of the computed eigenvalues
!>    REIG(1:K) + sqrt(-1)*IMEIG(1:K).
!>    The eigenvalues are determined as follows:
!>    If IMEIG(i) == 0, then the corresponding eigenvalue is
!>    real, LAMBDA(i) = REIG(i).
!>    If IMEIG(i)>0, then the corresponding complex
!>    conjugate pair of eigenvalues reads
!>    LAMBDA(i)   = REIG(i) + sqrt(-1)*IMAG(i)
!>    LAMBDA(i+1) = REIG(i) - sqrt(-1)*IMAG(i)
!>    That is, complex conjugate pairs have consequtive
!>    indices (i,i+1), with the positive imaginary part
!>    listed first.
!>    See the descriptions of K, REIG, Z.
!>
[out]Z
!>    Z (workspace/output) REAL(KIND=WP)  M-by-(N-1) array
!>    If JOBZ =='V' then
!>       Z contains real Ritz vectors as follows:
!>       If IMEIG(i)=0, then Z(:,i) is an eigenvector of
!>       the i-th Ritz value.
!>       If IMEIG(i) > 0 (and IMEIG(i+1) < 0) then
!>       [Z(:,i) Z(:,i+1)] span an invariant subspace and
!>       the Ritz values extracted from this subspace are
!>       REIG(i) + sqrt(-1)*IMEIG(i) and
!>       REIG(i) - sqrt(-1)*IMEIG(i).
!>       The corresponding eigenvectors are
!>       Z(:,i) + sqrt(-1)*Z(:,i+1) and
!>       Z(:,i) - sqrt(-1)*Z(:,i+1), respectively.
!>    If JOBZ == 'F', then the above descriptions hold for
!>    the columns of Z*V, where the columns of V are the
!>    eigenvectors of the K-by-K Rayleigh quotient, and Z is
!>    orthonormal. The columns of V are similarly structured:
!>    If IMEIG(i) == 0 then Z*V(:,i) is an eigenvector, and if
!>    IMEIG(i) > 0 then Z*V(:,i)+sqrt(-1)*Z*V(:,i+1) and
!>                      Z*V(:,i)-sqrt(-1)*Z*V(:,i+1)
!>    are the eigenvectors of LAMBDA(i), LAMBDA(i+1).
!>    See the descriptions of REIG, IMEIG, X and V.
!>
[in]LDZ
!>    LDZ (input) INTEGER , LDZ >= M
!>    The leading dimension of the array Z.
!>
[out]RES
!>    RES (output) REAL(KIND=WP) (N-1)-by-1 array
!>    RES(1:K) contains the residuals for the K computed
!>    Ritz pairs.
!>    If LAMBDA(i) is real, then
!>       RES(i) = || A * Z(:,i) - LAMBDA(i)*Z(:,i))||_2.
!>    If [LAMBDA(i), LAMBDA(i+1)] is a complex conjugate pair
!>    then
!>    RES(i)=RES(i+1) = || A * Z(:,i:i+1) - Z(:,i:i+1) *B||_F
!>    where B = [ real(LAMBDA(i)) imag(LAMBDA(i)) ]
!>              [-imag(LAMBDA(i)) real(LAMBDA(i)) ].
!>    It holds that
!>    RES(i)   = || A*ZC(:,i)   - LAMBDA(i)  *ZC(:,i)   ||_2
!>    RES(i+1) = || A*ZC(:,i+1) - LAMBDA(i+1)*ZC(:,i+1) ||_2
!>    where ZC(:,i)   =  Z(:,i) + sqrt(-1)*Z(:,i+1)
!>          ZC(:,i+1) =  Z(:,i) - sqrt(-1)*Z(:,i+1)
!>    See the description of Z.
!>
[out]B
!>    B (output) REAL(KIND=WP)  MIN(M,N)-by-(N-1) array.
!>    IF JOBF =='R', B(1:N,1:K) contains A*U(:,1:K), and can
!>    be used for computing the refined vectors; see further
!>    details in the provided references.
!>    If JOBF == 'E', B(1:N,1;K) contains
!>    A*U(:,1:K)*W(1:K,1:K), which are the vectors from the
!>    Exact DMD, up to scaling by the inverse eigenvalues.
!>    In both cases, the content of B can be lifted to the
!>    original dimension of the input data by pre-multiplying
!>    with the Q factor from the initial QR factorization.
!>    Here A denotes a compression of the underlying operator.
!>    See the descriptions of F and X.
!>    If JOBF =='N', then B is not referenced.
!>
[in]LDB
!>    LDB (input) INTEGER, LDB >= MIN(M,N)
!>    The leading dimension of the array B.
!>
[out]V
!>    V (workspace/output) REAL(KIND=WP) (N-1)-by-(N-1) array
!>    On exit, V(1:K,1:K) contains the K eigenvectors of
!>    the Rayleigh quotient. The eigenvectors of a complex
!>    conjugate pair of eigenvalues are returned in real form
!>    as explained in the description of Z. The Ritz vectors
!>    (returned in Z) are the product of X and V; see
!>    the descriptions of X and Z.
!>
[in]LDV
!>    LDV (input) INTEGER, LDV >= N-1
!>    The leading dimension of the array V.
!>
[out]S
!>    S (output) REAL(KIND=WP) (N-1)-by-(N-1) array
!>    The array S(1:K,1:K) is used for the matrix Rayleigh
!>    quotient. This content is overwritten during
!>    the eigenvalue decomposition by DGEEV.
!>    See the description of K.
!>
[in]LDS
!>    LDS (input) INTEGER, LDS >= N-1
!>    The leading dimension of the array S.
!>
[out]WORK
!>    WORK (workspace/output) REAL(KIND=WP) LWORK-by-1 array
!>    On exit,
!>    WORK(1:MIN(M,N)) contains the scalar factors of the
!>    elementary reflectors as returned by DGEQRF of the
!>    M-by-N input matrix F.
!>    WORK(MIN(M,N)+1:MIN(M,N)+N-1) contains the singular values of
!>    the input submatrix F(1:M,1:N-1).
!>    If the call to DGEDMDQ is only workspace query, then
!>    WORK(1) contains the minimal workspace length and
!>    WORK(2) is the optimal workspace length. Hence, the
!>    length of work is at least 2.
!>    See the description of LWORK.
!>
[in]LWORK
!>    LWORK (input) INTEGER
!>    The minimal length of the  workspace vector WORK.
!>    LWORK is calculated as follows:
!>    Let MLWQR  = N (minimal workspace for DGEQRF[M,N])
!>        MLWDMD = minimal workspace for DGEDMD (see the
!>                 description of LWORK in DGEDMD) for
!>                 snapshots of dimensions MIN(M,N)-by-(N-1)
!>        MLWMQR = N (minimal workspace for
!>                   DORMQR['L','N',M,N,N])
!>        MLWGQR = N (minimal workspace for DORGQR[M,N,N])
!>    Then
!>    LWORK = MAX(N+MLWQR, N+MLWDMD)
!>    is updated as follows:
!>       if   JOBZ == 'V' or JOBZ == 'F' THEN
!>            LWORK = MAX( LWORK, MIN(M,N)+N-1+MLWMQR )
!>       if   JOBQ == 'Q' THEN
!>            LWORK = MAX( LWORK, MIN(M,N)+N-1+MLWGQR)
!>    If on entry LWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace lengths for both WORK and
!>    IWORK. See the descriptions of WORK and IWORK.
!>
[out]IWORK
!>    IWORK (workspace/output) INTEGER LIWORK-by-1 array
!>    Workspace that is required only if WHTSVD equals
!>    2 , 3 or 4. (See the description of WHTSVD).
!>    If on entry LWORK =-1 or LIWORK=-1, then the
!>    minimal length of IWORK is computed and returned in
!>    IWORK(1). See the description of LIWORK.
!>
[in]LIWORK
!>    LIWORK (input) INTEGER
!>    The minimal length of the workspace vector IWORK.
!>    If WHTSVD == 1, then only IWORK(1) is used; LIWORK >=1
!>    Let M1=MIN(M,N), N1=N-1. Then
!>    If WHTSVD == 2, then LIWORK >= MAX(1,8*MIN(M1,N1))
!>    If WHTSVD == 3, then LIWORK >= MAX(1,M1+N1-1)
!>    If WHTSVD == 4, then LIWORK >= MAX(3,M1+3*N1)
!>    If on entry LIWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace lengths for both WORK and
!>    IWORK. See the descriptions of WORK and IWORK.
!>
[out]INFO
!>    INFO (output) INTEGER
!>    -i < 0 :: On entry, the i-th argument had an
!>              illegal value
!>       = 0 :: Successful return.
!>       = 1 :: Void input. Quick exit (M=0 or N=0).
!>       = 2 :: The SVD computation of X did not converge.
!>              Suggestion: Check the input data and/or
!>              repeat with different WHTSVD.
!>       = 3 :: The computation of the eigenvalues did not
!>              converge.
!>       = 4 :: If data scaling was requested on input and
!>              the procedure found inconsistency in the data
!>              such that for some column index i,
!>              X(:,i) = 0 but Y(:,i) /= 0, then Y(:,i) is set
!>              to zero if JOBS=='C'. The computation proceeds
!>              with original or modified data and warning
!>              flag is set with INFO=4.
!>
Author
Zlatko Drmac