1149 lines
47 KiB
Fortran
1149 lines
47 KiB
Fortran
!> \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 )
|
|
!
|
|
! -- 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
|
|
! ~~~~~~~~~~
|
|
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, &
|
|
SSUM, XSCL1, XSCL2
|
|
INTEGER :: i, j, IMINWR, INFO1, INFO2, &
|
|
LWRKEV, LWRSDD, LWRSVD, LWRSVJ, &
|
|
LWRSVQ, MLWORK, MWRKEV, MWRSDD, &
|
|
MWRSVD, MWRSVJ, MWRSVQ, NUMRNK, &
|
|
OLWORK, MLRWRK
|
|
LOGICAL :: BADXY, LQUERY, SCCOLX, SCCOLY, &
|
|
WNTEX, WNTREF, WNTRES, WNTVEC
|
|
CHARACTER :: JOBZL, T_OR_N
|
|
CHARACTER :: JSVOPT
|
|
!
|
|
! Local arrays
|
|
! ~~~~~~~~~~~~
|
|
REAL(KIND=WP) :: RDUMMY(2)
|
|
!
|
|
! External functions (BLAS and LAPACK)
|
|
! ~~~~~~~~~~~~~~~~~
|
|
REAL(KIND=WP) ZLANGE, DLAMCH, DZNRM2
|
|
EXTERNAL ZLANGE, DLAMCH, DZNRM2, IZAMAX
|
|
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
|
|
!............................................................
|
|
!
|
|
! Test the input arguments
|
|
!
|
|
WNTRES = LSAME(JOBR,'R')
|
|
SCCOLX = LSAME(JOBS,'S') .OR. LSAME(JOBS,'C')
|
|
SCCOLY = LSAME(JOBS,'Y')
|
|
WNTVEC = LSAME(JOBZ,'V')
|
|
WNTREF = LSAME(JOBF,'R')
|
|
WNTEX = LSAME(JOBF,'E')
|
|
INFO = 0
|
|
LQUERY = ( ( LZWORK == -1 ) .OR. ( LIWORK == -1 ) &
|
|
.OR. ( LRWORK == -1 ) )
|
|
!
|
|
IF ( .NOT. (SCCOLX .OR. SCCOLY .OR. &
|
|
LSAME(JOBS,'N')) ) THEN
|
|
INFO = -1
|
|
ELSE IF ( .NOT. (WNTVEC .OR. LSAME(JOBZ,'N') &
|
|
.OR. LSAME(JOBZ,'F')) ) THEN
|
|
INFO = -2
|
|
ELSE IF ( .NOT. (WNTRES .OR. LSAME(JOBR,'N')) .OR. &
|
|
( WNTRES .AND. (.NOT.WNTVEC) ) ) THEN
|
|
INFO = -3
|
|
ELSE IF ( .NOT. (WNTREF .OR. WNTEX .OR. &
|
|
LSAME(JOBF,'N') ) ) THEN
|
|
INFO = -4
|
|
ELSE IF ( .NOT.((WHTSVD == 1) .OR. (WHTSVD == 2) .OR. &
|
|
(WHTSVD == 3) .OR. (WHTSVD == 4) )) THEN
|
|
INFO = -5
|
|
ELSE IF ( M < 0 ) THEN
|
|
INFO = -6
|
|
ELSE IF ( ( N < 0 ) .OR. ( N > M ) ) THEN
|
|
INFO = -7
|
|
ELSE IF ( LDX < M ) THEN
|
|
INFO = -9
|
|
ELSE IF ( LDY < M ) THEN
|
|
INFO = -11
|
|
ELSE IF ( .NOT. (( NRNK == -2).OR.(NRNK == -1).OR. &
|
|
((NRNK >= 1).AND.(NRNK <=N ))) ) THEN
|
|
INFO = -12
|
|
ELSE IF ( ( TOL < ZERO ) .OR. ( TOL >= ONE ) ) THEN
|
|
INFO = -13
|
|
ELSE IF ( LDZ < M ) THEN
|
|
INFO = -17
|
|
ELSE IF ( (WNTREF .OR. WNTEX ) .AND. ( LDB < M ) ) THEN
|
|
INFO = -20
|
|
ELSE IF ( LDW < N ) THEN
|
|
INFO = -22
|
|
ELSE IF ( LDS < N ) THEN
|
|
INFO = -24
|
|
END IF
|
|
!
|
|
IF ( INFO == 0 ) THEN
|
|
! Compute the minimal and the optimal workspace
|
|
! requirements. Simulate running the code and
|
|
! determine minimal and optimal sizes of the
|
|
! workspace at any moment of the run.
|
|
IF ( N == 0 ) THEN
|
|
! Quick return. All output except K is void.
|
|
! INFO=1 signals the void input.
|
|
! In case of a workspace query, the default
|
|
! minimal workspace lengths are returned.
|
|
IF ( LQUERY ) THEN
|
|
IWORK(1) = 1
|
|
RWORK(1) = 1
|
|
ZWORK(1) = 2
|
|
ZWORK(2) = 2
|
|
ELSE
|
|
K = 0
|
|
END IF
|
|
INFO = 1
|
|
RETURN
|
|
END IF
|
|
|
|
IMINWR = 1
|
|
MLRWRK = MAX(1,N)
|
|
MLWORK = 2
|
|
OLWORK = 2
|
|
SELECT CASE ( WHTSVD )
|
|
CASE (1)
|
|
! The following is specified as the minimal
|
|
! length of WORK in the definition of ZGESVD:
|
|
! MWRSVD = MAX(1,2*MIN(M,N)+MAX(M,N))
|
|
MWRSVD = MAX(1,2*MIN(M,N)+MAX(M,N))
|
|
MLWORK = MAX(MLWORK,MWRSVD)
|
|
MLRWRK = MAX(MLRWRK,N + 5*MIN(M,N))
|
|
IF ( LQUERY ) THEN
|
|
CALL ZGESVD( 'O', 'S', M, N, X, LDX, RWORK, &
|
|
B, LDB, W, LDW, ZWORK, -1, RDUMMY, INFO1 )
|
|
LWRSVD = INT( ZWORK(1) )
|
|
OLWORK = MAX(OLWORK,LWRSVD)
|
|
END IF
|
|
CASE (2)
|
|
! The following is specified as the minimal
|
|
! length of WORK in the definition of ZGESDD:
|
|
! MWRSDD = 2*min(M,N)*min(M,N)+2*min(M,N)+max(M,N).
|
|
! RWORK length: 5*MIN(M,N)*MIN(M,N)+7*MIN(M,N)
|
|
! In LAPACK 3.10.1 RWORK is defined differently.
|
|
! Below we take max over the two versions.
|
|
! IMINWR = 8*MIN(M,N)
|
|
MWRSDD = 2*MIN(M,N)*MIN(M,N)+2*MIN(M,N)+MAX(M,N)
|
|
MLWORK = MAX(MLWORK,MWRSDD)
|
|
IMINWR = 8*MIN(M,N)
|
|
MLRWRK = MAX( MLRWRK, N + &
|
|
MAX( 5*MIN(M,N)*MIN(M,N)+7*MIN(M,N), &
|
|
5*MIN(M,N)*MIN(M,N)+5*MIN(M,N), &
|
|
2*MAX(M,N)*MIN(M,N)+ &
|
|
2*MIN(M,N)*MIN(M,N)+MIN(M,N) ) )
|
|
IF ( LQUERY ) THEN
|
|
CALL ZGESDD( 'O', M, N, X, LDX, RWORK, B,LDB,&
|
|
W, LDW, ZWORK, -1, RDUMMY, IWORK, INFO1 )
|
|
LWRSDD = MAX( MWRSDD,INT( ZWORK(1) ))
|
|
! Possible bug in ZGESDD optimal workspace size.
|
|
OLWORK = MAX(OLWORK,LWRSDD)
|
|
END IF
|
|
CASE (3)
|
|
CALL ZGESVDQ( 'H', 'P', 'N', 'R', 'R', M, N, &
|
|
X, LDX, RWORK, Z, LDZ, W, LDW, NUMRNK, &
|
|
IWORK, -1, ZWORK, -1, RDUMMY, -1, INFO1 )
|
|
IMINWR = IWORK(1)
|
|
MWRSVQ = INT(ZWORK(2))
|
|
MLWORK = MAX(MLWORK,MWRSVQ)
|
|
MLRWRK = MAX(MLRWRK,N + INT(RDUMMY(1)))
|
|
IF ( LQUERY ) THEN
|
|
LWRSVQ = INT(ZWORK(1))
|
|
OLWORK = MAX(OLWORK,LWRSVQ)
|
|
END IF
|
|
CASE (4)
|
|
JSVOPT = 'J'
|
|
CALL ZGEJSV( 'F', 'U', JSVOPT, 'R', 'N', 'P', M, &
|
|
N, X, LDX, RWORK, Z, LDZ, W, LDW, &
|
|
ZWORK, -1, RDUMMY, -1, IWORK, INFO1 )
|
|
IMINWR = IWORK(1)
|
|
MWRSVJ = INT(ZWORK(2))
|
|
MLWORK = MAX(MLWORK,MWRSVJ)
|
|
MLRWRK = MAX(MLRWRK,N + MAX(7,INT(RDUMMY(1))))
|
|
IF ( LQUERY ) THEN
|
|
LWRSVJ = INT(ZWORK(1))
|
|
OLWORK = MAX(OLWORK,LWRSVJ)
|
|
END IF
|
|
END SELECT
|
|
IF ( WNTVEC .OR. WNTEX .OR. LSAME(JOBZ,'F') ) THEN
|
|
JOBZL = 'V'
|
|
ELSE
|
|
JOBZL = 'N'
|
|
END IF
|
|
! Workspace calculation to the ZGEEV call
|
|
MWRKEV = MAX( 1, 2*N )
|
|
MLWORK = MAX(MLWORK,MWRKEV)
|
|
MLRWRK = MAX(MLRWRK,N+2*N)
|
|
IF ( LQUERY ) THEN
|
|
CALL ZGEEV( 'N', JOBZL, N, S, LDS, EIGS, &
|
|
W, LDW, W, LDW, ZWORK, -1, RWORK, INFO1 )
|
|
LWRKEV = INT(ZWORK(1))
|
|
OLWORK = MAX( OLWORK, LWRKEV )
|
|
END IF
|
|
!
|
|
IF ( LIWORK < IMINWR .AND. (.NOT.LQUERY) ) INFO = -30
|
|
IF ( LRWORK < MLRWRK .AND. (.NOT.LQUERY) ) INFO = -28
|
|
IF ( LZWORK < MLWORK .AND. (.NOT.LQUERY) ) INFO = -26
|
|
|
|
END IF
|
|
!
|
|
IF( INFO /= 0 ) THEN
|
|
CALL XERBLA( 'ZGEDMD', -INFO )
|
|
RETURN
|
|
ELSE IF ( LQUERY ) THEN
|
|
! Return minimal and optimal workspace sizes
|
|
IWORK(1) = IMINWR
|
|
RWORK(1) = MLRWRK
|
|
ZWORK(1) = MLWORK
|
|
ZWORK(2) = OLWORK
|
|
RETURN
|
|
END IF
|
|
!............................................................
|
|
!
|
|
OFL = DLAMCH('O')
|
|
SMALL = DLAMCH('S')
|
|
BADXY = .FALSE.
|
|
!
|
|
! <1> Optional scaling of the snapshots (columns of X, Y)
|
|
! ==========================================================
|
|
IF ( SCCOLX ) THEN
|
|
! The columns of X will be normalized.
|
|
! To prevent overflows, the column norms of X are
|
|
! carefully computed using ZLASSQ.
|
|
K = 0
|
|
DO i = 1, N
|
|
!WORK(i) = DZNRM2( M, X(1,i), 1 )
|
|
SSUM = ONE
|
|
SCALE = ZERO
|
|
CALL ZLASSQ( M, X(1,i), 1, SCALE, SSUM )
|
|
IF ( DISNAN(SCALE) .OR. DISNAN(SSUM) ) THEN
|
|
K = 0
|
|
INFO = -8
|
|
CALL XERBLA('ZGEDMD',-INFO)
|
|
END IF
|
|
IF ( (SCALE /= ZERO) .AND. (SSUM /= ZERO) ) THEN
|
|
ROOTSC = SQRT(SSUM)
|
|
IF ( SCALE .GE. (OFL / ROOTSC) ) THEN
|
|
! Norm of X(:,i) overflows. First, X(:,i)
|
|
! is scaled by
|
|
! ( ONE / ROOTSC ) / SCALE = 1/||X(:,i)||_2.
|
|
! Next, the norm of X(:,i) is stored without
|
|
! overflow as RWORK(i) = - SCALE * (ROOTSC/M),
|
|
! the minus sign indicating the 1/M factor.
|
|
! Scaling is performed without overflow, and
|
|
! underflow may occur in the smallest entries
|
|
! of X(:,i). The relative backward and forward
|
|
! errors are small in the ell_2 norm.
|
|
CALL ZLASCL( 'G', 0, 0, SCALE, ONE/ROOTSC, &
|
|
M, 1, X(1,i), LDX, INFO2 )
|
|
RWORK(i) = - SCALE * ( ROOTSC / DBLE(M) )
|
|
ELSE
|
|
! X(:,i) will be scaled to unit 2-norm
|
|
RWORK(i) = SCALE * ROOTSC
|
|
CALL ZLASCL( 'G',0, 0, RWORK(i), ONE, M, 1, &
|
|
X(1,i), LDX, INFO2 ) ! LAPACK CALL
|
|
! X(1:M,i) = (ONE/RWORK(i)) * X(1:M,i) ! INTRINSIC
|
|
END IF
|
|
ELSE
|
|
RWORK(i) = ZERO
|
|
K = K + 1
|
|
END IF
|
|
END DO
|
|
IF ( K == N ) THEN
|
|
! All columns of X are zero. Return error code -8.
|
|
! (the 8th input variable had an illegal value)
|
|
K = 0
|
|
INFO = -8
|
|
CALL XERBLA('ZGEDMD',-INFO)
|
|
RETURN
|
|
END IF
|
|
DO i = 1, N
|
|
! Now, apply the same scaling to the columns of Y.
|
|
IF ( RWORK(i) > ZERO ) THEN
|
|
CALL ZDSCAL( M, ONE/RWORK(i), Y(1,i), 1 ) ! BLAS CALL
|
|
! Y(1:M,i) = (ONE/RWORK(i)) * Y(1:M,i) ! INTRINSIC
|
|
ELSE IF ( RWORK(i) < ZERO ) THEN
|
|
CALL ZLASCL( 'G', 0, 0, -RWORK(i), &
|
|
ONE/DBLE(M), M, 1, Y(1,i), LDY, INFO2 ) ! LAPACK CALL
|
|
ELSE IF ( ABS(Y(IZAMAX(M, Y(1,i),1),i )) &
|
|
/= ZERO ) THEN
|
|
! X(:,i) is zero vector. For consistency,
|
|
! Y(:,i) should also be zero. If Y(:,i) is not
|
|
! zero, then the data might be inconsistent or
|
|
! corrupted. If JOBS == 'C', Y(:,i) is set to
|
|
! zero and a warning flag is raised.
|
|
! The computation continues but the
|
|
! situation will be reported in the output.
|
|
BADXY = .TRUE.
|
|
IF ( LSAME(JOBS,'C')) &
|
|
CALL ZDSCAL( M, ZERO, Y(1,i), 1 ) ! BLAS CALL
|
|
END IF
|
|
END DO
|
|
END IF
|
|
!
|
|
IF ( SCCOLY ) THEN
|
|
! The columns of Y will be normalized.
|
|
! To prevent overflows, the column norms of Y are
|
|
! carefully computed using ZLASSQ.
|
|
DO i = 1, N
|
|
!RWORK(i) = DZNRM2( M, Y(1,i), 1 )
|
|
SSUM = ONE
|
|
SCALE = ZERO
|
|
CALL ZLASSQ( M, Y(1,i), 1, SCALE, SSUM )
|
|
IF ( DISNAN(SCALE) .OR. DISNAN(SSUM) ) THEN
|
|
K = 0
|
|
INFO = -10
|
|
CALL XERBLA('ZGEDMD',-INFO)
|
|
END IF
|
|
IF ( SCALE /= ZERO .AND. (SSUM /= ZERO) ) THEN
|
|
ROOTSC = SQRT(SSUM)
|
|
IF ( SCALE .GE. (OFL / ROOTSC) ) THEN
|
|
! Norm of Y(:,i) overflows. First, Y(:,i)
|
|
! is scaled by
|
|
! ( ONE / ROOTSC ) / SCALE = 1/||Y(:,i)||_2.
|
|
! Next, the norm of Y(:,i) is stored without
|
|
! overflow as RWORK(i) = - SCALE * (ROOTSC/M),
|
|
! the minus sign indicating the 1/M factor.
|
|
! Scaling is performed without overflow, and
|
|
! underflow may occur in the smallest entries
|
|
! of Y(:,i). The relative backward and forward
|
|
! errors are small in the ell_2 norm.
|
|
CALL ZLASCL( 'G', 0, 0, SCALE, ONE/ROOTSC, &
|
|
M, 1, Y(1,i), LDY, INFO2 )
|
|
RWORK(i) = - SCALE * ( ROOTSC / DBLE(M) )
|
|
ELSE
|
|
! Y(:,i) will be scaled to unit 2-norm
|
|
RWORK(i) = SCALE * ROOTSC
|
|
CALL ZLASCL( 'G',0, 0, RWORK(i), ONE, M, 1, &
|
|
Y(1,i), LDY, INFO2 ) ! LAPACK CALL
|
|
! Y(1:M,i) = (ONE/RWORK(i)) * Y(1:M,i) ! INTRINSIC
|
|
END IF
|
|
ELSE
|
|
RWORK(i) = ZERO
|
|
END IF
|
|
END DO
|
|
DO i = 1, N
|
|
! Now, apply the same scaling to the columns of X.
|
|
IF ( RWORK(i) > ZERO ) THEN
|
|
CALL ZDSCAL( M, ONE/RWORK(i), X(1,i), 1 ) ! BLAS CALL
|
|
! X(1:M,i) = (ONE/RWORK(i)) * X(1:M,i) ! INTRINSIC
|
|
ELSE IF ( RWORK(i) < ZERO ) THEN
|
|
CALL ZLASCL( 'G', 0, 0, -RWORK(i), &
|
|
ONE/DBLE(M), M, 1, X(1,i), LDX, INFO2 ) ! LAPACK CALL
|
|
ELSE IF ( ABS(X(IZAMAX(M, X(1,i),1),i )) &
|
|
/= ZERO ) THEN
|
|
! Y(:,i) is zero vector. If X(:,i) is not
|
|
! zero, then a warning flag is raised.
|
|
! The computation continues but the
|
|
! situation will be reported in the output.
|
|
BADXY = .TRUE.
|
|
END IF
|
|
END DO
|
|
END IF
|
|
!
|
|
! <2> SVD of the data snapshot matrix X.
|
|
! =====================================
|
|
! The left singular vectors are stored in the array X.
|
|
! The right singular vectors are in the array W.
|
|
! The array W will later on contain the eigenvectors
|
|
! of a Rayleigh quotient.
|
|
NUMRNK = N
|
|
SELECT CASE ( WHTSVD )
|
|
CASE (1)
|
|
CALL ZGESVD( 'O', 'S', M, N, X, LDX, RWORK, B, &
|
|
LDB, W, LDW, ZWORK, LZWORK, RWORK(N+1), INFO1 ) ! LAPACK CALL
|
|
T_OR_N = 'C'
|
|
CASE (2)
|
|
CALL ZGESDD( 'O', M, N, X, LDX, RWORK, B, LDB, W, &
|
|
LDW, ZWORK, LZWORK, RWORK(N+1), IWORK, INFO1 ) ! LAPACK CALL
|
|
T_OR_N = 'C'
|
|
CASE (3)
|
|
CALL ZGESVDQ( 'H', 'P', 'N', 'R', 'R', M, N, &
|
|
X, LDX, RWORK, Z, LDZ, W, LDW, &
|
|
NUMRNK, IWORK, LIWORK, ZWORK, &
|
|
LZWORK, RWORK(N+1), LRWORK-N, INFO1) ! LAPACK CALL
|
|
CALL ZLACPY( 'A', M, NUMRNK, Z, LDZ, X, LDX ) ! LAPACK CALL
|
|
T_OR_N = 'C'
|
|
CASE (4)
|
|
CALL ZGEJSV( 'F', 'U', JSVOPT, 'R', 'N', 'P', M, &
|
|
N, X, LDX, RWORK, Z, LDZ, W, LDW, &
|
|
ZWORK, LZWORK, RWORK(N+1), LRWORK-N, IWORK, INFO1 ) ! LAPACK CALL
|
|
CALL ZLACPY( 'A', M, N, Z, LDZ, X, LDX ) ! LAPACK CALL
|
|
T_OR_N = 'N'
|
|
XSCL1 = RWORK(N+1)
|
|
XSCL2 = RWORK(N+2)
|
|
IF ( XSCL1 /= XSCL2 ) THEN
|
|
! This is an exceptional situation. If the
|
|
! data matrices are not scaled and the
|
|
! largest singular value of X overflows.
|
|
! In that case ZGEJSV can return the SVD
|
|
! in scaled form. The scaling factor can be used
|
|
! to rescale the data (X and Y).
|
|
CALL ZLASCL( 'G', 0, 0, XSCL1, XSCL2, M, N, Y, LDY, INFO2 )
|
|
END IF
|
|
END SELECT
|
|
!
|
|
IF ( INFO1 > 0 ) THEN
|
|
! The SVD selected subroutine did not converge.
|
|
! Return with an error code.
|
|
INFO = 2
|
|
RETURN
|
|
END IF
|
|
!
|
|
IF ( RWORK(1) == ZERO ) THEN
|
|
! The largest computed singular value of (scaled)
|
|
! X is zero. Return error code -8
|
|
! (the 8th input variable had an illegal value).
|
|
K = 0
|
|
INFO = -8
|
|
CALL XERBLA('ZGEDMD',-INFO)
|
|
RETURN
|
|
END IF
|
|
!
|
|
!<3> Determine the numerical rank of the data
|
|
! snapshots matrix X. This depends on the
|
|
! parameters NRNK and TOL.
|
|
|
|
SELECT CASE ( NRNK )
|
|
CASE ( -1 )
|
|
K = 1
|
|
DO i = 2, NUMRNK
|
|
IF ( ( RWORK(i) <= RWORK(1)*TOL ) .OR. &
|
|
( RWORK(i) <= SMALL ) ) EXIT
|
|
K = K + 1
|
|
END DO
|
|
CASE ( -2 )
|
|
K = 1
|
|
DO i = 1, NUMRNK-1
|
|
IF ( ( RWORK(i+1) <= RWORK(i)*TOL ) .OR. &
|
|
( RWORK(i) <= SMALL ) ) EXIT
|
|
K = K + 1
|
|
END DO
|
|
CASE DEFAULT
|
|
K = 1
|
|
DO i = 2, NRNK
|
|
IF ( RWORK(i) <= SMALL ) EXIT
|
|
K = K + 1
|
|
END DO
|
|
END SELECT
|
|
! Now, U = X(1:M,1:K) is the SVD/POD basis for the
|
|
! snapshot data in the input matrix X.
|
|
|
|
!<4> Compute the Rayleigh quotient S = U^H * A * U.
|
|
! Depending on the requested outputs, the computation
|
|
! is organized to compute additional auxiliary
|
|
! matrices (for the residuals and refinements).
|
|
!
|
|
! In all formulas below, we need V_k*Sigma_k^(-1)
|
|
! where either V_k is in W(1:N,1:K), or V_k^H is in
|
|
! W(1:K,1:N). Here Sigma_k=diag(WORK(1:K)).
|
|
IF ( LSAME(T_OR_N, 'N') ) THEN
|
|
DO i = 1, K
|
|
CALL ZDSCAL( N, ONE/RWORK(i), W(1,i), 1 ) ! BLAS CALL
|
|
! W(1:N,i) = (ONE/RWORK(i)) * W(1:N,i) ! INTRINSIC
|
|
END DO
|
|
ELSE
|
|
! This non-unit stride access is due to the fact
|
|
! that ZGESVD, ZGESVDQ and ZGESDD return the
|
|
! adjoint matrix of the right singular vectors.
|
|
!DO i = 1, K
|
|
! CALL ZDSCAL( N, ONE/RWORK(i), W(i,1), LDW ) ! BLAS CALL
|
|
! ! W(i,1:N) = (ONE/RWORK(i)) * W(i,1:N) ! INTRINSIC
|
|
!END DO
|
|
DO i = 1, K
|
|
RWORK(N+i) = ONE/RWORK(i)
|
|
END DO
|
|
DO j = 1, N
|
|
DO i = 1, K
|
|
W(i,j) = CMPLX(RWORK(N+i),ZERO,KIND=WP)*W(i,j)
|
|
END DO
|
|
END DO
|
|
END IF
|
|
!
|
|
IF ( WNTREF ) THEN
|
|
!
|
|
! Need A*U(:,1:K)=Y*V_k*inv(diag(WORK(1:K)))
|
|
! for computing the refined Ritz vectors
|
|
! (optionally, outside ZGEDMD).
|
|
CALL ZGEMM( 'N', T_OR_N, M, K, N, ZONE, Y, LDY, W, &
|
|
LDW, ZZERO, Z, LDZ ) ! BLAS CALL
|
|
! Z(1:M,1:K)=MATMUL(Y(1:M,1:N),TRANSPOSE(CONJG(W(1:K,1:N)))) ! INTRINSIC, for T_OR_N=='C'
|
|
! Z(1:M,1:K)=MATMUL(Y(1:M,1:N),W(1:N,1:K)) ! INTRINSIC, for T_OR_N=='N'
|
|
!
|
|
! At this point Z contains
|
|
! A * U(:,1:K) = Y * V_k * Sigma_k^(-1), and
|
|
! this is needed for computing the residuals.
|
|
! This matrix is returned in the array B and
|
|
! it can be used to compute refined Ritz vectors.
|
|
CALL ZLACPY( 'A', M, K, Z, LDZ, B, LDB ) ! BLAS CALL
|
|
! B(1:M,1:K) = Z(1:M,1:K) ! INTRINSIC
|
|
|
|
CALL ZGEMM( 'C', 'N', K, K, M, ZONE, X, LDX, Z, &
|
|
LDZ, ZZERO, S, LDS ) ! BLAS CALL
|
|
! S(1:K,1:K) = MATMUL(TRANSPOSE(CONJG(X(1:M,1:K))),Z(1:M,1:K)) ! INTRINSIC
|
|
! At this point S = U^H * A * U is the Rayleigh quotient.
|
|
ELSE
|
|
! A * U(:,1:K) is not explicitly needed and the
|
|
! computation is organized differently. The Rayleigh
|
|
! quotient is computed more efficiently.
|
|
CALL ZGEMM( 'C', 'N', K, N, M, ZONE, X, LDX, Y, LDY, &
|
|
ZZERO, Z, LDZ ) ! BLAS CALL
|
|
! Z(1:K,1:N) = MATMUL( TRANSPOSE(CONJG(X(1:M,1:K))), Y(1:M,1:N) ) ! INTRINSIC
|
|
!
|
|
CALL ZGEMM( 'N', T_OR_N, K, K, N, ZONE, Z, LDZ, W, &
|
|
LDW, ZZERO, S, LDS ) ! BLAS CALL
|
|
! S(1:K,1:K) = MATMUL(Z(1:K,1:N),TRANSPOSE(CONJG(W(1:K,1:N)))) ! INTRINSIC, for T_OR_N=='T'
|
|
! S(1:K,1:K) = MATMUL(Z(1:K,1:N),(W(1:N,1:K))) ! INTRINSIC, for T_OR_N=='N'
|
|
! At this point S = U^H * A * U is the Rayleigh quotient.
|
|
! If the residuals are requested, save scaled V_k into Z.
|
|
! Recall that V_k or V_k^H is stored in W.
|
|
IF ( WNTRES .OR. WNTEX ) THEN
|
|
IF ( LSAME(T_OR_N, 'N') ) THEN
|
|
CALL ZLACPY( 'A', N, K, W, LDW, Z, LDZ )
|
|
ELSE
|
|
CALL ZLACPY( 'A', K, N, W, LDW, Z, LDZ )
|
|
END IF
|
|
END IF
|
|
END IF
|
|
!
|
|
!<5> Compute the Ritz values and (if requested) the
|
|
! right eigenvectors of the Rayleigh quotient.
|
|
!
|
|
CALL ZGEEV( 'N', JOBZL, K, S, LDS, EIGS, W, LDW, &
|
|
W, LDW, ZWORK, LZWORK, RWORK(N+1), INFO1 ) ! LAPACK CALL
|
|
!
|
|
! W(1:K,1:K) contains the eigenvectors of the Rayleigh
|
|
! quotient. See the description of Z.
|
|
! Also, see the description of ZGEEV.
|
|
IF ( INFO1 > 0 ) THEN
|
|
! ZGEEV failed to compute the eigenvalues and
|
|
! eigenvectors of the Rayleigh quotient.
|
|
INFO = 3
|
|
RETURN
|
|
END IF
|
|
!
|
|
! <6> Compute the eigenvectors (if requested) and,
|
|
! the residuals (if requested).
|
|
!
|
|
IF ( WNTVEC .OR. WNTEX ) THEN
|
|
IF ( WNTRES ) THEN
|
|
IF ( WNTREF ) THEN
|
|
! Here, if the refinement is requested, we have
|
|
! A*U(:,1:K) already computed and stored in Z.
|
|
! For the residuals, need Y = A * U(:,1;K) * W.
|
|
CALL ZGEMM( 'N', 'N', M, K, K, ZONE, Z, LDZ, W, &
|
|
LDW, ZZERO, Y, LDY ) ! BLAS CALL
|
|
! Y(1:M,1:K) = Z(1:M,1:K) * W(1:K,1:K) ! INTRINSIC
|
|
! This frees Z; Y contains A * U(:,1:K) * W.
|
|
ELSE
|
|
! Compute S = V_k * Sigma_k^(-1) * W, where
|
|
! V_k * Sigma_k^(-1) (or its adjoint) is stored in Z
|
|
CALL ZGEMM( T_OR_N, 'N', N, K, K, ZONE, Z, LDZ, &
|
|
W, LDW, ZZERO, S, LDS )
|
|
! Then, compute Z = Y * S =
|
|
! = Y * V_k * Sigma_k^(-1) * W(1:K,1:K) =
|
|
! = A * U(:,1:K) * W(1:K,1:K)
|
|
CALL ZGEMM( 'N', 'N', M, K, N, ZONE, Y, LDY, S, &
|
|
LDS, ZZERO, Z, LDZ )
|
|
! Save a copy of Z into Y and free Z for holding
|
|
! the Ritz vectors.
|
|
CALL ZLACPY( 'A', M, K, Z, LDZ, Y, LDY )
|
|
IF ( WNTEX ) CALL ZLACPY( 'A', M, K, Z, LDZ, B, LDB )
|
|
END IF
|
|
ELSE IF ( WNTEX ) THEN
|
|
! Compute S = V_k * Sigma_k^(-1) * W, where
|
|
! V_k * Sigma_k^(-1) is stored in Z
|
|
CALL ZGEMM( T_OR_N, 'N', N, K, K, ZONE, Z, LDZ, &
|
|
W, LDW, ZZERO, S, LDS )
|
|
! Then, compute Z = Y * S =
|
|
! = Y * V_k * Sigma_k^(-1) * W(1:K,1:K) =
|
|
! = A * U(:,1:K) * W(1:K,1:K)
|
|
CALL ZGEMM( 'N', 'N', M, K, N, ZONE, Y, LDY, S, &
|
|
LDS, ZZERO, B, LDB )
|
|
! The above call replaces the following two calls
|
|
! that were used in the developing-testing phase.
|
|
! CALL ZGEMM( 'N', 'N', M, K, N, ZONE, Y, LDY, S, &
|
|
! LDS, ZZERO, Z, LDZ)
|
|
! Save a copy of Z into B and free Z for holding
|
|
! the Ritz vectors.
|
|
! CALL ZLACPY( 'A', M, K, Z, LDZ, B, LDB )
|
|
END IF
|
|
!
|
|
! Compute the Ritz vectors
|
|
IF ( WNTVEC ) CALL ZGEMM( 'N', 'N', M, K, K, ZONE, X, LDX, W, LDW, &
|
|
ZZERO, Z, LDZ ) ! BLAS CALL
|
|
! Z(1:M,1:K) = MATMUL(X(1:M,1:K), W(1:K,1:K)) ! INTRINSIC
|
|
!
|
|
IF ( WNTRES ) THEN
|
|
DO i = 1, K
|
|
CALL ZAXPY( M, -EIGS(i), Z(1,i), 1, Y(1,i), 1 ) ! BLAS CALL
|
|
! Y(1:M,i) = Y(1:M,i) - EIGS(i) * Z(1:M,i) ! INTRINSIC
|
|
RES(i) = DZNRM2( M, Y(1,i), 1 ) ! BLAS CALL
|
|
END DO
|
|
END IF
|
|
END IF
|
|
!
|
|
IF ( WHTSVD == 4 ) THEN
|
|
RWORK(N+1) = XSCL1
|
|
RWORK(N+2) = XSCL2
|
|
END IF
|
|
!
|
|
! Successful exit.
|
|
IF ( .NOT. BADXY ) THEN
|
|
INFO = 0
|
|
ELSE
|
|
! A warning on possible data inconsistency.
|
|
! This should be a rare event.
|
|
INFO = 4
|
|
END IF
|
|
!............................................................
|
|
RETURN
|
|
! ......
|
|
END SUBROUTINE ZGEDMD
|
|
|