cchkdmd.f90

Go to the documentation of this file.

!   This is a test program for checking the implementations of
!   the implementations of the following subroutines
!
!   CGEDMD,  for computation of the
!            Dynamic Mode Decomposition (DMD)
!   CGEDMDQ, 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.
!............................................................
 
!............................................................
!............................................................
!
      PROGRAM dmd_test
 
      use iso_fortran_env
      IMPLICIT NONE
      integer, parameter :: WP = real32
!............................................................
      REAL(KIND=wp), PARAMETER ::  one = 1.0_wp
      REAL(KIND=wp), PARAMETER :: zero = 0.0_wp
 
      COMPLEX(KIND=WP), PARAMETER ::  CONE = ( 1.0_wp, 0.0_wp )
      COMPLEX(KIND=WP), PARAMETER :: CZERO = ( 0.0_wp, 0.0_wp )
!............................................................
      REAL(KIND=wp), ALLOCATABLE, DIMENSION(:)   :: res, &
                     res1, resex, singvx, singvqx, work
      INTEGER      , ALLOCATABLE, DIMENSION(:)   ::   IWORK
      REAL(KIND=wp) :: wdummy(2)
      INTEGER       :: IDUMMY(4), ISEED(4)
      REAL(KIND=wp) :: anorm, cond, condl, condr, eps,       &
                       tol, tol2, svdiff, tmp, tmp_au,       &
                       tmp_fqr, tmp_rez, tmp_rezq,  tmp_xw, &
                       tmp_ex
!............................................................
      COMPLEX(KIND=WP) :: CMAX
      INTEGER :: LCWORK
      COMPLEX(KIND=WP), ALLOCATABLE, DIMENSION(:,:) ::  A, AC,  &
                                 au, f, f0, f1, s, w,  &
                                 x, x0, y, y0, y1, z, z1
      COMPLEX(KIND=WP), ALLOCATABLE, DIMENSION(:)   ::  CDA, CDR, &
                                       cdl, ceigs, ceigsa, cwork
      COMPLEX(KIND=WP) ::  CDUMMY(22), CDUM2X2(2,2)
!............................................................
      INTEGER :: K, KQ, LDF, LDS, LDA, LDAU, LDW, LDX, LDY,  &
                 ldz, liwork, lwork, m, n, lloop, nrnk
      INTEGER :: i, iJOBREF, iJOBZ, iSCALE, INFO, j,     &
                 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 :: GRADE, JOBREF, JOBZ, PIVTNG, RSIGN,   &
                   scale, resids, wantq, wantr
      LOGICAL :: TEST_QRDMD
 
!..... external subroutines (BLAS and LAPACK)
      EXTERNAL caxpy, cgeev, cgemm, cgemv, clascl
!.....external subroutines DMD package
!     subroutines under test
      EXTERNAL cgedmd, cgedmdq
!..... external functions (BLAS and LAPACK)
      EXTERNAL         scnrm2, slamch
      REAL(KIND=wp) :: scnrm2, slamch
      EXTERNAL         clange
      REAL(KIND=wp) :: clange
      EXTERNAL icamax
      INTEGER  ICAMAX
      EXTERNAL lsame
      LOGICAL  LSAME
 
      INTRINSIC abs, int, min, max, sign
!............................................................
 
 
      WRITE(*,*) 'COMPLEX CODE TESTING'
 
      ! The test is always in pairs : ( CGEDMD and CGEDMDQ)
      ! because the test includes comparing the results (in pairs).
