1 SUBROUTINE SGGES_F95( A, B, ALPHAR, ALPHAI, BETA, VSL, VSR, SELECT, SDIM, INFO ) 2! 3! -- LAPACK95 interface driver routine (version 3.0) -- 4! UNI-C, Denmark; Univ. of Tennessee, USA; NAG Ltd., UK 5! September, 2000 6! 7! .. USE STATEMENTS .. 8 USE LA_PRECISION, ONLY: WP => SP 9 USE LA_AUXMOD, ONLY: ERINFO, LSAME 10 USE F77_LAPACK, ONLY: GGES_F77 => LA_GGES 11! .. IMPLICIT STATEMENT .. 12 IMPLICIT NONE 13! .. SCALAR ARGUMENTS .. 14 INTEGER, INTENT(OUT), OPTIONAL :: INFO 15 INTEGER, INTENT(OUT), OPTIONAL :: SDIM 16! .. ARRAY ARGUMENTS .. 17 REAL(WP), INTENT(INOUT) :: A(:,:), B(:,:) 18 REAL(WP), INTENT(OUT) :: ALPHAR(:), ALPHAI(:), BETA(:) 19 REAL(WP), INTENT(OUT), OPTIONAL, TARGET :: VSL(:,:), VSR(:,:) 20! .. FUNCTIONAL ARGUMENTS .. 21 INTERFACE 22 LOGICAL FUNCTION SELECT(ALPHAR, ALPHAI, BETA) 23 USE LA_PRECISION, ONLY: WP => SP 24 REAL(WP), INTENT(IN) :: ALPHAR, ALPHAI, BETA 25 END FUNCTION SELECT 26 END INTERFACE 27 OPTIONAL :: SELECT 28!---------------------------------------------------------------------- 29! 30! Purpose 31! ======= 32! 33! LA_GGES computes for a pair of n by n real or complex matrices 34! (A, B) the (generalized) real or complex Schur form, the generalized 35! eigenvalues in the form of scalar pairs (alpha,beta), and, optionally, 36! the left and/or right Schur vectors. 37! If A and B are real then the real-Schur form is computed, 38! otherwise the complex-Schur form is computed. The real-Schur form is a 39! pair of real matrices (S,T) such that 1) S has block upper triangular 40! form, with 1 by 1 and 2 by 2 blocks along the main diagonal, 2) T has 41! upper triangular form with nonnegative elements on the main diagonal, 42! and 3) S = Q^T*A*Z and T = Q^T*B*Z, where Q and Z are orthogonal 43! matrices. The 2 by 2 blocks of S are "standardized" by making the 44! corresponding elements of T have the form 45! [ a 0 ] 46! [ 0 b ] 47! The complex-Schur form is a pair of matrices (S,T) such that 1) S has 48! upper triangular form, 2) T has upper triangular form with nonnegative 49! elements on the main diagonal, and 3) S = Q^H*A*Z and T = Q^H*B*Z, 50! where Q and Z are unitary matrices. 51! In both cases the columns of Q and Z are called, respectively, 52! the left and right (generalized) Schur vectors. 53! A generalized eigenvalue of the pair (A,B) is, roughly speaking, a 54! scalar of the form lambda = alpha/beta such that the matrix 55! A -lambda*B is singular. It is usually represented as the pair 56! (alpha, beta), as there is a reasonable interpretation of the case 57! beta = 0 (even if alpha = 0). 58! 59! ========= 60! 61! SUBROUTINE LA_GGES( A, B, <alpha>, BETA, VSL=vsl, & 62! VSR=vsr, SELECT=select, SDIM=sdim, INFO=info ) 63! <type>(<wp>), INTENT(INOUT) :: A(:,:), B(:,:) 64! <type>(<wp>), INTENT(OUT) :: <alpha(:)>, BETA(:) 65! <type>(<wp>), INTENT(OUT), OPTIONAL :: VSL(:,:), VSR(:,:) 66! INTERFACE 67! LOGICAL FUNCTION SELECT(<alpha(j)> , BETA(j)) 68! <type>(<wp>), INTENT(IN) :: <alpha(j)> , BETA(j) 69! END FUNCTION SELECT 70! END INTERFACE 71! OPTIONAL :: SELECT 72! INTEGER, INTENT(OUT), OPTIONAL :: SDIM 73! INTEGER, INTENT(OUT), OPTIONAL :: INFO 74! where 75! <type> ::= REAL | COMPLEX 76! <wp> ::= KIND(1.0) | KIND(1.0D0) 77! <alpha> ::= ALPHAR, ALPHAI | ALPHA 78! <alpha(:)> ::= ALPHAR(:), ALPHAI(:) | ALPHA(:) 79! <alpha(j)> ::= ALPHAR(j) , ALPHAI(j) | ALPHA(j) 80! 81! Arguments 82! ========= 83! 84! A (input/output) REAL or COMPLEX square array, shape (:,:). 85! On entry, the matrix A. 86! On exit, the matrix S. 87! B (input/output) REAL or COMPLEX square array, shape (:,:) with 88! size(B,1) = size(A,1). 89! On entry, the matrix B. 90! On exit, the matrix T . 91! <alpha> (output) REAL or COMPLEX array, shape (:) with size(alpha) = 92! size(A,1). 93! The values of alpha. 94! <alpha(:)> ::= ALPHAR(:), ALPHAI(:) | ALPHA(:), 95! where 96! ALPHAR(:), ALPHAI(:) are of REAL type (for the real and 97! imaginary parts) and ALPHA(:) is of COMPLEX type. 98! BETA (output) REAL or COMPLEX array, shape (:) with size(BETA) = 99! size(A,1). 100! The values of beta. 101! Note: The generalized eigenvalues of the pair (A,B) are the 102! scalars lambda(j) = alpha(j)/beta(j). These quotients may 103! easily over- or underflow, and beta(j) may even be zero. Thus, 104! the user should avoid computing them naively. 105! Note: If A and B are real then complex eigenvalues occur in 106! complex conjugate pairs. Each pair is stored consecutively. 107! Thus a complex conjugate pair is given by 108! lambda(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j) 109! lambda(j+1) = (ALPHAR(j+1) + i*ALPHAI(j+1))/BETA(j+1) 110! where 111! ALPHAI(j)/BETA(j) = -( ALPHAI(j+1)/BETA(j+1)) 112! VSL Optional (output) REAL or COMPLEX square array, shape (:,:) 113! with size(VSL,1) = size(A,1). 114! The left Schur vectors. 115! VSR Optional (output) REAL or COMPLEX square array, shape (:,:) 116! with size(VSR,1) = size(A,1). 117! The right Schur vectors. 118! SELECT Optional (input) LOGICAL FUNCTION 119! LOGICAL FUNCTION SELECT( <alpha(j)> , BETA(j)) 120! <type>(<wp>), INTENT(IN) :: <alpha(j)> , BETA(j) 121! where 122! <type> ::= REAL | COMPLEX 123! <wp> ::= KIND(1.0) | KIND(1.0D0) 124! <alpha(j)> ::= ALPHAR(j) , ALPHAI(j) | ALPHA(j) 125! 1. SELECT must be declared as EXTERNAL or as an explicit 126! interface in the calling (sub)program. 127! 2. SELECT is called by LA_GGES for every computed eigenvalue 128! (<alpha(j)> , BETA(j)) (but only once for a complex conjugate 129! pair when A and B are real). It is used to select the 130! eigenvalues that will be ordered to the top left of the Schur 131! form. The eigenvalue (<alpha(j)>, BETA(j)) is selected if 132! SELECT(<alpha(j)>, BETA(j)) has the value .TRUE. 133! 3. A selected complex eigenvalue may no longer satisfy 134! SELECT(<alpha(j)>, BETA(j)) = .TRUE. after ordering, since 135! ordering may change the value of complex eigenvalues 136! (especially if the eigenvalue is ill-conditioned); in this case 137! INFO is set to size(A,1) + 2 (see INFO below). 138! Note: Select must be present if SDIM is desired. 139! SDIM Optional (output) INTEGER. 140! The number of eigenvalues (after sorting) for which SELECT = 141! .TRUE. (If A and B are real, then complex conjugate pairs for 142! which SELECT = .TRUE. for either eigenvalue count as 2). 143! INFO Optional (output) INTEGER. 144! = 0: successful exit. 145! < 0: if INFO = -i, the i-th argument had an illegal value. 146! > 0: if INFO = i, and i is 147! <= n: the QZ iteration failed. The matrix pair (A,B) has not 148! been reduced to Schur form, but (<alpha(j)>,BETA(j)) 149! should be correct for j = INFO + 1,..., n. 150! = n+1: another part of the algorithm failed. 151! = n+2: after reordering, roundoff changed values of some 152! complex eigenvalues so that leading eigenvalues in the 153! Schur form no longer satisfy SELECT = .TRUE. This can be 154! caused by ordinary roundoff or underflow due to scaling. 155! = n+3: the reordering failed. 156! If INFO is not present and an error occurs, then the program 157! is terminated with an error message. 158!----------------------------------------------------------------------- 159! .. LOCAL PARAMETERS .. 160 CHARACTER(LEN=7), PARAMETER :: SRNAME = 'LA_GGES' 161! .. LOCAL SCALARS .. 162 CHARACTER(LEN=1) :: LJOBVSL, LJOBVSR, LSORT 163 INTEGER, SAVE :: LWORK = 0 164 INTEGER :: N, LINFO, LDA, LDB, ISTAT, S1VSL, S2VSL, S1VSR, S2VSR, & 165 & SALPHAR, SALPHAI, SBETA 166 INTEGER :: LSDIM 167! .. LOCAL ARRAYS .. 168 REAL(WP), TARGET :: LLVSL(1,1), LLVSR(1,1), WORKMIN(1) 169 REAL(WP), POINTER :: WORK(:) 170 LOGICAL, POINTER :: BWORK(:) 171 LOGICAL, TARGET :: LLBWORK(1) 172! .. INTRINSIC FUNCTIONS .. 173 INTRINSIC MAX, PRESENT, SIZE 174! .. EXECUTABLE STATEMENTS .. 175 LINFO = 0; ISTAT = 0; N = SIZE(A,1); LDA = MAX(1,N); LDB = MAX(1,SIZE(B,1)) 176 SALPHAR = SIZE(ALPHAR); SALPHAI = SIZE(ALPHAI) 177 IF (PRESENT (SELECT)) THEN 178 LSORT = 'S'; ELSE ; LSORT = 'N'; ENDIF 179 SBETA = SIZE(BETA) 180 IF( PRESENT(VSL) )THEN; S1VSL = SIZE(VSL,1); S2VSL = SIZE(VSL,2); LJOBVSL = 'V' 181 ELSE; S1VSL = 1; S2VSL = 1; LJOBVSL = 'N'; END IF 182 IF( PRESENT(VSR) )THEN; S1VSR = SIZE(VSR,1); S2VSR = SIZE(VSR,2); LJOBVSR = 'V' 183 ELSE; S1VSR = 1; S2VSR = 1; LJOBVSR = 'N'; END IF 184! .. TEST THE ARGUMENTS 185 IF( N < 0 .OR. SIZE(A,2) /= N )THEN; LINFO = -1 186 ELSE IF( SIZE(B,1) /= N .OR. SIZE(B,2) /= N )THEN; LINFO = -2 187 ELSE IF( SALPHAR /= N )THEN; LINFO = -3 188 ELSE IF( SALPHAI /= N )THEN; LINFO = -4 189 ELSE IF( SBETA /= N )THEN; LINFO = -5 190 ELSE IF( PRESENT(VSL) .AND. ( S1VSL /= N .OR. S2VSL /= N ) )THEN; LINFO = -6 191 ELSE IF( PRESENT(VSR) .AND. ( S1VSR /= N .OR. S2VSR /= N ) )THEN; LINFO = -7 192 ELSE IF( N >= 0 )THEN 193 IF(LSAME(LSORT,'S')) THEN; ALLOCATE(BWORK(N),STAT=ISTAT) 194 IF (ISTAT /= 0) THEN ; LINFO=-100; GOTO 500; ENDIF 195 ELSE; BWORK => LLBWORK; END IF 196 197 198! .. DETERMINE THE WORKSPACE .. 199! .. QUERING THE SIZE OF WORKSPACE .. 200 LWORK = -1 201 IF( PRESENT(VSL) )THEN 202 IF( PRESENT(VSR) )THEN 203 CALL GGES_F77( LJOBVSL, LJOBVSR, LSORT, SELECT, N, A, LDA, B, LDB, & 204& LSDIM, ALPHAR, ALPHAI, BETA, VSL, MAX(1,S1VSL), VSR, MAX(1,S1VSR), & 205& WORKMIN, LWORK, BWORK, LINFO ) 206 ELSE 207 CALL GGES_F77( LJOBVSL, LJOBVSR, LSORT, SELECT, N, A, LDA, B, LDB, & 208& LSDIM, ALPHAR, ALPHAI, BETA, VSL, MAX(1,S1VSL), LLVSR, MAX(1,S1VSR), & 209& WORKMIN, LWORK, BWORK, LINFO ) 210 ENDIF 211 ELSE 212 IF( PRESENT(VSR) )THEN 213 CALL GGES_F77( LJOBVSL, LJOBVSR, LSORT, SELECT, N, A, LDA, B, LDB, & 214& LSDIM, ALPHAR, ALPHAI, BETA, LLVSL, MAX(1,S1VSL), VSR, MAX(1,S1VSR), & 215& WORKMIN, LWORK, BWORK, LINFO ) 216 ELSE 217 CALL GGES_F77( LJOBVSL, LJOBVSR, LSORT, SELECT, N, A, LDA, B, LDB, & 218& LSDIM, ALPHAR, ALPHAI, BETA, LLVSL, MAX(1,S1VSL), LLVSR, MAX(1,S1VSR), & 219& WORKMIN, LWORK, BWORK, LINFO ) 220 ENDIF 221 ENDIF 222 LWORK = WORKMIN(1) 223 224 ALLOCATE(WORK(LWORK), STAT=ISTAT) 225 IF (ISTAT /= 0) THEN ; LINFO=-100; GOTO 400; ENDIF 226 227 IF( PRESENT(VSL) )THEN 228 IF( PRESENT(VSR) )THEN 229 CALL GGES_F77( LJOBVSL, LJOBVSR, LSORT, SELECT, N, A, LDA, B, LDB, & 230& LSDIM, ALPHAR, ALPHAI, BETA, VSL, MAX(1,S1VSL), VSR, MAX(1,S1VSR), & 231& WORK, LWORK, BWORK, LINFO ) 232 ELSE 233 CALL GGES_F77( LJOBVSL, LJOBVSR, LSORT, SELECT, N, A, LDA, B, LDB, & 234& LSDIM, ALPHAR, ALPHAI, BETA, VSL, MAX(1,S1VSL), LLVSR, MAX(1,S1VSR), & 235& WORK, LWORK, BWORK, LINFO ) 236 ENDIF 237 ELSE 238 IF( PRESENT(VSR) )THEN 239 CALL GGES_F77( LJOBVSL, LJOBVSR, LSORT, SELECT, N, A, LDA, B, LDB, & 240& LSDIM, ALPHAR, ALPHAI, BETA, LLVSL, MAX(1,S1VSL), VSR, MAX(1,S1VSR), & 241& WORK, LWORK, BWORK, LINFO ) 242 ELSE 243 CALL GGES_F77( LJOBVSL, LJOBVSR, LSORT, SELECT, N, A, LDA, B, LDB, & 244& LSDIM, ALPHAR, ALPHAI, BETA, LLVSL, MAX(1,S1VSL), LLVSR, MAX(1,S1VSR), & 245& WORK, LWORK, BWORK, LINFO ) 246 ENDIF 247 ENDIF 248 IF (PRESENT(SDIM)) SDIM = LSDIM 249 250 DEALLOCATE(WORK) 251400 IF(LSAME(LSORT,'S')) DEALLOCATE(BWORK) 252 ENDIF 253500 CALL ERINFO(LINFO,SRNAME,INFO,ISTAT) 254END SUBROUTINE SGGES_F95 255