zgedmd.f90

Go to the documentation of this file.

 
!
!  =========== DOCUMENTATION ===========
!
!  Definition:
!  ===========
!
!      SUBROUTINE ZGEDMD( JOBS, JOBZ, JOBR, JOBF,  WHTSVD,   &
!                         M, N, X, LDX, Y, LDY, NRNK, TOL,   &
!                         K, EIGS, Z, LDZ, RES, B,    LDB,   &
!                         W, LDW,  S, LDS, ZWORK,  LZWORK,   &
!                         RWORK, LRWORK, IWORK, LIWORK, INFO )
!......
!      USE, INTRINSIC :: iso_fortran_env, only: real64
!      IMPLICIT NONE
!      INTEGER, PARAMETER :: WP = real64
!
!......
!      Scalar arguments
!      CHARACTER, INTENT(IN)   :: JOBS,   JOBZ,  JOBR,  JOBF
!      INTEGER,   INTENT(IN)   :: WHTSVD, M, N,   LDX,  LDY, &
!                                 NRNK, LDZ, LDB, LDW,  LDS, &
!                                 LIWORK, LRWORK, LZWORK
!      INTEGER,       INTENT(OUT)  :: K, INFO
!      REAL(KIND=WP), INTENT(IN)   ::    TOL
!      Array arguments
!      COMPLEX(KIND=WP), INTENT(INOUT) :: X(LDX,*), Y(LDY,*)
!      COMPLEX(KIND=WP), INTENT(OUT)   :: Z(LDZ,*), B(LDB,*), &
!                                         W(LDW,*), S(LDS,*)
!      COMPLEX(KIND=WP), INTENT(OUT)   :: EIGS(*)
!      COMPLEX(KIND=WP), INTENT(OUT)   :: ZWORK(*)
!      REAL(KIND=WP),    INTENT(OUT)   :: RES(*)
!      REAL(KIND=WP),    INTENT(OUT)   :: RWORK(*)
!      INTEGER,          INTENT(OUT)   :: IWORK(*)
!
!............................................................
!     =============
!............................................................
!     ================
!......................................................................
!     ================================
!......................................................................
!     ==============================
!............................................................
!     Arguments
!     =========
!
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!
!  Authors:
!  ========
!
!
!
!.............................................................
!.............................................................
      SUBROUTINE zgedmd( JOBS, JOBZ, JOBR, JOBF,  WHTSVD,   &
                         M, N, X, LDX, Y, LDY, NRNK, TOL,   &
                         K, EIGS, Z, LDZ, RES, B,    LDB,   &
                         W, LDW,  S, LDS, ZWORK,  LZWORK,   &
                         RWORK, LRWORK, IWORK, LIWORK, INFO )
!
!  -- LAPACK driver routine                                           --
!
!  -- LAPACK is a software package provided by University of          --
!  -- Tennessee, University of California Berkeley, University of     --
!  -- Colorado Denver and NAG Ltd..                                   --
!
!.....
      use, INTRINSIC :: iso_fortran_env, only: real64
      IMPLICIT NONE
      INTEGER, PARAMETER :: WP = real64
!
!     Scalar arguments
!     ~~~~~~~~~~~~~~~~
      CHARACTER, INTENT(IN)   :: JOBS,   JOBZ,  JOBR,  JOBF
      INTEGER,   INTENT(IN)   :: WHTSVD, M, N,   LDX,  LDY, &
                                 NRNK, LDZ, LDB, LDW,  LDS, &
                                 LIWORK, LRWORK, LZWORK
      INTEGER,       INTENT(OUT)  :: K, INFO
      REAL(KIND=wp), INTENT(IN)   ::    tol
!
!     Array arguments
!     ~~~~~~~~~~~~~~~
      COMPLEX(KIND=WP), INTENT(INOUT) :: X(LDX,*), Y(LDY,*)
      COMPLEX(KIND=WP), INTENT(OUT)   :: Z(LDZ,*), B(LDB,*), &
                                         W(LDW,*), S(LDS,*)
      COMPLEX(KIND=WP), INTENT(OUT)   :: EIGS(*)
      COMPLEX(KIND=WP), INTENT(OUT)   :: ZWORK(*)
      REAL(KIND=wp),    INTENT(OUT)   :: res(*)
      REAL(KIND=wp),    INTENT(OUT)   :: rwork(*)
      INTEGER,          INTENT(OUT)   :: IWORK(*)
