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