LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
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 "Koopman Operator-Based Forecasting
for Nonstationary Processes from Near-Term, Limited
Observational Data" Contract No: W31P4Q-21-C-0007
- DARPA PAI project "Physics-Informed Machine Learning
Methodologies" Contract No: HR0011-18-9-0033
- DARPA MoDyL project "A Data-Driven, Operator-Theoretic
Framework for Space-Time Analysis of Process Dynamics"
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] 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. ``` [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. ```

Definition at line 550 of file zgedmdq.f90.

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