!
!     Parameters
!     ~~~~~~~~~~
      REAL(KIND=wp),    PARAMETER ::  one = 1.0_wp
      REAL(KIND=wp),    PARAMETER :: zero = 0.0_wp
      COMPLEX(KIND=WP), PARAMETER ::  ZONE = ( 1.0_wp, 0.0_wp )
      COMPLEX(KIND=WP), PARAMETER :: ZZERO = ( 0.0_wp, 0.0_wp )
!
!     Local scalars
!     ~~~~~~~~~~~~~
      REAL(KIND=wp) :: ofl,   rootsc, scale,  small,    &
                       ssum,  xscl1,  xscl2
      INTEGER       ::  i,  j,  IMINWR,  INFO1, INFO2,  &
                        LWRKEV, LWRSDD, LWRSVD, LWRSVJ, &
                        lwrsvq, mlwork, mwrkev, mwrsdd, &
                        mwrsvd, mwrsvj, mwrsvq, numrnk, &
                        olwork, mlrwrk
      LOGICAL       ::  BADXY, LQUERY, SCCOLX, SCCOLY,  &
                        wntex, wntref, wntres, wntvec
      CHARACTER     ::  JOBZL, T_OR_N
      CHARACTER     ::  JSVOPT
!
!     Local arrays
!     ~~~~~~~~~~~~
      REAL(KIND=wp) :: rdummy(2)
!
!     External functions (BLAS and LAPACK)
!     ~~~~~~~~~~~~~~~~~
      REAL(KIND=wp) zlange, dlamch, dznrm2
      EXTERNAL      zlange, dlamch, dznrm2, izamax
      INTEGER                               IZAMAX
      LOGICAL       DISNAN, LSAME
      EXTERNAL      disnan, lsame
!
!     External subroutines (BLAS and LAPACK)
!     ~~~~~~~~~~~~~~~~~~~~
      EXTERNAL      zaxpy,  zgemm,  zdscal
      EXTERNAL      zgeev,  zgejsv, zgesdd, zgesvd, zgesvdq, &
                    zlacpy, zlascl, zlassq, xerbla
!
!     Intrinsic functions
!     ~~~~~~~~~~~~~~~~~~~
      INTRINSIC     dble, int, max, sqrt
!............................................................
!
!    Test the input arguments
!
      wntres = lsame(jobr,'R')
      sccolx = lsame(jobs,'S') .OR. lsame(jobs,'C')
      sccoly = lsame(jobs,'Y')
      wntvec = lsame(jobz,'V')
      wntref = lsame(jobf,'R')
      wntex  = lsame(jobf,'E')
      info   = 0
      lquery = ( ( lzwork == -1 ) .OR. ( liwork == -1 ) &
                                  .OR. ( lrwork == -1 ) )
!
      IF ( .NOT. (sccolx .OR. sccoly .OR. &
                                  lsame(jobs,'N')) )   THEN
          info = -1
      ELSE IF ( .NOT. (wntvec .OR. lsame(jobz,'N')        &
                              .OR. lsame(jobz,'F')) )  THEN
          info = -2
      ELSE IF ( .NOT. (wntres .OR. lsame(jobr,'N')) .OR.  &
                ( wntres .AND. (.NOT.wntvec) ) )       THEN
          info = -3
      ELSE IF ( .NOT. (wntref .OR. wntex .OR.             &
                lsame(jobf,'N') ) )                    THEN
          info = -4
      ELSE IF ( .NOT.((whtsvd == 1) .OR. (whtsvd == 2) .OR.  &
                      (whtsvd == 3) .OR. (whtsvd == 4) )) THEN
          info = -5
      ELSE IF ( m < 0 )   THEN
          info = -6
      ELSE IF ( ( n < 0 ) .OR. ( n > m ) ) THEN
          info = -7
      ELSE IF ( ldx < m ) THEN
          info = -9
      ELSE IF ( ldy < m ) THEN
          info = -11
      ELSE IF ( .NOT. (( nrnk == -2).OR.(nrnk == -1).OR. &
                ((nrnk >= 1).AND.(nrnk <=n ))) )      THEN
          info = -12
      ELSE IF ( ( tol < zero ) .OR. ( tol >= one ) )  THEN
          info = -13
      ELSE IF ( ldz < m ) THEN
          info = -17
      ELSE IF ( (wntref .OR. wntex ) .AND. ( ldb < m ) ) THEN
          info = -20
      ELSE IF ( ldw < n ) THEN
          info = -22
      ELSE IF ( lds < n ) THEN
          info = -24
      END IF
!
      IF ( info == 0 ) THEN
          ! Compute the minimal and the optimal workspace
          ! requirements. Simulate running the code and
          ! determine minimal and optimal sizes of the
          ! workspace at any moment of the run.
         IF ( n == 0 ) THEN
             ! Quick return. All output except K is void.
             ! INFO=1 signals the void input.
             ! In case of a workspace query, the default
             ! minimal workspace lengths are returned.
            IF ( lquery ) THEN
                iwork(1) = 1
                rwork(1) = 1
                zwork(1) = 2
                zwork(2) = 2
            ELSE
               k   =  0
            END IF
            info = 1
            RETURN
         END IF
 
         iminwr = 1
         mlrwrk = max(1,n)
         mlwork = 2
         olwork = 2
         SELECT CASE ( whtsvd )
         CASE (1)
             ! The following is specified as the minimal
             ! length of WORK in the definition of ZGESVD:
             ! MWRSVD = MAX(1,2*MIN(M,N)+MAX(M,N))
             mwrsvd = max(1,2*min(m,n)+max(m,n))
             mlwork = max(mlwork,mwrsvd)
             mlrwrk = max(mlrwrk,n + 5*min(m,n))
             IF ( lquery ) THEN
                CALL zgesvd( 'O', 'S', m, n, x, ldx, rwork, &
                     b, ldb, w, ldw, zwork, -1, rdummy, info1 )
                lwrsvd = int( zwork(1) )
                olwork = max(olwork,lwrsvd)
             END IF
         CASE (2)
             ! The following is specified as the minimal
             ! length of WORK in the definition of ZGESDD:
             ! MWRSDD = 2*min(M,N)*min(M,N)+2*min(M,N)+max(M,N).
             ! RWORK length: 5*MIN(M,N)*MIN(M,N)+7*MIN(M,N)
             ! In LAPACK 3.10.1 RWORK is defined differently.
             ! Below we take max over the two versions.
             ! IMINWR = 8*MIN(M,N)
             mwrsdd = 2*min(m,n)*min(m,n)+2*min(m,n)+max(m,n)
             mlwork = max(mlwork,mwrsdd)
             iminwr = 8*min(m,n)
             mlrwrk = max( mlrwrk,  n +                    &
                      max( 5*min(m,n)*min(m,n)+7*min(m,n), &
                           5*min(m,n)*min(m,n)+5*min(m,n), &
                           2*max(m,n)*min(m,n)+            &
                           2*min(m,n)*min(m,n)+min(m,n) ) )
             IF ( lquery ) THEN
                CALL zgesdd( 'O', m, n, x, ldx, rwork, b,ldb,&
                     w, ldw, zwork, -1, rdummy, iwork, info1 )
                lwrsdd = max( mwrsdd,int( zwork(1) ))
                ! Possible bug in ZGESDD optimal workspace size.
                olwork = max(olwork,lwrsdd)
             END IF
         CASE (3)
             CALL zgesvdq( 'H', 'P', 'N', 'R', 'R', m, n, &
                  x, ldx, rwork, z, ldz, w, ldw, numrnk,  &
                  iwork, -1, zwork, -1, rdummy, -1, info1 )
             iminwr = iwork(1)
             mwrsvq = int(zwork(2))
             mlwork = max(mlwork,mwrsvq)
             mlrwrk  = max(mlrwrk,n + int(rdummy(1)))
             IF ( lquery ) THEN
                lwrsvq = int(zwork(1))
                olwork = max(olwork,lwrsvq)
             END IF
         CASE (4)
             jsvopt = 'J'
             CALL zgejsv( 'F', 'U', jsvopt, 'R', 'N', 'P', m, &
                   n, x, ldx, rwork, z, ldz, w, ldw,       &
                   zwork, -1, rdummy, -1, iwork, info1 )
             iminwr = iwork(1)
             mwrsvj = int(zwork(2))
             mlwork = max(mlwork,mwrsvj)
             mlrwrk = max(mlrwrk,n + max(7,int(rdummy(1))))
             IF ( lquery ) THEN
                lwrsvj = int(zwork(1))
                olwork = max(olwork,lwrsvj)
             END IF
         END SELECT
         IF ( wntvec .OR. wntex .OR. lsame(jobz,'F') ) THEN
             jobzl = 'V'
         ELSE
             jobzl = 'N'
         END IF
         ! Workspace calculation to the ZGEEV call
         mwrkev = max( 1, 2*n )
         mlwork = max(mlwork,mwrkev)
         mlrwrk = max(mlrwrk,n+2*n)
         IF ( lquery ) THEN
             CALL zgeev( 'N', jobzl, n, s, lds, eigs, &
              w, ldw, w, ldw, zwork, -1, rwork, info1 )
                lwrkev = int(zwork(1))
                olwork = max( olwork, lwrkev )
         END IF
!
         IF ( liwork < iminwr .AND. (.NOT.lquery) ) info = -30
         IF ( lrwork < mlrwrk .AND. (.NOT.lquery) ) info = -28
         IF ( lzwork < mlwork .AND. (.NOT.lquery) ) info = -26
 
      END IF
!
      IF( info /= 0 ) THEN
         CALL xerbla( 'ZGEDMD', -info )
         RETURN
      ELSE IF ( lquery ) THEN
!     Return minimal and optimal workspace sizes
          iwork(1) = iminwr
          rwork(1) = mlrwrk
          zwork(1) = mlwork
          zwork(2) = olwork
          RETURN
      END IF
!............................................................
!
      ofl   = dlamch('O')
      small = dlamch('S')
      badxy = .false.
!
!     <1> Optional scaling of the snapshots (columns of X, Y)
!     ==========================================================
      IF ( sccolx ) THEN
          ! The columns of X will be normalized.
          ! To prevent overflows, the column norms of X are
          ! carefully computed using ZLASSQ.
          k = 0
          DO i = 1, n
            !WORK(i) = DZNRM2( M, X(1,i), 1 )
            ssum  = one
            scale = zero
            CALL zlassq( m, x(1,i), 1, scale, ssum )
            IF ( disnan(scale) .OR. disnan(ssum) ) THEN
                k    =  0
                info = -8
                CALL xerbla('ZGEDMD',-info)
            END IF
            IF ( (scale /= zero) .AND. (ssum /= zero) ) THEN
               rootsc = sqrt(ssum)
               IF ( scale .GE. (ofl / rootsc) ) THEN
!                 Norm of X(:,i) overflows. First, X(:,i)
!                 is scaled by
!                 ( ONE / ROOTSC ) / SCALE = 1/||X(:,i)||_2.
!                 Next, the norm of X(:,i) is stored without
!                 overflow as RWORK(i) = - SCALE * (ROOTSC/M),
!                 the minus sign indicating the 1/M factor.
!                 Scaling is performed without overflow, and
!                 underflow may occur in the smallest entries
!                 of X(:,i). The relative backward and forward
!                 errors are small in the ell_2 norm.
                  CALL zlascl( 'G', 0, 0, scale, one/rootsc, &
                               m, 1, x(1,i), ldx, info2 )
                  rwork(i) = - scale * ( rootsc / dble(m) )
               ELSE
!                 X(:,i) will be scaled to unit 2-norm
                  rwork(i) =   scale * rootsc
                  CALL zlascl( 'G',0, 0, rwork(i), one, m, 1, &
                               x(1,i), ldx, info2 )      ! LAPACK CALL
!                 X(1:M,i) = (ONE/RWORK(i)) * X(1:M,i)   ! INTRINSIC
               END IF
            ELSE
               rwork(i) = zero
               k = k + 1
            END IF
          END DO
          IF ( k == n ) THEN
             ! All columns of X are zero. Return error code -8.
             ! (the 8th input variable had an illegal value)
             k = 0
             info = -8
             CALL xerbla('ZGEDMD',-info)
             RETURN
          END IF
          DO i = 1, n
