1055 lines
45 KiB
Fortran
1055 lines
45 KiB
Fortran
SUBROUTINE DGEDMD( JOBS, JOBZ, JOBR, JOBF, WHTSVD, &
|
|
M, N, X, LDX, Y, LDY, NRNK, TOL, &
|
|
K, REIG, IMEIG, Z, LDZ, RES, &
|
|
B, LDB, W, LDW, S, LDS, &
|
|
WORK, LWORK, IWORK, LIWORK, INFO )
|
|
! March 2023
|
|
!.....
|
|
USE iso_fortran_env
|
|
IMPLICIT NONE
|
|
INTEGER, PARAMETER :: WP = real64
|
|
!.....
|
|
! Scalar arguments
|
|
CHARACTER, INTENT(IN) :: JOBS, JOBZ, JOBR, JOBF
|
|
INTEGER, INTENT(IN) :: WHTSVD, M, N, LDX, LDY, &
|
|
NRNK, LDZ, LDB, LDW, LDS, &
|
|
LWORK, LIWORK
|
|
INTEGER, INTENT(OUT) :: K, INFO
|
|
REAL(KIND=WP), INTENT(IN) :: TOL
|
|
! Array arguments
|
|
REAL(KIND=WP), INTENT(INOUT) :: X(LDX,*), Y(LDY,*)
|
|
REAL(KIND=WP), INTENT(OUT) :: Z(LDZ,*), B(LDB,*), &
|
|
W(LDW,*), S(LDS,*)
|
|
REAL(KIND=WP), INTENT(OUT) :: REIG(*), IMEIG(*), &
|
|
RES(*)
|
|
REAL(KIND=WP), INTENT(OUT) :: WORK(*)
|
|
INTEGER, INTENT(OUT) :: IWORK(*)
|
|
!............................................................
|
|
! Purpose
|
|
! =======
|
|
! DGEDMD 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, DGEDMD 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, DGEDMD 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 :: DGESVD (the QR SVD algorithm)
|
|
! 2 :: DGESDD (the Divide and Conquer algorithm; if enough
|
|
! workspace available, this is the fastest option)
|
|
! 3 :: DGESVDQ (the preconditioned QR SVD ; this and 4
|
|
! are the most accurate options)
|
|
! 4 :: DGEJSV (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) REAL(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) REAL(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.
|
|
!.....
|
|
! REIG (output) REAL(KIND=WP) N-by-1 array
|
|
! The leading K (K<=N) entries of REIG contain
|
|
! the real parts of the computed eigenvalues
|
|
! REIG(1:K) + sqrt(-1)*IMEIG(1:K).
|
|
! See the descriptions of K, IMEIG, and Z.
|
|
!.....
|
|
! IMEIG (output) REAL(KIND=WP) N-by-1 array
|
|
! The leading K (K<=N) entries of IMEIG contain
|
|
! the imaginary parts of the computed eigenvalues
|
|
! REIG(1:K) + sqrt(-1)*IMEIG(1:K).
|
|
! The eigenvalues are determined as follows:
|
|
! If IMEIG(i) == 0, then the corresponding eigenvalue is
|
|
! real, LAMBDA(i) = REIG(i).
|
|
! If IMEIG(i)>0, then the corresponding complex
|
|
! conjugate pair of eigenvalues reads
|
|
! LAMBDA(i) = REIG(i) + sqrt(-1)*IMAG(i)
|
|
! LAMBDA(i+1) = REIG(i) - sqrt(-1)*IMAG(i)
|
|
! That is, complex conjugate pairs have consecutive
|
|
! indices (i,i+1), with the positive imaginary part
|
|
! listed first.
|
|
! See the descriptions of K, REIG, and Z.
|
|
!.....
|
|
! Z (workspace/output) REAL(KIND=WP) M-by-N array
|
|
! If JOBZ =='V' then
|
|
! Z contains real Ritz vectors as follows:
|
|
! If IMEIG(i)=0, then Z(:,i) is an eigenvector of
|
|
! the i-th Ritz value; ||Z(:,i)||_2=1.
|
|
! If IMEIG(i) > 0 (and IMEIG(i+1) < 0) then
|
|
! [Z(:,i) Z(:,i+1)] span an invariant subspace and
|
|
! the Ritz values extracted from this subspace are
|
|
! REIG(i) + sqrt(-1)*IMEIG(i) and
|
|
! REIG(i) - sqrt(-1)*IMEIG(i).
|
|
! The corresponding eigenvectors are
|
|
! Z(:,i) + sqrt(-1)*Z(:,i+1) and
|
|
! Z(:,i) - sqrt(-1)*Z(:,i+1), respectively.
|
|
! || Z(:,i:i+1)||_F = 1.
|
|
! If JOBZ == 'F', then the above descriptions hold for
|
|
! the columns of X(:,1:K)*W(1:K,1:K), where the columns
|
|
! of W(1:k,1:K) are the computed eigenvectors of the
|
|
! K-by-K Rayleigh quotient. The columns of W(1:K,1:K)
|
|
! are similarly structured: If IMEIG(i) == 0 then
|
|
! X(:,1:K)*W(:,i) is an eigenvector, and if IMEIG(i)>0
|
|
! then X(:,1:K)*W(:,i)+sqrt(-1)*X(:,1:K)*W(:,i+1) and
|
|
! X(:,1:K)*W(:,i)-sqrt(-1)*X(:,1:K)*W(:,i+1)
|
|
! are the eigenvectors of LAMBDA(i), LAMBDA(i+1).
|
|
! See the descriptions of REIG, IMEIG, 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.
|
|
! If LAMBDA(i) is real, then
|
|
! RES(i) = || A * Z(:,i) - LAMBDA(i)*Z(:,i))||_2.
|
|
! If [LAMBDA(i), LAMBDA(i+1)] is a complex conjugate pair
|
|
! then
|
|
! RES(i)=RES(i+1) = || A * Z(:,i:i+1) - Z(:,i:i+1) *B||_F
|
|
! where B = [ real(LAMBDA(i)) imag(LAMBDA(i)) ]
|
|
! [-imag(LAMBDA(i)) real(LAMBDA(i)) ].
|
|
! It holds that
|
|
! RES(i) = || A*ZC(:,i) - LAMBDA(i) *ZC(:,i) ||_2
|
|
! RES(i+1) = || A*ZC(:,i+1) - LAMBDA(i+1)*ZC(:,i+1) ||_2
|
|
! where ZC(:,i) = Z(:,i) + sqrt(-1)*Z(:,i+1)
|
|
! ZC(:,i+1) = Z(:,i) - sqrt(-1)*Z(:,i+1)
|
|
! See the description of REIG, IMEIG and Z.
|
|
!.....
|
|
! B (output) REAL(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) REAL(KIND=WP) N-by-N array
|
|
! On exit, W(1:K,1:K) contains the K computed
|
|
! eigenvectors of the matrix Rayleigh quotient (real and
|
|
! imaginary parts for each complex conjugate pair of the
|
|
! eigenvalues). 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) REAL(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 DGEEV.
|
|
! See the description of K.
|
|
!.....
|
|
! LDS (input) INTEGER, LDS >= N
|
|
! The leading dimension of the array S.
|
|
!.....
|
|
! WORK (workspace/output) REAL(KIND=WP) LWORK-by-1 array
|
|
! On exit, WORK(1:N) contains the singular values of
|
|
! X (for JOBS=='N') or column scaled X (JOBS=='S', 'C').
|
|
! If WHTSVD==4, then WORK(N+1) and WORK(N+2) contain
|
|
! scaling factor WORK(N+2)/WORK(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 DGEDMD is only workspace query, then
|
|
! WORK(1) contains the minimal workspace length and
|
|
! WORK(2) is the optimal workspace length. Hence, the
|
|
! leng of work is at least 2.
|
|
! See the description of LWORK.
|
|
!.....
|
|
! LWORK (input) INTEGER
|
|
! The minimal length of the workspace vector WORK.
|
|
! LWORK is calculated as follows:
|
|
! If WHTSVD == 1 ::
|
|
! If JOBZ == 'V', then
|
|
! LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,4*N)).
|
|
! If JOBZ == 'N' then
|
|
! LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,3*N)).
|
|
! Here LWORK_SVD = MAX(1,3*N+M,5*N) is the minimal
|
|
! workspace length of DGESVD.
|
|
! If WHTSVD == 2 ::
|
|
! If JOBZ == 'V', then
|
|
! LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,4*N))
|
|
! If JOBZ == 'N', then
|
|
! LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,3*N))
|
|
! Here LWORK_SVD = MAX(M, 5*N*N+4*N)+3*N*N is the
|
|
! minimal workspace length of DGESDD.
|
|
! If WHTSVD == 3 ::
|
|
! If JOBZ == 'V', then
|
|
! LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,4*N))
|
|
! If JOBZ == 'N', then
|
|
! LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,3*N))
|
|
! Here LWORK_SVD = N+M+MAX(3*N+1,
|
|
! MAX(1,3*N+M,5*N),MAX(1,N))
|
|
! is the minimal workspace length of DGESVDQ.
|
|
! If WHTSVD == 4 ::
|
|
! If JOBZ == 'V', then
|
|
! LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,4*N))
|
|
! If JOBZ == 'N', then
|
|
! LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,3*N))
|
|
! Here LWORK_SVD = MAX(7,2*M+N,6*N+2*N*N) is the
|
|
! minimal workspace length of DGEJSV.
|
|
! The above expressions are not simplified in order to
|
|
! make the usage of WORK more transparent, and for
|
|
! easier checking. In any case, LWORK >= 2.
|
|
! If on entry LWORK = -1, then a workspace query is
|
|
! assumed and the procedure only computes the minimal
|
|
! and the optimal workspace lengths for both WORK and
|
|
! IWORK. See the descriptions of WORK and IWORK.
|
|
!.....
|
|
! 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 both WORK and
|
|
! IWORK. See the descriptions of WORK 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
|
|
|
|
! Local scalars
|
|
! ~~~~~~~~~~~~~
|
|
REAL(KIND=WP) :: OFL, ROOTSC, SCALE, SMALL, &
|
|
SSUM, XSCL1, XSCL2
|
|
INTEGER :: i, j, IMINWR, INFO1, INFO2, &
|
|
LWRKEV, LWRSDD, LWRSVD, &
|
|
LWRSVQ, MLWORK, MWRKEV, MWRSDD, &
|
|
MWRSVD, MWRSVJ, MWRSVQ, NUMRNK, &
|
|
OLWORK
|
|
LOGICAL :: BADXY, LQUERY, SCCOLX, SCCOLY, &
|
|
WNTEX, WNTREF, WNTRES, WNTVEC
|
|
CHARACTER :: JOBZL, T_OR_N
|
|
CHARACTER :: JSVOPT
|
|
|
|
! Local arrays
|
|
! ~~~~~~~~~~~~
|
|
REAL(KIND=WP) :: AB(2,2), RDUMMY(2), RDUMMY2(2)
|
|
! External functions (BLAS and LAPACK)
|
|
! ~~~~~~~~~~~~~~~~~
|
|
REAL(KIND=WP) DLANGE, DLAMCH, DNRM2
|
|
EXTERNAL DLANGE, DLAMCH, DNRM2, IDAMAX
|
|
INTEGER IDAMAX
|
|
LOGICAL DISNAN, LSAME
|
|
EXTERNAL DISNAN, LSAME
|
|
|
|
! External subroutines (BLAS and LAPACK)
|
|
! ~~~~~~~~~~~~~~~~~~~~
|
|
EXTERNAL DAXPY, DGEMM, DSCAL
|
|
EXTERNAL DGEEV, DGEJSV, DGESDD, DGESVD, DGESVDQ, &
|
|
DLACPY, DLASCL, DLASSQ, 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 = ( ( LWORK == -1 ) .OR. ( LIWORK == -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 = -18
|
|
ELSE IF ( (WNTREF .OR. WNTEX ) .AND. ( LDB < M ) ) THEN
|
|
INFO = -21
|
|
ELSE IF ( LDW < N ) THEN
|
|
INFO = -23
|
|
ELSE IF ( LDS < N ) THEN
|
|
INFO = -25
|
|
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
|
|
WORK(1) = 2
|
|
WORK(2) = 2
|
|
ELSE
|
|
K = 0
|
|
END IF
|
|
INFO = 1
|
|
RETURN
|
|
END IF
|
|
MLWORK = MAX(2,N)
|
|
OLWORK = MAX(2,N)
|
|
IMINWR = 1
|
|
SELECT CASE ( WHTSVD )
|
|
CASE (1)
|
|
! The following is specified as the minimal
|
|
! length of WORK in the definition of DGESVD:
|
|
! MWRSVD = MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N))
|
|
MWRSVD = MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N))
|
|
MLWORK = MAX(MLWORK,N + MWRSVD)
|
|
IF ( LQUERY ) THEN
|
|
CALL DGESVD( 'O', 'S', M, N, X, LDX, WORK, &
|
|
B, LDB, W, LDW, RDUMMY, -1, INFO1 )
|
|
LWRSVD = MAX( MWRSVD, INT( RDUMMY(1) ) )
|
|
OLWORK = MAX(OLWORK,N + LWRSVD)
|
|
END IF
|
|
CASE (2)
|
|
! The following is specified as the minimal
|
|
! length of WORK in the definition of DGESDD:
|
|
! MWRSDD = 3*MIN(M,N)*MIN(M,N) +
|
|
! MAX( MAX(M,N),5*MIN(M,N)*MIN(M,N)+4*MIN(M,N) )
|
|
! IMINWR = 8*MIN(M,N)
|
|
MWRSDD = 3*MIN(M,N)*MIN(M,N) + &
|
|
MAX( MAX(M,N),5*MIN(M,N)*MIN(M,N)+4*MIN(M,N) )
|
|
MLWORK = MAX(MLWORK,N + MWRSDD)
|
|
IMINWR = 8*MIN(M,N)
|
|
IF ( LQUERY ) THEN
|
|
CALL DGESDD( 'O', M, N, X, LDX, WORK, B, &
|
|
LDB, W, LDW, RDUMMY, -1, IWORK, INFO1 )
|
|
LWRSDD = MAX( MWRSDD, INT( RDUMMY(1) ) )
|
|
OLWORK = MAX(OLWORK,N + LWRSDD)
|
|
END IF
|
|
CASE (3)
|
|
!LWQP3 = 3*N+1
|
|
!LWORQ = MAX(N, 1)
|
|
!MWRSVD = MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N))
|
|
!MWRSVQ = N + MAX( LWQP3, MWRSVD, LWORQ ) + MAX(M,2)
|
|
!MLWORK = N + MWRSVQ
|
|
!IMINWR = M+N-1
|
|
CALL DGESVDQ( 'H', 'P', 'N', 'R', 'R', M, N, &
|
|
X, LDX, WORK, Z, LDZ, W, LDW, &
|
|
NUMRNK, IWORK, LIWORK, RDUMMY, &
|
|
-1, RDUMMY2, -1, INFO1 )
|
|
IMINWR = IWORK(1)
|
|
MWRSVQ = INT(RDUMMY(2))
|
|
MLWORK = MAX(MLWORK,N+MWRSVQ+INT(RDUMMY2(1)))
|
|
IF ( LQUERY ) THEN
|
|
LWRSVQ = MAX( MWRSVQ, INT(RDUMMY(1)) )
|
|
OLWORK = MAX(OLWORK,N+LWRSVQ+INT(RDUMMY2(1)))
|
|
END IF
|
|
CASE (4)
|
|
JSVOPT = 'J'
|
|
!MWRSVJ = MAX( 7, 2*M+N, 6*N+2*N*N ) ! for JSVOPT='V'
|
|
MWRSVJ = MAX( 7, 2*M+N, 4*N+N*N, 2*N+N*N+6 )
|
|
MLWORK = MAX(MLWORK,N+MWRSVJ)
|
|
IMINWR = MAX( 3, M+3*N )
|
|
IF ( LQUERY ) THEN
|
|
OLWORK = MAX(OLWORK,N+MWRSVJ)
|
|
END IF
|
|
END SELECT
|
|
IF ( WNTVEC .OR. WNTEX .OR. LSAME(JOBZ,'F') ) THEN
|
|
JOBZL = 'V'
|
|
ELSE
|
|
JOBZL = 'N'
|
|
END IF
|
|
! Workspace calculation to the DGEEV call
|
|
IF ( LSAME(JOBZL,'V') ) THEN
|
|
MWRKEV = MAX( 1, 4*N )
|
|
ELSE
|
|
MWRKEV = MAX( 1, 3*N )
|
|
END IF
|
|
MLWORK = MAX(MLWORK,N+MWRKEV)
|
|
IF ( LQUERY ) THEN
|
|
CALL DGEEV( 'N', JOBZL, N, S, LDS, REIG, &
|
|
IMEIG, W, LDW, W, LDW, RDUMMY, -1, INFO1 )
|
|
LWRKEV = MAX( MWRKEV, INT(RDUMMY(1)) )
|
|
OLWORK = MAX( OLWORK, N+LWRKEV )
|
|
END IF
|
|
!
|
|
IF ( LIWORK < IMINWR .AND. (.NOT.LQUERY) ) INFO = -29
|
|
IF ( LWORK < MLWORK .AND. (.NOT.LQUERY) ) INFO = -27
|
|
END IF
|
|
!
|
|
IF( INFO /= 0 ) THEN
|
|
CALL XERBLA( 'DGEDMD', -INFO )
|
|
RETURN
|
|
ELSE IF ( LQUERY ) THEN
|
|
! Return minimal and optimal workspace sizes
|
|
IWORK(1) = IMINWR
|
|
WORK(1) = MLWORK
|
|
WORK(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 DLASSQ.
|
|
K = 0
|
|
DO i = 1, N
|
|
!WORK(i) = DNRM2( M, X(1,i), 1 )
|
|
SCALE = ZERO
|
|
CALL DLASSQ( M, X(1,i), 1, SCALE, SSUM )
|
|
IF ( DISNAN(SCALE) .OR. DISNAN(SSUM) ) THEN
|
|
K = 0
|
|
INFO = -8
|
|
CALL XERBLA('DGEDMD',-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 WORK(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 DLASCL( 'G', 0, 0, SCALE, ONE/ROOTSC, &
|
|
M, 1, X(1,i), M, INFO2 )
|
|
WORK(i) = - SCALE * ( ROOTSC / DBLE(M) )
|
|
ELSE
|
|
! X(:,i) will be scaled to unit 2-norm
|
|
WORK(i) = SCALE * ROOTSC
|
|
CALL DLASCL( 'G',0, 0, WORK(i), ONE, M, 1, &
|
|
X(1,i), M, INFO2 ) ! LAPACK CALL
|
|
! X(1:M,i) = (ONE/WORK(i)) * X(1:M,i) ! INTRINSIC
|
|
END IF
|
|
ELSE
|
|
WORK(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('DGEDMD',-INFO)
|
|
RETURN
|
|
END IF
|
|
DO i = 1, N
|
|
! Now, apply the same scaling to the columns of Y.
|
|
IF ( WORK(i) > ZERO ) THEN
|
|
CALL DSCAL( M, ONE/WORK(i), Y(1,i), 1 ) ! BLAS CALL
|
|
! Y(1:M,i) = (ONE/WORK(i)) * Y(1:M,i) ! INTRINSIC
|
|
ELSE IF ( WORK(i) < ZERO ) THEN
|
|
CALL DLASCL( 'G', 0, 0, -WORK(i), &
|
|
ONE/DBLE(M), M, 1, Y(1,i), M, INFO2 ) ! LAPACK CALL
|
|
ELSE IF ( Y(IDAMAX(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 DSCAL( 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 DLASSQ.
|
|
DO i = 1, N
|
|
!WORK(i) = DNRM2( M, Y(1,i), 1 )
|
|
SCALE = ZERO
|
|
CALL DLASSQ( M, Y(1,i), 1, SCALE, SSUM )
|
|
IF ( DISNAN(SCALE) .OR. DISNAN(SSUM) ) THEN
|
|
K = 0
|
|
INFO = -10
|
|
CALL XERBLA('DGEDMD',-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 WORK(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 DLASCL( 'G', 0, 0, SCALE, ONE/ROOTSC, &
|
|
M, 1, Y(1,i), M, INFO2 )
|
|
WORK(i) = - SCALE * ( ROOTSC / DBLE(M) )
|
|
ELSE
|
|
! X(:,i) will be scaled to unit 2-norm
|
|
WORK(i) = SCALE * ROOTSC
|
|
CALL DLASCL( 'G',0, 0, WORK(i), ONE, M, 1, &
|
|
Y(1,i), M, INFO2 ) ! LAPACK CALL
|
|
! Y(1:M,i) = (ONE/WORK(i)) * Y(1:M,i) ! INTRINSIC
|
|
END IF
|
|
ELSE
|
|
WORK(i) = ZERO
|
|
END IF
|
|
END DO
|
|
DO i = 1, N
|
|
! Now, apply the same scaling to the columns of X.
|
|
IF ( WORK(i) > ZERO ) THEN
|
|
CALL DSCAL( M, ONE/WORK(i), X(1,i), 1 ) ! BLAS CALL
|
|
! X(1:M,i) = (ONE/WORK(i)) * X(1:M,i) ! INTRINSIC
|
|
ELSE IF ( WORK(i) < ZERO ) THEN
|
|
CALL DLASCL( 'G', 0, 0, -WORK(i), &
|
|
ONE/DBLE(M), M, 1, X(1,i), M, INFO2 ) ! LAPACK CALL
|
|
ELSE IF ( X(IDAMAX(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 DGESVD( 'O', 'S', M, N, X, LDX, WORK, B, &
|
|
LDB, W, LDW, WORK(N+1), LWORK-N, INFO1 ) ! LAPACK CALL
|
|
T_OR_N = 'T'
|
|
CASE (2)
|
|
CALL DGESDD( 'O', M, N, X, LDX, WORK, B, LDB, W, &
|
|
LDW, WORK(N+1), LWORK-N, IWORK, INFO1 ) ! LAPACK CALL
|
|
T_OR_N = 'T'
|
|
CASE (3)
|
|
CALL DGESVDQ( 'H', 'P', 'N', 'R', 'R', M, N, &
|
|
X, LDX, WORK, Z, LDZ, W, LDW, &
|
|
NUMRNK, IWORK, LIWORK, WORK(N+MAX(2,M)+1),&
|
|
LWORK-N-MAX(2,M), WORK(N+1), MAX(2,M), INFO1) ! LAPACK CALL
|
|
CALL DLACPY( 'A', M, NUMRNK, Z, LDZ, X, LDX ) ! LAPACK CALL
|
|
T_OR_N = 'T'
|
|
CASE (4)
|
|
CALL DGEJSV( 'F', 'U', JSVOPT, 'N', 'N', 'P', M, &
|
|
N, X, LDX, WORK, Z, LDZ, W, LDW, &
|
|
WORK(N+1), LWORK-N, IWORK, INFO1 ) ! LAPACK CALL
|
|
CALL DLACPY( 'A', M, N, Z, LDZ, X, LDX ) ! LAPACK CALL
|
|
T_OR_N = 'N'
|
|
XSCL1 = WORK(N+1)
|
|
XSCL2 = WORK(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 DGEJSV can return the SVD
|
|
! in scaled form. The scaling factor can be used
|
|
! to rescale the data (X and Y).
|
|
CALL DLASCL( '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 ( WORK(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('DGEDMD',-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 ( ( WORK(i) <= WORK(1)*TOL ) .OR. &
|
|
( WORK(i) <= SMALL ) ) EXIT
|
|
K = K + 1
|
|
END DO
|
|
CASE ( -2 )
|
|
K = 1
|
|
DO i = 1, NUMRNK-1
|
|
IF ( ( WORK(i+1) <= WORK(i)*TOL ) .OR. &
|
|
( WORK(i) <= SMALL ) ) EXIT
|
|
K = K + 1
|
|
END DO
|
|
CASE DEFAULT
|
|
K = 1
|
|
DO i = 2, NRNK
|
|
IF ( WORK(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^T * 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^T 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 DSCAL( N, ONE/WORK(i), W(1,i), 1 ) ! BLAS CALL
|
|
! W(1:N,i) = (ONE/WORK(i)) * W(1:N,i) ! INTRINSIC
|
|
END DO
|
|
ELSE
|
|
! This non-unit stride access is due to the fact
|
|
! that DGESVD, DGESVDQ and DGESDD return the
|
|
! transposed matrix of the right singular vectors.
|
|
!DO i = 1, K
|
|
! CALL DSCAL( N, ONE/WORK(i), W(i,1), LDW ) ! BLAS CALL
|
|
! ! W(i,1:N) = (ONE/WORK(i)) * W(i,1:N) ! INTRINSIC
|
|
!END DO
|
|
DO i = 1, K
|
|
WORK(N+i) = ONE/WORK(i)
|
|
END DO
|
|
DO j = 1, N
|
|
DO i = 1, K
|
|
W(i,j) = (WORK(N+i))*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 DGEDMD).
|
|
CALL DGEMM( 'N', T_OR_N, M, K, N, ONE, Y, LDY, W, &
|
|
LDW, ZERO, Z, LDZ ) ! BLAS CALL
|
|
! Z(1:M,1:K)=MATMUL(Y(1:M,1:N),TRANSPOSE(W(1:K,1:N))) ! INTRINSIC, for T_OR_N=='T'
|
|
! 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 DLACPY( 'A', M, K, Z, LDZ, B, LDB ) ! BLAS CALL
|
|
! B(1:M,1:K) = Z(1:M,1:K) ! INTRINSIC
|
|
|
|
CALL DGEMM( 'T', 'N', K, K, M, ONE, X, LDX, Z, &
|
|
LDZ, ZERO, S, LDS ) ! BLAS CALL
|
|
! S(1:K,1:K) = MATMUL(TANSPOSE(X(1:M,1:K)),Z(1:M,1:K)) ! INTRINSIC
|
|
! At this point S = U^T * 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 DGEMM( 'T', 'N', K, N, M, ONE, X, LDX, Y, LDY, &
|
|
ZERO, Z, LDZ ) ! BLAS CALL
|
|
! Z(1:K,1:N) = MATMUL( TRANSPOSE(X(1:M,1:K)), Y(1:M,1:N) ) ! INTRINSIC
|
|
! In the two DGEMM calls here, can use K for LDZ.
|
|
CALL DGEMM( 'N', T_OR_N, K, K, N, ONE, Z, LDZ, W, &
|
|
LDW, ZERO, S, LDS ) ! BLAS CALL
|
|
! S(1:K,1:K) = MATMUL(Z(1:K,1:N),TRANSPOSE(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^T * A * U is the Rayleigh quotient.
|
|
! If the residuals are requested, save scaled V_k into Z.
|
|
! Recall that V_k or V_k^T is stored in W.
|
|
IF ( WNTRES .OR. WNTEX ) THEN
|
|
IF ( LSAME(T_OR_N, 'N') ) THEN
|
|
CALL DLACPY( 'A', N, K, W, LDW, Z, LDZ )
|
|
ELSE
|
|
CALL DLACPY( '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 DGEEV( 'N', JOBZL, K, S, LDS, REIG, IMEIG, W, &
|
|
LDW, W, LDW, WORK(N+1), LWORK-N, INFO1 ) ! LAPACK CALL
|
|
!
|
|
! W(1:K,1:K) contains the eigenvectors of the Rayleigh
|
|
! quotient. Even in the case of complex spectrum, all
|
|
! computation is done in real arithmetic. REIG and
|
|
! IMEIG are the real and the imaginary parts of the
|
|
! eigenvalues, so that the spectrum is given as
|
|
! REIG(:) + sqrt(-1)*IMEIG(:). Complex conjugate pairs
|
|
! are listed at consecutive positions. For such a
|
|
! complex conjugate pair of the eigenvalues, the
|
|
! corresponding eigenvectors are also a complex
|
|
! conjugate pair with the real and imaginary parts
|
|
! stored column-wise in W at the corresponding
|
|
! consecutive column indices. See the description of Z.
|
|
! Also, see the description of DGEEV.
|
|
IF ( INFO1 > 0 ) THEN
|
|
! DGEEV 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 DGEMM( 'N', 'N', M, K, K, ONE, Z, LDZ, W, &
|
|
LDW, ZERO, 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) is stored in Z
|
|
CALL DGEMM( T_OR_N, 'N', N, K, K, ONE, Z, LDZ, &
|
|
W, LDW, ZERO, 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 DGEMM( 'N', 'N', M, K, N, ONE, Y, LDY, S, &
|
|
LDS, ZERO, Z, LDZ)
|
|
! Save a copy of Z into Y and free Z for holding
|
|
! the Ritz vectors.
|
|
CALL DLACPY( 'A', M, K, Z, LDZ, Y, LDY )
|
|
IF ( WNTEX ) CALL DLACPY( '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 DGEMM( T_OR_N, 'N', N, K, K, ONE, Z, LDZ, &
|
|
W, LDW, ZERO, 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 DGEMM( 'N', 'N', M, K, N, ONE, Y, LDY, S, &
|
|
LDS, ZERO, B, LDB )
|
|
! The above call replaces the following two calls
|
|
! that were used in the developing-testing phase.
|
|
! CALL DGEMM( 'N', 'N', M, K, N, ONE, Y, LDY, S, &
|
|
! LDS, ZERO, Z, LDZ)
|
|
! Save a copy of Z into B and free Z for holding
|
|
! the Ritz vectors.
|
|
! CALL DLACPY( 'A', M, K, Z, LDZ, B, LDB )
|
|
END IF
|
|
!
|
|
! Compute the real form of the Ritz vectors
|
|
IF ( WNTVEC ) CALL DGEMM( 'N', 'N', M, K, K, ONE, X, LDX, W, LDW, &
|
|
ZERO, Z, LDZ ) ! BLAS CALL
|
|
! Z(1:M,1:K) = MATMUL(X(1:M,1:K), W(1:K,1:K)) ! INTRINSIC
|
|
!
|
|
IF ( WNTRES ) THEN
|
|
i = 1
|
|
DO WHILE ( i <= K )
|
|
IF ( IMEIG(i) == ZERO ) THEN
|
|
! have a real eigenvalue with real eigenvector
|
|
CALL DAXPY( M, -REIG(i), Z(1,i), 1, Y(1,i), 1 ) ! BLAS CALL
|
|
! Y(1:M,i) = Y(1:M,i) - REIG(i) * Z(1:M,i) ! INTRINSIC
|
|
RES(i) = DNRM2( M, Y(1,i), 1) ! BLAS CALL
|
|
i = i + 1
|
|
ELSE
|
|
! Have a complex conjugate pair
|
|
! REIG(i) +- sqrt(-1)*IMEIG(i).
|
|
! Since all computation is done in real
|
|
! arithmetic, the formula for the residual
|
|
! is recast for real representation of the
|
|
! complex conjugate eigenpair. See the
|
|
! description of RES.
|
|
AB(1,1) = REIG(i)
|
|
AB(2,1) = -IMEIG(i)
|
|
AB(1,2) = IMEIG(i)
|
|
AB(2,2) = REIG(i)
|
|
CALL DGEMM( 'N', 'N', M, 2, 2, -ONE, Z(1,i), &
|
|
LDZ, AB, 2, ONE, Y(1,i), LDY ) ! BLAS CALL
|
|
! Y(1:M,i:i+1) = Y(1:M,i:i+1) - Z(1:M,i:i+1) * AB ! INTRINSIC
|
|
RES(i) = DLANGE( 'F', M, 2, Y(1,i), LDY, &
|
|
WORK(N+1) ) ! LAPACK CALL
|
|
RES(i+1) = RES(i)
|
|
i = i + 2
|
|
END IF
|
|
END DO
|
|
END IF
|
|
END IF
|
|
!
|
|
IF ( WHTSVD == 4 ) THEN
|
|
WORK(N+1) = XSCL1
|
|
WORK(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 DGEDMD
|
|
|