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

◆ zgedmd()

subroutine zgedmd ( 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,
complex(kind=wp), dimension(ldx,*), intent(inout) x,
integer, intent(in) ldx,
complex(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,
complex(kind=wp), dimension(*), intent(out) eigs,
complex(kind=wp), dimension(ldz,*), intent(out) z,
integer, intent(in) ldz,
real(kind=wp), dimension(*), intent(out) res,
complex(kind=wp), dimension(ldb,*), intent(out) b,
integer, intent(in) ldb,
complex(kind=wp), dimension(ldw,*), intent(out) w,
integer, intent(in) ldw,
complex(kind=wp), dimension(lds,*), intent(out) s,
integer, intent(in) lds,
complex(kind=wp), dimension(*), intent(out) zwork,
integer, intent(in) lzwork,
real(kind=wp), dimension(*), intent(out) rwork,
integer, intent(in) lrwork,
integer, dimension(*), intent(out) iwork,
integer, intent(in) liwork,
integer, intent(out) info )

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

Purpose:
!>    ZGEDMD computes the Dynamic Mode Decomposition (DMD) for
!>    a pair of data snapshot matrices. For the input matrices
!>    X and Y such that Y = A*X with an unaccessible matrix
!>    A, ZGEDMD computes a certain number of Ritz pairs of A using
!>    the standard Rayleigh-Ritz extraction from a subspace of
!>    range(X) that is determined using the leading left singular
!>    vectors of X. Optionally, ZGEDMD returns the residuals
!>    of the computed Ritz pairs, the information needed for
!>    a refinement of the Ritz vectors, or the eigenvectors of
!>    the Exact DMD.
!>    For further details see the references listed
!>    below. For more details of the implementation see [3].
!>
References:
!>    [1] P. Schmid: Dynamic mode decomposition of numerical
!>        and experimental data,
!>        Journal of Fluid Mechanics 656, 5-28, 2010.
!>    [2] Z. Drmac, I. Mezic, R. Mohr: Data driven modal
!>        decompositions: analysis and enhancements,
!>        SIAM J. on Sci. Comp. 40 (4), A2253-A2285, 2018.
!>    [3] Z. Drmac: A LAPACK implementation of the Dynamic
!>        Mode Decomposition I. Technical report. AIMDyn Inc.
!>        and LAPACK Working Note 298.
!>    [4] J. Tu, C. W. Rowley, D. M. Luchtenburg, S. L.
!>        Brunton, N. Kutz: On Dynamic Mode Decomposition:
!>        Theory and Applications, Journal of Computational
!>        Dynamics 1(2), 391 -421, 2014.
!>
Developed and supported by:
!>    Developed and coded by Zlatko Drmac, Faculty of Science,
!>    University of Zagreb;  drmac@math.hr
!>    In cooperation with
!>    AIMdyn Inc., Santa Barbara, CA.
!>    and supported by
!>    - DARPA SBIR project  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) 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 :: ZGESVD (the QR SVD algorithm)
!>    2 :: ZGESDD (the Divide and Conquer algorithm; if enough
!>         workspace available, this is the fastest option)
!>    3 :: ZGESVDQ (the preconditioned QR SVD  ; this and 4
!>         are the most accurate options)
!>    4 :: ZGEJSV (the preconditioned Jacobi SVD; this and 3
!>         are the most accurate options)
!>    For the four methods above, a significant difference in
!>    the accuracy of small singular values is possible if
!>    the snapshots vary in norm so that X is severely
!>    ill-conditioned. If small (smaller than EPS*||X||)
!>    singular values are of interest and JOBS=='N',  then
!>    the options (3, 4) give the most accurate results, where
!>    the option 4 is slightly better and with stronger
!>    theoretical background.
!>    If JOBS=='S', i.e. the columns of X will be normalized,
!>    then all methods give nearly equally accurate results.
!>
[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) COMPLEX(KIND=WP) M-by-N array
!>    > On entry, X contains the data snapshot matrix X. It is
!>    assumed that the column norms of X are in the range of
!>    the normalized floating point numbers.
!>    < On exit, the leading K columns of X contain a POD basis,
!>    i.e. the leading K left singular vectors of the input
!>    data matrix X, U(:,1:K). All N columns of X contain all
!>    left singular vectors of the input matrix X.
!>    See the descriptions of K, Z and W.
!>
[in]LDX
!>    LDX (input) INTEGER, LDX >= M
!>    The leading dimension of the array X.
!>
[in,out]Y
!>    Y (input/workspace/output) COMPLEX(KIND=WP) M-by-N array
!>    > On entry, Y contains the data snapshot matrix Y
!>    < On exit,
!>    If JOBR == 'R', the leading K columns of Y  contain
!>    the residual vectors for the computed Ritz pairs.
!>    See the description of RES.
!>    If JOBR == 'N', Y contains the original input data,
!>                    scaled according to the value of JOBS.
!>
[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]EIGS
!>    EIGS (output) COMPLEX(KIND=WP) N-by-1 array
!>    The leading K (K<=N) entries of EIGS contain
!>    the computed eigenvalues (Ritz values).
!>    See the descriptions of K, and Z.
!>
[out]Z
!>    Z (workspace/output) COMPLEX(KIND=WP)  M-by-N array
!>    If JOBZ =='V' then Z contains the  Ritz vectors.  Z(:,i)
!>    is an eigenvector of the i-th Ritz value; ||Z(:,i)||_2=1.
!>    If JOBZ == 'F', then the Z(:,i)'s are given implicitly as
!>    the columns of X(:,1:K)*W(1:K,1:K), i.e. X(:,1:K)*W(:,i)
!>    is an eigenvector corresponding to EIGS(i). The columns
!>    of W(1:k,1:K) are the computed eigenvectors of the
!>    K-by-K Rayleigh quotient.
!>    See the descriptions of EIGS, X and W.
!>
[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,
!>    RES(i) = || A * Z(:,i) - EIGS(i)*Z(:,i))||_2.
!>    See the description of EIGS and Z.
!>
[out]B
!>    B (output) COMPLEX(KIND=WP)  M-by-N array.
!>    IF JOBF =='R', B(1:M,1:K) contains A*U(:,1:K), and can
!>    be used for computing the refined vectors; see further
!>    details in the provided references.
!>    If JOBF == 'E', B(1:M,1:K) contains
!>    A*U(:,1:K)*W(1:K,1:K), which are the vectors from the
!>    Exact DMD, up to scaling by the inverse eigenvalues.
!>    If JOBF =='N', then B is not referenced.
!>    See the descriptions of X, W, K.
!>
[in]LDB
!>    LDB (input) INTEGER, LDB >= M
!>    The leading dimension of the array B.
!>
[out]W
!>    W (workspace/output) COMPLEX(KIND=WP) N-by-N array
!>    On exit, W(1:K,1:K) contains the K computed
!>    eigenvectors of the matrix Rayleigh quotient.
!>    The Ritz vectors (returned in Z) are the
!>    product of X (containing a POD basis for the input
!>    matrix X) and W. See the descriptions of K, S, X and Z.
!>    W is also used as a workspace to temporarily store the
!>    right singular vectors of X.
!>
[in]LDW
!>    LDW (input) INTEGER, LDW >= N
!>    The leading dimension of the array W.
!>
[out]S
!>    S (workspace/output) COMPLEX(KIND=WP) N-by-N array
!>    The array S(1:K,1:K) is used for the matrix Rayleigh
!>    quotient. This content is overwritten during
!>    the eigenvalue decomposition by ZGEEV.
!>    See the description of K.
!>
[in]LDS
!>    LDS (input) INTEGER, LDS >= N
!>    The leading dimension of the array S.
!>
[out]ZWORK
!>    ZWORK (workspace/output) COMPLEX(KIND=WP) LZWORK-by-1 array
!>    ZWORK is used as complex workspace in the complex SVD, as
!>    specified by WHTSVD (1,2, 3 or 4) and for ZGEEV for computing
!>    the eigenvalues of a Rayleigh quotient.
!>    If the call to ZGEDMD is only workspace query, then
!>    ZWORK(1) contains the minimal complex workspace length and
!>    ZWORK(2) is the optimal complex workspace length.
!>    Hence, the length of work is at least 2.
!>    See the description of LZWORK.
!>
[in]LZWORK
!>    LZWORK (input) INTEGER
!>    The minimal length of the workspace vector ZWORK.
!>    LZWORK is calculated as MAX(LZWORK_SVD, LZWORK_ZGEEV),
!>    where LZWORK_ZGEEV = MAX( 1, 2*N )  and the minimal
!>    LZWORK_SVD is calculated as follows
!>    If WHTSVD == 1 :: ZGESVD ::
!>       LZWORK_SVD = MAX(1,2*MIN(M,N)+MAX(M,N))
!>    If WHTSVD == 2 :: ZGESDD ::
!>       LZWORK_SVD = 2*MIN(M,N)*MIN(M,N)+2*MIN(M,N)+MAX(M,N)
!>    If WHTSVD == 3 :: ZGESVDQ ::
!>       LZWORK_SVD = obtainable by a query
!>    If WHTSVD == 4 :: ZGEJSV ::
!>       LZWORK_SVD = obtainable by a query
!>    If on entry LZWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    and the optimal workspace lengths and returns them in
!>    LZWORK(1) and LZWORK(2), respectively.
!>
[out]RWORK
!>    RWORK (workspace/output) REAL(KIND=WP) LRWORK-by-1 array
!>    On exit, RWORK(1:N) contains the singular values of
!>    X (for JOBS=='N') or column scaled X (JOBS=='S', 'C').
!>    If WHTSVD==4, then RWORK(N+1) and RWORK(N+2) contain
!>    scaling factor RWORK(N+2)/RWORK(N+1) used to scale X
!>    and Y to avoid overflow in the SVD of X.
!>    This may be of interest if the scaling option is off
!>    and as many as possible smallest eigenvalues are
!>    desired to the highest feasible accuracy.
!>    If the call to ZGEDMD is only workspace query, then
!>    RWORK(1) contains the minimal workspace length.
!>    See the description of LRWORK.
!>
[in]LRWORK
!>    LRWORK (input) INTEGER
!>    The minimal length of the workspace vector RWORK.
!>    LRWORK is calculated as follows:
!>    LRWORK = MAX(1, N+LRWORK_SVD,N+LRWORK_ZGEEV), where
!>    LRWORK_ZGEEV = MAX(1,2*N) and RWORK_SVD is the real workspace
!>    for the SVD subroutine determined by the input parameter
!>    WHTSVD.
!>    If WHTSVD == 1 :: ZGESVD ::
!>       LRWORK_SVD = 5*MIN(M,N)
!>    If WHTSVD == 2 :: ZGESDD ::
!>       LRWORK_SVD =  MAX(5*MIN(M,N)*MIN(M,N)+7*MIN(M,N),
!>       2*MAX(M,N)*MIN(M,N)+2*MIN(M,N)*MIN(M,N)+MIN(M,N) ) )
!>    If WHTSVD == 3 :: ZGESVDQ ::
!>       LRWORK_SVD = obtainable by a query
!>    If WHTSVD == 4 :: ZGEJSV ::
!>       LRWORK_SVD = obtainable by a query
!>    If on entry LRWORK = -1, then a workspace query is
!>    assumed and the procedure only computes the minimal
!>    real workspace length and returns it in RWORK(1).
!>
[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  ZWORK, RWORK and
!>    IWORK. See the descriptions of ZWORK, RWORK 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 496 of file zgedmd.f90.

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