LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sgedmd.f90
Go to the documentation of this file.
1
2!
3! =========== DOCUMENTATION ===========
4!
5! Definition:
6! ===========
7!
8! SUBROUTINE SGEDMD( JOBS, JOBZ, JOBR, JOBF, WHTSVD, &
9! M, N, X, LDX, Y, LDY, NRNK, TOL, &
10! K, REIG, IMEIG, Z, LDZ, RES, &
11! B, LDB, W, LDW, S, LDS, &
12! WORK, LWORK, IWORK, LIWORK, INFO )
13!.....
14! USE, INTRINSIC :: iso_fortran_env, only: real32
15! IMPLICIT NONE
16! INTEGER, PARAMETER :: WP = real32
17!.....
18! Scalar arguments
19! CHARACTER, INTENT(IN) :: JOBS, JOBZ, JOBR, JOBF
20! INTEGER, INTENT(IN) :: WHTSVD, M, N, LDX, LDY, &
21! NRNK, LDZ, LDB, LDW, LDS, &
22! LWORK, LIWORK
23! INTEGER, INTENT(OUT) :: K, INFO
24! REAL(KIND=WP), INTENT(IN) :: TOL
25! Array arguments
26! REAL(KIND=WP), INTENT(INOUT) :: X(LDX,*), Y(LDY,*)
27! REAL(KIND=WP), INTENT(OUT) :: Z(LDZ,*), B(LDB,*), &
28! W(LDW,*), S(LDS,*)
29! REAL(KIND=WP), INTENT(OUT) :: REIG(*), IMEIG(*), &
30! RES(*)
31! REAL(KIND=WP), INTENT(OUT) :: WORK(*)
32! INTEGER, INTENT(OUT) :: IWORK(*)
33!
34!............................................................
36! =============
51!............................................................
53! ================
69!......................................................................
71! ================================
91!......................................................................
93! ==============================
99!============================================================
100! Arguments
101! =========
102!
121!.....
138!.....
149!.....
163!.....
187!.....
193!.....
200!.....
213!.....
219!.....
231!.....
237!.....
262!.....
269!.....
279!.....
288!.....
307!.....
335!.....
341!.....
361!.....
374!.....
380!.....
393!.....
399!.....
408!.....
414!.....
432!.....
475!.....
485!.....
499!.....
520!
521! Authors:
522! ========
523!
525!
527!
528!.............................................................
529!.............................................................
530 SUBROUTINE sgedmd( JOBS, JOBZ, JOBR, JOBF, WHTSVD, &
531 M, N, X, LDX, Y, LDY, NRNK, TOL, &
532 K, REIG, IMEIG, Z, LDZ, RES, &
533 B, LDB, W, LDW, S, LDS, &
534 WORK, LWORK, IWORK, LIWORK, INFO )
535!
536! -- LAPACK driver routine --
537!
538! -- LAPACK is a software package provided by University of --
539! -- Tennessee, University of California Berkeley, University of --
540! -- Colorado Denver and NAG Ltd.. --
541!
542!.....
543 use, INTRINSIC :: iso_fortran_env, only: real32
544 IMPLICIT NONE
545 INTEGER, PARAMETER :: WP = real32
546!
547! Scalar arguments
548! ~~~~~~~~~~~~~~~~
549 CHARACTER, INTENT(IN) :: JOBS, JOBZ, JOBR, JOBF
550 INTEGER, INTENT(IN) :: WHTSVD, M, N, LDX, LDY, &
551 NRNK, LDZ, LDB, LDW, LDS, &
552 LWORK, LIWORK
553 INTEGER, INTENT(OUT) :: K, INFO
554 REAL(KIND=wp), INTENT(IN) :: tol
555!
556! Array arguments
557! ~~~~~~~~~~~~~~~
558 REAL(KIND=wp), INTENT(INOUT) :: x(ldx,*), y(ldy,*)
559 REAL(KIND=wp), INTENT(OUT) :: z(ldz,*), b(ldb,*), &
560 w(ldw,*), s(lds,*)
561 REAL(KIND=wp), INTENT(OUT) :: reig(*), imeig(*), &
562 res(*)
563 REAL(KIND=wp), INTENT(OUT) :: work(*)
564 INTEGER, INTENT(OUT) :: IWORK(*)
565!
566! Parameters
567! ~~~~~~~~~~
568 REAL(KIND=wp), PARAMETER :: one = 1.0_wp
569 REAL(KIND=wp), PARAMETER :: zero = 0.0_wp
570!
571! Local scalars
572! ~~~~~~~~~~~~~
573 REAL(KIND=wp) :: ofl, rootsc, scale, small, &
574 ssum, xscl1, xscl2
575 INTEGER :: i, j, IMINWR, INFO1, INFO2, &
576 LWRKEV, LWRSDD, LWRSVD, &
577 lwrsvq, mlwork, mwrkev, mwrsdd, &
578 mwrsvd, mwrsvj, mwrsvq, numrnk, &
579 olwork
580 LOGICAL :: BADXY, LQUERY, SCCOLX, SCCOLY, &
581 wntex, wntref, wntres, wntvec
582 CHARACTER :: JOBZL, T_OR_N
583 CHARACTER :: JSVOPT
584!
585! Local arrays
586! ~~~~~~~~~~~~
587 REAL(KIND=wp) :: ab(2,2), rdummy(2), rdummy2(2)
588!
589! External functions (BLAS and LAPACK)
590! ~~~~~~~~~~~~~~~~~
591 REAL(KIND=wp) slange, slamch, snrm2
592 EXTERNAL slange, slamch, snrm2, isamax
593 INTEGER ISAMAX
594 LOGICAL SISNAN, LSAME
595 EXTERNAL sisnan, lsame
596!
597! External subroutines (BLAS and LAPACK)
598! ~~~~~~~~~~~~~~~~~~~~
599 EXTERNAL saxpy, sgemm, sscal
600 EXTERNAL sgeev, sgejsv, sgesdd, sgesvd, sgesvdq, &
602!
603! Intrinsic functions
604! ~~~~~~~~~~~~~~~~~~~
605 INTRINSIC int, float, max, sqrt
606!............................................................
607!
608! Test the input arguments
609!
610 wntres = lsame(jobr,'R')
611 sccolx = lsame(jobs,'S') .OR. lsame(jobs,'C')
612 sccoly = lsame(jobs,'Y')
613 wntvec = lsame(jobz,'V')
614 wntref = lsame(jobf,'R')
615 wntex = lsame(jobf,'E')
616 info = 0
617 lquery = ( ( lwork == -1 ) .OR. ( liwork == -1 ) )
618!
619 IF ( .NOT. (sccolx .OR. sccoly .OR. &
620 lsame(jobs,'N')) ) THEN
621 info = -1
622 ELSE IF ( .NOT. (wntvec .OR. lsame(jobz,'N') &
623 .OR. lsame(jobz,'F')) ) THEN
624 info = -2
625 ELSE IF ( .NOT. (wntres .OR. lsame(jobr,'N')) .OR. &
626 ( wntres .AND. (.NOT.wntvec) ) ) THEN
627 info = -3
628 ELSE IF ( .NOT. (wntref .OR. wntex .OR. &
629 lsame(jobf,'N') ) ) THEN
630 info = -4
631 ELSE IF ( .NOT.((whtsvd == 1) .OR. (whtsvd == 2) .OR. &
632 (whtsvd == 3) .OR. (whtsvd == 4) )) THEN
633 info = -5
634 ELSE IF ( m < 0 ) THEN
635 info = -6
636 ELSE IF ( ( n < 0 ) .OR. ( n > m ) ) THEN
637 info = -7
638 ELSE IF ( ldx < m ) THEN
639 info = -9
640 ELSE IF ( ldy < m ) THEN
641 info = -11
642 ELSE IF ( .NOT. (( nrnk == -2).OR.(nrnk == -1).OR. &
643 ((nrnk >= 1).AND.(nrnk <=n ))) ) THEN
644 info = -12
645 ELSE IF ( ( tol < zero ) .OR. ( tol >= one ) ) THEN
646 info = -13
647 ELSE IF ( ldz < m ) THEN
648 info = -18
649 ELSE IF ( (wntref .OR. wntex ) .AND. ( ldb < m ) ) THEN
650 info = -21
651 ELSE IF ( ldw < n ) THEN
652 info = -23
653 ELSE IF ( lds < n ) THEN
654 info = -25
655 END IF
656!
657 IF ( info == 0 ) THEN
658 ! Compute the minimal and the optimal workspace
659 ! requirements. Simulate running the code and
660 ! determine minimal and optimal sizes of the
661 ! workspace at any moment of the run.
662 IF ( n == 0 ) THEN
663 ! Quick return. All output except K is void.
664 ! INFO=1 signals the void input.
665 ! In case of a workspace query, the default
666 ! minimal workspace lengths are returned.
667 IF ( lquery ) THEN
668 iwork(1) = 1
669 work(1) = 2
670 work(2) = 2
671 ELSE
672 k = 0
673 END IF
674 info = 1
675 RETURN
676 END IF
677 mlwork = max(2,n)
678 olwork = max(2,n)
679 iminwr = 1
680 SELECT CASE ( whtsvd )
681 CASE (1)
682 ! The following is specified as the minimal
683 ! length of WORK in the definition of SGESVD:
684 ! MWRSVD = MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N))
685 mwrsvd = max(1,3*min(m,n)+max(m,n),5*min(m,n))
686 mlwork = max(mlwork,n + mwrsvd)
687 IF ( lquery ) THEN
688 CALL sgesvd( 'O', 'S', m, n, x, ldx, work, &
689 b, ldb, w, ldw, rdummy, -1, info1 )
690 lwrsvd = max( mwrsvd, int( rdummy(1) ) )
691 olwork = max(olwork,n + lwrsvd)
692 END IF
693 CASE (2)
694 ! The following is specified as the minimal
695 ! length of WORK in the definition of SGESDD:
696 ! MWRSDD = 3*MIN(M,N)*MIN(M,N) +
697 ! MAX( MAX(M,N),5*MIN(M,N)*MIN(M,N)+4*MIN(M,N) )
698 ! IMINWR = 8*MIN(M,N)
699 mwrsdd = 3*min(m,n)*min(m,n) + &
700 max( max(m,n),5*min(m,n)*min(m,n)+4*min(m,n) )
701 mlwork = max(mlwork,n + mwrsdd)
702 iminwr = 8*min(m,n)
703 IF ( lquery ) THEN
704 CALL sgesdd( 'O', m, n, x, ldx, work, b, &
705 ldb, w, ldw, rdummy, -1, iwork, info1 )
706 lwrsdd = max( mwrsdd, int( rdummy(1) ) )
707 olwork = max(olwork,n + lwrsdd)
708 END IF
709 CASE (3)
710 !LWQP3 = 3*N+1
711 !LWORQ = MAX(N, 1)
712 !MWRSVD = MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N))
713 !MWRSVQ = N + MAX( LWQP3, MWRSVD, LWORQ )+ MAX(M,2)
714 !MLWORK = N + MWRSVQ
715 !IMINWR = M+N-1
716 CALL sgesvdq( 'H', 'P', 'N', 'R', 'R', m, n, &
717 x, ldx, work, z, ldz, w, ldw, &
718 numrnk, iwork, -1, rdummy, &
719 -1, rdummy2, -1, info1 )
720 iminwr = iwork(1)
721 mwrsvq = int(rdummy(2))
722 mlwork = max(mlwork,n+mwrsvq+int(rdummy2(1)))
723 IF ( lquery ) THEN
724 lwrsvq = int(rdummy(1))
725 olwork = max(olwork,n+lwrsvq+int(rdummy2(1)))
726 END IF
727 CASE (4)
728 jsvopt = 'J'
729 !MWRSVJ = MAX( 7, 2*M+N, 6*N+2*N*N )! for JSVOPT='V'
730 mwrsvj = max( 7, 2*m+n, 4*n+n*n, 2*n+n*n+6 )
731 mlwork = max(mlwork,n+mwrsvj)
732 iminwr = max( 3, m+3*n )
733 IF ( lquery ) THEN
734 olwork = max(olwork,n+mwrsvj)
735 END IF
736 END SELECT
737 IF ( wntvec .OR. wntex .OR. lsame(jobz,'F') ) THEN
738 jobzl = 'V'
739 ELSE
740 jobzl = 'N'
741 END IF
742 ! Workspace calculation to the SGEEV call
743 IF ( lsame(jobzl,'V') ) THEN
744 mwrkev = max( 1, 4*n )
745 ELSE
746 mwrkev = max( 1, 3*n )
747 END IF
748 mlwork = max(mlwork,n+mwrkev)
749 IF ( lquery ) THEN
750 CALL sgeev( 'N', jobzl, n, s, lds, reig, &
751 imeig, w, ldw, w, ldw, rdummy, -1, info1 )
752 lwrkev = max( mwrkev, int(rdummy(1)) )
753 olwork = max( olwork, n+lwrkev )
754 END IF
755!
756 IF ( liwork < iminwr .AND. (.NOT.lquery) ) info = -29
757 IF ( lwork < mlwork .AND. (.NOT.lquery) ) info = -27
758 END IF
759!
760 IF( info /= 0 ) THEN
761 CALL xerbla( 'SGEDMD', -info )
762 RETURN
763 ELSE IF ( lquery ) THEN
764! Return minimal and optimal workspace sizes
765 iwork(1) = iminwr
766 work(1) = real(mlwork)
767 work(2) = real(olwork)
768 RETURN
769 END IF
770!............................................................
771!
772 ofl = slamch('O')
773 small = slamch('S')
774 badxy = .false.
775!
776! <1> Optional scaling of the snapshots (columns of X, Y)
777! ==========================================================
778 IF ( sccolx ) THEN
779 ! The columns of X will be normalized.
780 ! To prevent overflows, the column norms of X are
781 ! carefully computed using SLASSQ.
782 k = 0
783 DO i = 1, n
784 !WORK(i) = DNRM2( M, X(1,i), 1 )
785 ssum = one
786 scale = zero
787 CALL slassq( m, x(1,i), 1, scale, ssum )
788 IF ( sisnan(scale) .OR. sisnan(ssum) ) THEN
789 k = 0
790 info = -8
791 CALL xerbla('SGEDMD',-info)
792 END IF
793 IF ( (scale /= zero) .AND. (ssum /= zero) ) THEN
794 rootsc = sqrt(ssum)
795 IF ( scale .GE. (ofl / rootsc) ) THEN
796! Norm of X(:,i) overflows. First, X(:,i)
797! is scaled by
798! ( ONE / ROOTSC ) / SCALE = 1/||X(:,i)||_2.
799! Next, the norm of X(:,i) is stored without
800! overflow as WORK(i) = - SCALE * (ROOTSC/M),
801! the minus sign indicating the 1/M factor.
802! Scaling is performed without overflow, and
803! underflow may occur in the smallest entries
804! of X(:,i). The relative backward and forward
805! errors are small in the ell_2 norm.
806 CALL slascl( 'G', 0, 0, scale, one/rootsc, &
807 m, 1, x(1,i), m, info2 )
808 work(i) = - scale * ( rootsc / float(m) )
809 ELSE
810! X(:,i) will be scaled to unit 2-norm
811 work(i) = scale * rootsc
812 CALL slascl( 'G',0, 0, work(i), one, m, 1, &
813 x(1,i), m, info2 ) ! LAPACK CALL
814! X(1:M,i) = (ONE/WORK(i)) * X(1:M,i) ! INTRINSIC
815 END IF
816 ELSE
817 work(i) = zero
818 k = k + 1
819 END IF
820 END DO
821 IF ( k == n ) THEN
822 ! All columns of X are zero. Return error code -8.
823 ! (the 8th input variable had an illegal value)
824 k = 0
825 info = -8
826 CALL xerbla('SGEDMD',-info)
827 RETURN
828 END IF
829 DO i = 1, n
830! Now, apply the same scaling to the columns of Y.
831 IF ( work(i) > zero ) THEN
832 CALL sscal( m, one/work(i), y(1,i), 1 ) ! BLAS CALL
833! Y(1:M,i) = (ONE/WORK(i)) * Y(1:M,i) ! INTRINSIC
834 ELSE IF ( work(i) < zero ) THEN
835 CALL slascl( 'G', 0, 0, -work(i), &
836 one/float(m), m, 1, y(1,i), m, info2 ) ! LAPACK CALL
837 ELSE IF ( y(isamax(m, y(1,i),1),i ) &
838 /= zero ) THEN
839! X(:,i) is zero vector. For consistency,
840! Y(:,i) should also be zero. If Y(:,i) is not
841! zero, then the data might be inconsistent or
842! corrupted. If JOBS == 'C', Y(:,i) is set to
843! zero and a warning flag is raised.
844! The computation continues but the
845! situation will be reported in the output.
846 badxy = .true.
847 IF ( lsame(jobs,'C')) &
848 CALL sscal( m, zero, y(1,i), 1 ) ! BLAS CALL
849 END IF
850 END DO
851 END IF
852 !
853 IF ( sccoly ) THEN
854 ! The columns of Y will be normalized.
855 ! To prevent overflows, the column norms of Y are
856 ! carefully computed using SLASSQ.
857 DO i = 1, n
858 !WORK(i) = DNRM2( M, Y(1,i), 1 )
859 ssum = one
860 scale = zero
861 CALL slassq( m, y(1,i), 1, scale, ssum )
862 IF ( sisnan(scale) .OR. sisnan(ssum) ) THEN
863 k = 0
864 info = -10
865 CALL xerbla('SGEDMD',-info)
866 END IF
867 IF ( scale /= zero .AND. (ssum /= zero) ) THEN
868 rootsc = sqrt(ssum)
869 IF ( scale .GE. (ofl / rootsc) ) THEN
870! Norm of Y(:,i) overflows. First, Y(:,i)
871! is scaled by
872! ( ONE / ROOTSC ) / SCALE = 1/||Y(:,i)||_2.
873! Next, the norm of Y(:,i) is stored without
874! overflow as WORK(i) = - SCALE * (ROOTSC/M),
875! the minus sign indicating the 1/M factor.
876! Scaling is performed without overflow, and
877! underflow may occur in the smallest entries
878! of Y(:,i). The relative backward and forward
879! errors are small in the ell_2 norm.
880 CALL slascl( 'G', 0, 0, scale, one/rootsc, &
881 m, 1, y(1,i), m, info2 )
882 work(i) = - scale * ( rootsc / float(m) )
883 ELSE
884! X(:,i) will be scaled to unit 2-norm
885 work(i) = scale * rootsc
886 CALL slascl( 'G',0, 0, work(i), one, m, 1, &
887 y(1,i), m, info2 ) ! LAPACK CALL
888! Y(1:M,i) = (ONE/WORK(i)) * Y(1:M,i) ! INTRINSIC
889 END IF
890 ELSE
891 work(i) = zero
892 END IF
893 END DO
894 DO i = 1, n
895! Now, apply the same scaling to the columns of X.
896 IF ( work(i) > zero ) THEN
897 CALL sscal( m, one/work(i), x(1,i), 1 ) ! BLAS CALL
898! X(1:M,i) = (ONE/WORK(i)) * X(1:M,i) ! INTRINSIC
899 ELSE IF ( work(i) < zero ) THEN
900 CALL slascl( 'G', 0, 0, -work(i), &
901 one/float(m), m, 1, x(1,i), m, info2 ) ! LAPACK CALL
902 ELSE IF ( x(isamax(m, x(1,i),1),i ) &
903 /= zero ) THEN
904! Y(:,i) is zero vector. If X(:,i) is not
905! zero, then a warning flag is raised.
906! The computation continues but the
907! situation will be reported in the output.
908 badxy = .true.
909 END IF
910 END DO
911 END IF
912!
913! <2> SVD of the data snapshot matrix X.
914! =====================================
915! The left singular vectors are stored in the array X.
916! The right singular vectors are in the array W.
917! The array W will later on contain the eigenvectors
918! of a Rayleigh quotient.
919 numrnk = n
920 SELECT CASE ( whtsvd )
921 CASE (1)
922 CALL sgesvd( 'O', 'S', m, n, x, ldx, work, b, &
923 ldb, w, ldw, work(n+1), lwork-n, info1 ) ! LAPACK CALL
924 t_or_n = 'T'
925 CASE (2)
926 CALL sgesdd( 'O', m, n, x, ldx, work, b, ldb, w, &
927 ldw, work(n+1), lwork-n, iwork, info1 ) ! LAPACK CALL
928 t_or_n = 'T'
929 CASE (3)
930 CALL sgesvdq( 'H', 'P', 'N', 'R', 'R', m, n, &
931 x, ldx, work, z, ldz, w, ldw, &
932 numrnk, iwork, liwork, work(n+max(2,m)+1),&
933 lwork-n-max(2,m), work(n+1), max(2,m), info1) ! LAPACK CALL
934 CALL slacpy( 'A', m, numrnk, z, ldz, x, ldx ) ! LAPACK CALL
935 t_or_n = 'T'
936 CASE (4)
937 CALL sgejsv( 'F', 'U', jsvopt, 'N', 'N', 'P', m, &
938 n, x, ldx, work, z, ldz, w, ldw, &
939 work(n+1), lwork-n, iwork, info1 ) ! LAPACK CALL
940 CALL slacpy( 'A', m, n, z, ldz, x, ldx ) ! LAPACK CALL
941 t_or_n = 'N'
942 xscl1 = work(n+1)
943 xscl2 = work(n+2)
944 IF ( xscl1 /= xscl2 ) THEN
945 ! This is an exceptional situation. If the
946 ! data matrices are not scaled and the
947 ! largest singular value of X overflows.
948 ! In that case SGEJSV can return the SVD
949 ! in scaled form. The scaling factor can be used
950 ! to rescale the data (X and Y).
951 CALL slascl( 'G', 0, 0, xscl1, xscl2, m, n, y, ldy, info2 )
952 END IF
953 END SELECT
954!
955 IF ( info1 > 0 ) THEN
956 ! The SVD selected subroutine did not converge.
957 ! Return with an error code.
958 info = 2
959 RETURN
960 END IF
961!
962 IF ( work(1) == zero ) THEN
963 ! The largest computed singular value of (scaled)
964 ! X is zero. Return error code -8
965 ! (the 8th input variable had an illegal value).
966 k = 0
967 info = -8
968 CALL xerbla('SGEDMD',-info)
969 RETURN
970 END IF
971!
972
973 ! snapshots matrix X. This depends on the
974 ! parameters NRNK and TOL.
975
976 SELECT CASE ( nrnk )
977 CASE ( -1 )
978 k = 1
979 DO i = 2, numrnk
980 IF ( ( work(i) <= work(1)*tol ) .OR. &
981 ( work(i) <= small ) ) EXIT
982 k = k + 1
983 END DO
984 CASE ( -2 )
985 k = 1
986 DO i = 1, numrnk-1
987 IF ( ( work(i+1) <= work(i)*tol ) .OR. &
988 ( work(i) <= small ) ) EXIT
989 k = k + 1
990 END DO
991 CASE DEFAULT
992 k = 1
993 DO i = 2, nrnk
994 IF ( work(i) <= small ) EXIT
995 k = k + 1
996 END DO
997 END SELECT
998 ! Now, U = X(1:M,1:K) is the SVD/POD basis for the
999 ! snapshot data in the input matrix X.
1000
1001
1002 ! Depending on the requested outputs, the computation
1003 ! is organized to compute additional auxiliary
1004 ! matrices (for the residuals and refinements).
1005 !
1006 ! In all formulas below, we need V_k*Sigma_k^(-1)
1007 ! where either V_k is in W(1:N,1:K), or V_k^T is in
1008 ! W(1:K,1:N). Here Sigma_k=diag(WORK(1:K)).
1009 IF ( lsame(t_or_n, 'N') ) THEN
1010 DO i = 1, k
1011 CALL sscal( n, one/work(i), w(1,i), 1 ) ! BLAS CALL
1012 ! W(1:N,i) = (ONE/WORK(i)) * W(1:N,i) ! INTRINSIC
1013 END DO
1014 ELSE
1015 ! This non-unit stride access is due to the fact
1016 ! that SGESVD, SGESVDQ and SGESDD return the
1017 ! transposed matrix of the right singular vectors.
1018 !DO i = 1, K
1019 ! CALL SSCAL( N, ONE/WORK(i), W(i,1), LDW ) ! BLAS CALL
1020 ! ! W(i,1:N) = (ONE/WORK(i)) * W(i,1:N) ! INTRINSIC
1021 !END DO
1022 DO i = 1, k
1023 work(n+i) = one/work(i)
1024 END DO
1025 DO j = 1, n
1026 DO i = 1, k
1027 w(i,j) = (work(n+i))*w(i,j)
1028 END DO
1029 END DO
1030 END IF
1031!
1032 IF ( wntref ) THEN
1033 !
1034 ! Need A*U(:,1:K)=Y*V_k*inv(diag(WORK(1:K)))
1035 ! for computing the refined Ritz vectors
1036 ! (optionally, outside SGEDMD).
1037 CALL sgemm( 'N', t_or_n, m, k, n, one, y, ldy, w, &
1038 ldw, zero, z, ldz ) ! BLAS CALL
1039 ! Z(1:M,1:K)=MATMUL(Y(1:M,1:N),TRANSPOSE(W(1:K,1:N))) ! INTRINSIC, for T_OR_N=='T'
1040 ! Z(1:M,1:K)=MATMUL(Y(1:M,1:N),W(1:N,1:K)) ! INTRINSIC, for T_OR_N=='N'
1041 !
1042 ! At this point Z contains
1043 ! A * U(:,1:K) = Y * V_k * Sigma_k^(-1), and
1044 ! this is needed for computing the residuals.
1045 ! This matrix is returned in the array B and
1046 ! it can be used to compute refined Ritz vectors.
1047 CALL slacpy( 'A', m, k, z, ldz, b, ldb ) ! BLAS CALL
1048 ! B(1:M,1:K) = Z(1:M,1:K) ! INTRINSIC
1049
1050 CALL sgemm( 'T', 'N', k, k, m, one, x, ldx, z, &
1051 ldz, zero, s, lds ) ! BLAS CALL
1052 ! S(1:K,1:K) = MATMUL(TANSPOSE(X(1:M,1:K)),Z(1:M,1:K)) ! INTRINSIC
1053 ! At this point S = U^T * A * U is the Rayleigh quotient.
1054 ELSE
1055 ! A * U(:,1:K) is not explicitly needed and the
1056 ! computation is organized differently. The Rayleigh
1057 ! quotient is computed more efficiently.
1058 CALL sgemm( 'T', 'N', k, n, m, one, x, ldx, y, ldy, &
1059 zero, z, ldz ) ! BLAS CALL
1060 ! Z(1:K,1:N) = MATMUL( TRANSPOSE(X(1:M,1:K)), Y(1:M,1:N) ) ! INTRINSIC
1061 ! In the two SGEMM calls here, can use K for LDZ
1062 CALL sgemm( 'N', t_or_n, k, k, n, one, z, ldz, w, &
1063 ldw, zero, s, lds ) ! BLAS CALL
1064 ! S(1:K,1:K) = MATMUL(Z(1:K,1:N),TRANSPOSE(W(1:K,1:N))) ! INTRINSIC, for T_OR_N=='T'
1065 ! S(1:K,1:K) = MATMUL(Z(1:K,1:N),(W(1:N,1:K))) ! INTRINSIC, for T_OR_N=='N'
1066 ! At this point S = U^T * A * U is the Rayleigh quotient.
1067 ! If the residuals are requested, save scaled V_k into Z.
1068 ! Recall that V_k or V_k^T is stored in W.
1069 IF ( wntres .OR. wntex ) THEN
1070 IF ( lsame(t_or_n, 'N') ) THEN
1071 CALL slacpy( 'A', n, k, w, ldw, z, ldz )
1072 ELSE
1073 CALL slacpy( 'A', k, n, w, ldw, z, ldz )
1074 END IF
1075 END IF
1076 END IF
1077!
1078
1079 ! right eigenvectors of the Rayleigh quotient.
1080 !
1081 CALL sgeev( 'N', jobzl, k, s, lds, reig, imeig, w, &
1082 ldw, w, ldw, work(n+1), lwork-n, info1 ) ! LAPACK CALL
1083 !
1084 ! W(1:K,1:K) contains the eigenvectors of the Rayleigh
1085 ! quotient. Even in the case of complex spectrum, all
1086 ! computation is done in real arithmetic. REIG and
1087 ! IMEIG are the real and the imaginary parts of the
1088 ! eigenvalues, so that the spectrum is given as
1089 ! REIG(:) + sqrt(-1)*IMEIG(:). Complex conjugate pairs
1090 ! are listed at consecutive positions. For such a
1091 ! complex conjugate pair of the eigenvalues, the
1092 ! corresponding eigenvectors are also a complex
1093 ! conjugate pair with the real and imaginary parts
1094 ! stored column-wise in W at the corresponding
1095 ! consecutive column indices. See the description of Z.
1096 ! Also, see the description of SGEEV.
1097 IF ( info1 > 0 ) THEN
1098 ! SGEEV failed to compute the eigenvalues and
1099 ! eigenvectors of the Rayleigh quotient.
1100 info = 3
1101 RETURN
1102 END IF
1103!
1104 ! <6> Compute the eigenvectors (if requested) and,
1105 ! the residuals (if requested).
1106 !
1107 IF ( wntvec .OR. wntex ) THEN
1108 IF ( wntres ) THEN
1109 IF ( wntref ) THEN
1110 ! Here, if the refinement is requested, we have
1111 ! A*U(:,1:K) already computed and stored in Z.
1112 ! For the residuals, need Y = A * U(:,1;K) * W.
1113 CALL sgemm( 'N', 'N', m, k, k, one, z, ldz, w, &
1114 ldw, zero, y, ldy ) ! BLAS CALL
1115 ! Y(1:M,1:K) = Z(1:M,1:K) * W(1:K,1:K) ! INTRINSIC
1116 ! This frees Z; Y contains A * U(:,1:K) * W.
1117 ELSE
1118 ! Compute S = V_k * Sigma_k^(-1) * W, where
1119 ! V_k * Sigma_k^(-1) is stored in Z
1120 CALL sgemm( t_or_n, 'N', n, k, k, one, z, ldz, &
1121 w, ldw, zero, s, lds )
1122 ! Then, compute Z = Y * S =
1123 ! = Y * V_k * Sigma_k^(-1) * W(1:K,1:K) =
1124 ! = A * U(:,1:K) * W(1:K,1:K)
1125 CALL sgemm( 'N', 'N', m, k, n, one, y, ldy, s, &
1126 lds, zero, z, ldz )
1127 ! Save a copy of Z into Y and free Z for holding
1128 ! the Ritz vectors.
1129 CALL slacpy( 'A', m, k, z, ldz, y, ldy )
1130 IF ( wntex ) CALL slacpy( 'A', m, k, z, ldz, b, ldb )
1131 END IF
1132 ELSE IF ( wntex ) THEN
1133 ! Compute S = V_k * Sigma_k^(-1) * W, where
1134 ! V_k * Sigma_k^(-1) is stored in Z
1135 CALL sgemm( t_or_n, 'N', n, k, k, one, z, ldz, &
1136 w, ldw, zero, s, lds )
1137 ! Then, compute Z = Y * S =
1138 ! = Y * V_k * Sigma_k^(-1) * W(1:K,1:K) =
1139 ! = A * U(:,1:K) * W(1:K,1:K)
1140 CALL sgemm( 'N', 'N', m, k, n, one, y, ldy, s, &
1141 lds, zero, b, ldb )
1142 ! The above call replaces the following two calls
1143 ! that were used in the developing-testing phase.
1144 ! CALL SGEMM( 'N', 'N', M, K, N, ONE, Y, LDY, S, &
1145 ! LDS, ZERO, Z, LDZ)
1146 ! Save a copy of Z into B and free Z for holding
1147 ! the Ritz vectors.
1148 ! CALL SLACPY( 'A', M, K, Z, LDZ, B, LDB )
1149 END IF
1150!
1151 ! Compute the real form of the Ritz vectors
1152 IF ( wntvec ) CALL sgemm( 'N', 'N', m, k, k, one, x, ldx, w, ldw, &
1153 zero, z, ldz ) ! BLAS CALL
1154 ! Z(1:M,1:K) = MATMUL(X(1:M,1:K), W(1:K,1:K)) ! INTRINSIC
1155!
1156 IF ( wntres ) THEN
1157 i = 1
1158 DO WHILE ( i <= k )
1159 IF ( imeig(i) == zero ) THEN
1160 ! have a real eigenvalue with real eigenvector
1161 CALL saxpy( m, -reig(i), z(1,i), 1, y(1,i), 1 ) ! BLAS CALL
1162 ! Y(1:M,i) = Y(1:M,i) - REIG(i) * Z(1:M,i) ! INTRINSIC
1163 res(i) = snrm2( m, y(1,i), 1 ) ! BLAS CALL
1164 i = i + 1
1165 ELSE
1166 ! Have a complex conjugate pair
1167 ! REIG(i) +- sqrt(-1)*IMEIG(i).
1168 ! Since all computation is done in real
1169 ! arithmetic, the formula for the residual
1170 ! is recast for real representation of the
1171 ! complex conjugate eigenpair. See the
1172 ! description of RES.
1173 ab(1,1) = reig(i)
1174 ab(2,1) = -imeig(i)
1175 ab(1,2) = imeig(i)
1176 ab(2,2) = reig(i)
1177 CALL sgemm( 'N', 'N', m, 2, 2, -one, z(1,i), &
1178 ldz, ab, 2, one, y(1,i), ldy ) ! BLAS CALL
1179 ! Y(1:M,i:i+1) = Y(1:M,i:i+1) - Z(1:M,i:i+1) * AB ! INTRINSIC
1180 res(i) = slange( 'F', m, 2, y(1,i), ldy, &
1181 work(n+1) ) ! LAPACK CALL
1182 res(i+1) = res(i)
1183 i = i + 2
1184 END IF
1185 END DO
1186 END IF
1187 END IF
1188!
1189 IF ( whtsvd == 4 ) THEN
1190 work(n+1) = xscl1
1191 work(n+2) = xscl2
1192 END IF
1193!
1194! Successful exit.
1195 IF ( .NOT. badxy ) THEN
1196 info = 0
1197 ELSE
1198 ! A warning on possible data inconsistency.
1199 ! This should be a rare event.
1200 info = 4
1201 END IF
1202!............................................................
1203 RETURN
1204! ......
1205 END SUBROUTINE sgedmd
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
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
subroutine sgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)
SGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition sgeev.f:191
subroutine sgejsv(joba, jobu, jobv, jobr, jobt, jobp, m, n, a, lda, sva, u, ldu, v, ldv, work, lwork, iwork, info)
SGEJSV
Definition sgejsv.f:474
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine sgesdd(jobz, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, iwork, info)
SGESDD
Definition sgesdd.f:211
subroutine sgesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, info)
SGESVD computes the singular value decomposition (SVD) for GE matrices
Definition sgesvd.f:210
subroutine sgesvdq(joba, jobp, jobr, jobu, jobv, m, n, a, lda, s, u, ldu, v, ldv, numrank, iwork, liwork, work, lwork, rwork, lrwork, info)
SGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE...
Definition sgesvdq.f:413
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:101
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition slange.f:112
subroutine slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition slascl.f:142
subroutine slassq(n, x, incx, scale, sumsq)
SLASSQ updates a sum of squares represented in scaled form.
Definition slassq.f90:122
real(wp) function snrm2(n, x, incx)
SNRM2
Definition snrm2.f90:89
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79