dgedmd.f90

Go to the documentation of this file.

 
!
!  =========== DOCUMENTATION ===========
!
!  Definition:
!  ===========
!
!     SUBROUTINE DGEDMD( JOBS, JOBZ, JOBR, JOBF,  WHTSVD,  &
!                        M, N, X, LDX, Y, LDY, NRNK, TOL,  &
!                        K, REIG,  IMEIG,   Z, LDZ,  RES,  &
!                        B, LDB, W,  LDW,   S, LDS,        &
!                        WORK, LWORK, IWORK, LIWORK, INFO )
!.....
!     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, &
!                                LWORK,  LIWORK
!     INTEGER,   INTENT(OUT)  :: K, INFO
!     REAL(KIND=WP), INTENT(IN)  :: TOL
!     Array arguments
!     REAL(KIND=WP), INTENT(INOUT) :: X(LDX,*), Y(LDY,*)
!     REAL(KIND=WP), INTENT(OUT)   :: Z(LDZ,*), B(LDB,*), &
!                                     W(LDW,*), S(LDS,*)
!     REAL(KIND=WP), INTENT(OUT)   :: REIG(*),  IMEIG(*), &
!                                     RES(*)
!     REAL(KIND=WP), INTENT(OUT)   :: WORK(*)
!     INTEGER,       INTENT(OUT)   :: IWORK(*)
!
!............................................................
!     =============
!............................................................
!     ================
!......................................................................
!     ================================
!......................................................................
!     ==============================
!......................................................................
!     Arguments
!     =========
!
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!
!  Authors:
!  ========
!
!
!
!.............................................................
!.............................................................
      SUBROUTINE dgedmd( JOBS, JOBZ, JOBR, JOBF,  WHTSVD,  &
                         M, N, X, LDX, Y, LDY, NRNK, TOL,  &
                         K, REIG,  IMEIG,   Z, LDZ,  RES,  &
                         B, LDB, W,  LDW,   S, LDS,        &
                         WORK, LWORK, IWORK, LIWORK, INFO )
!
!  -- 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, &
                                 LWORK,  LIWORK
      INTEGER,   INTENT(OUT)  :: K, INFO
      REAL(KIND=wp), INTENT(IN)  :: tol
!
!     Array arguments
!     ~~~~~~~~~~~~~~~
      REAL(KIND=wp), INTENT(INOUT) :: x(ldx,*), y(ldy,*)
      REAL(KIND=wp), INTENT(OUT)   :: z(ldz,*), b(ldb,*), &
                                      w(ldw,*), s(lds,*)
      REAL(KIND=wp), INTENT(OUT)   :: reig(*),  imeig(*), &
                                      res(*)
      REAL(KIND=wp), INTENT(OUT)   :: work(*)
      INTEGER,       INTENT(OUT)   :: IWORK(*)
!
!     Parameters
!     ~~~~~~~~~~
      REAL(KIND=wp), PARAMETER ::  one = 1.0_wp
      REAL(KIND=wp), PARAMETER :: zero = 0.0_wp
!
!     Local scalars
!     ~~~~~~~~~~~~~
      REAL(KIND=wp) :: ofl,    rootsc, scale,  small,  &
                       ssum,   xscl1,  xscl2
      INTEGER       :: i,   j, IMINWR,  INFO1, INFO2,  &
                       LWRKEV, LWRSDD, LWRSVD,         &
                       lwrsvq, mlwork, mwrkev, mwrsdd, &
                       mwrsvd, mwrsvj, mwrsvq, numrnk, &
                       olwork
      LOGICAL       :: BADXY,  LQUERY, SCCOLX, SCCOLY, &
                       wntex,  wntref, wntres, wntvec
      CHARACTER     :: JOBZL,  T_OR_N
      CHARACTER     :: JSVOPT
!
!     Local arrays
!     ~~~~~~~~~~~~
      REAL(KIND=wp) :: ab(2,2), rdummy(2), rdummy2(2)
