OSchip
/
OpenBLAS
Not watched
1
0
Fork 0
Code
Releases 66 Wiki evaluate Activity
Issues 0 Pull Requests 0 Datasets Model Cloudbrain HPC
Browse Source

Merge pull request #4392 from martin-frbg/lapack959 

Fix issues related to the ?GEDMD functions (Reference-LAPACK PR 959)
tags/v0.3.26
Martin Kroeker GitHub 2 years ago
parent
05dde8ef04 c6fa921027
commit
dba404055d
No known key found for this signature in database GPG Key ID: 4AEE18F83AFDEB23
6 changed files with 2108 additions and 1472 deletions
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. +514
    -358
      lapack-netlib/SRC/cgedmd.f90
  2. +536
    -384
      lapack-netlib/SRC/dgedmd.f90
  3. +534
    -382
      lapack-netlib/SRC/sgedmd.f90
  4. +500
    -348
      lapack-netlib/SRC/zgedmd.f90
  5. +12
    -0
      lapack-netlib/TESTING/CMakeLists.txt
  6. +12
    -0
      lapack-netlib/TESTING/EIG/CMakeLists.txt

+ 514
- 358
lapack-netlib/SRC/cgedmd.f90 View File

@@ -1,22 +1,526 @@
!> \brief \b CGEDMD computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.
!
! =========== DOCUMENTATION ===========
!
! Definition:
! ===========
!
! 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 )
!.....
! USE iso_fortran_env
! IMPLICIT NONE
! INTEGER, PARAMETER :: WP = real32
!
!.....
! Scalar arguments
! CHARACTER, INTENT(IN) :: JOBS, JOBZ, JOBR, JOBF
! INTEGER, INTENT(IN) :: WHTSVD, M, N, LDX, LDY, &
! NRNK, LDZ, LDB, LDW, LDS, &
! LIWORK, LRWORK, LZWORK
! INTEGER, INTENT(OUT) :: K, INFO
! REAL(KIND=WP), INTENT(IN) :: TOL
! Array arguments
! COMPLEX(KIND=WP), INTENT(INOUT) :: X(LDX,*), Y(LDY,*)
! COMPLEX(KIND=WP), INTENT(OUT) :: Z(LDZ,*), B(LDB,*), &
! W(LDW,*), S(LDS,*)
! COMPLEX(KIND=WP), INTENT(OUT) :: EIGS(*)
! COMPLEX(KIND=WP), INTENT(OUT) :: ZWORK(*)
! REAL(KIND=WP), INTENT(OUT) :: RES(*)
! REAL(KIND=WP), INTENT(OUT) :: RWORK(*)
! INTEGER, INTENT(OUT) :: IWORK(*)
!
!............................................................
!> \par Purpose:
! =============
!> \verbatim
!> 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].
!> \endverbatim
!............................................................
!> \par References:
! ================
!> \verbatim
!> [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.
!> \endverbatim
!......................................................................
!> \par Developed and supported by:
! ================================
!> \verbatim
!> 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
!> \endverbatim
!......................................................................
!> \par Distribution Statement A:
! ==============================
!> \verbatim
!> Approved for Public Release, Distribution Unlimited.
!> Cleared by DARPA on September 29, 2022
!> \endverbatim
!......................................................................
! Arguments
! =========
!
!> \param[in] JOBS
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[in] JOBZ
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[in] JOBR
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[in] JOBF
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[in] WHTSVD
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[in] M
!> \verbatim
!> M (input) INTEGER, M>= 0
!> The state space dimension (the row dimension of X, Y).
!> \endverbatim
!.....
!> \param[in] N
!> \verbatim
!> N (input) INTEGER, 0 <= N <= M
!> The number of data snapshot pairs
!> (the number of columns of X and Y).
!> \endverbatim
!.....
!> \param[in,out] X
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[in] LDX
!> \verbatim
!> LDX (input) INTEGER, LDX >= M
!> The leading dimension of the array X.
!> \endverbatim
!.....
!> \param[in,out] Y
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[in] LDY
!> \verbatim
!> LDY (input) INTEGER , LDY >= M
!> The leading dimension of the array Y.
!> \endverbatim
!.....
!> \param[in] NRNK
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[in] TOL
!> \verbatim
!> TOL (input) REAL(KIND=WP), 0 <= TOL < 1
!> The tolerance for truncating small singular values.
!> See the description of NRNK.
!> \endverbatim
!.....
!> \param[out] K
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[out] EIGS
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[out] Z
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[in] LDZ
!> \verbatim
!> LDZ (input) INTEGER , LDZ >= M
!> The leading dimension of the array Z.
!> \endverbatim
!.....
!> \param[out] RES
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[out] B
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[in] LDB
!> \verbatim
!> LDB (input) INTEGER, LDB >= M
!> The leading dimension of the array B.
!> \endverbatim
!.....
!> \param[out] W
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[in] LDW
!> \verbatim
!> LDW (input) INTEGER, LDW >= N
!> The leading dimension of the array W.
!> \endverbatim
!.....
!> \param[out] S
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[in] LDS
!> \verbatim
!> LDS (input) INTEGER, LDS >= N
!> The leading dimension of the array S.
!> \endverbatim
!.....
!> \param[out] ZWORK
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[in] LZWORK
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[out] RWORK
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[in] LRWORK
!> \verbatim
!> 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).
!> \endverbatim
!.....
!> \param[out] IWORK
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[in] LIWORK
!> \verbatim
!> 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.
!> \endverbatim
!.....
!> \param[out] INFO
!> \verbatim
!> 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.
!> \endverbatim
!
! Authors:
! ========
!
!> \author Zlatko Drmac
!
!> \ingroup gedmd
!
!.............................................................
!.............................................................
SUBROUTINE CGEDMD( JOBS, JOBZ, JOBR, JOBF, WHTSVD, & SUBROUTINE CGEDMD( JOBS, JOBZ, JOBR, JOBF, WHTSVD, &
M, N, X, LDX, Y, LDY, NRNK, TOL, & M, N, X, LDX, Y, LDY, NRNK, TOL, &
K, EIGS, Z, LDZ, RES, B, LDB, & K, EIGS, Z, LDZ, RES, B, LDB, &
W, LDW, S, LDS, ZWORK, LZWORK, & W, LDW, S, LDS, ZWORK, LZWORK, &
RWORK, LRWORK, IWORK, LIWORK, INFO ) RWORK, LRWORK, IWORK, LIWORK, INFO )
! March 2023
!
! -- LAPACK driver routine --
!
! -- LAPACK is a software package provided by University of --
! -- Tennessee, University of California Berkeley, University of --
! -- Colorado Denver and NAG Ltd.. --
!
!..... !.....
USE iso_fortran_env USE iso_fortran_env
IMPLICIT NONE IMPLICIT NONE
INTEGER, PARAMETER :: WP = real32 INTEGER, PARAMETER :: WP = real32
!.....
!
! Scalar arguments ! Scalar arguments
! ~~~~~~~~~~~~~~~~
CHARACTER, INTENT(IN) :: JOBS, JOBZ, JOBR, JOBF CHARACTER, INTENT(IN) :: JOBS, JOBZ, JOBR, JOBF
INTEGER, INTENT(IN) :: WHTSVD, M, N, LDX, LDY, & INTEGER, INTENT(IN) :: WHTSVD, M, N, LDX, LDY, &
NRNK, LDZ, LDB, LDW, LDS, & NRNK, LDZ, LDB, LDW, LDS, &
LIWORK, LRWORK, LZWORK LIWORK, LRWORK, LZWORK
INTEGER, INTENT(OUT) :: K, INFO INTEGER, INTENT(OUT) :: K, INFO
REAL(KIND=WP), INTENT(IN) :: TOL REAL(KIND=WP), INTENT(IN) :: TOL
!
! Array arguments ! Array arguments
! ~~~~~~~~~~~~~~~
COMPLEX(KIND=WP), INTENT(INOUT) :: X(LDX,*), Y(LDY,*) COMPLEX(KIND=WP), INTENT(INOUT) :: X(LDX,*), Y(LDY,*)
COMPLEX(KIND=WP), INTENT(OUT) :: Z(LDZ,*), B(LDB,*), & COMPLEX(KIND=WP), INTENT(OUT) :: Z(LDZ,*), B(LDB,*), &
W(LDW,*), S(LDS,*) W(LDW,*), S(LDS,*)
@@ -25,364 +529,14 @@
REAL(KIND=WP), INTENT(OUT) :: RES(*) REAL(KIND=WP), INTENT(OUT) :: RES(*)
REAL(KIND=WP), INTENT(OUT) :: RWORK(*) REAL(KIND=WP), INTENT(OUT) :: RWORK(*)
INTEGER, INTENT(OUT) :: IWORK(*) INTEGER, INTENT(OUT) :: IWORK(*)
!............................................................
! 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 "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
!============================================================
! Distribution Statement A:
! Approved for Public Release, Distribution Unlimited.
! Cleared by DARPA on September 29, 2022
!============================================================
!......................................................................
! Arguments
! =========
! 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 (input) CHARACTER*1
! Determines whether the eigenvectors (Koopman modes) will
! be computed.
! 'V' :: The eigenvectors (Koopman modes) will be computed
! and returned in the matrix Z.
! See the description of Z.
! 'F' :: The eigenvectors (Koopman modes) will be returned
! in factored form as the product X(:,1:K)*W, where X
! contains a POD basis (leading left singular vectors
! of the data matrix X) and W contains the eigenvectors
! of the corresponding Rayleigh quotient.
! See the descriptions of K, X, W, Z.
! 'N' :: The eigenvectors are not computed.
!.....
! JOBR (input) CHARACTER*1
! Determines whether to compute the residuals.
! 'R' :: The residuals for the computed eigenpairs will be
! computed and stored in the array RES.
! See the description of RES.
! For this option to be legal, JOBZ must be 'V'.
! 'N' :: The residuals are not computed.
!.....
! JOBF (input) CHARACTER*1
! Specifies whether to store information needed for post-
! processing (e.g. computing refined Ritz vectors)
! 'R' :: The matrix needed for the refinement of the Ritz
! vectors is computed and stored in the array B.
! See the description of B.
! 'E' :: The unscaled eigenvectors of the Exact DMD are
! computed and returned in the array B. See the
! description of B.
! 'N' :: No eigenvector refinement data is computed.
!.....
! WHTSVD (input) INTEGER, WHSTVD in { 1, 2, 3, 4 }
! Allows for a selection of the SVD algorithm from the
! LAPACK library.
! 1 :: 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 (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 (input) INTEGER , LDY >= M
! The leading dimension of the array Y.
!.....
! NRNK (input) INTEGER
! Determines the mode how to compute the numerical rank,
! i.e. how to truncate small singular values of the input
! matrix X. On input, if
! NRNK = -1 :: i-th singular value sigma(i) is truncated
! if sigma(i) <= TOL*sigma(1)
! This option is recommended.
! NRNK = -2 :: i-th singular value sigma(i) is truncated
! if sigma(i) <= TOL*sigma(i-1)
! This option is included for R&D purposes.
! It requires highly accurate SVD, which
! may not be feasible.
! The numerical rank can be enforced by using positive
! value of NRNK as follows:
! 0 < NRNK <= N :: at most NRNK largest singular values
! will be used. If the number of the computed nonzero
! singular values is less than NRNK, then only those
! nonzero values will be used and the actually used
! dimension is less than NRNK. The actual number of
! the nonzero singular values is returned in the variable
! K. See the descriptions of TOL and K.
!.....
! TOL (input) REAL(KIND=WP), 0 <= TOL < 1
! The tolerance for truncating small singular values.
! See the description of NRNK.
!.....
! K (output) INTEGER, 0 <= K <= N
! The dimension of the POD basis for the data snapshot
! matrix X and the number of the computed Ritz pairs.
! The value of K is determined according to the rule set
! by the parameters NRNK and TOL.
! See the descriptions of NRNK and TOL.
!.....
! 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 (input) INTEGER , LDZ >= M
! The leading dimension of the array Z.
!.....
! RES (output) REAL(KIND=WP) N-by-1 array
! RES(1:K) contains the residuals for the K computed
! Ritz pairs,
! 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 (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 (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 (input) INTEGER, LDS >= N
! The leading dimension of the array S.
!.....
! 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 (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 (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 (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 (workspace/output) INTEGER LIWORK-by-1 array
! Workspace that is required only if WHTSVD equals
! 2 , 3 or 4. (See the description of WHTSVD).
! If on entry LWORK =-1 or LIWORK=-1, then the
! minimal length of IWORK is computed and returned in
! IWORK(1). See the description of LIWORK.
!.....
! LIWORK (input) INTEGER
! The minimal length of the workspace vector IWORK.
! If WHTSVD == 1, then only IWORK(1) is used; LIWORK >=1
! If WHTSVD == 2, then LIWORK >= MAX(1,8*MIN(M,N))
! If WHTSVD == 3, then LIWORK >= MAX(1,M+N-1)
! If WHTSVD == 4, then LIWORK >= MAX(3,M+3*N)
! If on entry LIWORK = -1, then a workspace query is
! assumed and the procedure only computes the minimal
! and the optimal workspace lengths for ZWORK, RWORK and
! IWORK. See the descriptions of ZWORK, RWORK and IWORK.
!.....
! INFO (output) INTEGER
! -i < 0 :: On entry, the i-th argument had an
! illegal value
! = 0 :: Successful return.
! = 1 :: Void input. Quick exit (M=0 or N=0).
! = 2 :: The SVD computation of X did not converge.
! Suggestion: Check the input data and/or
! repeat with different WHTSVD.
! = 3 :: The computation of the eigenvalues did not
! converge.
! = 4 :: If data scaling was requested on input and
! the procedure found inconsistency in the data
! such that for some column index i,
! X(:,i) = 0 but Y(:,i) /= 0, then Y(:,i) is set
! to zero if JOBS=='C'. The computation proceeds
! with original or modified data and warning
! flag is set with INFO=4.
!.............................................................
!.............................................................
! Parameters ! Parameters
! ~~~~~~~~~~ ! ~~~~~~~~~~
REAL(KIND=WP), PARAMETER :: ONE = 1.0_WP REAL(KIND=WP), PARAMETER :: ONE = 1.0_WP
REAL(KIND=WP), PARAMETER :: ZERO = 0.0_WP REAL(KIND=WP), PARAMETER :: ZERO = 0.0_WP
COMPLEX(KIND=WP), PARAMETER :: ZONE = ( 1.0_WP, 0.0_WP ) COMPLEX(KIND=WP), PARAMETER :: ZONE = ( 1.0_WP, 0.0_WP )
COMPLEX(KIND=WP), PARAMETER :: ZZERO = ( 0.0_WP, 0.0_WP ) COMPLEX(KIND=WP), PARAMETER :: ZZERO = ( 0.0_WP, 0.0_WP )
!
! Local scalars ! Local scalars
! ~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~
REAL(KIND=WP) :: OFL, ROOTSC, SCALE, SMALL, & REAL(KIND=WP) :: OFL, ROOTSC, SCALE, SMALL, &
@@ -400,7 +554,7 @@
! Local arrays ! Local arrays
! ~~~~~~~~~~~~ ! ~~~~~~~~~~~~
REAL(KIND=WP) :: RDUMMY(2) REAL(KIND=WP) :: RDUMMY(2)
!
! External functions (BLAS and LAPACK) ! External functions (BLAS and LAPACK)
! ~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~
REAL(KIND=WP) CLANGE, SLAMCH, SCNRM2 REAL(KIND=WP) CLANGE, SLAMCH, SCNRM2
@@ -408,13 +562,13 @@
INTEGER ICAMAX INTEGER ICAMAX
LOGICAL SISNAN, LSAME LOGICAL SISNAN, LSAME
EXTERNAL SISNAN, LSAME EXTERNAL SISNAN, LSAME
!
! External subroutines (BLAS and LAPACK) ! External subroutines (BLAS and LAPACK)
! ~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~
EXTERNAL CAXPY, CGEMM, CSSCAL EXTERNAL CAXPY, CGEMM, CSSCAL
EXTERNAL CGEEV, CGEJSV, CGESDD, CGESVD, CGESVDQ, & EXTERNAL CGEEV, CGEJSV, CGESDD, CGESVD, CGESVDQ, &
CLACPY, CLASCL, CLASSQ, XERBLA CLACPY, CLASCL, CLASSQ, XERBLA
!
! Intrinsic functions ! Intrinsic functions
! ~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~
INTRINSIC FLOAT, INT, MAX, SQRT INTRINSIC FLOAT, INT, MAX, SQRT
@@ -607,7 +761,8 @@
K = 0 K = 0
DO i = 1, N DO i = 1, N
!WORK(i) = SCNRM2( M, X(1,i), 1 ) !WORK(i) = SCNRM2( M, X(1,i), 1 )
SCALE = ZERO
SSUM = ONE
SCALE = ZERO
CALL CLASSQ( M, X(1,i), 1, SCALE, SSUM ) CALL CLASSQ( M, X(1,i), 1, SCALE, SSUM )
IF ( SISNAN(SCALE) .OR. SISNAN(SSUM) ) THEN IF ( SISNAN(SCALE) .OR. SISNAN(SSUM) ) THEN
K = 0 K = 0
@@ -680,7 +835,8 @@
! carefully computed using CLASSQ. ! carefully computed using CLASSQ.
DO i = 1, N DO i = 1, N
!RWORK(i) = SCNRM2( M, Y(1,i), 1 ) !RWORK(i) = SCNRM2( M, Y(1,i), 1 )
SCALE = ZERO
SSUM = ONE
SCALE = ZERO
CALL CLASSQ( M, Y(1,i), 1, SCALE, SSUM ) CALL CLASSQ( M, Y(1,i), 1, SCALE, SSUM )
IF ( SISNAN(SCALE) .OR. SISNAN(SSUM) ) THEN IF ( SISNAN(SCALE) .OR. SISNAN(SSUM) ) THEN
K = 0 K = 0


+ 536
- 384
lapack-netlib/SRC/dgedmd.f90
File diff suppressed because it is too large
View File


+ 534
- 382
lapack-netlib/SRC/sgedmd.f90
File diff suppressed because it is too large
View File


+ 500
- 348
lapack-netlib/SRC/zgedmd.f90 View File

@@ -1,389 +1,539 @@
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 )
! March 2023
!.....
USE iso_fortran_env
IMPLICIT NONE
INTEGER, PARAMETER :: WP = real64
!.....
! Scalar arguments
CHARACTER, INTENT(IN) :: JOBS, JOBZ, JOBR, JOBF
INTEGER, INTENT(IN) :: WHTSVD, M, N, LDX, LDY, &
NRNK, LDZ, LDB, LDW, LDS, &
LIWORK, LRWORK, LZWORK
INTEGER, INTENT(OUT) :: K, INFO
REAL(KIND=WP), INTENT(IN) :: TOL
! Array arguments
COMPLEX(KIND=WP), INTENT(INOUT) :: X(LDX,*), Y(LDY,*)
COMPLEX(KIND=WP), INTENT(OUT) :: Z(LDZ,*), B(LDB,*), &
W(LDW,*), S(LDS,*)
COMPLEX(KIND=WP), INTENT(OUT) :: EIGS(*)
COMPLEX(KIND=WP), INTENT(OUT) :: ZWORK(*)
REAL(KIND=WP), INTENT(OUT) :: RES(*)
REAL(KIND=WP), INTENT(OUT) :: RWORK(*)
INTEGER, INTENT(OUT) :: IWORK(*)
!............................................................
! 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.
!> \brief \b ZGEDMD computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.
!
! =========== DOCUMENTATION ===========
!
! Definition:
! ===========
! !
! 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 )
!......
! USE iso_fortran_env
! IMPLICIT NONE
! INTEGER, PARAMETER :: WP = real64
!
!......
! Scalar arguments
! CHARACTER, INTENT(IN) :: JOBS, JOBZ, JOBR, JOBF
! INTEGER, INTENT(IN) :: WHTSVD, M, N, LDX, LDY, &
! NRNK, LDZ, LDB, LDW, LDS, &
! LIWORK, LRWORK, LZWORK
! INTEGER, INTENT(OUT) :: K, INFO
! REAL(KIND=WP), INTENT(IN) :: TOL
! Array arguments
! COMPLEX(KIND=WP), INTENT(INOUT) :: X(LDX,*), Y(LDY,*)
! COMPLEX(KIND=WP), INTENT(OUT) :: Z(LDZ,*), B(LDB,*), &
! W(LDW,*), S(LDS,*)
! COMPLEX(KIND=WP), INTENT(OUT) :: EIGS(*)
! COMPLEX(KIND=WP), INTENT(OUT) :: ZWORK(*)
! REAL(KIND=WP), INTENT(OUT) :: RES(*)
! REAL(KIND=WP), INTENT(OUT) :: RWORK(*)
! INTEGER, INTENT(OUT) :: IWORK(*)
!
!............................................................
!> \par Purpose:
! =============
!> \verbatim
!> 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].
!> \endverbatim
!............................................................
!> \par References:
! ================
!> \verbatim
!> [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.
!> \endverbatim
!...................................................................... !......................................................................
! 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
!============================================================
! Distribution Statement A:
! Approved for Public Release, Distribution Unlimited.
! Cleared by DARPA on September 29, 2022
!============================================================
!> \par Developed and supported by:
! ================================
!> \verbatim
!> 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
!> \endverbatim
!......................................................................
!> \par Distribution Statement A:
! ==============================
!> \verbatim
!> Approved for Public Release, Distribution Unlimited.
!> Cleared by DARPA on September 29, 2022
!> \endverbatim
!............................................................ !............................................................
! Arguments ! Arguments
! ========= ! =========
! 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.
!
!> \param[in] JOBS
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! 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.
!> \param[in] JOBZ
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! 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.
!> \param[in] JOBR
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! 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.
!> \param[in] JOBF
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! 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.
!> \param[in] WHTSVD
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! M (input) INTEGER, M>= 0
! The state space dimension (the row dimension of X, Y).
!> \param[in] M
!> \verbatim
!> M (input) INTEGER, M>= 0
!> The state space dimension (the row dimension of X, Y).
!> \endverbatim
!..... !.....
! N (input) INTEGER, 0 <= N <= M
! The number of data snapshot pairs
! (the number of columns of X and Y).
!> \param[in] N
!> \verbatim
!> N (input) INTEGER, 0 <= N <= M
!> The number of data snapshot pairs
!> (the number of columns of X and Y).
!> \endverbatim
!..... !.....
! 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.
!> \param[in] LDX
!> \verbatim
!> 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.
!> LDX (input) INTEGER, LDX >= M
!> The leading dimension of the array X.
!> \endverbatim
!..... !.....
! 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.
!> \param[in,out] Y
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! LDY (input) INTEGER , LDY >= M
! The leading dimension of the array Y.
!> \param[in] LDY
!> \verbatim
!> LDY (input) INTEGER , LDY >= M
!> The leading dimension of the array Y.
!> \endverbatim
!..... !.....
! 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.
!> \param[in] NRNK
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! TOL (input) REAL(KIND=WP), 0 <= TOL < 1
! The tolerance for truncating small singular values.
! See the description of NRNK.
!> \param[in] TOL
!> \verbatim
!> TOL (input) REAL(KIND=WP), 0 <= TOL < 1
!> The tolerance for truncating small singular values.
!> See the description of NRNK.
!> \endverbatim
!..... !.....
! 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.
!> \param[out] K
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! 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.
!> \param[out] EIGS
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! 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.
!> \param[out] Z
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! LDZ (input) INTEGER , LDZ >= M
! The leading dimension of the array Z.
!> \param[in] LDZ
!> \verbatim
!> LDZ (input) INTEGER , LDZ >= M
!> The leading dimension of the array Z.
!> \endverbatim
!..... !.....
! 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.
!> \param[out] RES
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! 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.
!> \param[out] B
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! LDB (input) INTEGER, LDB >= M
! The leading dimension of the array B.
!> \param[in] LDB
!> \verbatim
!> LDB (input) INTEGER, LDB >= M
!> The leading dimension of the array B.
!> \endverbatim
!..... !.....
! 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.
!> \param[out] W
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! LDW (input) INTEGER, LDW >= N
! The leading dimension of the array W.
!> \param[in] LDW
!> \verbatim
!> LDW (input) INTEGER, LDW >= N
!> The leading dimension of the array W.
!> \endverbatim
!..... !.....
! 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.
!> \param[out] S
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! LDS (input) INTEGER, LDS >= N
! The leading dimension of the array S.
!> \param[in] LDS
!> \verbatim
!> LDS (input) INTEGER, LDS >= N
!> The leading dimension of the array S.
!> \endverbatim
!..... !.....
! 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.
!> \param[out] ZWORK
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! 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.
!> \param[in] LZWORK
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! 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.
!> \param[out] RWORK
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! 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).
!> \param[in] LRWORK
!> \verbatim
!> 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).
!> \endverbatim
!..... !.....
! 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.
!> \param[out] IWORK
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! 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.
!> \param[in] LIWORK
!> \verbatim
!> 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.
!> \endverbatim
!..... !.....
! 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.
!> \param[out] INFO
!> \verbatim
!> 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.
!> \endverbatim
!
! Authors:
! ========
!
!> \author Zlatko Drmac
!
!> \ingroup gedmd
!
!............................................................. !.............................................................
!............................................................. !.............................................................
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 )
!
! -- LAPACK driver routine --
!
! -- LAPACK is a software package provided by University of --
! -- Tennessee, University of California Berkeley, University of --
! -- Colorado Denver and NAG Ltd.. --
!
!.....
USE iso_fortran_env
IMPLICIT NONE
INTEGER, PARAMETER :: WP = real64
!
! Scalar arguments
! ~~~~~~~~~~~~~~~~
CHARACTER, INTENT(IN) :: JOBS, JOBZ, JOBR, JOBF
INTEGER, INTENT(IN) :: WHTSVD, M, N, LDX, LDY, &
NRNK, LDZ, LDB, LDW, LDS, &
LIWORK, LRWORK, LZWORK
INTEGER, INTENT(OUT) :: K, INFO
REAL(KIND=WP), INTENT(IN) :: TOL
!
! Array arguments
! ~~~~~~~~~~~~~~~
COMPLEX(KIND=WP), INTENT(INOUT) :: X(LDX,*), Y(LDY,*)
COMPLEX(KIND=WP), INTENT(OUT) :: Z(LDZ,*), B(LDB,*), &
W(LDW,*), S(LDS,*)
COMPLEX(KIND=WP), INTENT(OUT) :: EIGS(*)
COMPLEX(KIND=WP), INTENT(OUT) :: ZWORK(*)
REAL(KIND=WP), INTENT(OUT) :: RES(*)
REAL(KIND=WP), INTENT(OUT) :: RWORK(*)
INTEGER, INTENT(OUT) :: IWORK(*)
!
! Parameters ! Parameters
! ~~~~~~~~~~ ! ~~~~~~~~~~
REAL(KIND=WP), PARAMETER :: ONE = 1.0_WP REAL(KIND=WP), PARAMETER :: ONE = 1.0_WP
REAL(KIND=WP), PARAMETER :: ZERO = 0.0_WP REAL(KIND=WP), PARAMETER :: ZERO = 0.0_WP
COMPLEX(KIND=WP), PARAMETER :: ZONE = ( 1.0_WP, 0.0_WP ) COMPLEX(KIND=WP), PARAMETER :: ZONE = ( 1.0_WP, 0.0_WP )
COMPLEX(KIND=WP), PARAMETER :: ZZERO = ( 0.0_WP, 0.0_WP ) COMPLEX(KIND=WP), PARAMETER :: ZZERO = ( 0.0_WP, 0.0_WP )
!
! Local scalars ! Local scalars
! ~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~
REAL(KIND=WP) :: OFL, ROOTSC, SCALE, SMALL, & REAL(KIND=WP) :: OFL, ROOTSC, SCALE, SMALL, &
@@ -401,7 +551,7 @@
! Local arrays ! Local arrays
! ~~~~~~~~~~~~ ! ~~~~~~~~~~~~
REAL(KIND=WP) :: RDUMMY(2) REAL(KIND=WP) :: RDUMMY(2)
!
! External functions (BLAS and LAPACK) ! External functions (BLAS and LAPACK)
! ~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~
REAL(KIND=WP) ZLANGE, DLAMCH, DZNRM2 REAL(KIND=WP) ZLANGE, DLAMCH, DZNRM2
@@ -409,13 +559,13 @@
INTEGER IZAMAX INTEGER IZAMAX
LOGICAL DISNAN, LSAME LOGICAL DISNAN, LSAME
EXTERNAL DISNAN, LSAME EXTERNAL DISNAN, LSAME
!
! External subroutines (BLAS and LAPACK) ! External subroutines (BLAS and LAPACK)
! ~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~
EXTERNAL ZAXPY, ZGEMM, ZDSCAL EXTERNAL ZAXPY, ZGEMM, ZDSCAL
EXTERNAL ZGEEV, ZGEJSV, ZGESDD, ZGESVD, ZGESVDQ, & EXTERNAL ZGEEV, ZGEJSV, ZGESDD, ZGESVD, ZGESVDQ, &
ZLACPY, ZLASCL, ZLASSQ, XERBLA ZLACPY, ZLASCL, ZLASSQ, XERBLA
!
! Intrinsic functions ! Intrinsic functions
! ~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~
INTRINSIC DBLE, INT, MAX, SQRT INTRINSIC DBLE, INT, MAX, SQRT
@@ -608,7 +758,8 @@
K = 0 K = 0
DO i = 1, N DO i = 1, N
!WORK(i) = DZNRM2( M, X(1,i), 1 ) !WORK(i) = DZNRM2( M, X(1,i), 1 )
SCALE = ZERO
SSUM = ONE
SCALE = ZERO
CALL ZLASSQ( M, X(1,i), 1, SCALE, SSUM ) CALL ZLASSQ( M, X(1,i), 1, SCALE, SSUM )
IF ( DISNAN(SCALE) .OR. DISNAN(SSUM) ) THEN IF ( DISNAN(SCALE) .OR. DISNAN(SSUM) ) THEN
K = 0 K = 0
@@ -681,7 +832,8 @@
! carefully computed using ZLASSQ. ! carefully computed using ZLASSQ.
DO i = 1, N DO i = 1, N
!RWORK(i) = DZNRM2( M, Y(1,i), 1 ) !RWORK(i) = DZNRM2( M, Y(1,i), 1 )
SCALE = ZERO
SSUM = ONE
SCALE = ZERO
CALL ZLASSQ( M, Y(1,i), 1, SCALE, SSUM ) CALL ZLASSQ( M, Y(1,i), 1, SCALE, SSUM )
IF ( DISNAN(SCALE) .OR. DISNAN(SSUM) ) THEN IF ( DISNAN(SCALE) .OR. DISNAN(SSUM) ) THEN
K = 0 K = 0


+ 12
- 0
lapack-netlib/TESTING/CMakeLists.txt View File

@@ -54,6 +54,9 @@ add_lapack_test(sgqr.out gqr.in xeigtsts)
add_lapack_test(sgsv.out gsv.in xeigtsts) add_lapack_test(sgsv.out gsv.in xeigtsts)
add_lapack_test(scsd.out csd.in xeigtsts) add_lapack_test(scsd.out csd.in xeigtsts)
add_lapack_test(slse.out lse.in xeigtsts) add_lapack_test(slse.out lse.in xeigtsts)
#
# ======== SINGLE DMD EIG TESTS ===========================
add_lapack_test(sdmd.out sdmd.in xdmdeigtsts)
endif() endif()


if(BUILD_DOUBLE) if(BUILD_DOUBLE)
@@ -85,6 +88,9 @@ add_lapack_test(dgqr.out gqr.in xeigtstd)
add_lapack_test(dgsv.out gsv.in xeigtstd) add_lapack_test(dgsv.out gsv.in xeigtstd)
add_lapack_test(dcsd.out csd.in xeigtstd) add_lapack_test(dcsd.out csd.in xeigtstd)
add_lapack_test(dlse.out lse.in xeigtstd) add_lapack_test(dlse.out lse.in xeigtstd)
#
# ======== DOUBLE DMD EIG TESTS ===========================
add_lapack_test(ddmd.out ddmd.in xdmdeigtstd)
endif() endif()


if(BUILD_COMPLEX) if(BUILD_COMPLEX)
@@ -114,6 +120,9 @@ add_lapack_test(cgqr.out gqr.in xeigtstc)
add_lapack_test(cgsv.out gsv.in xeigtstc) add_lapack_test(cgsv.out gsv.in xeigtstc)
add_lapack_test(ccsd.out csd.in xeigtstc) add_lapack_test(ccsd.out csd.in xeigtstc)
add_lapack_test(clse.out lse.in xeigtstc) add_lapack_test(clse.out lse.in xeigtstc)
#
# ======== COMPLEX DMD EIG TESTS ===========================
add_lapack_test(cdmd.out cdmd.in xdmdeigtstc)
endif() endif()


if(BUILD_COMPLEX16) if(BUILD_COMPLEX16)
@@ -145,6 +154,9 @@ add_lapack_test(zgqr.out gqr.in xeigtstz)
add_lapack_test(zgsv.out gsv.in xeigtstz) add_lapack_test(zgsv.out gsv.in xeigtstz)
add_lapack_test(zcsd.out csd.in xeigtstz) add_lapack_test(zcsd.out csd.in xeigtstz)
add_lapack_test(zlse.out lse.in xeigtstz) add_lapack_test(zlse.out lse.in xeigtstz)
#
# ======== COMPLEX16 DMD EIG TESTS ===========================
add_lapack_test(zdmd.out zdmd.in xdmdeigtstz)
endif() endif()






+ 12
- 0
lapack-netlib/TESTING/EIG/CMakeLists.txt View File

@@ -42,6 +42,8 @@ set(SEIGTST schkee.F
sort03.f ssbt21.f ssgt01.f sslect.f sspt21.f sstt21.f sort03.f ssbt21.f ssgt01.f sslect.f sspt21.f sstt21.f
sstt22.f ssyl01.f ssyt21.f ssyt22.f) sstt22.f ssyl01.f ssyt21.f ssyt22.f)


set(SDMDEIGTST schkdmd.f90)

set(CEIGTST cchkee.F set(CEIGTST cchkee.F
cbdt01.f cbdt02.f cbdt03.f cbdt05.f cbdt01.f cbdt02.f cbdt03.f cbdt05.f
cchkbb.f cchkbd.f cchkbk.f cchkbl.f cchkec.f cchkbb.f cchkbd.f cchkbk.f cchkbl.f cchkec.f
@@ -59,6 +61,8 @@ set(CEIGTST cchkee.F
csgt01.f cslect.f csyl01.f csgt01.f cslect.f csyl01.f
cstt21.f cstt22.f cunt01.f cunt03.f) cstt21.f cstt22.f cunt01.f cunt03.f)


set(CDMDEIGTST cchkdmd.f90)

set(DZIGTST dlafts.f dlahd2.f dlasum.f dlatb9.f dstech.f dstect.f set(DZIGTST dlafts.f dlahd2.f dlasum.f dlatb9.f dstech.f dstect.f
dsvdch.f dsvdct.f dsxt1.f) dsvdch.f dsvdct.f dsxt1.f)


@@ -79,6 +83,8 @@ set(DEIGTST dchkee.F
dort03.f dsbt21.f dsgt01.f dslect.f dspt21.f dstt21.f dort03.f dsbt21.f dsgt01.f dslect.f dspt21.f dstt21.f
dstt22.f dsyl01.f dsyt21.f dsyt22.f) dstt22.f dsyl01.f dsyt21.f dsyt22.f)


set(DDMDEIGTST dchkdmd.f90)

set(ZEIGTST zchkee.F set(ZEIGTST zchkee.F
zbdt01.f zbdt02.f zbdt03.f zbdt05.f zbdt01.f zbdt02.f zbdt03.f zbdt05.f
zchkbb.f zchkbd.f zchkbk.f zchkbl.f zchkec.f zchkbb.f zchkbd.f zchkbk.f zchkbl.f zchkec.f
@@ -96,6 +102,8 @@ set(ZEIGTST zchkee.F
zsgt01.f zslect.f zsyl01.f zsgt01.f zslect.f zsyl01.f
zstt21.f zstt22.f zunt01.f zunt03.f) zstt21.f zstt22.f zunt01.f zunt03.f)


set(ZDMDEIGTST zchkdmd.f90)

macro(add_eig_executable name) macro(add_eig_executable name)
add_executable(${name} ${ARGN}) add_executable(${name} ${ARGN})
target_link_libraries(${name} openblas${SUFFIX64_UNDERSCORE}) target_link_libraries(${name} openblas${SUFFIX64_UNDERSCORE})
@@ -104,16 +112,20 @@ endmacro()


if(BUILD_SINGLE) if(BUILD_SINGLE)
add_eig_executable(xeigtsts ${SEIGTST} ${SCIGTST} ${AEIGTST}) add_eig_executable(xeigtsts ${SEIGTST} ${SCIGTST} ${AEIGTST})
add_eig_executable(xdmdeigtsts ${SDMDEIGTST})
endif() endif()


if(BUILD_COMPLEX) if(BUILD_COMPLEX)
add_eig_executable(xeigtstc ${CEIGTST} ${SCIGTST} ${AEIGTST}) add_eig_executable(xeigtstc ${CEIGTST} ${SCIGTST} ${AEIGTST})
add_eig_executable(xdmdeigtstc ${CDMDEIGTST})
endif() endif()


if(BUILD_DOUBLE) if(BUILD_DOUBLE)
add_eig_executable(xeigtstd ${DEIGTST} ${DZIGTST} ${AEIGTST}) add_eig_executable(xeigtstd ${DEIGTST} ${DZIGTST} ${AEIGTST})
add_eig_executable(xdmdeigtstd ${DDMDEIGTST})
endif() endif()


if(BUILD_COMPLEX16) if(BUILD_COMPLEX16)
add_eig_executable(xeigtstz ${ZEIGTST} ${DZIGTST} ${AEIGTST}) add_eig_executable(xeigtstz ${ZEIGTST} ${DZIGTST} ${AEIGTST})
add_eig_executable(xdmdeigtstz ${ZDMDEIGTST})
endif() endif()

Write Preview
Loading…
Cancel
Save