Fix issues related to ?GEDMD (Reference-LAPACK PR 959)
This commit is contained in:
Martin Kroeker
GitHub
parent
05dde8ef04
commit
201f22f49a
|
|
@ -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, &
|
||||
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
|
||||
!
|
||||
! -- 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 = 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,*)
|
||||
|
|
@ -25,364 +529,14 @@
|
|||
REAL(KIND=WP), INTENT(OUT) :: RES(*)
|
||||
REAL(KIND=WP), INTENT(OUT) :: RWORK(*)
|
||||
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
|
||||
! ~~~~~~~~~~
|
||||
REAL(KIND=WP), PARAMETER :: ONE = 1.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 :: ZZERO = ( 0.0_WP, 0.0_WP )
|
||||
|
||||
!
|
||||
! Local scalars
|
||||
! ~~~~~~~~~~~~~
|
||||
REAL(KIND=WP) :: OFL, ROOTSC, SCALE, SMALL, &
|
||||
|
|
@ -400,7 +554,7 @@
|
|||
! Local arrays
|
||||
! ~~~~~~~~~~~~
|
||||
REAL(KIND=WP) :: RDUMMY(2)
|
||||
|
||||
!
|
||||
! External functions (BLAS and LAPACK)
|
||||
! ~~~~~~~~~~~~~~~~~
|
||||
REAL(KIND=WP) CLANGE, SLAMCH, SCNRM2
|
||||
|
|
@ -408,13 +562,13 @@
|
|||
INTEGER ICAMAX
|
||||
LOGICAL SISNAN, LSAME
|
||||
EXTERNAL SISNAN, LSAME
|
||||
|
||||
!
|
||||
! External subroutines (BLAS and LAPACK)
|
||||
! ~~~~~~~~~~~~~~~~~~~~
|
||||
EXTERNAL CAXPY, CGEMM, CSSCAL
|
||||
EXTERNAL CGEEV, CGEJSV, CGESDD, CGESVD, CGESVDQ, &
|
||||
CLACPY, CLASCL, CLASSQ, XERBLA
|
||||
|
||||
!
|
||||
! Intrinsic functions
|
||||
! ~~~~~~~~~~~~~~~~~~~
|
||||
INTRINSIC FLOAT, INT, MAX, SQRT
|
||||
|
|
@ -607,7 +761,8 @@
|
|||
K = 0
|
||||
DO i = 1, N
|
||||
!WORK(i) = SCNRM2( M, X(1,i), 1 )
|
||||
SCALE = ZERO
|
||||
SSUM = ONE
|
||||
SCALE = ZERO
|
||||
CALL CLASSQ( M, X(1,i), 1, SCALE, SSUM )
|
||||
IF ( SISNAN(SCALE) .OR. SISNAN(SSUM) ) THEN
|
||||
K = 0
|
||||
|
|
@ -680,7 +835,8 @@
|
|||
! carefully computed using CLASSQ.
|
||||
DO i = 1, N
|
||||
!RWORK(i) = SCNRM2( M, Y(1,i), 1 )
|
||||
SCALE = ZERO
|
||||
SSUM = ONE
|
||||
SCALE = ZERO
|
||||
CALL CLASSQ( M, Y(1,i), 1, SCALE, SSUM )
|
||||
IF ( SISNAN(SCALE) .OR. SISNAN(SSUM) ) THEN
|
||||
K = 0
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
|
|
@ -1,23 +1,523 @@
|
|||
!> \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
|
||||
!......................................................................
|
||||
!> \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 :: 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
|
||||
!.....
|
||||
!> \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] 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.
|
||||
!> \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 ZGEEV.
|
||||
!> 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 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
|
||||
!.....
|
||||
!> \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
|
||||
!.....
|
||||
!> \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
|
||||
!.....
|
||||
!> \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
|
||||
!.....
|
||||
!> \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 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
|
||||
!
|
||||
! -- 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,*)
|
||||
|
|
@ -26,364 +526,14 @@
|
|||
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.
|
||||
!
|
||||
!......................................................................
|
||||
! 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 :: ZGESVD (the QR SVD algorithm)
|
||||
! 2 :: ZGESDD (the Divide and Conquer algorithm; if enough
|
||||
! workspace available, this is the fastest option)
|
||||
! 3 :: ZGESVDQ (the preconditioned QR SVD ; this and 4
|
||||
! are the most accurate options)
|
||||
! 4 :: ZGEJSV (the preconditioned Jacobi SVD; this and 3
|
||||
! are the most accurate options)
|
||||
! For the four methods above, a significant difference in
|
||||
! the accuracy of small singular values is possible if
|
||||
! the snapshots vary in norm so that X is severely
|
||||
! ill-conditioned. If small (smaller than EPS*||X||)
|
||||
! singular values are of interest and JOBS=='N', then
|
||||
! the options (3, 4) give the most accurate results, where
|
||||
! the option 4 is slightly better and with stronger
|
||||
! theoretical background.
|
||||
! If JOBS=='S', i.e. the columns of X will be normalized,
|
||||
! then all methods give nearly equally accurate results.
|
||||
!.....
|
||||
! M (input) INTEGER, M>= 0
|
||||
! The state space dimension (the row dimension of X, Y).
|
||||
!.....
|
||||
! N (input) INTEGER, 0 <= N <= M
|
||||
! The number of data snapshot pairs
|
||||
! (the number of columns of X and Y).
|
||||
!.....
|
||||
! 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 ZGEEV.
|
||||
! 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 ZGEEV for computing
|
||||
! the eigenvalues of a Rayleigh quotient.
|
||||
! If the call to ZGEDMD is only workspace query, then
|
||||
! ZWORK(1) contains the minimal complex workspace length and
|
||||
! ZWORK(2) is the optimal complex workspace length.
|
||||
! Hence, the length of work is at least 2.
|
||||
! See the description of LZWORK.
|
||||
!.....
|
||||
! LZWORK (input) INTEGER
|
||||
! The minimal length of the workspace vector ZWORK.
|
||||
! LZWORK is calculated as MAX(LZWORK_SVD, LZWORK_ZGEEV),
|
||||
! where LZWORK_ZGEEV = MAX( 1, 2*N ) and the minimal
|
||||
! LZWORK_SVD is calculated as follows
|
||||
! If WHTSVD == 1 :: ZGESVD ::
|
||||
! LZWORK_SVD = MAX(1,2*MIN(M,N)+MAX(M,N))
|
||||
! If WHTSVD == 2 :: ZGESDD ::
|
||||
! LZWORK_SVD = 2*MIN(M,N)*MIN(M,N)+2*MIN(M,N)+MAX(M,N)
|
||||
! If WHTSVD == 3 :: ZGESVDQ ::
|
||||
! LZWORK_SVD = obtainable by a query
|
||||
! If WHTSVD == 4 :: ZGEJSV ::
|
||||
! LZWORK_SVD = obtainable by a query
|
||||
! If on entry LZWORK = -1, then a workspace query is
|
||||
! assumed and the procedure only computes the minimal
|
||||
! and the optimal workspace lengths and returns them in
|
||||
! LZWORK(1) and LZWORK(2), respectively.
|
||||
!.....
|
||||
! RWORK (workspace/output) REAL(KIND=WP) LRWORK-by-1 array
|
||||
! On exit, RWORK(1:N) contains the singular values of
|
||||
! X (for JOBS=='N') or column scaled X (JOBS=='S', 'C').
|
||||
! If WHTSVD==4, then RWORK(N+1) and RWORK(N+2) contain
|
||||
! scaling factor RWORK(N+2)/RWORK(N+1) used to scale X
|
||||
! and Y to avoid overflow in the SVD of X.
|
||||
! This may be of interest if the scaling option is off
|
||||
! and as many as possible smallest eigenvalues are
|
||||
! desired to the highest feasible accuracy.
|
||||
! If the call to ZGEDMD is only workspace query, then
|
||||
! RWORK(1) contains the minimal workspace length.
|
||||
! See the description of LRWORK.
|
||||
!.....
|
||||
! LRWORK (input) INTEGER
|
||||
! The minimal length of the workspace vector RWORK.
|
||||
! LRWORK is calculated as follows:
|
||||
! LRWORK = MAX(1, N+LRWORK_SVD,N+LRWORK_ZGEEV), where
|
||||
! LRWORK_ZGEEV = MAX(1,2*N) and RWORK_SVD is the real workspace
|
||||
! for the SVD subroutine determined by the input parameter
|
||||
! WHTSVD.
|
||||
! If WHTSVD == 1 :: ZGESVD ::
|
||||
! LRWORK_SVD = 5*MIN(M,N)
|
||||
! If WHTSVD == 2 :: ZGESDD ::
|
||||
! LRWORK_SVD = MAX(5*MIN(M,N)*MIN(M,N)+7*MIN(M,N),
|
||||
! 2*MAX(M,N)*MIN(M,N)+2*MIN(M,N)*MIN(M,N)+MIN(M,N) ) )
|
||||
! If WHTSVD == 3 :: ZGESVDQ ::
|
||||
! LRWORK_SVD = obtainable by a query
|
||||
! If WHTSVD == 4 :: ZGEJSV ::
|
||||
! LRWORK_SVD = obtainable by a query
|
||||
! If on entry LRWORK = -1, then a workspace query is
|
||||
! assumed and the procedure only computes the minimal
|
||||
! real workspace length and returns it in RWORK(1).
|
||||
!.....
|
||||
! IWORK (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
|
||||
! ~~~~~~~~~~
|
||||
REAL(KIND=WP), PARAMETER :: ONE = 1.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 :: ZZERO = ( 0.0_WP, 0.0_WP )
|
||||
|
||||
!
|
||||
! Local scalars
|
||||
! ~~~~~~~~~~~~~
|
||||
REAL(KIND=WP) :: OFL, ROOTSC, SCALE, SMALL, &
|
||||
|
|
@ -401,7 +551,7 @@
|
|||
! Local arrays
|
||||
! ~~~~~~~~~~~~
|
||||
REAL(KIND=WP) :: RDUMMY(2)
|
||||
|
||||
!
|
||||
! External functions (BLAS and LAPACK)
|
||||
! ~~~~~~~~~~~~~~~~~
|
||||
REAL(KIND=WP) ZLANGE, DLAMCH, DZNRM2
|
||||
|
|
@ -409,13 +559,13 @@
|
|||
INTEGER IZAMAX
|
||||
LOGICAL DISNAN, LSAME
|
||||
EXTERNAL DISNAN, LSAME
|
||||
|
||||
!
|
||||
! External subroutines (BLAS and LAPACK)
|
||||
! ~~~~~~~~~~~~~~~~~~~~
|
||||
EXTERNAL ZAXPY, ZGEMM, ZDSCAL
|
||||
EXTERNAL ZGEEV, ZGEJSV, ZGESDD, ZGESVD, ZGESVDQ, &
|
||||
ZLACPY, ZLASCL, ZLASSQ, XERBLA
|
||||
|
||||
!
|
||||
! Intrinsic functions
|
||||
! ~~~~~~~~~~~~~~~~~~~
|
||||
INTRINSIC DBLE, INT, MAX, SQRT
|
||||
|
|
@ -608,7 +758,8 @@
|
|||
K = 0
|
||||
DO i = 1, N
|
||||
!WORK(i) = DZNRM2( M, X(1,i), 1 )
|
||||
SCALE = ZERO
|
||||
SSUM = ONE
|
||||
SCALE = ZERO
|
||||
CALL ZLASSQ( M, X(1,i), 1, SCALE, SSUM )
|
||||
IF ( DISNAN(SCALE) .OR. DISNAN(SSUM) ) THEN
|
||||
K = 0
|
||||
|
|
@ -681,7 +832,8 @@
|
|||
! carefully computed using ZLASSQ.
|
||||
DO i = 1, N
|
||||
!RWORK(i) = DZNRM2( M, Y(1,i), 1 )
|
||||
SCALE = ZERO
|
||||
SSUM = ONE
|
||||
SCALE = ZERO
|
||||
CALL ZLASSQ( M, Y(1,i), 1, SCALE, SSUM )
|
||||
IF ( DISNAN(SCALE) .OR. DISNAN(SSUM) ) THEN
|
||||
K = 0
|
||||
|
|
|
|||
Loading…
Reference in New Issue