!
!     External functions (BLAS and LAPACK)
!     ~~~~~~~~~~~~~~~~~
      REAL(KIND=wp) dlange, dlamch, dnrm2
      EXTERNAL      dlange, dlamch, dnrm2, idamax
      INTEGER       IDAMAX
      LOGICAL       DISNAN, LSAME
      EXTERNAL      disnan, lsame
!
!     External subroutines (BLAS and LAPACK)
!     ~~~~~~~~~~~~~~~~~~~~
      EXTERNAL      daxpy,  dgemm,  dscal
      EXTERNAL      dgeev,  dgejsv, dgesdd, dgesvd, dgesvdq, &
                    dlacpy, dlascl, dlassq, xerbla
!
!     Intrinsic functions
!     ~~~~~~~~~~~~~~~~~~~
      INTRINSIC     dble, int, max, sqrt
!............................................................
!
!    Test the input arguments
!
      wntres = lsame(jobr,'R')
      sccolx = lsame(jobs,'S') .OR. lsame(jobs,'C')
      sccoly = lsame(jobs,'Y')
      wntvec = lsame(jobz,'V')
      wntref = lsame(jobf,'R')
      wntex  = lsame(jobf,'E')
      info   = 0
      lquery = ( ( lwork == -1 ) .OR. ( liwork == -1 ) )
!
      IF ( .NOT. (sccolx .OR. sccoly .OR. &
                                  lsame(jobs,'N')) )   THEN
          info = -1
      ELSE IF ( .NOT. (wntvec .OR. lsame(jobz,'N')        &
                              .OR. lsame(jobz,'F')) )  THEN
          info = -2
      ELSE IF ( .NOT. (wntres .OR. lsame(jobr,'N')) .OR.  &
                ( wntres .AND. (.NOT.wntvec) ) )       THEN
          info = -3
      ELSE IF ( .NOT. (wntref .OR. wntex .OR.             &
                lsame(jobf,'N') ) )                    THEN
          info = -4
      ELSE IF ( .NOT.((whtsvd == 1) .OR. (whtsvd == 2) .OR.  &
                      (whtsvd == 3) .OR. (whtsvd == 4) )) THEN
          info = -5
      ELSE IF ( m < 0 )   THEN
          info = -6
      ELSE IF ( ( n < 0 ) .OR. ( n > m ) ) THEN
          info = -7
      ELSE IF ( ldx < m ) THEN
          info = -9
      ELSE IF ( ldy < m ) THEN
          info = -11
      ELSE IF ( .NOT. (( nrnk == -2).OR.(nrnk == -1).OR. &
                ((nrnk >= 1).AND.(nrnk <=n ))) )      THEN
          info = -12
      ELSE IF ( ( tol < zero ) .OR. ( tol >= one ) )  THEN
          info = -13
      ELSE IF ( ldz < m ) THEN
          info = -18
      ELSE IF ( (wntref .OR. wntex ) .AND. ( ldb < m ) ) THEN
          info = -21
      ELSE IF ( ldw < n ) THEN
          info = -23
      ELSE IF ( lds < n ) THEN
          info = -25
      END IF
!
      IF ( info == 0 ) THEN
          ! Compute the minimal and the optimal workspace
          ! requirements. Simulate running the code and
          ! determine minimal and optimal sizes of the
          ! workspace at any moment of the run.
         IF ( n == 0 ) THEN
             ! Quick return. All output except K is void.
             ! INFO=1 signals the void input.
             ! In case of a workspace query, the default
             ! minimal workspace lengths are returned.
            IF ( lquery ) THEN
                iwork(1) = 1
                work(1)  = 2
                work(2)  = 2
            ELSE
               k = 0
            END IF
            info = 1
            RETURN
         END IF
         mlwork = max(2,n)
         olwork = max(2,n)
         iminwr = 1
         SELECT CASE ( whtsvd )
         CASE (1)
             ! The following is specified as the minimal
             ! length of WORK in the definition of DGESVD:
             ! MWRSVD = MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N))
             mwrsvd = max(1,3*min(m,n)+max(m,n),5*min(m,n))
             mlwork = max(mlwork,n + mwrsvd)
             IF ( lquery ) THEN
                CALL dgesvd( 'O', 'S', m, n, x, ldx, work, &
                           b, ldb, w, ldw, rdummy, -1, info1 )
                lwrsvd = max( mwrsvd, int( rdummy(1) ) )
                olwork = max(olwork,n + lwrsvd)
             END IF
         CASE (2)
             ! The following is specified as the minimal
             ! length of WORK in the definition of DGESDD:
             ! MWRSDD = 3*MIN(M,N)*MIN(M,N) +
             ! MAX( MAX(M,N),5*MIN(M,N)*MIN(M,N)+4*MIN(M,N) )
             ! IMINWR = 8*MIN(M,N)
             mwrsdd = 3*min(m,n)*min(m,n) +                &
              max( max(m,n),5*min(m,n)*min(m,n)+4*min(m,n) )
             mlwork = max(mlwork,n + mwrsdd)
             iminwr = 8*min(m,n)
             IF ( lquery ) THEN
                CALL dgesdd( 'O', m, n, x, ldx, work, b,     &
                     ldb, w, ldw, rdummy, -1, iwork, info1 )
                lwrsdd = max( mwrsdd, int( rdummy(1) ) )
                olwork = max(olwork,n + lwrsdd)
             END IF
         CASE (3)
             !LWQP3 = 3*N+1
             !LWORQ = MAX(N, 1)
             !MWRSVD = MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N))
             !MWRSVQ = N + MAX( LWQP3, MWRSVD, LWORQ ) + MAX(M,2)
             !MLWORK = N +  MWRSVQ
             !IMINWR = M+N-1
             CALL dgesvdq( 'H', 'P', 'N', 'R', 'R', m, n, &
                             x, ldx, work, z, ldz, w, ldw,   &
                             numrnk, iwork, liwork, rdummy,  &
                             -1, rdummy2, -1, info1 )
             iminwr = iwork(1)
             mwrsvq = int(rdummy(2))
             mlwork = max(mlwork,n+mwrsvq+int(rdummy2(1)))
             IF ( lquery ) THEN
                lwrsvq = max( mwrsvq, int(rdummy(1)) )
                olwork = max(olwork,n+lwrsvq+int(rdummy2(1)))
             END IF
         CASE (4)
             jsvopt = 'J'
             !MWRSVJ = MAX( 7, 2*M+N, 6*N+2*N*N ) ! for JSVOPT='V'
             mwrsvj = max( 7, 2*m+n, 4*n+n*n, 2*n+n*n+6 )
             mlwork = max(mlwork,n+mwrsvj)
             iminwr = max( 3, m+3*n )
             IF ( lquery ) THEN
                olwork =  max(olwork,n+mwrsvj)
             END IF
         END SELECT
         IF ( wntvec .OR. wntex .OR. lsame(jobz,'F') ) THEN
             jobzl = 'V'
         ELSE
             jobzl = 'N'
         END IF
         ! Workspace calculation to the DGEEV call
         IF ( lsame(jobzl,'V') ) THEN
             mwrkev = max( 1, 4*n )
         ELSE
             mwrkev = max( 1, 3*n )
         END IF
         mlwork = max(mlwork,n+mwrkev)
         IF ( lquery ) THEN
                CALL dgeev( 'N', jobzl, n, s, lds, reig, &
                    imeig, w, ldw, w, ldw, rdummy, -1, info1 )
                lwrkev = max( mwrkev, int(rdummy(1)) )
                olwork = max( olwork, n+lwrkev )
         END IF
!
         IF ( liwork < iminwr .AND. (.NOT.lquery) ) info = -29
         IF (  lwork < mlwork .AND. (.NOT.lquery) ) info = -27
      END IF
!
      IF( info /= 0 ) THEN
         CALL xerbla( 'DGEDMD', -info )
         RETURN
      ELSE IF ( lquery ) THEN
!     Return minimal and optimal workspace sizes
          iwork(1) = iminwr
          work(1)  = mlwork
          work(2)  = olwork
          RETURN
      END IF
!............................................................
!
      ofl   = dlamch('O')
      small = dlamch('S')
      badxy = .false.
!
!     <1> Optional scaling of the snapshots (columns of X, Y)
!     ==========================================================
      IF ( sccolx ) THEN
          ! The columns of X will be normalized.
          ! To prevent overflows, the column norms of X are
          ! carefully computed using DLASSQ.
          k = 0
          DO i = 1, n
            !WORK(i) = DNRM2( M, X(1,i), 1 )
            ssum  = one
            scale = zero
            CALL dlassq( m, x(1,i), 1, scale, ssum )
            IF ( disnan(scale) .OR. disnan(ssum) ) THEN
                k    =  0
                info = -8
                CALL xerbla('DGEDMD',-info)
            END IF
            IF ( (scale /= zero) .AND. (ssum /= zero) ) THEN
               rootsc = sqrt(ssum)
               IF ( scale .GE. (ofl / rootsc) ) THEN
!                 Norm of X(:,i) overflows. First, X(:,i)
!                 is scaled by
!                 ( ONE / ROOTSC ) / SCALE = 1/||X(:,i)||_2.
!                 Next, the norm of X(:,i) is stored without
!                 overflow as WORK(i) = - SCALE * (ROOTSC/M),
!                 the minus sign indicating the 1/M factor.
!                 Scaling is performed without overflow, and
!                 underflow may occur in the smallest entries
!                 of X(:,i). The relative backward and forward
!                 errors are small in the ell_2 norm.
                  CALL dlascl( 'G', 0, 0, scale, one/rootsc, &
                               m, 1, x(1,i), m, info2 )
                  work(i) = - scale * ( rootsc / dble(m) )
               ELSE
!                 X(:,i) will be scaled to unit 2-norm
                  work(i) =   scale * rootsc
                  CALL dlascl( 'G',0, 0, work(i), one, m, 1, &
                               x(1,i), m, info2 )              ! LAPACK CALL
!                 X(1:M,i) = (ONE/WORK(i)) * X(1:M,i)          ! INTRINSIC
               END IF
            ELSE
               work(i) = zero
               k = k + 1
            END IF
          END DO
          IF ( k == n ) THEN
          ! All columns of X are zero. Return error code -8.
          ! (the 8th input variable had an illegal value)
          k = 0
          info = -8
          CALL xerbla('DGEDMD',-info)
          RETURN
          END IF
          DO i = 1, n
!           Now, apply the same scaling to the columns of Y.
            IF ( work(i) >  zero ) THEN
                CALL dscal( m, one/work(i), y(1,i), 1 )  ! BLAS CALL
!               Y(1:M,i) = (ONE/WORK(i)) * Y(1:M,i)      ! INTRINSIC
            ELSE IF ( work(i) < zero ) THEN
                CALL dlascl( 'G', 0, 0, -work(i),          &
                     one/dble(m), m, 1, y(1,i), m, info2 ) ! LAPACK CALL
            ELSE IF ( y(idamax(m, y(1,i),1),i )  &
                                            /= zero ) THEN
!               X(:,i) is zero vector. For consistency,
!               Y(:,i) should also be zero. If Y(:,i) is not
!               zero, then the data might be inconsistent or
!               corrupted. If JOBS == 'C', Y(:,i) is set to
!               zero and a warning flag is raised.
!               The computation continues but the
!               situation will be reported in the output.
                badxy = .true.
                IF ( lsame(jobs,'C')) &
                CALL dscal( m, zero, y(1,i), 1 )  ! BLAS CALL
            END IF
          END DO
      END IF
  !
      IF ( sccoly ) THEN
          ! The columns of Y will be normalized.
          ! To prevent overflows, the column norms of Y are
          ! carefully computed using DLASSQ.
          DO i = 1, n
            !WORK(i) = DNRM2( M, Y(1,i), 1 )
            ssum  = one
            scale = zero
            CALL dlassq( m, y(1,i), 1, scale, ssum )
            IF ( disnan(scale) .OR. disnan(ssum) ) THEN
                k    =  0
                info = -10
                CALL xerbla('DGEDMD',-info)
            END IF
            IF ( scale /= zero  .AND. (ssum /= zero) ) THEN
               rootsc = sqrt(ssum)
               IF ( scale .GE. (ofl / rootsc) ) THEN
!                 Norm of Y(:,i) overflows. First, Y(:,i)
!                 is scaled by
!                 ( ONE / ROOTSC ) / SCALE = 1/||Y(:,i)||_2.
!                 Next, the norm of Y(:,i) is stored without
!                 overflow as WORK(i) = - SCALE * (ROOTSC/M),
!                 the minus sign indicating the 1/M factor.
!                 Scaling is performed without overflow, and
!                 underflow may occur in the smallest entries
!                 of Y(:,i). The relative backward and forward
!                 errors are small in the ell_2 norm.
                  CALL dlascl( 'G', 0, 0, scale, one/rootsc, &
                               m, 1, y(1,i), m, info2 )
                  work(i) = - scale * ( rootsc / dble(m) )
               ELSE
!                 X(:,i) will be scaled to unit 2-norm
                  work(i) =   scale * rootsc
                  CALL dlascl( 'G',0, 0, work(i), one, m, 1, &
                               y(1,i), m, info2 )              ! LAPACK CALL
!                 Y(1:M,i) = (ONE/WORK(i)) * Y(1:M,i)          ! INTRINSIC
               END IF
            ELSE
               work(i) = zero
            END IF
         END DO
         DO i = 1, n
!           Now, apply the same scaling to the columns of X.
            IF ( work(i) >  zero ) THEN
                CALL dscal( m, one/work(i), x(1,i), 1 )  ! BLAS CALL
!               X(1:M,i) = (ONE/WORK(i)) * X(1:M,i)      ! INTRINSIC
            ELSE IF ( work(i) < zero ) THEN
                CALL dlascl( 'G', 0, 0, -work(i),          &
                     one/dble(m), m, 1, x(1,i), m, info2 ) ! LAPACK CALL
            ELSE IF ( x(idamax(m, x(1,i),1),i )  &
                                           /= zero ) THEN
!               Y(:,i) is zero vector.  If X(:,i) is not
!               zero, then a warning flag is raised.
!               The computation continues but the
!               situation will be reported in the output.
                badxy = .true.
            END IF
         END DO
       END IF
!
!     <2> SVD of the data snapshot matrix X.
!     =====================================
!     The left singular vectors are stored in the array X.
!     The right singular vectors are in the array W.
!     The array W will later on contain the eigenvectors
!     of a Rayleigh quotient.
      numrnk = n
      SELECT CASE ( whtsvd )
         CASE (1)
             CALL dgesvd( 'O', 'S', m, n, x, ldx, work, b, &
                  ldb, w, ldw, work(n+1), lwork-n, info1 ) ! LAPACK CALL
             t_or_n = 'T'
         CASE (2)
            CALL dgesdd( 'O', m, n, x, ldx, work, b, ldb, w, &
                 ldw, work(n+1), lwork-n, iwork, info1 )   ! LAPACK CALL
            t_or_n = 'T'
         CASE (3)
              CALL dgesvdq( 'H', 'P', 'N', 'R', 'R', m, n, &
                   x, ldx, work, z, ldz, w, ldw, &
                   numrnk, iwork, liwork, work(n+max(2,m)+1),&
                   lwork-n-max(2,m), work(n+1), max(2,m), info1)     ! LAPACK CALL
              CALL dlacpy( 'A', m, numrnk, z, ldz, x, ldx )   ! LAPACK CALL
         t_or_n = 'T'
         CASE (4)
              CALL dgejsv( 'F', 'U', jsvopt, 'N', 'N', 'P', m, &
                   n, x, ldx, work, z, ldz, w, ldw, &
                   work(n+1), lwork-n, iwork, info1 )    ! LAPACK CALL
              CALL dlacpy( 'A', m, n, z, ldz, x, ldx )   ! LAPACK CALL
              t_or_n = 'N'
              xscl1 = work(n+1)
              xscl2 = work(n+2)
              IF ( xscl1 /=  xscl2 ) THEN
                 ! This is an exceptional situation. If the
                 ! data matrices are not scaled and the
                 ! largest singular value of X overflows.
                 ! In that case DGEJSV can return the SVD
                 ! in scaled form. The scaling factor can be used
                 ! to rescale the data (X and Y).
                 CALL dlascl( 'G', 0, 0, xscl1, xscl2, m, n, y, ldy, info2  )
              END IF
      END SELECT
