LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zgedmdq.f90
Go to the documentation of this file.
1
2!
3! =========== DOCUMENTATION ===========
4!
5! Definition:
6! ===========
7!
8! SUBROUTINE ZGEDMDQ( JOBS, JOBZ, JOBR, JOBQ, JOBT, JOBF, &
9! WHTSVD, M, N, F, LDF, X, LDX, Y, &
10! LDY, NRNK, TOL, K, EIGS, &
11! Z, LDZ, RES, B, LDB, V, LDV, &
12! S, LDS, ZWORK, LZWORK, WORK, LWORK, &
13! IWORK, LIWORK, INFO )
14!.....
15! USE iso_fortran_env
16! IMPLICIT NONE
17! INTEGER, PARAMETER :: WP = real64
18!.....
19! Scalar arguments
20! CHARACTER, INTENT(IN) :: JOBS, JOBZ, JOBR, JOBQ, &
21! JOBT, JOBF
22! INTEGER, INTENT(IN) :: WHTSVD, M, N, LDF, LDX, &
23! LDY, NRNK, LDZ, LDB, LDV, &
24! LDS, LZWORK, LWORK, LIWORK
25! INTEGER, INTENT(OUT) :: INFO, K
26! REAL(KIND=WP), INTENT(IN) :: TOL
27! Array arguments
28! COMPLEX(KIND=WP), INTENT(INOUT) :: F(LDF,*)
29! COMPLEX(KIND=WP), INTENT(OUT) :: X(LDX,*), Y(LDY,*), &
30! Z(LDZ,*), B(LDB,*), &
31! V(LDV,*), S(LDS,*)
32! COMPLEX(KIND=WP), INTENT(OUT) :: EIGS(*)
33! COMPLEX(KIND=WP), INTENT(OUT) :: ZWORK(*)
34! REAL(KIND=WP), INTENT(OUT) :: RES(*)
35! REAL(KIND=WP), INTENT(OUT) :: WORK(*)
36! INTEGER, INTENT(OUT) :: IWORK(*)
37!............................................................
39! =============
55!............................................................
57! ================
73!......................................................................
75! ================================
95!......................................................................
97! ================================
103!============================================================
104! Arguments
105! =========
106!
127!.....
149!.....
160!.....
171!.....
182!.....
197!.....
220!.....
226!.....
234!.....
254!.....
260!.....
274!.....
280!.....
292!.....
298!.....
323!.....
330!.....
340!.....
348!.....
365!.....
371!.....
380!.....
397!.....
403!.....
413!.....
419!.....
428!.....
434!.....
448!.....
469!.....
482!.....
494!.....
504!.....
519!.....
540!
541! Authors:
542! ========
543!
545!
547!
548!.............................................................
549!.............................................................
550SUBROUTINE zgedmdq( JOBS, JOBZ, JOBR, JOBQ, JOBT, JOBF, &
551 WHTSVD, M, N, F, LDF, X, LDX, Y, &
552 LDY, NRNK, TOL, K, EIGS, &
553 Z, LDZ, RES, B, LDB, V, LDV, &
554 S, LDS, ZWORK, LZWORK, WORK, LWORK, &
555 IWORK, LIWORK, INFO )
556!
557! -- LAPACK driver routine --
558!
559! -- LAPACK is a software package provided by University of --
560! -- Tennessee, University of California Berkeley, University of --
561! -- Colorado Denver and NAG Ltd.. --
562!
563!.....
564 USE iso_fortran_env
565 IMPLICIT NONE
566 INTEGER, PARAMETER :: WP = real64
567!
568! Scalar arguments
569! ~~~~~~~~~~~~~~~~
570 CHARACTER, INTENT(IN) :: JOBS, JOBZ, JOBR, JOBQ, &
571 JOBT, JOBF
572 INTEGER, INTENT(IN) :: WHTSVD, M, N, LDF, LDX, &
573 LDY, NRNK, LDZ, LDB, LDV, &
574 LDS, LZWORK, LWORK, LIWORK
575 INTEGER, INTENT(OUT) :: INFO, K
576 REAL(KIND=wp), INTENT(IN) :: tol
577!
578! Array arguments
579! ~~~~~~~~~~~~~~~
580 COMPLEX(KIND=WP), INTENT(INOUT) :: F(LDF,*)
581 COMPLEX(KIND=WP), INTENT(OUT) :: X(LDX,*), Y(LDY,*), &
582 Z(LDZ,*), B(LDB,*), &
583 V(LDV,*), S(LDS,*)
584 COMPLEX(KIND=WP), INTENT(OUT) :: EIGS(*)
585 COMPLEX(KIND=WP), INTENT(OUT) :: ZWORK(*)
586 REAL(KIND=wp), INTENT(OUT) :: res(*)
587 REAL(KIND=wp), INTENT(OUT) :: work(*)
588 INTEGER, INTENT(OUT) :: IWORK(*)
589!
590! Parameters
591! ~~~~~~~~~~
592 REAL(KIND=wp), PARAMETER :: one = 1.0_wp
593 REAL(KIND=wp), PARAMETER :: zero = 0.0_wp
594! COMPLEX(KIND=WP), PARAMETER :: ZONE = ( 1.0_WP, 0.0_WP )
595 COMPLEX(KIND=WP), PARAMETER :: ZZERO = ( 0.0_wp, 0.0_wp )
596!
597! Local scalars
598! ~~~~~~~~~~~~~
599 INTEGER :: IMINWR, INFO1, MINMN, MLRWRK, &
600 MLWDMD, MLWGQR, MLWMQR, MLWORK, &
601 MLWQR, OLWDMD, OLWGQR, OLWMQR, &
602 OLWORK, OLWQR
603 LOGICAL :: LQUERY, SCCOLX, SCCOLY, WANTQ, &
604 WNTTRF, WNTRES, WNTVEC, WNTVCF, &
605 WNTVCQ, WNTREF, WNTEX
606 CHARACTER(LEN=1) :: JOBVL
607!
608! External functions (BLAS and LAPACK)
609! ~~~~~~~~~~~~~~~~~
610 LOGICAL LSAME
611 EXTERNAL lsame
612!
613! External subroutines (BLAS and LAPACK)
614! ~~~~~~~~~~~~~~~~~~~~
615 EXTERNAL zgedmd, zgeqrf, zlacpy, zlaset, zungqr, &
616 zunmqr, xerbla
617!
618! Intrinsic functions
619! ~~~~~~~~~~~~~~~~~~~
620 INTRINSIC max, min, int
621!..........................................................
622!
623! Test the input arguments
624 wntres = lsame(jobr,'R')
625 sccolx = lsame(jobs,'S') .OR. lsame( jobs, 'C' )
626 sccoly = lsame(jobs,'Y')
627 wntvec = lsame(jobz,'V')
628 wntvcf = lsame(jobz,'F')
629 wntvcq = lsame(jobz,'Q')
630 wntref = lsame(jobf,'R')
631 wntex = lsame(jobf,'E')
632 wantq = lsame(jobq,'Q')
633 wnttrf = lsame(jobt,'R')
634 minmn = min(m,n)
635 info = 0
636 lquery = ( (lzwork == -1) .OR. (lwork == -1) .OR. (liwork == -1) )
637!
638 IF ( .NOT. (sccolx .OR. sccoly .OR. &
639 lsame(jobs,'N')) ) THEN
640 info = -1
641 ELSE IF ( .NOT. (wntvec .OR. wntvcf .OR. wntvcq &
642 .OR. lsame(jobz,'N')) ) THEN
643 info = -2
644 ELSE IF ( .NOT. (wntres .OR. lsame(jobr,'N')) .OR. &
645 ( wntres .AND. lsame(jobz,'N') ) ) THEN
646 info = -3
647 ELSE IF ( .NOT. (wantq .OR. lsame(jobq,'N')) ) THEN
648 info = -4
649 ELSE IF ( .NOT. ( wnttrf .OR. lsame(jobt,'N') ) ) THEN
650 info = -5
651 ELSE IF ( .NOT. (wntref .OR. wntex .OR. &
652 lsame(jobf,'N') ) ) THEN
653 info = -6
654 ELSE IF ( .NOT. ((whtsvd == 1).OR.(whtsvd == 2).OR. &
655 (whtsvd == 3).OR.(whtsvd == 4)) ) THEN
656 info = -7
657 ELSE IF ( m < 0 ) THEN
658 info = -8
659 ELSE IF ( ( n < 0 ) .OR. ( n > m+1 ) ) THEN
660 info = -9
661 ELSE IF ( ldf < m ) THEN
662 info = -11
663 ELSE IF ( ldx < minmn ) THEN
664 info = -13
665 ELSE IF ( ldy < minmn ) THEN
666 info = -15
667 ELSE IF ( .NOT. (( nrnk == -2).OR.(nrnk == -1).OR. &
668 ((nrnk >= 1).AND.(nrnk <=n ))) ) THEN
669 info = -16
670 ELSE IF ( ( tol < zero ) .OR. ( tol >= one ) ) THEN
671 info = -17
672 ELSE IF ( ldz < m ) THEN
673 info = -21
674 ELSE IF ( (wntref.OR.wntex ).AND.( ldb < minmn ) ) THEN
675 info = -24
676 ELSE IF ( ldv < n-1 ) THEN
677 info = -26
678 ELSE IF ( lds < n-1 ) THEN
679 info = -28
680 END IF
681!
682 IF ( wntvec .OR. wntvcf .OR. wntvcq ) THEN
683 jobvl = 'V'
684 ELSE
685 jobvl = 'N'
686 END IF
687 IF ( info == 0 ) THEN
688 ! Compute the minimal and the optimal workspace
689 ! requirements. Simulate running the code and
690 ! determine minimal and optimal sizes of the
691 ! workspace at any moment of the run.
692 IF ( ( n == 0 ) .OR. ( n == 1 ) ) THEN
693 ! All output except K is void. INFO=1 signals
694 ! the void input. In case of a workspace query,
695 ! the minimal workspace lengths are returned.
696 IF ( lquery ) THEN
697 iwork(1) = 1
698 zwork(1) = 2
699 zwork(2) = 2
700 work(1) = 2
701 work(2) = 2
702 ELSE
703 k = 0
704 END IF
705 info = 1
706 RETURN
707 END IF
708
709 mlrwrk = 2
710 mlwork = 2
711 olwork = 2
712 iminwr = 1
713 mlwqr = max(1,n) ! Minimal workspace length for ZGEQRF.
714 mlwork = max(mlwork,minmn + mlwqr)
715
716 IF ( lquery ) THEN
717 CALL zgeqrf( m, n, f, ldf, zwork, zwork, -1, &
718 info1 )
719 olwqr = int(zwork(1))
720 olwork = max(olwork,minmn + olwqr)
721 END IF
722 CALL zgedmd( jobs, jobvl, jobr, jobf, whtsvd, minmn,&
723 n-1, x, ldx, y, ldy, nrnk, tol, k, &
724 eigs, z, ldz, res, b, ldb, v, ldv, &
725 s, lds, zwork, -1, work, -1, iwork,&
726 -1, info1 )
727 mlwdmd = int(zwork(1))
728 mlwork = max(mlwork, minmn + mlwdmd)
729 mlrwrk = max(mlrwrk, int(work(1)))
730 iminwr = max(iminwr, iwork(1))
731 IF ( lquery ) THEN
732 olwdmd = int(zwork(2))
733 olwork = max(olwork, minmn+olwdmd)
734 END IF
735 IF ( wntvec .OR. wntvcf ) THEN
736 mlwmqr = max(1,n)
737 mlwork = max(mlwork,minmn+mlwmqr)
738 IF ( lquery ) THEN
739 CALL zunmqr( 'L','N', m, n, minmn, f, ldf, &
740 zwork, z, ldz, zwork, -1, info1 )
741 olwmqr = int(zwork(1))
742 olwork = max(olwork,minmn+olwmqr)
743 END IF
744 END IF
745 IF ( wantq ) THEN
746 mlwgqr = max(1,n)
747 mlwork = max(mlwork,minmn+mlwgqr)
748 IF ( lquery ) THEN
749 CALL zungqr( m, minmn, minmn, f, ldf, zwork, &
750 zwork, -1, info1 )
751 olwgqr = int(zwork(1))
752 olwork = max(olwork,minmn+olwgqr)
753 END IF
754 END IF
755 IF ( liwork < iminwr .AND. (.NOT.lquery) ) info = -34
756 IF ( lwork < mlrwrk .AND. (.NOT.lquery) ) info = -32
757 IF ( lzwork < mlwork .AND. (.NOT.lquery) ) info = -30
758 END IF
759 IF( info /= 0 ) THEN
760 CALL xerbla( 'ZGEDMDQ', -info )
761 RETURN
762 ELSE IF ( lquery ) THEN
763! Return minimal and optimal workspace sizes
764 iwork(1) = iminwr
765 zwork(1) = mlwork
766 zwork(2) = olwork
767 work(1) = mlrwrk
768 work(2) = mlrwrk
769 RETURN
770 END IF
771!.....
772! Initial QR factorization that is used to represent the
773! snapshots as elements of lower dimensional subspace.
774! For large scale computation with M >> N, at this place
775! one can use an out of core QRF.
776!
777 CALL zgeqrf( m, n, f, ldf, zwork, &
778 zwork(minmn+1), lzwork-minmn, info1 )
779!
780! Define X and Y as the snapshots representations in the
781! orthogonal basis computed in the QR factorization.
782! X corresponds to the leading N-1 and Y to the trailing
783! N-1 snapshots.
784 CALL zlaset( 'L', minmn, n-1, zzero, zzero, x, ldx )
785 CALL zlacpy( 'U', minmn, n-1, f, ldf, x, ldx )
786 CALL zlacpy( 'A', minmn, n-1, f(1,2), ldf, y, ldy )
787 IF ( m >= 3 ) THEN
788 CALL zlaset( 'L', minmn-2, n-2, zzero, zzero, &
789 y(3,1), ldy )
790 END IF
791!
792! Compute the DMD of the projected snapshot pairs (X,Y)
793 CALL zgedmd( jobs, jobvl, jobr, jobf, whtsvd, minmn, &
794 n-1, x, ldx, y, ldy, nrnk, tol, k, &
795 eigs, z, ldz, res, b, ldb, v, ldv, &
796 s, lds, zwork(minmn+1), lzwork-minmn, &
797 work, lwork, iwork, liwork, info1 )
798 IF ( info1 == 2 .OR. info1 == 3 ) THEN
799 ! Return with error code. See ZGEDMD for details.
800 info = info1
801 RETURN
802 ELSE
803 info = info1
804 END IF
805!
806! The Ritz vectors (Koopman modes) can be explicitly
807! formed or returned in factored form.
808 IF ( wntvec ) THEN
809 ! Compute the eigenvectors explicitly.
810 IF ( m > minmn ) CALL zlaset( 'A', m-minmn, k, zzero, &
811 zzero, z(minmn+1,1), ldz )
812 CALL zunmqr( 'L','N', m, k, minmn, f, ldf, zwork, z, &
813 ldz, zwork(minmn+1), lzwork-minmn, info1 )
814 ELSE IF ( wntvcf ) THEN
815 ! Return the Ritz vectors (eigenvectors) in factored
816 ! form Z*V, where Z contains orthonormal matrix (the
817 ! product of Q from the initial QR factorization and
818 ! the SVD/POD_basis returned by ZGEDMD in X) and the
819 ! second factor (the eigenvectors of the Rayleigh
820 ! quotient) is in the array V, as returned by ZGEDMD.
821 CALL zlacpy( 'A', n, k, x, ldx, z, ldz )
822 IF ( m > n ) CALL zlaset( 'A', m-n, k, zzero, zzero, &
823 z(n+1,1), ldz )
824 CALL zunmqr( 'L','N', m, k, minmn, f, ldf, zwork, z, &
825 ldz, zwork(minmn+1), lzwork-minmn, info1 )
826 END IF
827!
828! Some optional output variables:
829!
830! The upper triangular factor R in the initial QR
831! factorization is optionally returned in the array Y.
832! This is useful if this call to ZGEDMDQ is to be
833! followed by a streaming DMD that is implemented in a
834! QR compressed form.
835 IF ( wnttrf ) THEN ! Return the upper triangular R in Y
836 CALL zlaset( 'A', minmn, n, zzero, zzero, y, ldy )
837 CALL zlacpy( 'U', minmn, n, f, ldf, y, ldy )
838 END IF
839!
840! The orthonormal/unitary factor Q in the initial QR
841! factorization is optionally returned in the array F.
842! Same as with the triangular factor above, this is
843! useful in a streaming DMD.
844 IF ( wantq ) THEN ! Q overwrites F
845 CALL zungqr( m, minmn, minmn, f, ldf, zwork, &
846 zwork(minmn+1), lzwork-minmn, info1 )
847 END IF
848!
849 RETURN
850!
851 END SUBROUTINE zgedmdq
852
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zgedmdq
subroutine zgedmdq(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)
ZGEDMDQ computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.
Definition zgedmdq.f90:556
zgedmd
subroutine zgedmd(jobs, jobz, jobr, jobf, whtsvd, m, n, x, ldx, y, ldy, nrnk, tol, k, eigs, z, ldz, res, b, ldb, w, ldw, s, lds, zwork, lzwork, rwork, lrwork, iwork, liwork, info)
ZGEDMD computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.
Definition zgedmd.f90:498
zgeqrf
subroutine zgeqrf(m, n, a, lda, tau, work, lwork, info)
ZGEQRF
Definition zgeqrf.f:146
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
zlaset
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:106
zungqr
subroutine zungqr(m, n, k, a, lda, tau, work, lwork, info)
ZUNGQR
Definition zungqr.f:128
zunmqr
subroutine zunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMQR
Definition zunmqr.f:167