LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sgedmdq.f90
Go to the documentation of this file.
1
2!
3! =========== DOCUMENTATION ===========
4!
5! Definition:
6! ===========
7!
8! SUBROUTINE SGEDMDQ( 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 = real32
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 sgedmdq( 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 = real32
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 sgemm
637 EXTERNAL sgedmd, sgeqrf, slacpy, slaset, sorgqr, &
638 sormqr, 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. lsame(jobs,'N')) ) THEN
661 info = -1
662 ELSE IF ( .NOT. (wntvec .OR. wntvcf .OR. wntvcq &
663 .OR. lsame(jobz,'N')) ) THEN
664 info = -2
665 ELSE IF ( .NOT. (wntres .OR. lsame(jobr,'N')) .OR. &
666 ( wntres .AND. lsame(jobz,'N') ) ) THEN
667 info = -3
668 ELSE IF ( .NOT. (wantq .OR. lsame(jobq,'N')) ) THEN
669 info = -4
670 ELSE IF ( .NOT. ( wnttrf .OR. lsame(jobt,'N') ) ) THEN
671 info = -5
672 ELSE IF ( .NOT. (wntref .OR. wntex .OR. &
673 lsame(jobf,'N') ) ) THEN
674 info = -6
675 ELSE IF ( .NOT. ((whtsvd == 1).OR.(whtsvd == 2).OR. &
676 (whtsvd == 3).OR.(whtsvd == 4)) ) THEN
677 info = -7
678 ELSE IF ( m < 0 ) THEN
679 info = -8
680 ELSE IF ( ( n < 0 ) .OR. ( n > m+1 ) ) THEN
681 info = -9
682 ELSE IF ( ldf < m ) THEN
683 info = -11
684 ELSE IF ( ldx < minmn ) THEN
685 info = -13
686 ELSE IF ( ldy < minmn ) THEN
687 info = -15
688 ELSE IF ( .NOT. (( nrnk == -2).OR.(nrnk == -1).OR. &
689 ((nrnk >= 1).AND.(nrnk <=n ))) ) THEN
690 info = -16
691 ELSE IF ( ( tol < zero ) .OR. ( tol >= one ) ) THEN
692 info = -17
693 ELSE IF ( ldz < m ) THEN
694 info = -22
695 ELSE IF ( (wntref.OR.wntex ).AND.( ldb < minmn ) ) THEN
696 info = -25
697 ELSE IF ( ldv < n-1 ) THEN
698 info = -27
699 ELSE IF ( lds < n-1 ) THEN
700 info = -29
701 END IF
702!
703 IF ( wntvec .OR. wntvcf ) THEN
704 jobvl = 'V'
705 ELSE
706 jobvl = 'N'
707 END IF
708 IF ( info == 0 ) THEN
709 ! Compute the minimal and the optimal workspace
710 ! requirements. Simulate running the code and
711 ! determine minimal and optimal sizes of the
712 ! workspace at any moment of the run.
713 IF ( ( n == 0 ) .OR. ( n == 1 ) ) THEN
714 ! All output except K is void. INFO=1 signals
715 ! the void input. In case of a workspace query,
716 ! the minimal workspace lengths are returned.
717 IF ( lquery ) THEN
718 iwork(1) = 1
719 work(1) = 2
720 work(2) = 2
721 ELSE
722 k = 0
723 END IF
724 info = 1
725 RETURN
726 END IF
727 mlwqr = max(1,n) ! Minimal workspace length for SGEQRF.
728 mlwork = min(m,n) + mlwqr
729 IF ( lquery ) THEN
730 CALL sgeqrf( m, n, f, ldf, work, rdummy, -1, &
731 info1 )
732 olwqr = int(rdummy(1))
733 olwork = min(m,n) + olwqr
734 END IF
735 CALL sgedmd( jobs, jobvl, jobr, jobf, whtsvd, minmn,&
736 n-1, x, ldx, y, ldy, nrnk, tol, k, &
737 reig, imeig, z, ldz, res, b, ldb, &
738 v, ldv, s, lds, work, -1, iwork, &
739 liwork, info1 )
740 mlwdmd = int(work(1))
741 mlwork = max(mlwork, minmn + mlwdmd)
742 iminwr = iwork(1)
743 IF ( lquery ) THEN
744 olwdmd = int(work(2))
745 olwork = max(olwork, minmn+olwdmd)
746 END IF
747 IF ( wntvec .OR. wntvcf ) THEN
748 mlwmqr = max(1,n)
749 mlwork = max(mlwork,minmn+n-1+mlwmqr)
750 IF ( lquery ) THEN
751 CALL sormqr( 'L','N', m, n, minmn, f, ldf, &
752 work, z, ldz, work, -1, info1 )
753 olwmqr = int(work(1))
754 olwork = max(olwork,minmn+n-1+olwmqr)
755 END IF
756 END IF
757 IF ( wantq ) THEN
758 mlwgqr = n
759 mlwork = max(mlwork,minmn+n-1+mlwgqr)
760 IF ( lquery ) THEN
761 CALL sorgqr( m, minmn, minmn, f, ldf, work, &
762 work, -1, info1 )
763 olwgqr = int(work(1))
764 olwork = max(olwork,minmn+n-1+olwgqr)
765 END IF
766 END IF
767 iminwr = max( 1, iminwr )
768 mlwork = max( 2, mlwork )
769 IF ( lwork < mlwork .AND. (.NOT.lquery) ) info = -31
770 IF ( liwork < iminwr .AND. (.NOT.lquery) ) info = -33
771 END IF
772 IF( info /= 0 ) THEN
773 CALL xerbla( 'SGEDMDQ', -info )
774 RETURN
775 ELSE IF ( lquery ) THEN
776! Return minimal and optimal workspace sizes
777 iwork(1) = iminwr
778 work(1) = mlwork
779 work(2) = olwork
780 RETURN
781 END IF
782!.....
783! Initial QR factorization that is used to represent the
784! snapshots as elements of lower dimensional subspace.
785! For large scale computation with M >>N , at this place
786! one can use an out of core QRF.
787!
788 CALL sgeqrf( m, n, f, ldf, work, &
789 work(minmn+1), lwork-minmn, info1 )
790!
791! Define X and Y as the snapshots representations in the
792! orthogonal basis computed in the QR factorization.
793! X corresponds to the leading N-1 and Y to the trailing
794! N-1 snapshots.
795 CALL slaset( 'L', minmn, n-1, zero, zero, x, ldx )
796 CALL slacpy( 'U', minmn, n-1, f, ldf, x, ldx )
797 CALL slacpy( 'A', minmn, n-1, f(1,2), ldf, y, ldy )
798 IF ( m >= 3 ) THEN
799 CALL slaset( 'L', minmn-2, n-2, zero, zero, &
800 y(3,1), ldy )
801 END IF
802!
803! Compute the DMD of the projected snapshot pairs (X,Y)
804 CALL sgedmd( jobs, jobvl, jobr, jobf, whtsvd, minmn, &
805 n-1, x, ldx, y, ldy, nrnk, tol, k, &
806 reig, imeig, z, ldz, res, b, ldb, v, &
807 ldv, s, lds, work(minmn+1), lwork-minmn, iwork, &
808 liwork, info1 )
809 IF ( info1 == 2 .OR. info1 == 3 ) THEN
810 ! Return with error code.
811 info = info1
812 RETURN
813 ELSE
814 info = info1
815 END IF
816!
817! The Ritz vectors (Koopman modes) can be explicitly
818! formed or returned in factored form.
819 IF ( wntvec ) THEN
820 ! Compute the eigenvectors explicitly.
821 IF ( m > minmn ) CALL slaset( 'A', m-minmn, k, zero, &
822 zero, z(minmn+1,1), ldz )
823 CALL sormqr( 'L','N', m, k, minmn, f, ldf, work, z, &
824 ldz, work(minmn+n), lwork-(minmn+n-1), info1 )
825 ELSE IF ( wntvcf ) THEN
826 ! Return the Ritz vectors (eigenvectors) in factored
827 ! form Z*V, where Z contains orthonormal matrix (the
828 ! product of Q from the initial QR factorization and
829 ! the SVD/POD_basis returned by SGEDMD in X) and the
830 ! second factor (the eigenvectors of the Rayleigh
831 ! quotient) is in the array V, as returned by SGEDMD.
832 CALL slacpy( 'A', n, k, x, ldx, z, ldz )
833 IF ( m > n ) CALL slaset( 'A', m-n, k, zero, zero, &
834 z(n+1,1), ldz )
835 CALL sormqr( 'L','N', m, k, minmn, f, ldf, work, z, &
836 ldz, work(minmn+n), lwork-(minmn+n-1), info1 )
837 END IF
838!
839! Some optional output variables:
840!
841! The upper triangular factor in the initial QR
842! factorization is optionally returned in the array Y.
843! This is useful if this call to SGEDMDQ is to be
844! followed by a streaming DMD that is implemented in a
845! QR compressed form.
846 IF ( wnttrf ) THEN ! Return the upper triangular R in Y
847 CALL slaset( 'A', minmn, n, zero, zero, y, ldy )
848 CALL slacpy( 'U', minmn, n, f, ldf, y, ldy )
849 END IF
850!
851! The orthonormal/orthogonal factor in the initial QR
852! factorization is optionally returned in the array F.
853! Same as with the triangular factor above, this is
854! useful in a streaming DMD.
855 IF ( wantq ) THEN ! Q overwrites F
856 CALL sorgqr( m, minmn, minmn, f, ldf, work, &
857 work(minmn+n), lwork-(minmn+n-1), info1 )
858 END IF
859!
860 RETURN
861!
862 END SUBROUTINE sgedmdq
863
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
sgedmd
subroutine sgedmd(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)
SGEDMD computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.
Definition sgedmd.f90:535
sgedmdq
subroutine sgedmdq(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)
SGEDMDQ computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.
Definition sgedmdq.f90:576
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
sgeqrf
subroutine sgeqrf(m, n, a, lda, tau, work, lwork, info)
SGEQRF
Definition sgeqrf.f:146
slacpy
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
slaset
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:110
sorgqr
subroutine sorgqr(m, n, k, a, lda, tau, work, lwork, info)
SORGQR
Definition sorgqr.f:128
sormqr
subroutine sormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQR
Definition sormqr.f:168