!.....................................................................................
      ! This code by default performs tests on CGEDMDQ
      ! 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 = slamch( 'P' )  ! machine precision WP
 
      ! 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
      nfail_svdiff = 0
      nfail_total  = 0
      nfailq_total = 0
 
      DO lloop = 1, 4
 
      WRITE(*,*) 'L Loop Index = ', lloop
 
      ! Set the dimensions of the problem ...
      READ(*,*) m
      WRITE(*,*) 'M = ', m
      ! ... and the number of snapshots.
      READ(*,*) n
      WRITE(*,*) 'N = ', 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  = m
      ldy  = m
      ldw  = n
      ldz  = m
      ldau = m
      lds  = n
 
      tmp_xw  = zero
      tmp_au   = zero
      tmp_rez  = zero
      tmp_rezq = zero
      svdiff   = zero
      tmp_ex   = zero
 
      ALLOCATE( a(lda,m) )
      ALLOCATE( ac(lda,m) )
      ALLOCATE( f(ldf,n+1) )
      ALLOCATE( f0(ldf,n+1) )
      ALLOCATE( f1(ldf,n+1) )
      ALLOCATE( x(ldx,n) )
      ALLOCATE( x0(ldx,n) )
      ALLOCATE( y(ldy,n+1) )
      ALLOCATE( y0(ldy,n+1) )
      ALLOCATE( y1(ldy,n+1) )
      ALLOCATE( au(ldau,n) )
      ALLOCATE( w(ldw,n) )
      ALLOCATE( s(lds,n) )
      ALLOCATE( z(ldz,n) )
      ALLOCATE( z1(ldz,n) )
      ALLOCATE( res(n) )
      ALLOCATE( res1(n) )
      ALLOCATE( resex(n) )
      ALLOCATE( ceigs(n) )
      ALLOCATE( singvx(n) )
      ALLOCATE( singvqx(n) )
 
      tol  = 10*m*eps
      tol2 = 10*m*n*eps
 
!.............
 
      DO k_traj = 1, 2
      !  Number of intial conditions in the simulation/trajectories (1 or 2)
 
      cond   = 1.0d4
      cmax   = (1.0d1,1.0d1)
      rsign  = 'F'
      grade  = 'N'
      model  = 6
      condl  = 1.0d1
      moder  = 6
      condr  = 1.0d1
      pivtng = 'N'
      ! Loop over all parameter MODE values for CLATMR (+-1,..,+-6)
 
      DO mode = 1, 6
 
      ALLOCATE( iwork(2*m) )
      ALLOCATE( cda(m) )
      ALLOCATE( cdl(m) )
      ALLOCATE( cdr(m) )
 
      CALL clatmr( m, m, 'N', iseed, 'N', cda, mode, cond, &
                   cmax, rsign, grade, cdl, model,  condl, &
                   cdr, moder, condr, pivtng, iwork, m, m, &
                   zero, -one, 'N', a, lda, iwork(m+1), info )
      DEALLOCATE( cdr )
      DEALLOCATE( cdl )
      DEALLOCATE( cda )
      DEALLOCATE( iwork )
 
      lcwork = max(1,2*m)
      ALLOCATE( ceigsa(m) )
      ALLOCATE( cwork(lcwork) )
      ALLOCATE( work(2*m) )
      ac(1:m,1:m) = a(1:m,1:m)
      CALL cgeev( 'N','N', m, ac, lda, ceigsa, cdum2x2, 2, &
                  cdum2x2, 2, cwork, lcwork, work, info ) ! LAPACK CALL
      DEALLOCATE(work)
      DEALLOCATE(cwork)
 
      tmp = abs(ceigsa(icamax(m, ceigsa, 1))) ! The spectral radius of A
      ! Scale the matrix A to have unit spectral radius.
      CALL clascl( 'G',0, 0, tmp, one, m, m, &
                   a, lda, info )
      CALL clascl( 'G',0, 0, tmp, one, m, 1, &
                   ceigsa, m, info )
      anorm = clange( 'F', m, m, a, lda, wdummy )
 
      IF ( k_traj == 2 ) THEN
          ! generate data as two trajectories
          ! with two inital conditions
          CALL clarnv(2, iseed, m, f(1,1) )
          DO i = 1, n/2
             CALL cgemv( 'N', m, m, cone, a, lda, f(1,i), 1,  &
                  czero, f(1,i+1), 1 )
          END DO
          x0(1:m,1:n/2) = f(1:m,1:n/2)
          y0(1:m,1:n/2) = f(1:m,2:n/2+1)
 
          CALL clarnv(2, iseed, m, f(1,1) )
          DO i = 1, n-n/2
             CALL cgemv( 'N', m, m, cone, a, lda, f(1,i), 1,  &
                  czero, f(1,i+1), 1 )
          END DO
          x0(1:m,n/2+1:n) = f(1:m,1:n-n/2)
          y0(1:m,n/2+1:n) = f(1:m,2:n-n/2+1)
      ELSE
          CALL clarnv(2, iseed, m, f(1,1) )
          DO i = 1, n
             CALL cgemv( 'N', m, m, cone, a, m, f(1,i), 1,  &
                  czero, f(1,i+1), 1 )
          END DO
          f0(1:m,1:n+1) = f(1:m,1:n+1)
          x0(1:m,1:n) = f0(1:m,1:n)
          y0(1:m,1:n) = f0(1:m,2:n+1)
      END IF
 
      DEALLOCATE( ceigsa )
!........................................................................
 
      DO ijobz = 1, 4
 
          SELECT CASE ( ijobz )
          CASE(1)
              jobz   = 'V'
              resids = 'R'
          CASE(2)
              jobz   = 'V'
              resids = 'N'
          CASE(3)
              jobz   = 'F'
              resids = 'N'
          CASE(4)
              jobz   = 'N'
              resids = 'N'
          END SELECT
 
      DO ijobref = 1, 3
 
          SELECT CASE ( ijobref )
          CASE(1)
              jobref = 'R'
          CASE(2)
              jobref = 'E'
          CASE(3)
              jobref = 'N'
          END SELECT
 
      DO iscale = 1, 4
 
          SELECT CASE ( iscale )
          CASE(1)
              scale = 'S'
          CASE(2)
              scale = 'C'
          CASE(3)
              scale = 'Y'
          CASE(4)
              scale = 'N'
          END SELECT
 
      DO inrnk = -1, -2, -1
          nrnk   = inrnk
 
      DO iwhtsvd = 1,  3
         ! 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.
 
      ! CGEDMD is always tested and its results are also used for
      ! comparisons with CGEDMDQ.
 
      x(1:m,1:n) = x0(1:m,1:n)
      y(1:m,1:n) = y0(1:m,1:n)
 
      CALL cgedmd( scale, jobz, resids, jobref, whtsvd,  &
                m,  n, x, ldx, y, ldy, nrnk, tol,  &
                k, ceigs, z, ldz,  res,  &
                au, ldau, w,  ldw,   s, lds,        &
                cdummy, -1, wdummy, -1, idummy, -1, info )
 
      IF ( (info .EQ. 2) .OR. ( info .EQ. 3 ) &
                       .OR. ( info < 0 ) ) THEN
        WRITE(*,*) 'Call to CGEDMD workspace query failed. &
                   &Check the calling sequence and the code.'
        WRITE(*,*) 'The error code is ', info
        WRITE(*,*) 'The input parameters were ',      &
        scale, jobz, resids, jobref, whtsvd,          &
        m, n, ldx, ldy, nrnk, tol, ldz, ldau, ldw, lds
        stop
      ELSE
        !WRITE(*,*) '... done. Workspace length computed.'
      END IF
 
      lcwork = int(cdummy(lwminopt))
      ALLOCATE(cwork(lcwork))
      liwork = idummy(1)
      ALLOCATE(iwork(liwork))
      lwork = int(wdummy(1))
      ALLOCATE(work(lwork))
 
      CALL cgedmd( scale, jobz, resids, jobref, whtsvd,  &
                   m,  n, x, ldx, y, ldy, nrnk, tol,  &
                   k, ceigs, z, ldz,  res,  &
                   au, ldau, w,  ldw,   s, lds,        &
                   cwork, lcwork, work, lwork, iwork, liwork, info )
      IF ( info /= 0 ) THEN
           WRITE(*,*) 'Call to CGEDMD failed. &
           &Check the calling sequence and the code.'
           WRITE(*,*) 'The error code is ', info
           WRITE(*,*) 'The input parameters were ',&
           scale, jobz, resids, jobref, whtsvd, &
           m, n, ldx, ldy, nrnk, tol
           stop
      END IF
      singvx(1:n) = work(1:n)
 
      !...... CGEDMD check point
      IF ( lsame(jobz,'V')  ) THEN
          ! Check that Z = X*W, on return from CGEDMD
          ! This checks that the returned eigenvectors in Z are
          ! the product of the SVD'POD basis returned in X
          ! and the eigenvectors of the Rayleigh quotient
          ! returned in W
          CALL cgemm( 'N', 'N', m, k, k, cone, x, ldx, w, ldw, &
                      czero, z1, ldz )
          tmp = zero
          DO i = 1, k
             CALL caxpy( m, -cone, z(1,i), 1, z1(1,i), 1)
             tmp = max(tmp, scnrm2( m, z1(1,i), 1 ) )
          END DO
          tmp_xw = max(tmp_xw, tmp )
          IF ( tmp_xw <= tol ) THEN
              !WRITE(*,*) ' :) .... OK .........CGEDMD PASSED.'
          ELSE
              nfail_z_xv = nfail_z_xv + 1
              WRITE(*,*) ':( .................CGEDMD 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
      !...... CGEDMD 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 cgemm( 'N', 'N', m, k, m, cone, a, lda, x, ldx, &
                      czero, z1, ldz )
          tmp = zero
          DO i = 1, k
             CALL caxpy( m, -cone, au(1,i), 1, z1(1,i), 1)
             tmp = max( tmp, scnrm2( m, z1(1,i),1 ) * &
                     singvx(k)/(anorm*singvx(1)) )
          END DO
          tmp_au = max( tmp_au, tmp )
          IF ( tmp <= tol2 ) THEN
              !WRITE(*,*) ':) .... OK .........CGEDMD PASSED.'
          ELSE
              nfail_au = nfail_au + 1
              WRITE(*,*) ':( .................CGEDMD FAILED!', &
                  'Check the code for implementation errors.'
              WRITE(*,*) 'The input parameters were ',&
                 scale, jobz, resids, jobref, whtsvd, &
                 m, n, ldx, ldy, nrnk, tol2
          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 cgemm( 'N', 'N', m, k, m, cone, a, lda, au, ldau, czero, y1, ldy )
 
          DO i=1, k
             CALL caxpy( m, -ceigs(i), au(1,i), 1, y1(1,i), 1 )
             resex(i) = scnrm2( m, y1(1,i), 1) / scnrm2(m,au(1,i),1)
          END DO
      END IF
      !...... CGEDMD check point
 
      IF ( lsame(resids, 'R') ) THEN
          ! Compare the residuals returned by CGEDMD with the
          ! explicitly computed residuals using the matrix A.
          ! Compute explicitly Y1 = A*Z
          CALL cgemm( 'N', 'N', m, k, m, cone, a, lda, z, ldz, czero, y1, ldy )
          ! ... 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 CGEDMD,)
 
          DO i=1, k
                ! have a real eigenvalue with real eigenvector
                CALL caxpy( m, -ceigs(i), z(1,i), 1, y1(1,i), 1 )
                res1(i) = scnrm2( m, y1(1,i), 1)
          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
              !WRITE(*,*) ':) .... OK ..........CGEDMD PASSED.'
          ELSE
              nfail_rez = nfail_rez + 1
              WRITE(*,*) ':( ..................CGEDMD 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
 
      DEALLOCATE(cwork)
      DEALLOCATE(work)
      DEALLOCATE(iwork)
 
!.......................................................................................................
 
      IF ( k_traj == 1 ) THEN
 
          f(1:m,1:n+1) = f0(1:m,1:n+1)
          CALL cgedmdq( scale, jobz, resids, wantq, wantr, jobref, &
                    whtsvd, m, n+1, f, ldf,  x, ldx,  y, ldy,  &
                    nrnk,  tol, k, ceigs, z, ldz, res,  au,  &
                    ldau, w, ldw, s, lds, cdummy, -1,   &
                    wdummy,  -1, idummy, -1, info )
 
          lcwork = int(cdummy(lwminopt))
          ALLOCATE(cwork(lcwork))
          liwork = idummy(1)
          ALLOCATE(iwork(liwork))
          lwork = int(wdummy(1))
          ALLOCATE(work(lwork))
 
          CALL cgedmdq( scale, jobz, resids, wantq, wantr, jobref, &
                        whtsvd, m, n+1, f, ldf,  x, ldx,  y, ldy,  &
                        nrnk,  tol, kq, ceigs, z, ldz, res,  au,  &
                        ldau, w, ldw, s, lds, cwork, lcwork,   &
                        work,  lwork, iwork, liwork, info )
          IF ( info /= 0 ) THEN
                 WRITE(*,*) 'Call to CGEDMDQ failed. &
                 &Check the calling sequence and the code.'
                 WRITE(*,*) 'The error code is ', info
                 WRITE(*,*) 'The input parameters were ',&
                 scale, jobz, resids, wantq, wantr, whtsvd, &
                 m, n, ldx, ldy, nrnk, tol
                 stop
          END IF
          singvqx(1:n) =work(1:n)
 
          !..... ZGEDMDQ check point
 
          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 > tol2 ) THEN
               WRITE(*,*) 'FAILED! Something was wrong with the run.'
             nfail_svdiff = nfail_svdiff + 1
          END IF
          !..... CGEDMDQ check point
 
          !..... CGEDMDQ 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.
             f1(1:m,1:n+1) = f0(1:m,1:n+1)
             CALL cgemm( 'N', 'N', m, n+1, min(m,n+1), -cone, f, &
                         ldf, y, ldy, cone, f1, ldf )
             tmp_fqr = clange( 'F', m, n+1, f1, ldf, work ) / &
                   clange( 'F', m, n+1, f0,  ldf, work )
             IF ( tmp_fqr <= tol2 ) THEN
                !WRITE(*,*) ':) CGEDMDQ ........ PASSED.'
             ELSE
                WRITE(*,*) ':( CGEDMDQ ........ FAILED.'
                nfail_f_qr = nfail_f_qr + 1
             END IF
          END IF
          !..... ZGEDMDQ checkpoint
                 !..... ZGEDMDQ checkpoint
          IF ( lsame(resids, 'R') ) THEN
              ! Compare the residuals returned by ZGEDMDQ with the
              ! explicitly computed residuals using the matrix A.
              ! Compute explicitly Y1 = A*Z
              CALL cgemm( 'N', 'N', m, kq, m, cone, a, lda, z, ldz, czero, y1, ldy )
              ! ... 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 ZGEDMDQ)
              DO i = 1, kq
                    ! have a real eigenvalue with real eigenvector
                    CALL caxpy( m, -ceigs(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) = scnrm2( m, y1(1,i), 1)
              END DO
              tmp = zero
              DO i = 1, kq
              tmp = max( tmp, abs(res(i) - res1(i)) * &
                  singvqx(kq)/(anorm*singvqx(1)) )
              END DO
              tmp_rezq = max( tmp_rezq, tmp )
              IF ( tmp <= tol2 ) THEN
                  !WRITE(*,*) '.... OK ........ CGEDMDQ PASSED.'
              ELSE
                  nfail_rezq = nfail_rezq + 1
                  WRITE(*,*) '................ CGEDMDQ FAILED!', &
                      'Check the code for implementation errors.'
              END IF
          END IF
 
          DEALLOCATE(cwork)
          DEALLOCATE(work)
          DEALLOCATE(iwork)
 
      END IF
 
      END DO   ! LWMINOPT
      !write(*,*) 'LWMINOPT loop completed'
      END DO   ! iWHTSVD
      !write(*,*) 'WHTSVD loop completed'
      END DO   ! iNRNK  -2:-1
      !write(*,*) 'NRNK loop completed'
      END DO   ! iSCALE  1:4
      !write(*,*) 'SCALE loop completed'
      END DO
      !write(*,*) 'JOBREF loop completed'
      END DO   ! iJOBZ
      !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( z )
      DEALLOCATE( f )
      DEALLOCATE( f0 )
      DEALLOCATE( f1 )
      DEALLOCATE( x )
      DEALLOCATE( x0 )
      DEALLOCATE( y )
      DEALLOCATE( y0 )
      DEALLOCATE( y1 )
      DEALLOCATE( au )
      DEALLOCATE( w )
      DEALLOCATE( s )
      DEALLOCATE( z1 )
      DEALLOCATE( res )
      DEALLOCATE( res1 )
      DEALLOCATE( resex )
      DEALLOCATE( ceigs )
      DEALLOCATE( singvx )
      DEALLOCATE( singvqx )
 
      END DO ! LLOOP
 
      WRITE(*,*)
      WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
      WRITE(*,*) ' Test summary for CGEDMD :'
      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_xw
          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(*,*) '>>>> CGEDMD :: ALL TESTS PASSED.'
      ELSE
        WRITE(*,*) nfail_total, 'FAILURES!'
        WRITE(*,*) '>>>>>>>>>>>>>> CGEDMD :: TESTS FAILED. CHECK THE IMPLEMENTATION.'
      END IF
 
      WRITE(*,*)
      WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
      WRITE(*,*) ' Test summary for CGEDMDQ :'
      WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
      WRITE(*,*)
 
      IF ( nfail_svdiff == 0 ) THEN
        WRITE(*,*) '>>>> CGEDMD and CGEDMDQ computed singular &
           &values test PASSED.'
      ELSE
        WRITE(*,*) 'ZGEDMD and ZGEDMDQ 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(*,*) '>>>>>>> CGEDMDQ :: ALL TESTS PASSED.'
      ELSE
        WRITE(*,*) nfailq_total, 'FAILURES!'
        WRITE(*,*) '>>>>>>> CGEDMDQ :: TESTS FAILED. CHECK THE IMPLEMENTATION.'
      END IF
 
      WRITE(*,*)
      WRITE(*,*) 'Test completed.'
      stop
      END