!           Now, apply the same scaling to the columns of Y.
            IF ( rwork(i) >  zero ) THEN
                CALL zdscal( m, one/rwork(i), y(1,i), 1 )  ! BLAS CALL
!               Y(1:M,i) = (ONE/RWORK(i)) * Y(1:M,i)       ! INTRINSIC
            ELSE IF ( rwork(i) < zero ) THEN
                CALL zlascl( 'G', 0, 0, -rwork(i),          &
                     one/dble(m), m, 1, y(1,i), ldy, info2 ) ! LAPACK CALL
            ELSE IF ( abs(y(izamax(m, y(1,i),1),i ))  &
                                            /= zero ) THEN
!               X(:,i) is zero vector. For consistency,
!               Y(:,i) should also be zero. If Y(:,i) is not
!               zero, then the data might be inconsistent or
!               corrupted. If JOBS == 'C', Y(:,i) is set to
!               zero and a warning flag is raised.
!               The computation continues but the
!               situation will be reported in the output.
                badxy = .true.
                IF ( lsame(jobs,'C')) &
                CALL zdscal( m, zero, y(1,i), 1 )  ! BLAS CALL
            END IF
          END DO
      END IF
  !
      IF ( sccoly ) THEN
          ! The columns of Y will be normalized.
          ! To prevent overflows, the column norms of Y are
          ! carefully computed using ZLASSQ.
          DO i = 1, n
            !RWORK(i) = DZNRM2( M, Y(1,i), 1 )
            ssum  = one
            scale = zero
            CALL zlassq( m, y(1,i), 1, scale, ssum )
            IF ( disnan(scale) .OR. disnan(ssum) ) THEN
                k    =  0
                info = -10
                CALL xerbla('ZGEDMD',-info)
            END IF
            IF ( scale /= zero  .AND. (ssum /= zero) ) THEN
               rootsc = sqrt(ssum)
               IF ( scale .GE. (ofl / rootsc) ) THEN
!                 Norm of Y(:,i) overflows. First, Y(:,i)
!                 is scaled by
!                 ( ONE / ROOTSC ) / SCALE = 1/||Y(:,i)||_2.
!                 Next, the norm of Y(:,i) is stored without
!                 overflow as RWORK(i) = - SCALE * (ROOTSC/M),
!                 the minus sign indicating the 1/M factor.
!                 Scaling is performed without overflow, and
!                 underflow may occur in the smallest entries
!                 of Y(:,i). The relative backward and forward
!                 errors are small in the ell_2 norm.
                  CALL zlascl( 'G', 0, 0, scale, one/rootsc, &
                               m, 1, y(1,i), ldy, info2 )
                  rwork(i) = - scale * ( rootsc / dble(m) )
               ELSE
!                 Y(:,i) will be scaled to unit 2-norm
                  rwork(i) =   scale * rootsc
                  CALL zlascl( 'G',0, 0, rwork(i), one, m, 1, &
                               y(1,i), ldy, info2 )             ! LAPACK CALL
!                 Y(1:M,i) = (ONE/RWORK(i)) * Y(1:M,i)          ! INTRINSIC
               END IF
            ELSE
               rwork(i) = zero
            END IF
         END DO
         DO i = 1, n
!           Now, apply the same scaling to the columns of X.
            IF ( rwork(i) >  zero ) THEN
                CALL zdscal( m, one/rwork(i), x(1,i), 1 ) ! BLAS CALL
!               X(1:M,i) = (ONE/RWORK(i)) * X(1:M,i)      ! INTRINSIC
            ELSE IF ( rwork(i) < zero ) THEN
                CALL zlascl( 'G', 0, 0, -rwork(i),          &
                     one/dble(m), m, 1, x(1,i), ldx, info2 ) ! LAPACK CALL
            ELSE IF ( abs(x(izamax(m, x(1,i),1),i ))  &
                                           /= zero ) THEN
!               Y(:,i) is zero vector.  If X(:,i) is not
!               zero, then a warning flag is raised.
!               The computation continues but the
!               situation will be reported in the output.
                badxy = .true.
            END IF
         END DO
       END IF
