LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dgedmd()

subroutine dgedmd ( character, intent(in) jobs,
character, intent(in) jobz,
character, intent(in) jobr,
character, intent(in) jobf,
integer, intent(in) whtsvd,
integer, intent(in) m,
integer, intent(in) n,
real(kind=wp), dimension(ldx,*), intent(inout) x,
integer, intent(in) ldx,
real(kind=wp), dimension(ldy,*), intent(inout) y,
integer, intent(in) ldy,
integer, intent(in) nrnk,
real(kind=wp), intent(in) tol,
integer, intent(out) k,
real(kind=wp), dimension(*), intent(out) reig,
real(kind=wp), dimension(*), intent(out) imeig,
real(kind=wp), dimension(ldz,*), intent(out) z,
integer, intent(in) ldz,
real(kind=wp), dimension(*), intent(out) res,
real(kind=wp), dimension(ldb,*), intent(out) b,
integer, intent(in) ldb,
real(kind=wp), dimension(ldw,*), intent(out) w,
integer, intent(in) ldw,
real(kind=wp), dimension(lds,*), intent(out) s,
integer, intent(in) lds,
real(kind=wp), dimension(*), intent(out) work,
integer, intent(in) lwork,
integer, dimension(*), intent(out) iwork,
integer, intent(in) liwork,
integer, intent(out) info )

DGEDMD computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.

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  Contract No: W31P4Q-21-C-0007
!>    - DARPA PAI project  Contract No: HR0011-18-9-0033
!>    - DARPA MoDyL project 
!>    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
!>
Parameters
[in]JOBS
!>    JOBS (input) is 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.
!>
[in]JOBZ
!>    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.
!>
[in]JOBR
!>    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.
!>
[in]JOBF
!>    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.
!>
[in]WHTSVD
!>    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.
!>
[in]M
!>    M (input) INTEGER, M>= 0
!>    The state space dimension (the row dimension of X, Y).
!>
[in]N
!>    N (input) INTEGER, 0 <= N <= M
!>    The number of data snapshot pairs
!>    (the number of columns of X and Y).
!>
[in,out]X
!>    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.
!>
[in]LDX
!>    LDX (input) INTEGER, LDX >= M
!>    The leading dimension of the array X.
!>
[in,out]Y
!>    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.
!>
[in]LDY
!>    LDY (input) INTEGER , LDY >= M
!>    The leading dimension of the array Y.
!>
[in]NRNK
!>    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.
!>
[in]TOL
!>    TOL (input) REAL(KIND=WP), 0 <= TOL < 1
!>    The tolerance for truncating small singular values.
!>    See the description of NRNK.
!>
[out]K
!>    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.
!>
[out]REIG
!>    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.
!>
[out]IMEIG
!>    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.
!>
[out]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.
!>
[in]LDZ
!>    LDZ (input) INTEGER , LDZ >= M
!>    The leading dimension of the array Z.
!>
[out]RES
!>    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.
!>
[out]B
!>    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.
!>
[in]LDB
!>    LDB (input) INTEGER, LDB >= M
!>    The leading dimension of the array B.
!>
[out]W
!>    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.
!>
[in]LDW
!>    LDW (input) INTEGER, LDW >= N
!>    The leading dimension of the array W.
!>
[out]S
!>    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.
!>
[in]LDS
!>    LDS (input) INTEGER, LDS >= N
!>    The leading dimension of the array S.
!>
[out]WORK
!>    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.
!>
[in]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.
!>
[out]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.
!>
[in]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.
!>
[out]INFO
!>    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.
!>
Author
Zlatko Drmac

Definition at line 530 of file dgedmd.f90.

535!
536! -- LAPACK driver routine --
537!
538! -- LAPACK is a software package provided by University of --
539! -- Tennessee, University of California Berkeley, University of --
540! -- Colorado Denver and NAG Ltd.. --
541!
542!.....
543 use, INTRINSIC :: iso_fortran_env, only: real64
544 IMPLICIT NONE
545 INTEGER, PARAMETER :: WP = real64
546!
547! Scalar arguments
548! ~~~~~~~~~~~~~~~~
549 CHARACTER, INTENT(IN) :: JOBS, JOBZ, JOBR, JOBF
550 INTEGER, INTENT(IN) :: WHTSVD, M, N, LDX, LDY, &
551 nrnk, ldz, ldb, ldw, lds, &
552 lwork, liwork
553 INTEGER, INTENT(OUT) :: K, INFO
554 REAL(KIND=wp), INTENT(IN) :: tol
555!
556! Array arguments
557! ~~~~~~~~~~~~~~~
558 REAL(KIND=wp), INTENT(INOUT) :: x(ldx,*), y(ldy,*)
559 REAL(KIND=wp), INTENT(OUT) :: z(ldz,*), b(ldb,*), &
560 w(ldw,*), s(lds,*)
561 REAL(KIND=wp), INTENT(OUT) :: reig(*), imeig(*), &
562 res(*)
563 REAL(KIND=wp), INTENT(OUT) :: work(*)
564 INTEGER, INTENT(OUT) :: IWORK(*)
565!
566! Parameters
567! ~~~~~~~~~~
568 REAL(KIND=wp), PARAMETER :: one = 1.0_wp
569 REAL(KIND=wp), PARAMETER :: zero = 0.0_wp
570!
571! Local scalars
572! ~~~~~~~~~~~~~
573 REAL(KIND=wp) :: ofl, rootsc, scale, small, &
574 ssum, xscl1, xscl2
575 INTEGER :: i, j, IMINWR, INFO1, INFO2, &
576 lwrkev, lwrsdd, lwrsvd, &
577 lwrsvq, mlwork, mwrkev, mwrsdd, &
578 mwrsvd, mwrsvj, mwrsvq, numrnk, &
579 olwork
580 LOGICAL :: BADXY, LQUERY, SCCOLX, SCCOLY, &
581 wntex, wntref, wntres, wntvec
582 CHARACTER :: JOBZL, T_OR_N
583 CHARACTER :: JSVOPT
584!
585! Local arrays
586! ~~~~~~~~~~~~
587 REAL(KIND=wp) :: ab(2,2), rdummy(2), rdummy2(2)
588!
589! External functions (BLAS and LAPACK)
590! ~~~~~~~~~~~~~~~~~
591 REAL(KIND=wp) dlange, dlamch, dnrm2
592 EXTERNAL dlange, dlamch, dnrm2, idamax
593 INTEGER IDAMAX
594 LOGICAL DISNAN, LSAME
595 EXTERNAL disnan, lsame
596!
597! External subroutines (BLAS and LAPACK)
598! ~~~~~~~~~~~~~~~~~~~~
599 EXTERNAL daxpy, dgemm, dscal
600 EXTERNAL dgeev, dgejsv, dgesdd, dgesvd, dgesvdq, &
601 dlacpy, dlascl, dlassq, xerbla
602!
603! Intrinsic functions
604! ~~~~~~~~~~~~~~~~~~~
605 INTRINSIC dble, int, max, sqrt
606!............................................................
607!
608! Test the input arguments
609!
610 wntres = lsame(jobr,'R')
611 sccolx = lsame(jobs,'S') .OR. lsame(jobs,'C')
612 sccoly = lsame(jobs,'Y')
613 wntvec = lsame(jobz,'V')
614 wntref = lsame(jobf,'R')
615 wntex = lsame(jobf,'E')
616 info = 0
617 lquery = ( ( lwork == -1 ) .OR. ( liwork == -1 ) )
618!
619 IF ( .NOT. (sccolx .OR. sccoly .OR. &
620 lsame(jobs,'N')) ) THEN
621 info = -1
622 ELSE IF ( .NOT. (wntvec .OR. lsame(jobz,'N') &
623 .OR. lsame(jobz,'F')) ) THEN
624 info = -2
625 ELSE IF ( .NOT. (wntres .OR. lsame(jobr,'N')) .OR. &
626 ( wntres .AND. (.NOT.wntvec) ) ) THEN
627 info = -3
628 ELSE IF ( .NOT. (wntref .OR. wntex .OR. &
629 lsame(jobf,'N') ) ) THEN
630 info = -4
631 ELSE IF ( .NOT.((whtsvd == 1) .OR. (whtsvd == 2) .OR. &
632 (whtsvd == 3) .OR. (whtsvd == 4) )) THEN
633 info = -5
634 ELSE IF ( m < 0 ) THEN
635 info = -6
636 ELSE IF ( ( n < 0 ) .OR. ( n > m ) ) THEN
637 info = -7
638 ELSE IF ( ldx < m ) THEN
639 info = -9
640 ELSE IF ( ldy < m ) THEN
641 info = -11
642 ELSE IF ( .NOT. (( nrnk == -2).OR.(nrnk == -1).OR. &
643 ((nrnk >= 1).AND.(nrnk <=n ))) ) THEN
644 info = -12
645 ELSE IF ( ( tol < zero ) .OR. ( tol >= one ) ) THEN
646 info = -13
647 ELSE IF ( ldz < m ) THEN
648 info = -18
649 ELSE IF ( (wntref .OR. wntex ) .AND. ( ldb < m ) ) THEN
650 info = -21
651 ELSE IF ( ldw < n ) THEN
652 info = -23
653 ELSE IF ( lds < n ) THEN
654 info = -25
655 END IF
656!
657 IF ( info == 0 ) THEN
658 ! Compute the minimal and the optimal workspace
659 ! requirements. Simulate running the code and
660 ! determine minimal and optimal sizes of the
661 ! workspace at any moment of the run.
662 IF ( n == 0 ) THEN
663 ! Quick return. All output except K is void.
664 ! INFO=1 signals the void input.
665 ! In case of a workspace query, the default
666 ! minimal workspace lengths are returned.
667 IF ( lquery ) THEN
668 iwork(1) = 1
669 work(1) = 2
670 work(2) = 2
671 ELSE
672 k = 0
673 END IF
674 info = 1
675 RETURN
676 END IF
677 mlwork = max(2,n)
678 olwork = max(2,n)
679 iminwr = 1
680 SELECT CASE ( whtsvd )
681 CASE (1)
682 ! The following is specified as the minimal
683 ! length of WORK in the definition of DGESVD:
684 ! MWRSVD = MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N))
685 mwrsvd = max(1,3*min(m,n)+max(m,n),5*min(m,n))
686 mlwork = max(mlwork,n + mwrsvd)
687 IF ( lquery ) THEN
688 CALL dgesvd( 'O', 'S', m, n, x, ldx, work, &
689 b, ldb, w, ldw, rdummy, -1, info1 )
690 lwrsvd = max( mwrsvd, int( rdummy(1) ) )
691 olwork = max(olwork,n + lwrsvd)
692 END IF
693 CASE (2)
694 ! The following is specified as the minimal
695 ! length of WORK in the definition of DGESDD:
696 ! MWRSDD = 3*MIN(M,N)*MIN(M,N) +
697 ! MAX( MAX(M,N),5*MIN(M,N)*MIN(M,N)+4*MIN(M,N) )
698 ! IMINWR = 8*MIN(M,N)
699 mwrsdd = 3*min(m,n)*min(m,n) + &
700 max( max(m,n),5*min(m,n)*min(m,n)+4*min(m,n) )
701 mlwork = max(mlwork,n + mwrsdd)
702 iminwr = 8*min(m,n)
703 IF ( lquery ) THEN
704 CALL dgesdd( 'O', m, n, x, ldx, work, b, &
705 ldb, w, ldw, rdummy, -1, iwork, info1 )
706 lwrsdd = max( mwrsdd, int( rdummy(1) ) )
707 olwork = max(olwork,n + lwrsdd)
708 END IF
709 CASE (3)
710 !LWQP3 = 3*N+1
711 !LWORQ = MAX(N, 1)
712 !MWRSVD = MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N))
713 !MWRSVQ = N + MAX( LWQP3, MWRSVD, LWORQ ) + MAX(M,2)
714 !MLWORK = N + MWRSVQ
715 !IMINWR = M+N-1
716 CALL dgesvdq( 'H', 'P', 'N', 'R', 'R', m, n, &
717 x, ldx, work, z, ldz, w, ldw, &
718 numrnk, iwork, liwork, rdummy, &
719 -1, rdummy2, -1, info1 )
720 iminwr = iwork(1)
721 mwrsvq = int(rdummy(2))
722 mlwork = max(mlwork,n+mwrsvq+int(rdummy2(1)))
723 IF ( lquery ) THEN
724 lwrsvq = max( mwrsvq, int(rdummy(1)) )
725 olwork = max(olwork,n+lwrsvq+int(rdummy2(1)))
726 END IF
727 CASE (4)
728 jsvopt = 'J'
729 !MWRSVJ = MAX( 7, 2*M+N, 6*N+2*N*N ) ! for JSVOPT='V'
730 mwrsvj = max( 7, 2*m+n, 4*n+n*n, 2*n+n*n+6 )
731 mlwork = max(mlwork,n+mwrsvj)
732 iminwr = max( 3, m+3*n )
733 IF ( lquery ) THEN
734 olwork = max(olwork,n+mwrsvj)
735 END IF
736 END SELECT
737 IF ( wntvec .OR. wntex .OR. lsame(jobz,'F') ) THEN
738 jobzl = 'V'
739 ELSE
740 jobzl = 'N'
741 END IF
742 ! Workspace calculation to the DGEEV call
743 IF ( lsame(jobzl,'V') ) THEN
744 mwrkev = max( 1, 4*n )
745 ELSE
746 mwrkev = max( 1, 3*n )
747 END IF
748 mlwork = max(mlwork,n+mwrkev)
749 IF ( lquery ) THEN
750 CALL dgeev( 'N', jobzl, n, s, lds, reig, &
751 imeig, w, ldw, w, ldw, rdummy, -1, info1 )
752 lwrkev = max( mwrkev, int(rdummy(1)) )
753 olwork = max( olwork, n+lwrkev )
754 END IF
755!
756 IF ( liwork < iminwr .AND. (.NOT.lquery) ) info = -29
757 IF ( lwork < mlwork .AND. (.NOT.lquery) ) info = -27
758 END IF
759!
760 IF( info /= 0 ) THEN
761 CALL xerbla( 'DGEDMD', -info )
762 RETURN
763 ELSE IF ( lquery ) THEN
764! Return minimal and optimal workspace sizes
765 iwork(1) = iminwr
766 work(1) = mlwork
767 work(2) = olwork
768 RETURN
769 END IF
770!............................................................
771!
772 ofl = dlamch('O')
773 small = dlamch('S')
774 badxy = .false.
775!
776! <1> Optional scaling of the snapshots (columns of X, Y)
777! ==========================================================
778 IF ( sccolx ) THEN
779 ! The columns of X will be normalized.
780 ! To prevent overflows, the column norms of X are
781 ! carefully computed using DLASSQ.
782 k = 0
783 DO i = 1, n
784 !WORK(i) = DNRM2( M, X(1,i), 1 )
785 ssum = one
786 scale = zero
787 CALL dlassq( m, x(1,i), 1, scale, ssum )
788 IF ( disnan(scale) .OR. disnan(ssum) ) THEN
789 k = 0
790 info = -8
791 CALL xerbla('DGEDMD',-info)
792 END IF
793 IF ( (scale /= zero) .AND. (ssum /= zero) ) THEN
794 rootsc = sqrt(ssum)
795 IF ( scale .GE. (ofl / rootsc) ) THEN
796! Norm of X(:,i) overflows. First, X(:,i)
797! is scaled by
798! ( ONE / ROOTSC ) / SCALE = 1/||X(:,i)||_2.
799! Next, the norm of X(:,i) is stored without
800! overflow as WORK(i) = - SCALE * (ROOTSC/M),
801! the minus sign indicating the 1/M factor.
802! Scaling is performed without overflow, and
803! underflow may occur in the smallest entries
804! of X(:,i). The relative backward and forward
805! errors are small in the ell_2 norm.
806 CALL dlascl( 'G', 0, 0, scale, one/rootsc, &
807 m, 1, x(1,i), m, info2 )
808 work(i) = - scale * ( rootsc / dble(m) )
809 ELSE
810! X(:,i) will be scaled to unit 2-norm
811 work(i) = scale * rootsc
812 CALL dlascl( 'G',0, 0, work(i), one, m, 1, &
813 x(1,i), m, info2 ) ! LAPACK CALL
814! X(1:M,i) = (ONE/WORK(i)) * X(1:M,i) ! INTRINSIC
815 END IF
816 ELSE
817 work(i) = zero
818 k = k + 1
819 END IF
820 END DO
821 IF ( k == n ) THEN
822 ! All columns of X are zero. Return error code -8.
823 ! (the 8th input variable had an illegal value)
824 k = 0
825 info = -8
826 CALL xerbla('DGEDMD',-info)
827 RETURN
828 END IF
829 DO i = 1, n
830! Now, apply the same scaling to the columns of Y.
831 IF ( work(i) > zero ) THEN
832 CALL dscal( m, one/work(i), y(1,i), 1 ) ! BLAS CALL
833! Y(1:M,i) = (ONE/WORK(i)) * Y(1:M,i) ! INTRINSIC
834 ELSE IF ( work(i) < zero ) THEN
835 CALL dlascl( 'G', 0, 0, -work(i), &
836 one/dble(m), m, 1, y(1,i), m, info2 ) ! LAPACK CALL
837 ELSE IF ( y(idamax(m, y(1,i),1),i ) &
838 /= zero ) THEN
839! X(:,i) is zero vector. For consistency,
840! Y(:,i) should also be zero. If Y(:,i) is not
841! zero, then the data might be inconsistent or
842! corrupted. If JOBS == 'C', Y(:,i) is set to
843! zero and a warning flag is raised.
844! The computation continues but the
845! situation will be reported in the output.
846 badxy = .true.
847 IF ( lsame(jobs,'C')) &
848 CALL dscal( m, zero, y(1,i), 1 ) ! BLAS CALL
849 END IF
850 END DO
851 END IF
852 !
853 IF ( sccoly ) THEN
854 ! The columns of Y will be normalized.
855 ! To prevent overflows, the column norms of Y are
856 ! carefully computed using DLASSQ.
857 DO i = 1, n
858 !WORK(i) = DNRM2( M, Y(1,i), 1 )
859 ssum = one
860 scale = zero
861 CALL dlassq( m, y(1,i), 1, scale, ssum )
862 IF ( disnan(scale) .OR. disnan(ssum) ) THEN
863 k = 0
864 info = -10
865 CALL xerbla('DGEDMD',-info)
866 END IF
867 IF ( scale /= zero .AND. (ssum /= zero) ) THEN
868 rootsc = sqrt(ssum)
869 IF ( scale .GE. (ofl / rootsc) ) THEN
870! Norm of Y(:,i) overflows. First, Y(:,i)
871! is scaled by
872! ( ONE / ROOTSC ) / SCALE = 1/||Y(:,i)||_2.
873! Next, the norm of Y(:,i) is stored without
874! overflow as WORK(i) = - SCALE * (ROOTSC/M),
875! the minus sign indicating the 1/M factor.
876! Scaling is performed without overflow, and
877! underflow may occur in the smallest entries
878! of Y(:,i). The relative backward and forward
879! errors are small in the ell_2 norm.
880 CALL dlascl( 'G', 0, 0, scale, one/rootsc, &
881 m, 1, y(1,i), m, info2 )
882 work(i) = - scale * ( rootsc / dble(m) )
883 ELSE
884! X(:,i) will be scaled to unit 2-norm
885 work(i) = scale * rootsc
886 CALL dlascl( 'G',0, 0, work(i), one, m, 1, &
887 y(1,i), m, info2 ) ! LAPACK CALL
888! Y(1:M,i) = (ONE/WORK(i)) * Y(1:M,i) ! INTRINSIC
889 END IF
890 ELSE
891 work(i) = zero
892 END IF
893 END DO
894 DO i = 1, n
895! Now, apply the same scaling to the columns of X.
896 IF ( work(i) > zero ) THEN
897 CALL dscal( m, one/work(i), x(1,i), 1 ) ! BLAS CALL
898! X(1:M,i) = (ONE/WORK(i)) * X(1:M,i) ! INTRINSIC
899 ELSE IF ( work(i) < zero ) THEN
900 CALL dlascl( 'G', 0, 0, -work(i), &
901 one/dble(m), m, 1, x(1,i), m, info2 ) ! LAPACK CALL
902 ELSE IF ( x(idamax(m, x(1,i),1),i ) &
903 /= zero ) THEN
904! Y(:,i) is zero vector. If X(:,i) is not
905! zero, then a warning flag is raised.
906! The computation continues but the
907! situation will be reported in the output.
908 badxy = .true.
909 END IF
910 END DO
911 END IF
912!
913! <2> SVD of the data snapshot matrix X.
914! =====================================
915! The left singular vectors are stored in the array X.
916! The right singular vectors are in the array W.
917! The array W will later on contain the eigenvectors
918! of a Rayleigh quotient.
919 numrnk = n
920 SELECT CASE ( whtsvd )
921 CASE (1)
922 CALL dgesvd( 'O', 'S', m, n, x, ldx, work, b, &
923 ldb, w, ldw, work(n+1), lwork-n, info1 ) ! LAPACK CALL
924 t_or_n = 'T'
925 CASE (2)
926 CALL dgesdd( 'O', m, n, x, ldx, work, b, ldb, w, &
927 ldw, work(n+1), lwork-n, iwork, info1 ) ! LAPACK CALL
928 t_or_n = 'T'
929 CASE (3)
930 CALL dgesvdq( 'H', 'P', 'N', 'R', 'R', m, n, &
931 x, ldx, work, z, ldz, w, ldw, &
932 numrnk, iwork, liwork, work(n+max(2,m)+1),&
933 lwork-n-max(2,m), work(n+1), max(2,m), info1) ! LAPACK CALL
934 CALL dlacpy( 'A', m, numrnk, z, ldz, x, ldx ) ! LAPACK CALL
935 t_or_n = 'T'
936 CASE (4)
937 CALL dgejsv( 'F', 'U', jsvopt, 'N', 'N', 'P', m, &
938 n, x, ldx, work, z, ldz, w, ldw, &
939 work(n+1), lwork-n, iwork, info1 ) ! LAPACK CALL
940 CALL dlacpy( 'A', m, n, z, ldz, x, ldx ) ! LAPACK CALL
941 t_or_n = 'N'
942 xscl1 = work(n+1)
943 xscl2 = work(n+2)
944 IF ( xscl1 /= xscl2 ) THEN
945 ! This is an exceptional situation. If the
946 ! data matrices are not scaled and the
947 ! largest singular value of X overflows.
948 ! In that case DGEJSV can return the SVD
949 ! in scaled form. The scaling factor can be used
950 ! to rescale the data (X and Y).
951 CALL dlascl( 'G', 0, 0, xscl1, xscl2, m, n, y, ldy, info2 )
952 END IF
953 END SELECT
954!
955 IF ( info1 > 0 ) THEN
956 ! The SVD selected subroutine did not converge.
957 ! Return with an error code.
958 info = 2
959 RETURN
960 END IF
961!
962 IF ( work(1) == zero ) THEN
963 ! The largest computed singular value of (scaled)
964 ! X is zero. Return error code -8
965 ! (the 8th input variable had an illegal value).
966 k = 0
967 info = -8
968 CALL xerbla('DGEDMD',-info)
969 RETURN
970 END IF
971!
972
973 ! snapshots matrix X. This depends on the
974 ! parameters NRNK and TOL.
975
976 SELECT CASE ( nrnk )
977 CASE ( -1 )
978 k = 1
979 DO i = 2, numrnk
980 IF ( ( work(i) <= work(1)*tol ) .OR. &
981 ( work(i) <= small ) ) EXIT
982 k = k + 1
983 END DO
984 CASE ( -2 )
985 k = 1
986 DO i = 1, numrnk-1
987 IF ( ( work(i+1) <= work(i)*tol ) .OR. &
988 ( work(i) <= small ) ) EXIT
989 k = k + 1
990 END DO
991 CASE DEFAULT
992 k = 1
993 DO i = 2, nrnk
994 IF ( work(i) <= small ) EXIT
995 k = k + 1
996 END DO
997 END SELECT
998 ! Now, U = X(1:M,1:K) is the SVD/POD basis for the
999 ! snapshot data in the input matrix X.
1000
1001
1002 ! Depending on the requested outputs, the computation
1003 ! is organized to compute additional auxiliary
1004 ! matrices (for the residuals and refinements).
1005 !
1006 ! In all formulas below, we need V_k*Sigma_k^(-1)
1007 ! where either V_k is in W(1:N,1:K), or V_k^T is in
1008 ! W(1:K,1:N). Here Sigma_k=diag(WORK(1:K)).
1009 IF ( lsame(t_or_n, 'N') ) THEN
1010 DO i = 1, k
1011 CALL dscal( n, one/work(i), w(1,i), 1 ) ! BLAS CALL
1012 ! W(1:N,i) = (ONE/WORK(i)) * W(1:N,i) ! INTRINSIC
1013 END DO
1014 ELSE
1015 ! This non-unit stride access is due to the fact
1016 ! that DGESVD, DGESVDQ and DGESDD return the
1017 ! transposed matrix of the right singular vectors.
1018 !DO i = 1, K
1019 ! CALL DSCAL( N, ONE/WORK(i), W(i,1), LDW ) ! BLAS CALL
1020 ! ! W(i,1:N) = (ONE/WORK(i)) * W(i,1:N) ! INTRINSIC
1021 !END DO
1022 DO i = 1, k
1023 work(n+i) = one/work(i)
1024 END DO
1025 DO j = 1, n
1026 DO i = 1, k
1027 w(i,j) = (work(n+i))*w(i,j)
1028 END DO
1029 END DO
1030 END IF
1031!
1032 IF ( wntref ) THEN
1033 !
1034 ! Need A*U(:,1:K)=Y*V_k*inv(diag(WORK(1:K)))
1035 ! for computing the refined Ritz vectors
1036 ! (optionally, outside DGEDMD).
1037 CALL dgemm( 'N', t_or_n, m, k, n, one, y, ldy, w, &
1038 ldw, zero, z, ldz ) ! BLAS CALL
1039 ! Z(1:M,1:K)=MATMUL(Y(1:M,1:N),TRANSPOSE(W(1:K,1:N))) ! INTRINSIC, for T_OR_N=='T'
1040 ! Z(1:M,1:K)=MATMUL(Y(1:M,1:N),W(1:N,1:K)) ! INTRINSIC, for T_OR_N=='N'
1041 !
1042 ! At this point Z contains
1043 ! A * U(:,1:K) = Y * V_k * Sigma_k^(-1), and
1044 ! this is needed for computing the residuals.
1045 ! This matrix is returned in the array B and
1046 ! it can be used to compute refined Ritz vectors.
1047 CALL dlacpy( 'A', m, k, z, ldz, b, ldb ) ! BLAS CALL
1048 ! B(1:M,1:K) = Z(1:M,1:K) ! INTRINSIC
1049
1050 CALL dgemm( 'T', 'N', k, k, m, one, x, ldx, z, &
1051 ldz, zero, s, lds ) ! BLAS CALL
1052 ! S(1:K,1:K) = MATMUL(TANSPOSE(X(1:M,1:K)),Z(1:M,1:K)) ! INTRINSIC
1053 ! At this point S = U^T * A * U is the Rayleigh quotient.
1054 ELSE
1055 ! A * U(:,1:K) is not explicitly needed and the
1056 ! computation is organized differently. The Rayleigh
1057 ! quotient is computed more efficiently.
1058 CALL dgemm( 'T', 'N', k, n, m, one, x, ldx, y, ldy, &
1059 zero, z, ldz ) ! BLAS CALL
1060 ! Z(1:K,1:N) = MATMUL( TRANSPOSE(X(1:M,1:K)), Y(1:M,1:N) ) ! INTRINSIC
1061 ! In the two DGEMM calls here, can use K for LDZ.
1062 CALL dgemm( 'N', t_or_n, k, k, n, one, z, ldz, w, &
1063 ldw, zero, s, lds ) ! BLAS CALL
1064 ! S(1:K,1:K) = MATMUL(Z(1:K,1:N),TRANSPOSE(W(1:K,1:N))) ! INTRINSIC, for T_OR_N=='T'
1065 ! S(1:K,1:K) = MATMUL(Z(1:K,1:N),(W(1:N,1:K))) ! INTRINSIC, for T_OR_N=='N'
1066 ! At this point S = U^T * A * U is the Rayleigh quotient.
1067 ! If the residuals are requested, save scaled V_k into Z.
1068 ! Recall that V_k or V_k^T is stored in W.
1069 IF ( wntres .OR. wntex ) THEN
1070 IF ( lsame(t_or_n, 'N') ) THEN
1071 CALL dlacpy( 'A', n, k, w, ldw, z, ldz )
1072 ELSE
1073 CALL dlacpy( 'A', k, n, w, ldw, z, ldz )
1074 END IF
1075 END IF
1076 END IF
1077!
1078
1079 ! right eigenvectors of the Rayleigh quotient.
1080 !
1081 CALL dgeev( 'N', jobzl, k, s, lds, reig, imeig, w, &
1082 ldw, w, ldw, work(n+1), lwork-n, info1 ) ! LAPACK CALL
1083 !
1084 ! W(1:K,1:K) contains the eigenvectors of the Rayleigh
1085 ! quotient. Even in the case of complex spectrum, all
1086 ! computation is done in real arithmetic. REIG and
1087 ! IMEIG are the real and the imaginary parts of the
1088 ! eigenvalues, so that the spectrum is given as
1089 ! REIG(:) + sqrt(-1)*IMEIG(:). Complex conjugate pairs
1090 ! are listed at consecutive positions. For such a
1091 ! complex conjugate pair of the eigenvalues, the
1092 ! corresponding eigenvectors are also a complex
1093 ! conjugate pair with the real and imaginary parts
1094 ! stored column-wise in W at the corresponding
1095 ! consecutive column indices. See the description of Z.
1096 ! Also, see the description of DGEEV.
1097 IF ( info1 > 0 ) THEN
1098 ! DGEEV failed to compute the eigenvalues and
1099 ! eigenvectors of the Rayleigh quotient.
1100 info = 3
1101 RETURN
1102 END IF
1103!
1104 ! <6> Compute the eigenvectors (if requested) and,
1105 ! the residuals (if requested).
1106 !
1107 IF ( wntvec .OR. wntex ) THEN
1108 IF ( wntres ) THEN
1109 IF ( wntref ) THEN
1110 ! Here, if the refinement is requested, we have
1111 ! A*U(:,1:K) already computed and stored in Z.
1112 ! For the residuals, need Y = A * U(:,1;K) * W.
1113 CALL dgemm( 'N', 'N', m, k, k, one, z, ldz, w, &
1114 ldw, zero, y, ldy ) ! BLAS CALL
1115 ! Y(1:M,1:K) = Z(1:M,1:K) * W(1:K,1:K) ! INTRINSIC
1116 ! This frees Z; Y contains A * U(:,1:K) * W.
1117 ELSE
1118 ! Compute S = V_k * Sigma_k^(-1) * W, where
1119 ! V_k * Sigma_k^(-1) is stored in Z
1120 CALL dgemm( t_or_n, 'N', n, k, k, one, z, ldz, &
1121 w, ldw, zero, s, lds)
1122 ! Then, compute Z = Y * S =
1123 ! = Y * V_k * Sigma_k^(-1) * W(1:K,1:K) =
1124 ! = A * U(:,1:K) * W(1:K,1:K)
1125 CALL dgemm( 'N', 'N', m, k, n, one, y, ldy, s, &
1126 lds, zero, z, ldz)
1127 ! Save a copy of Z into Y and free Z for holding
1128 ! the Ritz vectors.
1129 CALL dlacpy( 'A', m, k, z, ldz, y, ldy )
1130 IF ( wntex ) CALL dlacpy( 'A', m, k, z, ldz, b, ldb )
1131 END IF
1132 ELSE IF ( wntex ) THEN
1133 ! Compute S = V_k * Sigma_k^(-1) * W, where
1134 ! V_k * Sigma_k^(-1) is stored in Z
1135 CALL dgemm( t_or_n, 'N', n, k, k, one, z, ldz, &
1136 w, ldw, zero, s, lds )
1137 ! Then, compute Z = Y * S =
1138 ! = Y * V_k * Sigma_k^(-1) * W(1:K,1:K) =
1139 ! = A * U(:,1:K) * W(1:K,1:K)
1140 CALL dgemm( 'N', 'N', m, k, n, one, y, ldy, s, &
1141 lds, zero, b, ldb )
1142 ! The above call replaces the following two calls
1143 ! that were used in the developing-testing phase.
1144 ! CALL DGEMM( 'N', 'N', M, K, N, ONE, Y, LDY, S, &
1145 ! LDS, ZERO, Z, LDZ)
1146 ! Save a copy of Z into B and free Z for holding
1147 ! the Ritz vectors.
1148 ! CALL DLACPY( 'A', M, K, Z, LDZ, B, LDB )
1149 END IF
1150!
1151 ! Compute the real form of the Ritz vectors
1152 IF ( wntvec ) CALL dgemm( 'N', 'N', m, k, k, one, x, ldx, w, ldw, &
1153 zero, z, ldz ) ! BLAS CALL
1154 ! Z(1:M,1:K) = MATMUL(X(1:M,1:K), W(1:K,1:K)) ! INTRINSIC
1155!
1156 IF ( wntres ) THEN
1157 i = 1
1158 DO WHILE ( i <= k )
1159 IF ( imeig(i) == zero ) THEN
1160 ! have a real eigenvalue with real eigenvector
1161 CALL daxpy( m, -reig(i), z(1,i), 1, y(1,i), 1 ) ! BLAS CALL
1162 ! Y(1:M,i) = Y(1:M,i) - REIG(i) * Z(1:M,i) ! INTRINSIC
1163 res(i) = dnrm2( m, y(1,i), 1) ! BLAS CALL
1164 i = i + 1
1165 ELSE
1166 ! Have a complex conjugate pair
1167 ! REIG(i) +- sqrt(-1)*IMEIG(i).
1168 ! Since all computation is done in real
1169 ! arithmetic, the formula for the residual
1170 ! is recast for real representation of the
1171 ! complex conjugate eigenpair. See the
1172 ! description of RES.
1173 ab(1,1) = reig(i)
1174 ab(2,1) = -imeig(i)
1175 ab(1,2) = imeig(i)
1176 ab(2,2) = reig(i)
1177 CALL dgemm( 'N', 'N', m, 2, 2, -one, z(1,i), &
1178 ldz, ab, 2, one, y(1,i), ldy ) ! BLAS CALL
1179 ! Y(1:M,i:i+1) = Y(1:M,i:i+1) - Z(1:M,i:i+1) * AB ! INTRINSIC
1180 res(i) = dlange( 'F', m, 2, y(1,i), ldy, &
1181 work(n+1) ) ! LAPACK CALL
1182 res(i+1) = res(i)
1183 i = i + 2
1184 END IF
1185 END DO
1186 END IF
1187 END IF
1188!
1189 IF ( whtsvd == 4 ) THEN
1190 work(n+1) = xscl1
1191 work(n+2) = xscl2
1192 END IF
1193!
1194! Successful exit.
1195 IF ( .NOT. badxy ) THEN
1196 info = 0
1197 ELSE
1198 ! A warning on possible data inconsistency.
1199 ! This should be a rare event.
1200 info = 4
1201 END IF
1202!............................................................
1203 RETURN
1204! ......
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dgeev
subroutine dgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)
DGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition dgeev.f:191
dgejsv
subroutine dgejsv(joba, jobu, jobv, jobr, jobt, jobp, m, n, a, lda, sva, u, ldu, v, ldv, work, lwork, iwork, info)
DGEJSV
Definition dgejsv.f:474
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
dgesdd
subroutine dgesdd(jobz, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, iwork, info)
DGESDD
Definition dgesdd.f:211
dgesvd
subroutine dgesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, info)
DGESVD computes the singular value decomposition (SVD) for GE matrices
Definition dgesvd.f:209
dgesvdq
subroutine dgesvdq(joba, jobp, jobr, jobu, jobv, m, n, a, lda, s, u, ldu, v, ldv, numrank, iwork, liwork, work, lwork, rwork, lrwork, info)
DGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE...
Definition dgesvdq.f:413
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
disnan
logical function disnan(din)
DISNAN tests input for NaN.
Definition disnan.f:57
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:101
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlange
double precision function dlange(norm, m, n, a, lda, work)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition dlange.f:112
dlascl
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:142
dlassq
subroutine dlassq(n, x, incx, scale, sumsq)
DLASSQ updates a sum of squares represented in scaled form.
Definition dlassq.f90:122
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dnrm2
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition dnrm2.f90:89
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
Here is the call graph for this function:
Here is the caller graph for this function: