dchkdmd.f90

Go to the documentation of this file.

!     This is a test program for checking the implementations of
!     the implementations of the following subroutines
!
!     DGEDMD  for computation of the
!             Dynamic Mode Decomposition (DMD)
!     DGEDMDQ for computation of a
!             QR factorization based compressed DMD
!
!     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.
!     ========================================================
!     How to run the code (compiler, link info)
!     ========================================================
!     Compile as FORTRAN 90 (or later) and link with BLAS and
!     LAPACK libraries.
!     NOTE: The code is developed and tested on top of the
!     Intel MKL library (versions 2022.0.3 and 2022.2.0),
!     using the Intel Fortran compiler.
!
!     For developers of the C++ implementation
!     ========================================================
!     See the LAPACK++ and Template Numerical Toolkit (TNT)
!
!     Note on a development of the GPU HP implementation
!     ========================================================
!     Work in progress. See CUDA, MAGMA, SLATE.
!     NOTE: The four SVD subroutines used in this code are
!     included as a part of R&D and for the completeness.
!     This was also an opportunity to test those SVD codes.
!     If the scaling option is used all four are essentially
!     equally good. For implementations on HP platforms,
!     one can use whichever SVD is available.
!...  .........................................................
!     NOTE:
!     When using the Intel MKL 2022.0.3 the subroutine xGESVDQ
!     (optionally used in xGEDMD) may cause access violation
!     error for x = S, D, C, Z, but only if called with the
!     work space query. (At least in our Windows 10 MSVS 2019.)
!     The problem can be mitigated by downloading the source
!     code of xGESVDQ from the LAPACK repository and use it
!     localy instead of the one in the MKL. This seems to
!     indicate that the problem is indeed in the MKL.
!     This problem did not appear whith Intel MKL 2022.2.0.
!
!     NOTE:
!     xGESDD seems to have a problem with workspace. In some
!     cases the length of the optimal workspace is returned
!     smaller than the minimal workspace, as specified in the
!     code. As a precaution, all optimal workspaces are
!     set as MAX(minimal, optimal).
!     Latest implementations of complex xGESDD have different
!     length of the real worksapce. We use max value over
!     two versions.
!............................................................
!............................................................
!
      PROGRAM dmd_test
      use iso_fortran_env, only: real64
      IMPLICIT NONE
      integer, parameter :: WP = real64
 
!............................................................
      REAL(KIND=wp), PARAMETER ::  one = 1.0_wp
      REAL(KIND=wp), PARAMETER :: zero = 0.0_wp
!............................................................
      REAL(KIND=wp), ALLOCATABLE, DIMENSION(:,:) ::          &
                     a, ac, eiga, lambda, lambdaq, f, f1, f2,&
                     z, z1, s, au, w, va, x, x0, y, y0, y1
      REAL(KIND=wp), ALLOCATABLE, DIMENSION(:)   ::          &
                     da, dl, dr, reig, reiga, reigq, ieig,   &
                     ieiga, ieigq,  res, res1, resex, singvx,&
                     singvqx, work
      INTEGER      , ALLOCATABLE, DIMENSION(:)   ::   IWORK
      REAL(KIND=wp) :: ab(2,2),   wdummy(2)
      INTEGER       :: IDUMMY(2), ISEED(4), RJOBDATA(8)
      REAL(KIND=wp) :: anorm, cond, condl, condr, dmax, eps, &
                       tol, tol2, svdiff, tmp, tmp_au,       &
                       tmp_fqr, tmp_rez, tmp_rezq,  tmp_zxw, &
                       tmp_ex, xnorm, ynorm
!............................................................
      INTEGER :: K, KQ, LDF, LDS, LDA, LDAU, LDW, LDX, LDY,  &
                 ldz, liwork, lwork, m, n, l, lloop, nrnk
      INTEGER :: i, iJOBREF, iJOBZ, iSCALE, INFO, j, KDIFF,  &
                 nfail, nfail_au, nfail_f_qr, nfail_rez,     &
                 nfail_rezq, nfail_svdiff, nfail_total, nfailq_total, &
                 nfail_z_xv, mode, model, moder, whtsvd
      INTEGER    iNRNK, iWHTSVD, K_TRAJ, LWMINOPT
      CHARACTER(LEN=1) GRADE, JOBREF, JOBZ, PIVTNG, RSIGN,   &
                       scale, resids, wantq, wantr
 
      LOGICAL  TEST_QRDMD
!..... external subroutines (BLAS and LAPACK)
      EXTERNAL daxpy,  dgeev, dgemm, dgemv, dlacpy, dlascl
      EXTERNAL dlarnv, dlatmr
!.....external subroutines DMD package, part 1
!     subroutines under test
      EXTERNAL dgedmd, dgedmdq
 
!..... external functions (BLAS and LAPACK)
      EXTERNAL         dlamch, dlange, dnrm2
      REAL(KIND=wp) :: dlamch, dlange, dnrm2
      EXTERNAL         lsame
      LOGICAL          LSAME
 
      INTRINSIC abs, int, min, max
!............................................................
 
      ! The test is always in pairs : ( DGEDMD and DGEDMDQ )
      ! because the test includes comparing the results (in pairs).
!.....................................................................................
      test_qrdmd = .true. ! This code by default performs tests on DGEDMDQ
                          ! Since the QR factorizations based algorithm is designed for
                          ! single trajectory data, only single trajectory tests will
                          ! be performed with xGEDMDQ.
      wantq = 'Q'
      wantr = 'R'
!.................................................................................
 
      eps = dlamch( 'P' )  ! machine precision DP
 
      ! Global counters of failures of some particular tests
      nfail      = 0
      nfail_rez  = 0
      nfail_rezq = 0
      nfail_z_xv = 0
      nfail_f_qr = 0
      nfail_au   = 0
      kdiff      = 0
      nfail_svdiff = 0
      nfail_total  = 0
      nfailq_total = 0
 
 
      DO lloop = 1, 4
 
      WRITE(*,*) 'L Loop Index = ', lloop
 
      ! Set the dimensions of the problem ...
      WRITE(*,*) 'M = '
      READ(*,*) m
      WRITE(*,*) m
      ! ... and the number of snapshots.
      WRITE(*,*) 'N = '
      READ(*,*) n
      WRITE(*,*) n
 
      ! ... Test the dimensions
      IF ( ( min(m,n) == 0 ) .OR. ( m < n )  ) THEN
          WRITE(*,*) 'Bad dimensions. Required: M >= N > 0.'
          stop
      END IF
!.............
      ! The seed inside the LLOOP so that each pass can be reproduced easily.
 
      iseed(1) = 4
      iseed(2) = 3
      iseed(3) = 2
      iseed(4) = 1
 
      lda  = m
      ldf  = m
      ldx  = max(m,n+1)
      ldy  = max(m,n+1)
      ldw  = n
      ldz  = m
      ldau = max(m,n+1)
      lds  = n
 
      tmp_zxw  = zero
      tmp_au   = zero
      tmp_rez  = zero
      tmp_rezq = zero
      svdiff   = zero
      tmp_ex   = zero
 
      !
      ! Test the subroutines on real data snapshots. All
      ! computation is done in real arithmetic, even when
      ! Koopman eigenvalues and modes are real.
      !
      ! Allocate memory space
      ALLOCATE( a(lda,m) )
      ALLOCATE( ac(lda,m) )
      ALLOCATE( da(m) )
      ALLOCATE( dl(m) )
      ALLOCATE( f(ldf,n+1) )
      ALLOCATE( f1(ldf,n+1) )
      ALLOCATE( f2(ldf,n+1) )
      ALLOCATE( x(ldx,n) )
      ALLOCATE( x0(ldx,n) )
      ALLOCATE( singvx(n) )
      ALLOCATE( singvqx(n) )
      ALLOCATE( y(ldy,n+1) )
      ALLOCATE( y0(ldy,n+1) )
      ALLOCATE( y1(m,n+1) )
      ALLOCATE( z(ldz,n) )
      ALLOCATE( z1(ldz,n) )
      ALLOCATE( res(n)  )
      ALLOCATE( res1(n) )
      ALLOCATE( resex(n) )
      ALLOCATE( reig(n) )
      ALLOCATE( ieig(n) )
      ALLOCATE( reigq(n) )
      ALLOCATE( ieigq(n) )
      ALLOCATE( reiga(m) )
      ALLOCATE( ieiga(m) )
      ALLOCATE( va(lda,m) )
      ALLOCATE( lambda(n,2) )
      ALLOCATE( lambdaq(n,2) )
      ALLOCATE( eiga(m,2) )
      ALLOCATE( w(ldw,n) )
      ALLOCATE( au(ldau,n) )
      ALLOCATE( s(n,n) )
 
      tol  = m*eps
      ! This mimics O(M*N)*EPS bound for accumulated roundoff error.
      ! The factor 10 is somewhat arbitrary.
      tol2 = 10*m*n*eps
 
!.............
 
      DO k_traj = 1, 2
      !  Number of intial conditions in the simulation/trajectories (1 or 2)
 
      cond = 1.0d8
      dmax = 1.0d2
      rsign = 'F'
      grade = 'N'
      model = 6
      condl = 1.0d2
      moder = 6
      condr = 1.0d2
      pivtng = 'N'
 
      ! Loop over all parameter MODE values for ZLATMR (+1,..,+6)
      DO mode = 1, 6
 
      ALLOCATE( iwork(2*m) )
      ALLOCATE(dr(n))
      CALL dlatmr( m, m, 'S', iseed, 'N', da, mode, cond, &
                   dmax, rsign, grade, dl, model,  condl, &
                   dr, moder, condr, pivtng, iwork, m, m, &
                   zero, -one, 'N', a, lda, iwork(m+1), info )
      DEALLOCATE(iwork)
      DEALLOCATE(dr)
 
      lwork = 4*m+1
      ALLOCATE(work(lwork))
      ac  = a
      CALL dgeev( 'N','V', m, ac, m, reiga, ieiga, va, m, &
                  va, m, work, lwork, info ) ! LAPACK CALL
      DEALLOCATE(work)
      tmp = zero
      DO i = 1, m
         eiga(i,1) = reiga(i)
         eiga(i,2) = ieiga(i)
         tmp = max( tmp, sqrt(reiga(i)**2+ieiga(i)**2))
      END DO
 
      ! Scale A to have the desirable spectral radius.
      CALL dlascl( 'G', 0, 0, tmp, one, m, m, a, m, info )
      CALL dlascl( 'G', 0, 0, tmp, one, m, 2, eiga, m, info )
 
      ! Compute the norm of A
      anorm = dlange( 'F', n, n, a, m, wdummy )
 
      IF ( k_traj == 2 ) THEN
          ! generate data with two inital conditions
      CALL dlarnv(2, iseed, m, f1(1,1) )
      f1(1:m,1) = 1.0e-10*f1(1:m,1)
      DO i = 1, n/2
         CALL dgemv( 'N', m, m, one, a, m, f1(1,i), 1, zero, &
              f1(1,i+1), 1 )
      END DO
      x0(1:m,1:n/2) = f1(1:m,1:n/2)
      y0(1:m,1:n/2) = f1(1:m,2:n/2+1)
 
      CALL dlarnv(2, iseed, m, f1(1,1) )
      DO i = 1, n-n/2
         CALL dgemv( 'N', m, m, one, a, m, f1(1,i), 1, zero, &
              f1(1,i+1), 1 )
      END DO
      x0(1:m,n/2+1:n) = f1(1:m,1:n-n/2)
      y0(1:m,n/2+1:n) = f1(1:m,2:n-n/2+1)
      ELSE
      CALL dlarnv(2, iseed, m, f(1,1) )
      DO i = 1, n
         CALL dgemv( 'N', m, m, one, a, m, f(1,i), 1, zero, &
              f(1,i+1), 1 )
      END DO
      x0(1:m,1:n) = f(1:m,1:n)
      y0(1:m,1:n) = f(1:m,2:n+1)
      END IF
 
      xnorm = dlange( 'F', m, n, x0, ldx, wdummy )
      ynorm = dlange( 'F', m, n, y0, ldx, wdummy )
!............................................................
 
      DO ijobz = 1, 4
 
          SELECT CASE ( ijobz )
          CASE(1)
              jobz   = 'V' ! Ritz vectors will be computed
              resids = 'R' ! Residuals will be computed
          CASE(2)
              jobz   = 'V'
              resids = 'N'
          CASE(3)
              jobz   = 'F' ! Ritz vectors in factored form
              resids = 'N'
          CASE(4)
              jobz   = 'N'
              resids = 'N'
          END SELECT
 
      DO ijobref = 1, 3
 
          SELECT CASE ( ijobref )
          CASE(1)
              jobref = 'R' ! Data for refined Ritz vectors
          CASE(2)
              jobref = 'E' ! Exact DMD vectors
          CASE(3)
              jobref = 'N'
          END SELECT
 
      DO iscale = 1, 4
 
          SELECT CASE ( iscale )
          CASE(1)
              scale = 'S' ! X data normalized
          CASE(2)
              scale = 'C' ! X normalized, consist. check
          CASE(3)
              scale = 'Y' ! Y data normalized
          CASE(4)
              scale = 'N'
          END SELECT
 
      DO inrnk = -1, -2, -1
          ! Two truncation strategies. The "-2" case for R&D
          ! purposes only - it uses possibly low accuracy small
          ! singular values, in which case the formulas used in
          ! the DMD are highly sensitive.
          nrnk   = inrnk
 
      DO iwhtsvd = 1, 4
          ! Check all four options to compute the POD basis
          ! via the SVD.
          whtsvd   = iwhtsvd
 
      DO lwminopt = 1, 2
          ! Workspace query for the minimal (1) and for the optimal
          ! (2) workspace lengths determined by workspace query.
 
      x(1:m,1:n) = x0(1:m,1:n)
      y(1:m,1:n) = y0(1:m,1:n)
 
      ! DGEDMD: Workspace query and workspace allocation
      CALL dgedmd( scale, jobz, resids, jobref, whtsvd, m, &
           n, x, ldx, y, ldy, nrnk, tol, k, reig, ieig, z, &
           ldz, res, au, ldau, w, ldw, s, lds, wdummy, -1, &
           idummy, -1, info )
 
      liwork = idummy(1)
      ALLOCATE( iwork(liwork) )
      lwork = int(wdummy(lwminopt))
      ALLOCATE( work(lwork) )
 
      ! DGEDMD test: CALL DGEDMD
      CALL dgedmd( scale, jobz, resids, jobref, whtsvd, m, &
           n, x, ldx, y, ldy, nrnk, tol, k, reig, ieig, z, &
           ldz, res, au, ldau, w, ldw, s, lds, work, lwork,&
           iwork, liwork, info )
 
      singvx(1:n) = work(1:n)
 
      !...... DGEDMD check point
      IF ( lsame(jobz,'V')  ) THEN
          ! Check that Z = X*W, on return from DGEDMD
          ! This checks that the returned aigenvectors in Z are
          ! the product of the SVD'POD basis returned in X
          ! and the eigenvectors of the rayleigh quotient
          ! returned in W
          CALL dgemm( 'N', 'N', m, k, k, one, x, ldx, w, ldw, &
                      zero, z1, ldz )
          tmp = zero
          DO i = 1, k
             CALL daxpy( m, -one, z(1,i), 1, z1(1,i), 1)
             tmp = max(tmp, dnrm2( m, z1(1,i), 1 ) )
          END DO
          tmp_zxw = max(tmp_zxw, tmp )
 
          IF ( tmp_zxw > 10*m*eps ) THEN
              nfail_z_xv = nfail_z_xv + 1
              WRITE(*,*) ':( .................DGEDMD FAILED!', &
                  'Check the code for implementation errors.'
              WRITE(*,*) 'The input parameters were ',&
                 scale, jobz, resids, jobref, whtsvd, &
                 m, n, ldx, ldy, nrnk, tol
          END IF
 
      END IF
 
      !...... DGEDMD check point
      IF ( lsame(jobref,'R') ) THEN
          ! The matrix A*U is returned for computing refined Ritz vectors.
          ! Check that A*U is computed correctly using the formula
          ! A*U = Y * V * inv(SIGMA). This depends on the
          ! accuracy in the computed singular values and vectors of X.
          ! See the paper for an error analysis.
          ! Note that the left singular vectors of the input matrix X
          ! are returned in the array X.
          CALL dgemm( 'N', 'N', m, k, m, one, a, lda, x, ldx, &
                     zero, z1, ldz )
          tmp = zero
          DO i = 1, k
              CALL daxpy( m, -one, au(1,i), 1, z1(1,i), 1)
              tmp = max( tmp, dnrm2( m, z1(1,i),1 ) * &
                       singvx(k)/(anorm*singvx(1)) )
          END DO
          tmp_au = max( tmp_au, tmp )
 
          IF ( tmp > tol2 ) THEN
              nfail_au = nfail_au + 1
              WRITE(*,*) ':( .................DGEDMD FAILED!', &
                  'Check the code for implementation errors.'
              WRITE(*,*) 'The input parameters were ',&
                 scale, jobz, resids, jobref, whtsvd, &
                 m, n, ldx, ldy, nrnk, tol
          END IF
 
      ELSEIF ( lsame(jobref,'E') ) THEN
      ! The unscaled vectors of the Exact DMD are computed.
      ! This option is included for the sake of completeness,
      ! for users who prefer the Exact DMD vectors. The
      ! returned vectors are in the real form, in the same way
      ! as the Ritz vectors. Here we just save the vectors
      ! and test them separately using a Matlab script.
 
       CALL dgemm( 'N', 'N', m, k, m, one, a, lda, au, ldau, zero, y1, m )
       i=1
       DO WHILE ( i <= k )
           IF ( ieig(i) == zero ) THEN
           ! have a real eigenvalue with real eigenvector
           CALL daxpy( m, -reig(i), au(1,i), 1, y1(1,i), 1 )
           resex(i) = dnrm2( m, y1(1,i), 1) / dnrm2(m,au(1,i),1)
           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) = -ieig(i)
           ab(1,2) =  ieig(i)
           ab(2,2) =  reig(i)
           CALL dgemm( 'N', 'N', m, 2, 2, -one, au(1,i), &
                       m, ab, 2, one, y1(1,i), m )
           resex(i)   = dlange( 'F', m, 2, y1(1,i), m, &
                        work )/ dlange( 'F', m, 2, au(1,i), m, &
                        work )
           resex(i+1) = resex(i)
           i = i + 2
           END IF
       END DO
 
      END IF
 
      !...... DGEDMD check point
      IF ( lsame(resids, 'R') ) THEN
          ! Compare the residuals returned by DGEDMD with the
          ! explicitly computed residuals using the matrix A.
          ! Compute explicitly Y1 = A*Z
          CALL dgemm( 'N', 'N', m, k, m, one, a, lda, z, ldz, zero, y1, m )
          ! ... and then A*Z(:,i) - LAMBDA(i)*Z(:,i), using the real forms
          ! of the invariant subspaces that correspond to complex conjugate
          ! pairs of eigencalues. (See the description of Z in DGEDMD,)
          i = 1
          DO WHILE ( i <= k )
              IF ( ieig(i) == zero ) THEN
                  ! have a real eigenvalue with real eigenvector
                  CALL daxpy( m, -reig(i), z(1,i), 1, y1(1,i), 1 )
                  res1(i) = dnrm2( m, y1(1,i), 1)
                  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) = -ieig(i)
                  ab(1,2) =  ieig(i)
                  ab(2,2) =  reig(i)
                  CALL dgemm( 'N', 'N', m, 2, 2, -one, z(1,i), &
                              m, ab, 2, one, y1(1,i), m )
                  res1(i)   = dlange( 'F', m, 2, y1(1,i), m, &
                                     work )
                  res1(i+1) = res1(i)
                  i = i + 2
              END IF
          END DO
          tmp = zero
          DO i = 1, k
              tmp = max( tmp, abs(res(i) - res1(i)) * &
                        singvx(k)/(anorm*singvx(1)) )
          END DO
          tmp_rez = max( tmp_rez, tmp )
 
          IF ( tmp > tol2 ) THEN
              nfail_rez = nfail_rez + 1
              WRITE(*,*) ':( ..................DGEDMD FAILED!', &
                  'Check the code for implementation errors.'
              WRITE(*,*) 'The input parameters were ',&
                 scale, jobz, resids, jobref, whtsvd, &
                 m, n, ldx, ldy, nrnk, tol
          END IF
 
          IF ( lsame(jobref,'E') ) THEN
              tmp = zero
              DO i = 1, k
                  tmp = max( tmp, abs(res1(i) - resex(i))/(res1(i)+resex(i)) )
              END DO
              tmp_ex = max(tmp_ex,tmp)
          END IF
 
      END IF
 
      !..... store the results for inspection
      DO i = 1, k
          lambda(i,1) = reig(i)
          lambda(i,2) = ieig(i)
      END DO
 
      DEALLOCATE(iwork)
      DEALLOCATE(work)
 
      !======================================================================
      !     Now test the DGEDMDQ
      !======================================================================
      IF ( test_qrdmd .AND. (k_traj == 1) ) THEN
          rjobdata(2) = 1
          f1 = f
 
          ! DGEDMDQ test: Workspace query and workspace allocation
          CALL dgedmdq( scale, jobz, resids, wantq, wantr, &
               jobref, whtsvd, m, n+1, f1, ldf, x, ldx, y, &
               ldy, nrnk, tol, kq, reigq, ieigq, z, ldz,   &
               res, au, ldau, w, ldw, s, lds, wdummy,      &
               -1, idummy, -1, info )
          liwork = idummy(1)
          ALLOCATE( iwork(liwork) )
          lwork = int(wdummy(lwminopt))
          ALLOCATE(work(lwork))
          ! DGEDMDQ test: CALL DGEDMDQ
          CALL dgedmdq( scale, jobz, resids, wantq, wantr, &
               jobref, whtsvd, m, n+1, f1, ldf, x, ldx, y, &
               ldy, nrnk, tol, kq, reigq, ieigq, z, ldz,   &
               res, au, ldau, w, ldw, s, lds,              &
               work, lwork, iwork, liwork, info )
 
          singvqx(1:kq) = work(min(m,n+1)+1: min(m,n+1)+kq)
 
          !..... DGEDMDQ check point
          IF ( kq /= k ) THEN
              kdiff = kdiff+1
          END IF
 
          tmp = zero
          DO i = 1, min(k, kq)
              tmp = max(tmp, abs(singvx(i)-singvqx(i)) / &
                                    singvx(1) )
          END DO
          svdiff = max( svdiff, tmp )
          IF ( tmp > m*n*eps ) THEN
              WRITE(*,*) 'FAILED! Something was wrong with the run.'
              nfail_svdiff = nfail_svdiff + 1
              DO j =1, 3
                  write(*,*) j, singvx(j), singvqx(j)
                  read(*,*)
              END DO
          END IF
 
          !..... DGEDMDQ check point
          IF ( lsame(wantq,'Q') .AND. lsame(wantr,'R') ) THEN
              ! Check that the QR factors are computed and returned
              ! as requested. The residual ||F-Q*R||_F / ||F||_F
              ! is compared to M*N*EPS.
              f2 = f
              CALL dgemm( 'N', 'N', m, n+1, min(m,n+1), -one, f1, &
                          ldf, y, ldy, one, f2, ldf )
              tmp_fqr = dlange( 'F', m, n+1, f2, ldf, work ) / &
                    dlange( 'F', m, n+1, f,  ldf, work )
              IF ( tmp_fqr > tol2 ) THEN
                  WRITE(*,*) 'FAILED! Something was wrong with the run.'
                  nfail_f_qr = nfail_f_qr + 1
              END IF
          END IF
 
          !..... DGEDMDQ check point
          IF ( lsame(resids, 'R') ) THEN
              ! Compare the residuals returned by DGEDMDQ with the
              ! explicitly computed residuals using the matrix A.
              ! Compute explicitly Y1 = A*Z
              CALL dgemm( 'N', 'N', m, kq, m, one, a, m, z, m, zero, y1, m )
              ! ... and then A*Z(:,i) - LAMBDA(i)*Z(:,i), using the real forms
              ! of the invariant subspaces that correspond to complex conjugate
              ! pairs of eigencalues. (See the description of Z in DGEDMDQ)
              i = 1
              DO WHILE ( i <= kq )
                  IF ( ieigq(i) == zero ) THEN
                      ! have a real eigenvalue with real eigenvector
                      CALL daxpy( m, -reigq(i), z(1,i), 1, y1(1,i), 1 )
                      ! Y(1:M,i) = Y(1:M,i) - REIG(i)*Z(1:M,i)
                      res1(i) = dnrm2( m, y1(1,i), 1)
                      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) =  reigq(i)
                     ab(2,1) = -ieigq(i)
                     ab(1,2) =  ieigq(i)
                     ab(2,2) =  reigq(i)
                     CALL dgemm( 'N', 'N', m, 2, 2, -one, z(1,i), &
                                 m, ab, 2, one, y1(1,i), m )             ! BLAS CALL
                     ! Y(1:M,i:i+1) = Y(1:M,i:i+1) - Z(1:M,i:i+1) * AB   ! INTRINSIC
                     res1(i)   = dlange( 'F', m, 2, y1(1,i), m, &
                                        work )                           ! LAPACK CALL
                     res1(i+1) = res1(i)
                     i = i + 2
                  END IF
              END DO
              tmp = zero
              DO i = 1, kq
                  tmp = max( tmp, abs(res(i) - res1(i)) * &
                      singvqx(k)/(anorm*singvqx(1)) )
              END DO
              tmp_rezq = max( tmp_rezq, tmp )
              IF ( tmp > tol2 ) THEN
                  nfail_rezq = nfail_rezq + 1
                  WRITE(*,*) '................ DGEDMDQ FAILED!', &
                      'Check the code for implementation errors.'
                  stop
              END IF
 
          END IF
 
          DO i = 1, kq
              lambdaq(i,1) = reigq(i)
              lambdaq(i,2) = ieigq(i)
          END DO
 
          DEALLOCATE(work)
          DEALLOCATE(iwork)
      END IF ! TEST_QRDMD
!======================================================================
 
      END DO ! LWMINOPT
      !write(*,*) 'LWMINOPT loop completed'
      END DO ! WHTSVD LOOP
      !write(*,*) 'WHTSVD loop completed'
      END DO ! NRNK LOOP
      !write(*,*) 'NRNK loop completed'
      END DO ! SCALE LOOP
      !write(*,*) 'SCALE loop completed'
      END DO ! JOBF LOOP
      !write(*,*) 'JOBREF loop completed'
      END DO ! JOBZ LOOP
      !write(*,*) 'JOBZ loop completed'
 
      END DO ! MODE -6:6
      !write(*,*) 'MODE loop completed'
      END DO ! 1 or 2 trajectories
      !write(*,*) 'trajectories  loop completed'
 
      DEALLOCATE(a)
      DEALLOCATE(ac)
      DEALLOCATE(da)
      DEALLOCATE(dl)
      DEALLOCATE(f)
      DEALLOCATE(f1)
      DEALLOCATE(f2)
      DEALLOCATE(x)
      DEALLOCATE(x0)
      DEALLOCATE(singvx)
      DEALLOCATE(singvqx)
      DEALLOCATE(y)
      DEALLOCATE(y0)
      DEALLOCATE(y1)
      DEALLOCATE(z)
      DEALLOCATE(z1)
      DEALLOCATE(res)
      DEALLOCATE(res1)
      DEALLOCATE(resex)
      DEALLOCATE(reig)
      DEALLOCATE(ieig)
      DEALLOCATE(reigq)
      DEALLOCATE(ieigq)
      DEALLOCATE(reiga)
      DEALLOCATE(ieiga)
      DEALLOCATE(va)
      DEALLOCATE(lambda)
      DEALLOCATE(lambdaq)
      DEALLOCATE(eiga)
      DEALLOCATE(w)
      DEALLOCATE(au)
      DEALLOCATE(s)
 
!............................................................
      !     Generate random M-by-M matrix A. Use DLATMR from
      END DO ! LLOOP
 
      WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
      WRITE(*,*) ' Test summary for DGEDMD :'
      WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
      WRITE(*,*)
      IF ( nfail_z_xv == 0 ) THEN
          WRITE(*,*) '>>>> Z - U*V test PASSED.'
      ELSE
          WRITE(*,*) 'Z - U*V test FAILED ', nfail_z_xv, ' time(s)'
          WRITE(*,*) 'Max error ||Z-U*V||_F was ', tmp_zxw
          nfail_total = nfail_total + nfail_z_xv
      END IF
      IF ( nfail_au == 0 ) THEN
          WRITE(*,*) '>>>> A*U test PASSED. '
      ELSE
          WRITE(*,*) 'A*U test FAILED ', nfail_au, ' time(s)'
          WRITE(*,*) 'Max A*U test adjusted error measure was ', tmp_au
          WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', eps
          nfail_total = nfail_total + nfail_au
      END IF
 
      IF ( nfail_rez == 0 ) THEN
          WRITE(*,*) '>>>> Rezidual computation test PASSED.'
      ELSE
          WRITE(*,*) 'Rezidual computation test FAILED ', nfail_rez, 'time(s)'
          WRITE(*,*) 'Max residual computing test adjusted error measure was ', tmp_rez
          WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', eps
          nfail_total = nfail_total + nfail_rez
      END IF
 
      IF ( nfail_total == 0 ) THEN
          WRITE(*,*) '>>>> DGEDMD :: ALL TESTS PASSED.'
      ELSE
          WRITE(*,*) nfail_total, 'FAILURES!'
          WRITE(*,*) '>>>>>>>>>>>>>> DGEDMD :: TESTS FAILED. CHECK THE IMPLEMENTATION.'
      END IF
 
      IF ( test_qrdmd ) THEN
      WRITE(*,*)
      WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
      WRITE(*,*) ' Test summary for DGEDMDQ :'
      WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
      WRITE(*,*)
 
      IF ( nfail_svdiff == 0 ) THEN
          WRITE(*,*) '>>>> DGEDMD and DGEDMDQ computed singular &
              &values test PASSED.'
      ELSE
          WRITE(*,*) 'DGEDMD and DGEDMDQ discrepancies in &
              &the singular values unacceptable ', &
              nfail_svdiff, ' times. Test FAILED.'
          WRITE(*,*) 'The maximal discrepancy in the singular values (relative to the norm) was ', svdiff
          WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', eps
          nfailq_total = nfailq_total + nfail_svdiff
      END IF
 
      IF ( nfail_f_qr == 0 ) THEN
          WRITE(*,*) '>>>> F - Q*R test PASSED.'
      ELSE
          WRITE(*,*) 'F - Q*R test FAILED ', nfail_f_qr, ' time(s)'
          WRITE(*,*) 'The largest relative residual was ', tmp_fqr
          WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', eps
          nfailq_total = nfailq_total + nfail_f_qr
      END IF
 
      IF ( nfail_rezq == 0 ) THEN
          WRITE(*,*) '>>>> Rezidual computation test PASSED.'
      ELSE
          WRITE(*,*) 'Rezidual computation test FAILED ', nfail_rezq, 'time(s)'
          WRITE(*,*) 'Max residual computing test adjusted error measure was ', tmp_rezq
          WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', eps
          nfailq_total = nfailq_total + nfail_rezq
      END IF
 
      IF ( nfailq_total == 0 ) THEN
          WRITE(*,*) '>>>>>>> DGEDMDQ :: ALL TESTS PASSED.'
      ELSE
         WRITE(*,*) nfailq_total, 'FAILURES!'
         WRITE(*,*) '>>>>>>> DGEDMDQ :: TESTS FAILED. CHECK THE IMPLEMENTATION.'
      END IF
 
      END IF
 
      WRITE(*,*)
      WRITE(*,*) 'Test completed.'
      stop
      END