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

◆ zgedmdq()

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

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

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
[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
!>    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.
!>
[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 :: 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.
!>
[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) 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.
!>
[in]LDF
!>    LDF (input) INTEGER, LDF >= M
!>    The leading dimension of the array F.
!>
[in,out]X
!>    X (workspace/output) COMPLEX(KIND=WP) MIN(M,N)-by-(N-1) array
!>    X is used as workspace to hold representations of the
!>    leading N-1 snapshots in the orthonormal basis computed
!>    in the QR factorization of F.
!>    On exit, the leading K columns of X contain the leading
!>    K left singular vectors of the above described content
!>    of X. To lift them to the space of the left singular
!>    vectors U(:,1:K) of the input data, pre-multiply with the
!>    Q factor from the initial QR factorization.
!>    See the descriptions of F, K, V  and Z.
!>
[in]LDX
!>    LDX (input) INTEGER, LDX >= N
!>    The leading dimension of the array X.
!>
[in,out]Y
!>    Y (workspace/output) COMPLEX(KIND=WP) MIN(M,N)-by-(N) array
!>    Y is used as workspace to hold representations of the
!>    trailing N-1 snapshots in the orthonormal basis computed
!>    in the QR factorization of F.
!>    On exit,
!>    If JOBT == 'R', Y contains the MIN(M,N)-by-N upper
!>    triangular factor from the QR factorization of the data
!>    snapshot matrix F.
!>
[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]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.
!>
[out]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.
!>
[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,
!>    RES(i) = || A * Z(:,i) - EIGS(i)*Z(:,i))||_2.
!>    See the description of EIGS and Z.
!>
[out]B
!>    B (output) COMPLEX(KIND=WP)  MIN(M,N)-by-(N-1) array.
!>    IF JOBF =='R', B(1:N,1:K) contains A*U(:,1:K), and can
!>    be used for computing the refined vectors; see further
!>    details in the provided references.
!>    If JOBF == 'E', B(1:N,1;K) contains
!>    A*U(:,1:K)*W(1:K,1:K), which are the vectors from the
!>    Exact DMD, up to scaling by the inverse eigenvalues.
!>    In both cases, the content of B can be lifted to the
!>    original dimension of the input data by pre-multiplying
!>    with the Q factor from the initial QR factorization.
!>    Here A denotes a compression of the underlying operator.
!>    See the descriptions of F and X.
!>    If JOBF =='N', then B is not referenced.
!>
[in]LDB
!>    LDB (input) INTEGER, LDB >= MIN(M,N)
!>    The leading dimension of the array B.
!>
[out]V
!>    V (workspace/output) COMPLEX(KIND=WP) (N-1)-by-(N-1) array
!>    On exit, V(1:K,1:K) V contains the K eigenvectors of
!>    the Rayleigh quotient. The Ritz vectors
!>    (returned in Z) are the product of Q from the initial QR
!>    factorization (see the description of F) X (see the
!>    description of X) and V.
!>
[in]LDV
!>    LDV (input) INTEGER, LDV >= N-1
!>    The leading dimension of the array V.
!>
[out]S
!>    S (output) COMPLEX(KIND=WP) (N-1)-by-(N-1) array
!>    The array S(1:K,1:K) is used for the matrix Rayleigh
!>    quotient. This content is overwritten during
!>    the eigenvalue decomposition by ZGEEV.
!>    See the description of K.
!>
[in]LDS
!>    LDS (input) INTEGER, LDS >= N-1
!>    The leading dimension of the array S.
!>
[out]ZWORK
!>    ZWORK (workspace/output) COMPLEX(KIND=WP) LWORK-by-1 array
!>    On exit,
!>    ZWORK(1:MIN(M,N)) contains the scalar factors of the
!>    elementary reflectors as returned by ZGEQRF of the
!>    M-by-N input matrix F.
!>    If the call to ZGEDMDQ is only workspace query, then
!>    ZWORK(1) contains the minimal complex workspace length and
!>    ZWORK(2) is the optimal complex workspace length.
!>    Hence, the length of work is at least 2.
!>    See the description of LZWORK.
!>
[in]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)
!>
[out]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.
!>
[in]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.
!>
[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

Definition at line 551 of file zgedmdq.f90.

557!
558! -- LAPACK driver routine --
559!
560! -- LAPACK is a software package provided by University of --
561! -- Tennessee, University of California Berkeley, University of --
562! -- Colorado Denver and NAG Ltd.. --
563!
564!.....
565 use, INTRINSIC :: iso_fortran_env, only: real64
566 IMPLICIT NONE
567 INTEGER, PARAMETER :: WP = real64
568!
569! Scalar arguments
570! ~~~~~~~~~~~~~~~~
571 CHARACTER, INTENT(IN) :: JOBS, JOBZ, JOBR, JOBQ, &
572 jobt, jobf
573 INTEGER, INTENT(IN) :: WHTSVD, M, N, LDF, LDX, &
574 ldy, nrnk, ldz, ldb, ldv, &
575 lds, lzwork, lwork, liwork
576 INTEGER, INTENT(OUT) :: INFO, K
577 REAL(KIND=wp), INTENT(IN) :: tol
578!
579! Array arguments
580! ~~~~~~~~~~~~~~~
581 COMPLEX(KIND=WP), INTENT(INOUT) :: F(LDF,*)
582 COMPLEX(KIND=WP), INTENT(OUT) :: X(LDX,*), Y(LDY,*), &
583 z(ldz,*), b(ldb,*), &
584 v(ldv,*), s(lds,*)
585 COMPLEX(KIND=WP), INTENT(OUT) :: EIGS(*)
586 COMPLEX(KIND=WP), INTENT(OUT) :: ZWORK(*)
587 REAL(KIND=wp), INTENT(OUT) :: res(*)
588 REAL(KIND=wp), INTENT(OUT) :: work(*)
589 INTEGER, INTENT(OUT) :: IWORK(*)
590!
591! Parameters
592! ~~~~~~~~~~
593 REAL(KIND=wp), PARAMETER :: one = 1.0_wp
594 REAL(KIND=wp), PARAMETER :: zero = 0.0_wp
595! COMPLEX(KIND=WP), PARAMETER :: ZONE = ( 1.0_WP, 0.0_WP )
596 COMPLEX(KIND=WP), PARAMETER :: ZZERO = ( 0.0_wp, 0.0_wp )
597!
598! Local scalars
599! ~~~~~~~~~~~~~
600 INTEGER :: IMINWR, INFO1, MINMN, MLRWRK, &
601 mlwdmd, mlwgqr, mlwmqr, mlwork, &
602 mlwqr, olwdmd, olwgqr, olwmqr, &
603 olwork, olwqr
604 LOGICAL :: LQUERY, SCCOLX, SCCOLY, WANTQ, &
605 wnttrf, wntres, wntvec, wntvcf, &
606 wntvcq, wntref, wntex
607 CHARACTER(LEN=1) :: JOBVL
608!
609! External functions (BLAS and LAPACK)
610! ~~~~~~~~~~~~~~~~~
611 LOGICAL LSAME
612 EXTERNAL lsame
613!
614! External subroutines (BLAS and LAPACK)
615! ~~~~~~~~~~~~~~~~~~~~
616 EXTERNAL zgedmd, zgeqrf, zlacpy, zlaset, zungqr, &
617 zunmqr, xerbla
618!
619! Intrinsic functions
620! ~~~~~~~~~~~~~~~~~~~
621 INTRINSIC max, min, int
622!..........................................................
623!
624! Test the input arguments
625 wntres = lsame(jobr,'R')
626 sccolx = lsame(jobs,'S') .OR. lsame( jobs, 'C' )
627 sccoly = lsame(jobs,'Y')
628 wntvec = lsame(jobz,'V')
629 wntvcf = lsame(jobz,'F')
630 wntvcq = lsame(jobz,'Q')
631 wntref = lsame(jobf,'R')
632 wntex = lsame(jobf,'E')
633 wantq = lsame(jobq,'Q')
634 wnttrf = lsame(jobt,'R')
635 minmn = min(m,n)
636 info = 0
637 lquery = ( (lzwork == -1) .OR. (lwork == -1) .OR. (liwork == -1) )
638!
639 IF ( .NOT. (sccolx .OR. sccoly .OR. &
640 lsame(jobs,'N')) ) THEN
641 info = -1
642 ELSE IF ( .NOT. (wntvec .OR. wntvcf .OR. wntvcq &
643 .OR. lsame(jobz,'N')) ) THEN
644 info = -2
645 ELSE IF ( .NOT. (wntres .OR. lsame(jobr,'N')) .OR. &
646 ( wntres .AND. lsame(jobz,'N') ) ) THEN
647 info = -3
648 ELSE IF ( .NOT. (wantq .OR. lsame(jobq,'N')) ) THEN
649 info = -4
650 ELSE IF ( .NOT. ( wnttrf .OR. lsame(jobt,'N') ) ) THEN
651 info = -5
652 ELSE IF ( .NOT. (wntref .OR. wntex .OR. &
653 lsame(jobf,'N') ) ) THEN
654 info = -6
655 ELSE IF ( .NOT. ((whtsvd == 1).OR.(whtsvd == 2).OR. &
656 (whtsvd == 3).OR.(whtsvd == 4)) ) THEN
657 info = -7
658 ELSE IF ( m < 0 ) THEN
659 info = -8
660 ELSE IF ( ( n < 0 ) .OR. ( n > m+1 ) ) THEN
661 info = -9
662 ELSE IF ( ldf < m ) THEN
663 info = -11
664 ELSE IF ( ldx < minmn ) THEN
665 info = -13
666 ELSE IF ( ldy < minmn ) THEN
667 info = -15
668 ELSE IF ( .NOT. (( nrnk == -2).OR.(nrnk == -1).OR. &
669 ((nrnk >= 1).AND.(nrnk <=n ))) ) THEN
670 info = -16
671 ELSE IF ( ( tol < zero ) .OR. ( tol >= one ) ) THEN
672 info = -17
673 ELSE IF ( ldz < m ) THEN
674 info = -21
675 ELSE IF ( (wntref.OR.wntex ).AND.( ldb < minmn ) ) THEN
676 info = -24
677 ELSE IF ( ldv < n-1 ) THEN
678 info = -26
679 ELSE IF ( lds < n-1 ) THEN
680 info = -28
681 END IF
682!
683 IF ( wntvec .OR. wntvcf .OR. wntvcq ) THEN
684 jobvl = 'V'
685 ELSE
686 jobvl = 'N'
687 END IF
688 IF ( info == 0 ) THEN
689 ! Compute the minimal and the optimal workspace
690 ! requirements. Simulate running the code and
691 ! determine minimal and optimal sizes of the
692 ! workspace at any moment of the run.
693 IF ( ( n == 0 ) .OR. ( n == 1 ) ) THEN
694 ! All output except K is void. INFO=1 signals
695 ! the void input. In case of a workspace query,
696 ! the minimal workspace lengths are returned.
697 IF ( lquery ) THEN
698 iwork(1) = 1
699 zwork(1) = 2
700 zwork(2) = 2
701 work(1) = 2
702 work(2) = 2
703 ELSE
704 k = 0
705 END IF
706 info = 1
707 RETURN
708 END IF
709
710 mlrwrk = 2
711 mlwork = 2
712 olwork = 2
713 iminwr = 1
714 mlwqr = max(1,n) ! Minimal workspace length for ZGEQRF.
715 mlwork = max(mlwork,minmn + mlwqr)
716
717 IF ( lquery ) THEN
718 CALL zgeqrf( m, n, f, ldf, zwork, zwork, -1, &
719 info1 )
720 olwqr = int(zwork(1))
721 olwork = max(olwork,minmn + olwqr)
722 END IF
723 CALL zgedmd( jobs, jobvl, jobr, jobf, whtsvd, minmn,&
724 n-1, x, ldx, y, ldy, nrnk, tol, k, &
725 eigs, z, ldz, res, b, ldb, v, ldv, &
726 s, lds, zwork, -1, work, -1, iwork,&
727 -1, info1 )
728 mlwdmd = int(zwork(1))
729 mlwork = max(mlwork, minmn + mlwdmd)
730 mlrwrk = max(mlrwrk, int(work(1)))
731 iminwr = max(iminwr, iwork(1))
732 IF ( lquery ) THEN
733 olwdmd = int(zwork(2))
734 olwork = max(olwork, minmn+olwdmd)
735 END IF
736 IF ( wntvec .OR. wntvcf ) THEN
737 mlwmqr = max(1,n)
738 mlwork = max(mlwork,minmn+mlwmqr)
739 IF ( lquery ) THEN
740 CALL zunmqr( 'L','N', m, n, minmn, f, ldf, &
741 zwork, z, ldz, zwork, -1, info1 )
742 olwmqr = int(zwork(1))
743 olwork = max(olwork,minmn+olwmqr)
744 END IF
745 END IF
746 IF ( wantq ) THEN
747 mlwgqr = max(1,n)
748 mlwork = max(mlwork,minmn+mlwgqr)
749 IF ( lquery ) THEN
750 CALL zungqr( m, minmn, minmn, f, ldf, zwork, &
751 zwork, -1, info1 )
752 olwgqr = int(zwork(1))
753 olwork = max(olwork,minmn+olwgqr)
754 END IF
755 END IF
756 IF ( liwork < iminwr .AND. (.NOT.lquery) ) info = -34
757 IF ( lwork < mlrwrk .AND. (.NOT.lquery) ) info = -32
758 IF ( lzwork < mlwork .AND. (.NOT.lquery) ) info = -30
759 END IF
760 IF( info /= 0 ) THEN
761 CALL xerbla( 'ZGEDMDQ', -info )
762 RETURN
763 ELSE IF ( lquery ) THEN
764! Return minimal and optimal workspace sizes
765 iwork(1) = iminwr
766 zwork(1) = mlwork
767 zwork(2) = olwork
768 work(1) = mlrwrk
769 work(2) = mlrwrk
770 RETURN
771 END IF
772!.....
773! Initial QR factorization that is used to represent the
774! snapshots as elements of lower dimensional subspace.
775! For large scale computation with M >> N, at this place
776! one can use an out of core QRF.
777!
778 CALL zgeqrf( m, n, f, ldf, zwork, &
779 zwork(minmn+1), lzwork-minmn, info1 )
780!
781! Define X and Y as the snapshots representations in the
782! orthogonal basis computed in the QR factorization.
783! X corresponds to the leading N-1 and Y to the trailing
784! N-1 snapshots.
785 CALL zlaset( 'L', minmn, n-1, zzero, zzero, x, ldx )
786 CALL zlacpy( 'U', minmn, n-1, f, ldf, x, ldx )
787 CALL zlacpy( 'A', minmn, n-1, f(1,2), ldf, y, ldy )
788 IF ( m >= 3 ) THEN
789 CALL zlaset( 'L', minmn-2, n-2, zzero, zzero, &
790 y(3,1), ldy )
791 END IF
792!
793! Compute the DMD of the projected snapshot pairs (X,Y)
794 CALL zgedmd( jobs, jobvl, jobr, jobf, whtsvd, minmn, &
795 n-1, x, ldx, y, ldy, nrnk, tol, k, &
796 eigs, z, ldz, res, b, ldb, v, ldv, &
797 s, lds, zwork(minmn+1), lzwork-minmn, &
798 work, lwork, iwork, liwork, info1 )
799 IF ( info1 == 2 .OR. info1 == 3 ) THEN
800 ! Return with error code. See ZGEDMD for details.
801 info = info1
802 RETURN
803 ELSE
804 info = info1
805 END IF
806!
807! The Ritz vectors (Koopman modes) can be explicitly
808! formed or returned in factored form.
809 IF ( wntvec ) THEN
810 ! Compute the eigenvectors explicitly.
811 IF ( m > minmn ) CALL zlaset( 'A', m-minmn, k, zzero, &
812 zzero, z(minmn+1,1), ldz )
813 CALL zunmqr( 'L','N', m, k, minmn, f, ldf, zwork, z, &
814 ldz, zwork(minmn+1), lzwork-minmn, info1 )
815 ELSE IF ( wntvcf ) THEN
816 ! Return the Ritz vectors (eigenvectors) in factored
817 ! form Z*V, where Z contains orthonormal matrix (the
818 ! product of Q from the initial QR factorization and
819 ! the SVD/POD_basis returned by ZGEDMD in X) and the
820 ! second factor (the eigenvectors of the Rayleigh
821 ! quotient) is in the array V, as returned by ZGEDMD.
822 CALL zlacpy( 'A', n, k, x, ldx, z, ldz )
823 IF ( m > n ) CALL zlaset( 'A', m-n, k, zzero, zzero, &
824 z(n+1,1), ldz )
825 CALL zunmqr( 'L','N', m, k, minmn, f, ldf, zwork, z, &
826 ldz, zwork(minmn+1), lzwork-minmn, info1 )
827 END IF
828!
829! Some optional output variables:
830!
831! The upper triangular factor R in the initial QR
832! factorization is optionally returned in the array Y.
833! This is useful if this call to ZGEDMDQ is to be
834! followed by a streaming DMD that is implemented in a
835! QR compressed form.
836 IF ( wnttrf ) THEN ! Return the upper triangular R in Y
837 CALL zlaset( 'A', minmn, n, zzero, zzero, y, ldy )
838 CALL zlacpy( 'U', minmn, n, f, ldf, y, ldy )
839 END IF
840!
841! The orthonormal/unitary factor Q in the initial QR
842! factorization is optionally returned in the array F.
843! Same as with the triangular factor above, this is
844! useful in a streaming DMD.
845 IF ( wantq ) THEN ! Q overwrites F
846 CALL zungqr( m, minmn, minmn, f, ldf, zwork, &
847 zwork(minmn+1), lzwork-minmn, info1 )
848 END IF
849!
850 RETURN
851!
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zgedmd
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.
Definition zgedmd.f90:501
zgeqrf
subroutine zgeqrf(m, n, a, lda, tau, work, lwork, info)
ZGEQRF
Definition zgeqrf.f:144
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:101
zlaset
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:104
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
zungqr
subroutine zungqr(m, n, k, a, lda, tau, work, lwork, info)
ZUNGQR
Definition zungqr.f:126
zunmqr
subroutine zunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMQR
Definition zunmqr.f:165
Here is the call graph for this function: