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