cgedmdq.f90

Go to the documentation of this file.

 
!
!  =========== DOCUMENTATION ===========
!
!  Definition:
!  ===========
!
!     SUBROUTINE CGEDMDQ( JOBS,  JOBZ, JOBR, JOBQ, JOBT, JOBF,   &
!                   WHTSVD,   M, N, F, LDF,  X, LDX,  Y,   &
!                   LDY,   NRNK,  TOL,   K,  EIGS,         &
!                   Z, LDZ, RES,  B,     LDB,   V, LDV,    & 
!                   S, LDS, ZWORK, LZWORK, WORK,  LWORK,   &
!                   IWORK, LIWORK, INFO )
!.....
!     USE                   iso_fortran_env
!     IMPLICIT NONE 
!     INTEGER, PARAMETER :: WP = real32
!.....      
!     Scalar arguments       
!     CHARACTER, INTENT(IN)  :: JOBS, JOBZ, JOBR, JOBQ,    &
!                               JOBT, JOBF
!     INTEGER,   INTENT(IN)  :: WHTSVD, M, N,   LDF, LDX,  &
!                               LDY, NRNK, LDZ, LDB, LDV,  &
!                               LDS, LZWORK,  LWORK, LIWORK
!     INTEGER,   INTENT(OUT) :: INFO,   K      
!     REAL(KIND=WP), INTENT(IN)    ::   TOL     
!     Array arguments      
!     COMPLEX(KIND=WP), INTENT(INOUT) :: F(LDF,*)
!     COMPLEX(KIND=WP), INTENT(OUT)   :: X(LDX,*), Y(LDY,*), &
!                                        Z(LDZ,*), B(LDB,*), &
!                                        V(LDV,*), 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)   :: WORK(*)  
!     INTEGER,       INTENT(OUT)   :: IWORK(*)
!
!............................................................
!     =============
!............................................................
!     ================
!......................................................................
!     ================================
!......................................................................
!     ================================
!......................................................................
!     Arguments
!     =========
!
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!.....
!
!  Authors:
!  ========
!
!
!
!.............................................................
!.............................................................
SUBROUTINE cgedmdq( JOBS,  JOBZ, JOBR, JOBQ, JOBT, JOBF,   &
                    WHTSVD,   M, N, F, LDF,  X, LDX,  Y,   &
                    LDY,   NRNK,  TOL,   K,  EIGS,         &
                    Z, LDZ, RES,  B,     LDB,   V, LDV,    & 
                    S, LDS, ZWORK, LZWORK, 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                   iso_fortran_env
      IMPLICIT NONE 
      INTEGER, PARAMETER :: WP = real32
!
!     Scalar arguments       
!     ~~~~~~~~~~~~~~~~
      CHARACTER, INTENT(IN)  :: JOBS, JOBZ, JOBR, JOBQ,    &
                                JOBT, JOBF
      INTEGER,   INTENT(IN)  :: WHTSVD, M, N,   LDF, LDX,  &
                                LDY, NRNK, LDZ, LDB, LDV,  &
                                LDS, LZWORK,  LWORK, LIWORK
      INTEGER,   INTENT(OUT) :: INFO,   K      
      REAL(KIND=wp), INTENT(IN)    ::   tol     
!
!     Array arguments      
!     ~~~~~~~~~~~~~~~
      COMPLEX(KIND=WP), INTENT(INOUT) :: F(LDF,*)
      COMPLEX(KIND=WP), INTENT(OUT)   :: X(LDX,*), Y(LDY,*), &
                                         Z(LDZ,*), B(LDB,*), &
                                         V(LDV,*), 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)   :: work(*)  
      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      
!     ~~~~~~~~~~~~~
      INTEGER           :: IMINWR, INFO1,  MINMN, MLRWRK,   &
                           MLWDMD, MLWGQR, MLWMQR, MLWORK,  & 
                           MLWQR,  OLWDMD, OLWGQR, OLWMQR,  &
                           OLWORK, OLWQR
      LOGICAL           :: LQUERY, SCCOLX, SCCOLY, WANTQ,  &
                           WNTTRF, WNTRES, WNTVEC, WNTVCF, &
                           WNTVCQ, WNTREF, WNTEX
      CHARACTER(LEN=1)  :: JOBVL
!      
!     External functions (BLAS and LAPACK)
!     ~~~~~~~~~~~~~~~~~
      LOGICAL       LSAME
      EXTERNAL      lsame 
!
!     External subroutines (BLAS and LAPACK)
!     ~~~~~~~~~~~~~~~~~~~~
      EXTERNAL      cgedmd, cgeqrf, clacpy, claset, cungqr, & 
                    cunmqr, xerbla
!
!     Intrinsic functions
!     ~~~~~~~~~~~~~~~~~~~
      INTRINSIC      max, min, int         
!..........................................................  
!
!     Test the input arguments    
      wntres = lsame(jobr,'R')
      sccolx = lsame(jobs,'S') .OR. lsame( jobs, 'C' )
      sccoly = lsame(jobs,'Y')
      wntvec = lsame(jobz,'V')
      wntvcf = lsame(jobz,'F')
      wntvcq = lsame(jobz,'Q')
      wntref = lsame(jobf,'R') 
      wntex  = lsame(jobf,'E')
      wantq  = lsame(jobq,'Q')
      wnttrf = lsame(jobt,'R')     
      minmn  = min(m,n)
      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. wntvcf .OR. wntvcq       &
                              .OR. lsame(jobz,'N')) ) THEN
          info = -2
      ELSE IF ( .NOT. (wntres .OR. lsame(jobr,'N')) .OR.    & 
          ( wntres .AND. lsame(jobz,'N') ) ) THEN
          info = -3
      ELSE IF ( .NOT. (wantq .OR. lsame(jobq,'N')) ) THEN
           info = -4                 
      ELSE IF ( .NOT. ( wnttrf .OR. lsame(jobt,'N') ) )  THEN
          info = -5
       ELSE IF ( .NOT. (wntref .OR. wntex .OR.             & 
                lsame(jobf,'N') ) )                     THEN
          info = -6    
      ELSE IF ( .NOT. ((whtsvd == 1).OR.(whtsvd == 2).OR.   &
                       (whtsvd == 3).OR.(whtsvd == 4)) ) THEN
          info = -7
      ELSE IF ( m < 0 ) THEN
          info = -8
      ELSE IF ( ( n < 0 ) .OR. ( n > m+1 ) ) THEN
          info = -9
      ELSE IF ( ldf < m ) THEN
          info = -11
      ELSE IF ( ldx < minmn ) THEN
          info = -13
      ELSE IF ( ldy < minmn ) THEN
          info = -15
      ELSE IF ( .NOT. (( nrnk == -2).OR.(nrnk == -1).OR.    & 
                       ((nrnk >= 1).AND.(nrnk <=n ))) )  THEN
          info = -16
      ELSE IF ( ( tol < zero ) .OR. ( tol >= one ) ) THEN
          info = -17
      ELSE IF ( ldz < m ) THEN
          info = -21
      ELSE IF ( (wntref.OR.wntex ).AND.( ldb < minmn ) ) THEN
          info = -24
      ELSE IF ( ldv < n-1 ) THEN
          info = -26
      ELSE IF ( lds < n-1 ) THEN
          info = -28
      END IF
!      
      IF ( wntvec .OR. wntvcf .OR. wntvcq ) THEN
          jobvl = 'V'
      ELSE
          jobvl = 'N'
      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 ) .OR. ( n == 1 ) ) THEN
             ! All output except K is void. INFO=1 signals
             ! the void input. In case of a workspace query,
             ! the 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     
         
         mlrwrk = 2
         mlwork = 2
         olwork = 2 
         iminwr = 1
         mlwqr  = max(1,n)  ! Minimal workspace length for CGEQRF.
         mlwork = max(mlwork,minmn + mlwqr) 
 
         IF ( lquery ) THEN 
             CALL cgeqrf( m, n, f, ldf, zwork, zwork, -1, &
                          info1 )
             olwqr  = int(zwork(1))
             olwork = max(olwork,minmn + olwqr)           
         END IF
         CALL cgedmd( jobs, jobvl, jobr, jobf, whtsvd, minmn,& 
                      n-1, x, ldx, y, ldy, nrnk, tol, k,     & 
                      eigs, z, ldz, res,  b, ldb, v, ldv,    & 
                      s, lds, zwork, lzwork, work, -1, iwork,&
                      liwork, info1 )
         mlwdmd = int(zwork(1))
         mlwork = max(mlwork, minmn + mlwdmd)
         mlrwrk = max(mlrwrk, int(work(1)))
         iminwr = max(iminwr, iwork(1))
         IF ( lquery ) THEN 
             olwdmd = int(zwork(2))
             olwork = max(olwork, minmn+olwdmd)
         END IF
         IF ( wntvec .OR. wntvcf ) THEN
            mlwmqr = max(1,n) 
            mlwork = max(mlwork, minmn+mlwmqr)
            IF ( lquery ) THEN
               CALL cunmqr( 'L','N', m, n, minmn, f, ldf,  & 
                            zwork, z, ldz, zwork, -1, info1 )
               olwmqr = int(zwork(1))
               olwork = max(olwork, minmn+olwmqr)
            END IF
         END IF  
         IF ( wantq ) THEN
            mlwgqr = max(1,n)
            mlwork = max(mlwork, minmn+mlwgqr)
            IF ( lquery ) THEN 
                CALL cungqr( m, minmn, minmn, f, ldf, zwork, &
                             zwork, -1, info1 )               
                olwgqr = int(zwork(1))
                olwork = max(olwork, minmn+olwgqr)
            END IF            
         END IF          
         IF ( liwork < iminwr .AND. (.NOT.lquery) ) info = -34
         IF ( lwork  < mlrwrk .AND. (.NOT.lquery) ) info = -32
         IF ( lzwork < mlwork .AND. (.NOT.lquery) ) info = -30
      END IF  
      IF( info /= 0 ) THEN
         CALL xerbla( 'CGEDMDQ', -info )
         RETURN
      ELSE IF ( lquery ) THEN
!     Return minimal and optimal workspace sizes
          iwork(1) = iminwr
          zwork(1) = mlwork
          zwork(2) = olwork
          work(1)  = mlrwrk
          work(2)  = mlrwrk
          RETURN
      END IF   
!.....    
!     Initial QR factorization that is used to represent the
!     snapshots as elements of lower dimensional subspace.
!     For large scale computation with M >>N , at this place 
!     one can use an out of core QRF.
!   
      CALL cgeqrf( m, n, f, ldf, zwork,                & 
                   zwork(minmn+1), lzwork-minmn, info1 )
!      
!     Define X and Y as the snapshots representations in the
!     orthogonal basis computed in the QR factorization.
!     X corresponds to the leading N-1 and Y to the trailing
!     N-1 snapshots.
      CALL claset( 'L', minmn, n-1, zzero,  zzero, x, ldx )
      CALL clacpy( 'U', minmn, n-1, f,      ldf, x, ldx )
      CALL clacpy( 'A', minmn, n-1, f(1,2), ldf, y, ldy )
      IF ( m >= 3 ) THEN
          CALL claset( 'L', minmn-2, n-2, zzero,  zzero, &
                       y(3,1), ldy )  
      END IF
!
!     Compute the DMD of the projected snapshot pairs (X,Y)   
      CALL cgedmd( jobs, jobvl, jobr, jobf, whtsvd, minmn, &
                  n-1,  x, ldx, y, ldy, nrnk,   tol, k,    &
                  eigs, z, ldz, res, b,  ldb,   v, ldv,    &
                  s, lds, zwork(minmn+1), lzwork-minmn,    & 
                  work,   lwork, iwork, liwork, info1 )
      IF ( info1 == 2 .OR. info1 == 3 ) THEN
          ! Return with error code. See CGEDMD for details.
          info = info1
          RETURN
      ELSE
          info = info1
      END IF    
!      
!     The Ritz vectors (Koopman modes) can be explicitly 
!     formed or returned in factored form.
      IF ( wntvec ) THEN
        ! Compute the eigenvectors explicitly.  
        IF ( m > minmn ) CALL claset( 'A', m-minmn, k, zzero, &
                                     zzero, z(minmn+1,1), ldz )
        CALL cunmqr( 'L','N', m, k, minmn, f, ldf, zwork, z,  &
             ldz, zwork(minmn+1), lzwork-minmn, info1 )
      ELSE IF ( wntvcf ) THEN   
        !   Return the Ritz vectors (eigenvectors) in factored
        !   form Z*V, where Z contains orthonormal matrix (the
        !   product of Q from the initial QR factorization and 
        !   the SVD/POD_basis returned by CGEDMD in X) and the 
        !   second factor (the eigenvectors of the Rayleigh 
        !   quotient) is in the array V, as returned by CGEDMD.
        CALL clacpy( 'A', n, k, x, ldx, z, ldz )
        IF ( m > n ) CALL claset( 'A', m-n, k, zzero, zzero, & 
                                 z(n+1,1), ldz )
        CALL cunmqr( 'L','N', m, k, minmn, f, ldf, zwork, z, &
                    ldz, zwork(minmn+1), lzwork-minmn, info1 )
      END IF
!     
!     Some optional output variables:
!
!     The upper triangular factor R in the initial QR 
!     factorization is optionally returned in the array Y.
!     This is useful if this call to CGEDMDQ is to be 
      
!     followed by a streaming DMD that is implemented in a 
!     QR compressed form.
      IF ( wnttrf ) THEN ! Return the upper triangular R in Y 
         CALL claset( 'A', minmn, n, zzero,  zzero, y, ldy )
         CALL clacpy( 'U', minmn, n, f, ldf,        y, ldy )
      END IF    
!
!     The orthonormal/unitary factor Q in the initial QR 
!     factorization is optionally returned in the array F. 
!     Same as with the triangular factor above, this is 
!     useful in a streaming DMD.
      IF ( wantq ) THEN                   ! Q overwrites F 
         CALL cungqr( m, minmn, minmn, f, ldf, zwork,     &
                      zwork(minmn+1), lzwork-minmn, info1 )  
      END IF
!      
      RETURN
!      
      END SUBROUTINE cgedmdq