!
!     <2> SVD of the data snapshot matrix X.
!     =====================================
!     The left singular vectors are stored in the array X.
!     The right singular vectors are in the array W.
!     The array W will later on contain the eigenvectors
!     of a Rayleigh quotient.
      numrnk = n
      SELECT CASE ( whtsvd )
         CASE (1)
             CALL zgesvd( 'O', 'S', m, n, x, ldx, rwork, b, &
                  ldb, w, ldw, zwork, lzwork,  rwork(n+1), info1 ) ! LAPACK CALL
             t_or_n = 'C'
         CASE (2)
            CALL zgesdd( 'O', m, n, x, ldx, rwork, b, ldb, w, &
                 ldw, zwork, lzwork, rwork(n+1), iwork, info1 )   ! LAPACK CALL
            t_or_n = 'C'
         CASE (3)
              CALL zgesvdq( 'H', 'P', 'N', 'R', 'R', m, n, &
                   x, ldx, rwork, z, ldz, w, ldw, &
                   numrnk, iwork, liwork, zwork,     &
                   lzwork, rwork(n+1), lrwork-n, info1)     ! LAPACK CALL
              CALL zlacpy( 'A', m, numrnk, z, ldz, x, ldx )   ! LAPACK CALL
         t_or_n = 'C'
         CASE (4)
              CALL zgejsv( 'F', 'U', jsvopt, 'R', 'N', 'P', m, &
                   n, x, ldx, rwork, z, ldz, w, ldw, &
                   zwork, lzwork, rwork(n+1), lrwork-n, iwork, info1 )    ! LAPACK CALL
              CALL zlacpy( 'A', m, n, z, ldz, x, ldx )   ! LAPACK CALL
              t_or_n = 'N'
              xscl1 = rwork(n+1)
              xscl2 = rwork(n+2)
              IF ( xscl1 /=  xscl2 ) THEN
                 ! This is an exceptional situation. If the
                 ! data matrices are not scaled and the
                 ! largest singular value of X overflows.
                 ! In that case ZGEJSV can return the SVD
                 ! in scaled form. The scaling factor can be used
                 ! to rescale the data (X and Y).
                 CALL zlascl( 'G', 0, 0, xscl1, xscl2, m, n, y, ldy, info2  )
              END IF
      END SELECT
!
      IF ( info1 > 0 ) THEN
         ! The SVD selected subroutine did not converge.
         ! Return with an error code.
         info = 2
         RETURN
      END IF
!
      IF ( rwork(1) == zero ) THEN
          ! The largest computed singular value of (scaled)
          ! X is zero. Return error code -8
          ! (the 8th input variable had an illegal value).
          k = 0
          info = -8
          CALL xerbla('ZGEDMD',-info)
          RETURN
      END IF
!
      
      !    snapshots matrix X. This depends on the
      !    parameters NRNK and TOL.
 
      SELECT CASE ( nrnk )
          CASE ( -1 )
               k = 1
               DO i = 2, numrnk
                 IF ( ( rwork(i) <= rwork(1)*tol ) .OR. &
                      ( rwork(i) <= small ) ) EXIT
                 k = k + 1
               END DO
          CASE ( -2 )
               k = 1
               DO i = 1, numrnk-1
                 IF ( ( rwork(i+1) <= rwork(i)*tol  ) .OR. &
                      ( rwork(i) <= small ) ) EXIT
                 k = k + 1
               END DO
          CASE DEFAULT
               k = 1
               DO i = 2, nrnk
                  IF ( rwork(i) <= small ) EXIT
                  k = k + 1
               END DO
          END SELECT
      !   Now, U = X(1:M,1:K) is the SVD/POD basis for the
      !   snapshot data in the input matrix X.
 
      
      !    Depending on the requested outputs, the computation
      !    is organized to compute additional auxiliary
      !    matrices (for the residuals and refinements).
      !
      !    In all formulas below, we need V_k*Sigma_k^(-1)
      !    where either V_k is in W(1:N,1:K), or V_k^H is in
      !    W(1:K,1:N). Here Sigma_k=diag(WORK(1:K)).
      IF ( lsame(t_or_n, 'N') ) THEN
          DO i = 1, k
           CALL zdscal( n, one/rwork(i), w(1,i), 1 )    ! BLAS CALL
           ! W(1:N,i) = (ONE/RWORK(i)) * W(1:N,i)      ! INTRINSIC
          END DO
      ELSE
          ! This non-unit stride access is due to the fact
          ! that ZGESVD, ZGESVDQ and ZGESDD return the
          ! adjoint matrix of the right singular vectors.
          !DO i = 1, K
          ! CALL ZDSCAL( N, ONE/RWORK(i), W(i,1), LDW )    ! BLAS CALL
          ! ! W(i,1:N) = (ONE/RWORK(i)) * W(i,1:N)      ! INTRINSIC
          !END DO
          DO i = 1, k
              rwork(n+i) = one/rwork(i)
          END DO
          DO j = 1, n
             DO i = 1, k
                 w(i,j) = cmplx(rwork(n+i),zero,kind=wp)*w(i,j)
             END DO
          END DO
      END IF
!
      IF ( wntref ) THEN
         !
         ! Need A*U(:,1:K)=Y*V_k*inv(diag(WORK(1:K)))
         ! for computing the refined Ritz vectors
         ! (optionally, outside ZGEDMD).
          CALL zgemm( 'N', t_or_n, m, k, n, zone, y, ldy, w, &
                      ldw, zzero, z, ldz )                       ! BLAS CALL
          ! Z(1:M,1:K)=MATMUL(Y(1:M,1:N),TRANSPOSE(CONJG(W(1:K,1:N)))) ! INTRINSIC, for T_OR_N=='C'
          ! Z(1:M,1:K)=MATMUL(Y(1:M,1:N),W(1:N,1:K))                   ! INTRINSIC, for T_OR_N=='N'
          !
          ! At this point Z contains
          ! A * U(:,1:K) = Y * V_k * Sigma_k^(-1), and
          ! this is needed for computing the residuals.
          ! This matrix is  returned in the array B and
          ! it can be used to compute refined Ritz vectors.
          CALL zlacpy( 'A', m, k, z, ldz, b, ldb )   ! BLAS CALL
          ! B(1:M,1:K) = Z(1:M,1:K)                  ! INTRINSIC
 
          CALL zgemm( 'C', 'N', k, k, m, zone, x, ldx, z, &
                      ldz, zzero, s, lds )                        ! BLAS CALL
          ! S(1:K,1:K) = MATMUL(TRANSPOSE(CONJG(X(1:M,1:K))),Z(1:M,1:K)) ! INTRINSIC
          ! At this point S = U^H * A * U is the Rayleigh quotient.
      ELSE
        ! A * U(:,1:K) is not explicitly needed and the
        ! computation is organized differently. The Rayleigh
        ! quotient is computed more efficiently.
        CALL zgemm( 'C', 'N', k, n, m, zone, x, ldx, y, ldy, &
                   zzero, z, ldz )                                         ! BLAS CALL
        ! Z(1:K,1:N) = MATMUL( TRANSPOSE(CONJG(X(1:M,1:K))), Y(1:M,1:N) )  ! INTRINSIC
        !
        CALL zgemm( 'N', t_or_n, k, k, n, zone, z, ldz, w, &
                    ldw, zzero, s, lds )                         ! BLAS CALL
        ! S(1:K,1:K) = MATMUL(Z(1:K,1:N),TRANSPOSE(CONJG(W(1:K,1:N)))) ! INTRINSIC, for T_OR_N=='T'
        ! S(1:K,1:K) = MATMUL(Z(1:K,1:N),(W(1:N,1:K)))                 ! INTRINSIC, for T_OR_N=='N'
        ! At this point S = U^H * A * U is the Rayleigh quotient.
        ! If the residuals are requested, save scaled V_k into Z.
        ! Recall that V_k or V_k^H is stored in W.
        IF ( wntres .OR. wntex ) THEN
          IF ( lsame(t_or_n, 'N') ) THEN
              CALL zlacpy( 'A', n, k, w, ldw, z, ldz )
          ELSE
              CALL zlacpy( 'A', k, n, w, ldw, z, ldz )
          END IF
        END IF
      END IF
!
      
      !   right eigenvectors of the Rayleigh quotient.
      !
      CALL zgeev( 'N', jobzl, k, s, lds, eigs, w, ldw, &
            w, ldw, zwork, lzwork, rwork(n+1), info1 )  ! LAPACK CALL
      !
      ! W(1:K,1:K) contains the eigenvectors of the Rayleigh
      ! quotient.  See the description of Z.
      ! Also, see the description of ZGEEV.
      IF ( info1 > 0 ) THEN
         ! ZGEEV failed to compute the eigenvalues and
         ! eigenvectors of the Rayleigh quotient.
         info = 3
         RETURN
      END IF
!
      ! <6> Compute the eigenvectors (if requested) and,
      ! the residuals (if requested).
      !
      IF ( wntvec .OR. wntex ) THEN
      IF ( wntres ) THEN
          IF ( wntref ) THEN
            ! Here, if the refinement is requested, we have
            ! A*U(:,1:K) already computed and stored in Z.
            ! For the residuals, need Y = A * U(:,1;K) * W.
            CALL zgemm( 'N', 'N', m, k, k, zone, z, ldz, w, &
                       ldw, zzero, y, ldy )               ! BLAS CALL
            ! Y(1:M,1:K) = Z(1:M,1:K) * W(1:K,1:K)        ! INTRINSIC
            ! This frees Z; Y contains A * U(:,1:K) * W.
          ELSE
            ! Compute S = V_k * Sigma_k^(-1) * W, where
            ! V_k * Sigma_k^(-1) (or its adjoint) is stored in Z
            CALL zgemm( t_or_n, 'N', n, k, k, zone, z, ldz, &
                       w, ldw, zzero, s, lds )
            ! Then, compute Z = Y * S =
            ! = Y * V_k * Sigma_k^(-1) * W(1:K,1:K) =
            ! = A * U(:,1:K) * W(1:K,1:K)
            CALL zgemm( 'N', 'N', m, k, n, zone, y, ldy, s, &
                       lds, zzero, z, ldz )
            ! Save a copy of Z into Y and free Z for holding
            ! the Ritz vectors.
            CALL zlacpy( 'A', m, k, z, ldz, y, ldy )
            IF ( wntex ) CALL zlacpy( 'A', m, k, z, ldz, b, ldb )
          END IF
      ELSE IF ( wntex ) THEN
          ! Compute S = V_k * Sigma_k^(-1) * W, where
            ! V_k * Sigma_k^(-1) is stored in Z
            CALL zgemm( t_or_n, 'N', n, k, k, zone, z, ldz, &
                       w, ldw, zzero, s, lds )
            ! Then, compute Z = Y * S =
            ! = Y * V_k * Sigma_k^(-1) * W(1:K,1:K) =
            ! = A * U(:,1:K) * W(1:K,1:K)
            CALL zgemm( 'N', 'N', m, k, n, zone, y, ldy, s, &
                       lds, zzero, b, ldb )
            ! The above call replaces the following two calls
            ! that were used in the developing-testing phase.
            ! CALL ZGEMM( 'N', 'N', M, K, N, ZONE, Y, LDY, S, &
            !           LDS, ZZERO, Z, LDZ)
            ! Save a copy of Z into B and free Z for holding
            ! the Ritz vectors.
            ! CALL ZLACPY( 'A', M, K, Z, LDZ, B, LDB )
      END IF
!
      ! Compute the Ritz vectors
      IF ( wntvec ) CALL zgemm( 'N', 'N', m, k, k, zone, x, ldx, w, ldw, &
                   zzero, z, ldz )                          ! BLAS CALL
      ! Z(1:M,1:K) = MATMUL(X(1:M,1:K), W(1:K,1:K))         ! INTRINSIC
!
      IF ( wntres ) THEN
         DO i = 1, k
            CALL zaxpy( m, -eigs(i), z(1,i), 1, y(1,i), 1 )       ! BLAS CALL
            ! Y(1:M,i) = Y(1:M,i) - EIGS(i) * Z(1:M,i)            ! INTRINSIC
            res(i) = dznrm2( m, y(1,i), 1 )                       ! BLAS CALL
         END DO
      END IF
      END IF
!
      IF ( whtsvd == 4 ) THEN
          rwork(n+1) = xscl1
          rwork(n+2) = xscl2
      END IF
!
!     Successful exit.
      IF ( .NOT. badxy ) THEN
         info = 0
      ELSE
         ! A warning on possible data inconsistency.
         ! This should be a rare event.
         info = 4
      END IF
!............................................................
      RETURN
!     ......
      END SUBROUTINE zgedmd