!
      IF ( info1 > 0 ) THEN
         ! The SVD selected subroutine did not converge.
         ! Return with an error code.
         info = 2
         RETURN
      END IF
!
      IF ( work(1) == zero ) THEN
          ! The largest computed singular value of (scaled)
          ! X is zero. Return error code -8
          ! (the 8th input variable had an illegal value).
          k = 0
          info = -8
          CALL xerbla('DGEDMD',-info)
          RETURN
      END IF
!
      
      !    snapshots matrix X. This depends on the
      !    parameters NRNK and TOL.
 
      SELECT CASE ( nrnk )
          CASE ( -1 )
               k = 1
               DO i = 2, numrnk
                 IF ( ( work(i) <= work(1)*tol ) .OR. &
                      ( work(i) <= small ) ) EXIT
                 k = k + 1
               END DO
          CASE ( -2 )
               k = 1
               DO i = 1, numrnk-1
                 IF ( ( work(i+1) <= work(i)*tol  ) .OR. &
                      ( work(i) <= small ) ) EXIT
                 k = k + 1
               END DO
          CASE DEFAULT
               k = 1
               DO i = 2, nrnk
                  IF ( work(i) <= small ) EXIT
                  k = k + 1
               END DO
          END SELECT
      !   Now, U = X(1:M,1:K) is the SVD/POD basis for the
      !   snapshot data in the input matrix X.
 
      
      !    Depending on the requested outputs, the computation
      !    is organized to compute additional auxiliary
      !    matrices (for the residuals and refinements).
      !
      !    In all formulas below, we need V_k*Sigma_k^(-1)
      !    where either V_k is in W(1:N,1:K), or V_k^T is in
      !    W(1:K,1:N). Here Sigma_k=diag(WORK(1:K)).
      IF ( lsame(t_or_n, 'N') ) THEN
          DO i = 1, k
           CALL dscal( n, one/work(i), w(1,i), 1 )    ! BLAS CALL
           ! W(1:N,i) = (ONE/WORK(i)) * W(1:N,i)      ! INTRINSIC
          END DO
      ELSE
          ! This non-unit stride access is due to the fact
          ! that DGESVD, DGESVDQ and DGESDD return the
          ! transposed matrix of the right singular vectors.
          !DO i = 1, K
          ! CALL DSCAL( N, ONE/WORK(i), W(i,1), LDW )    ! BLAS CALL
          ! ! W(i,1:N) = (ONE/WORK(i)) * W(i,1:N)      ! INTRINSIC
          !END DO
          DO i = 1, k
              work(n+i) = one/work(i)
          END DO
          DO j = 1, n
             DO i = 1, k
                 w(i,j) = (work(n+i))*w(i,j)
             END DO
          END DO
      END IF
!
      IF ( wntref ) THEN
         !
         ! Need A*U(:,1:K)=Y*V_k*inv(diag(WORK(1:K)))
         ! for computing the refined Ritz vectors
         ! (optionally, outside DGEDMD).
          CALL dgemm( 'N', t_or_n, m, k, n, one, y, ldy, w, &
                      ldw, zero, z, ldz )                        ! BLAS CALL
          ! Z(1:M,1:K)=MATMUL(Y(1:M,1:N),TRANSPOSE(W(1:K,1:N)))  ! INTRINSIC, for T_OR_N=='T'
          ! Z(1:M,1:K)=MATMUL(Y(1:M,1:N),W(1:N,1:K))             ! INTRINSIC, for T_OR_N=='N'
          !
          ! At this point Z contains
          ! A * U(:,1:K) = Y * V_k * Sigma_k^(-1), and
          ! this is needed for computing the residuals.
          ! This matrix is  returned in the array B and
          ! it can be used to compute refined Ritz vectors.
          CALL dlacpy( 'A', m, k, z, ldz, b, ldb )   ! BLAS CALL
          ! B(1:M,1:K) = Z(1:M,1:K)                  ! INTRINSIC
 
          CALL dgemm( 'T', 'N', k, k, m, one, x, ldx, z, &
                      ldz, zero, s, lds )                        ! BLAS CALL
          ! S(1:K,1:K) = MATMUL(TANSPOSE(X(1:M,1:K)),Z(1:M,1:K)) ! INTRINSIC
          ! At this point S = U^T * A * U is the Rayleigh quotient.
      ELSE
        ! A * U(:,1:K) is not explicitly needed and the
        ! computation is organized differently. The Rayleigh
        ! quotient is computed more efficiently.
        CALL dgemm( 'T', 'N', k, n, m, one, x, ldx, y, ldy, &
                   zero, z, ldz )                                   ! BLAS CALL
        ! Z(1:K,1:N) = MATMUL( TRANSPOSE(X(1:M,1:K)), Y(1:M,1:N) )  ! INTRINSIC
        ! In the two DGEMM calls here, can use K for LDZ.
        CALL dgemm( 'N', t_or_n, k, k, n, one, z, ldz, w, &
                    ldw, zero, s, lds )                         ! BLAS CALL
        ! S(1:K,1:K) = MATMUL(Z(1:K,1:N),TRANSPOSE(W(1:K,1:N))) ! INTRINSIC, for T_OR_N=='T'
        ! S(1:K,1:K) = MATMUL(Z(1:K,1:N),(W(1:N,1:K)))          ! INTRINSIC, for T_OR_N=='N'
        ! At this point S = U^T * A * U is the Rayleigh quotient.
        ! If the residuals are requested, save scaled V_k into Z.
        ! Recall that V_k or V_k^T is stored in W.
        IF ( wntres .OR. wntex ) THEN
          IF ( lsame(t_or_n, 'N') ) THEN
              CALL dlacpy( 'A', n, k, w, ldw, z, ldz )
          ELSE
              CALL dlacpy( 'A', k, n, w, ldw, z, ldz )
          END IF
        END IF
      END IF
!
      
      !   right eigenvectors of the Rayleigh quotient.
      !
      CALL dgeev( 'N', jobzl, k, s, lds, reig, imeig, w, &
                  ldw, w, ldw, work(n+1), lwork-n, info1 )   ! LAPACK CALL
      !
      ! W(1:K,1:K) contains the eigenvectors of the Rayleigh
      ! quotient. Even in the case of complex spectrum, all
      ! computation is done in real arithmetic. REIG and
      ! IMEIG are the real and the imaginary parts of the
      ! eigenvalues, so that the spectrum is given as
      ! REIG(:) + sqrt(-1)*IMEIG(:). Complex conjugate pairs
      ! are listed at consecutive positions. For such a
      ! complex conjugate pair of the eigenvalues, the
      ! corresponding eigenvectors are also a complex
      ! conjugate pair with the real and imaginary parts
      ! stored column-wise in W at the corresponding
      ! consecutive column indices. See the description of Z.
      ! Also, see the description of DGEEV.
      IF ( info1 > 0 ) THEN
         ! DGEEV failed to compute the eigenvalues and
         ! eigenvectors of the Rayleigh quotient.
         info = 3
         RETURN
      END IF
!
      ! <6> Compute the eigenvectors (if requested) and,
      ! the residuals (if requested).
      !
      IF ( wntvec .OR. wntex ) THEN
      IF ( wntres ) THEN
          IF ( wntref ) THEN
            ! Here, if the refinement is requested, we have
            ! A*U(:,1:K) already computed and stored in Z.
            ! For the residuals, need Y = A * U(:,1;K) * W.
            CALL dgemm( 'N', 'N', m, k, k, one, z, ldz, w, &
                       ldw, zero, y, ldy )               ! BLAS CALL
            ! Y(1:M,1:K) = Z(1:M,1:K) * W(1:K,1:K)       ! INTRINSIC
            ! This frees Z; Y contains A * U(:,1:K) * W.
          ELSE
            ! Compute S = V_k * Sigma_k^(-1) * W, where
            ! V_k * Sigma_k^(-1) is stored in Z
            CALL dgemm( t_or_n, 'N', n, k, k, one, z, ldz, &
                       w, ldw, zero, s, lds)
            ! Then, compute Z = Y * S =
            ! = Y * V_k * Sigma_k^(-1) * W(1:K,1:K) =
            ! = A * U(:,1:K) * W(1:K,1:K)
            CALL dgemm( 'N', 'N', m, k, n, one, y, ldy, s, &
                       lds, zero, z, ldz)
            ! Save a copy of Z into Y and free Z for holding
            ! the Ritz vectors.
            CALL dlacpy( 'A', m, k, z, ldz, y, ldy )
            IF ( wntex ) CALL dlacpy( 'A', m, k, z, ldz, b, ldb )
          END IF
      ELSE IF ( wntex ) THEN
          ! Compute S = V_k * Sigma_k^(-1) * W, where
            ! V_k * Sigma_k^(-1) is stored in Z
            CALL dgemm( t_or_n, 'N', n, k, k, one, z, ldz, &
                       w, ldw, zero, s, lds )
            ! Then, compute Z = Y * S =
            ! = Y * V_k * Sigma_k^(-1) * W(1:K,1:K) =
            ! = A * U(:,1:K) * W(1:K,1:K)
            CALL dgemm( 'N', 'N', m, k, n, one, y, ldy, s, &
                       lds, zero, b, ldb )
            ! The above call replaces the following two calls
            ! that were used in the developing-testing phase.
            ! CALL DGEMM( 'N', 'N', M, K, N, ONE, Y, LDY, S, &
            !           LDS, ZERO, Z, LDZ)
            ! Save a copy of Z into B and free Z for holding
            ! the Ritz vectors.
            ! CALL DLACPY( 'A', M, K, Z, LDZ, B, LDB )
      END IF
!
      ! Compute the real form of the Ritz vectors
      IF ( wntvec ) CALL dgemm( 'N', 'N', m, k, k, one, x, ldx, w, ldw, &
                   zero, z, ldz )                           ! BLAS CALL
      ! Z(1:M,1:K) = MATMUL(X(1:M,1:K), W(1:K,1:K))         ! INTRINSIC
!
      IF ( wntres ) THEN
         i = 1
         DO WHILE ( i <= k )
            IF ( imeig(i) == zero ) THEN
                ! have a real eigenvalue with real eigenvector
                CALL daxpy( m, -reig(i), z(1,i), 1, y(1,i), 1 )       ! BLAS CALL
                ! Y(1:M,i) = Y(1:M,i) - REIG(i) * Z(1:M,i)            ! INTRINSIC
                res(i) = dnrm2( m, y(1,i), 1)                         ! BLAS CALL
                i = i + 1
            ELSE
               ! Have a complex conjugate pair
               ! REIG(i) +- sqrt(-1)*IMEIG(i).
               ! Since all computation is done in real
               ! arithmetic, the formula for the residual
               ! is recast for real representation of the
               ! complex conjugate eigenpair. See the
               ! description of RES.
               ab(1,1) =  reig(i)
               ab(2,1) = -imeig(i)
               ab(1,2) =  imeig(i)
               ab(2,2) =  reig(i)
               CALL dgemm( 'N', 'N', m, 2, 2, -one, z(1,i), &
                           ldz, ab, 2, one, y(1,i), ldy )          ! BLAS CALL
               ! Y(1:M,i:i+1) = Y(1:M,i:i+1) - Z(1:M,i:i+1) * AB   ! INTRINSIC
               res(i)   = dlange( 'F', m, 2, y(1,i), ldy, &
                                  work(n+1) )                      ! LAPACK CALL
               res(i+1) = res(i)
               i = i + 2
            END IF
         END DO
      END IF
      END IF
!
      IF ( whtsvd == 4 ) THEN
          work(n+1) = xscl1
          work(n+2) = xscl2
      END IF
!
!     Successful exit.
      IF ( .NOT. badxy ) THEN
         info = 0
      ELSE
         ! A warning on possible data inconsistency.
         ! This should be a rare event.
         info = 4
      END IF
!............................................................
      RETURN
!     ......
      END SUBROUTINE dgedmd