1SUBROUTINE DGGSVD_F95( A, B, ALPHA, BETA, K, L, U, V, Q, IWORK, 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 => DP 9 USE LA_AUXMOD, ONLY: ERINFO 10 USE F77_LAPACK, ONLY: GGSVD_F77 => LA_GGSVD 11! .. IMPLICIT STATEMENT .. 12 IMPLICIT NONE 13! .. SCALAR ARGUMENTS .. 14 INTEGER, INTENT(OUT), OPTIONAL :: INFO, K, L 15! .. ARRAY ARGUMENTS .. 16 REAL(WP), INTENT(INOUT) :: A(:,:), B(:,:) 17 REAL(WP), INTENT(OUT) :: ALPHA(:), BETA(:) 18 REAL(WP), INTENT(OUT), OPTIONAL, TARGET :: U(:,:), V(:,:), Q(:,:) 19 INTEGER, INTENT(OUT), OPTIONAL, TARGET :: IWORK(:) 20!---------------------------------------------------------------------- 21! 22! Purpose 23! ======= 24! 25! LA_GGSVD computes the generalized singular values and, optionally, 26! the transformation matrices from the generalized singular value 27! decomposition (GSVD) of a real or complex matrix pair (A,B), where A 28! is m by n and B is p by n. The GSVD of (A,B) is written 29! A = U * SIGMA1(0, R)*Q^H , B = V * SIGMA2(0, R)*Q^H 30! where U , V and Q are orthogonal (unitary) matrices of dimensions m by m, 31! p by p and n by n, respectively. Let l be the rank of B and r the rank of 32! the (m + p) * n matrix ( A ) 33! ( B ) 34! , and let k = r-l. Then SIGMA1 and SIGMA2 are m*(k + l) and p * (k + l) 35! "diagonal" matrices, respectively, and R is a (k + l) * (k + l) 36! nonsingular triangular matrix. The detailed structure of SIGMA1 ,SIGMA2 37! and R depends on the sign of (m - k - l) as follows: 38! The case m-k-l>=0: 39! 40! k l 41! k ( I 0 ) 42! SIGMA1 = l ( 0 C ) 43! m-k-l ( 0 0 ) 44! 45! 46! k l 47! SIGMA2 = l ( 0 S ) 48! p - l ( 0 S ) 49! 50! 51! n-k-l k l 52! (0, R) = k ( 0 R11 R12 ) 53! l ( 0 0 R22 ) 54! 55! where C^2 + S^2 = I . We define 56! alpha(1)=alpha(2)=...=alpha(k) = 1, alpha(k+i)=c(i i), i=1,2,...,l 57! beta(1) = beta(2) = ... = beta(k) = 0, beta(k+i) = s(i i), i=1,2,...,l 58! 59! The case m-k-l < 0: 60! 61! k m-k k+l-m 62! SIGMA1 = k ( I 0 0 ) 63! m-k ( 0 C 0 ) 64! 65! 66! k m-k k+l-m 67! m-k ( 0 S 0 ) 68! SIGMA2 = k+l-m ( 0 0 I ) 69! p-l ( 0 0 0 ) 70! 71! 72! n-k-l k m-k k+l-m 73! k ( 0 R11 R12 R13 ) 74! (0,R) = m-k ( 0 0 R22 R23 ) 75! k+l-m ( 0 0 0 R33 ) 76! 77! where C^2 + S^2 = I . We define 78! alpha(1)=alpha(2)=...=alpha(k) = 1, alpha(k+i)=c(i i), i =1,2,...,m-k, 79! alpha(m+1) = alpha(m+2)=...= alpha(k+l) = 0 80! beta(1)=beta(2)= ... =beta(k)=0, beta(k+i)=s(i i), i=1,2,...,m-k, 81! beta(m+1) = beta(m+2) = ... = beta(k+l) = 1 82! 83! In both cases the generalized singular values of the pair (A,B) are the 84! ratios 85! sigma(i) = alpha(i)/beta(i), i = 1,2, ... ,k+l 86! 87! The first k singular values are infinite. The finite singular values 88! are real and nonnegative. 89! LA_GGSVD computes the real (nonnegative) scalars alpha(i), beta(i), 90! i=1,2,..., k+l , the matrix R, and, optionally, the transformation 91! matrices U , V and Q. 92! 93! ========= 94! 95! SUBROUTINE LA_GGSVD( A, B, ALPHA, BETA, K=k, L=l, & 96! U=u, V=v, Q=q, IWORK=iwork, INFO=info ) 97! <type>(<wp>), INTENT(INOUT) :: A(:,:), B(:,:) 98! REAL(<wp>), INTENT(OUT) :: ALPHA(:), BETA(:) 99! INTEGER, INTENT(OUT), OPTIONAL :: K, L 100! <type>(<wp>), INTENT(OUT), OPTIONAL :: U(:,:), V(:,:), Q(:,:) 101! INTEGER, INTENT(IN), OPTIONAL :: IWORK(:) 102! INTEGER, INTENT(OUT), OPTIONAL :: INFO 103! where 104! <type> ::= REAL | COMPLEX 105! <wp> ::= KIND(1.0) | KIND(1.0D0) 106! 107! Arguments 108! ========= 109! 110! A (input/output) REAL or COMPLEX array, shape (:,:) with 111! size(A,1) = m and size(A,2) = n. 112! On entry, the matrix A. 113! On exit, A contains the triangular matrix R, or part of R, as 114! follows: 115! If m-k-l >= 0, then R is stored in A(1:k+l,n-k-l+1:n). 116! If m-k-l < 0, then the matrix 117! ( R11 R12 R13 ) 118! ( 0 R22 R23 ) 119! is stored in A(1:m,n-k-l+1:n). 120! B (input/output) REAL or COMPLEX array, shape (:,:) with 121! size(B,1) = p and size(B,2) = n. 122! On entry, the matrix B. 123! On exit, if m-k-l < 0, then R33 is stored in 124! B(m-k+1:l,n+m-k-l+1:n). 125! ALPHA (output) REAL array, shape (:) with size(ALPHA) = n 126! The real scalars alpha(i) , i = 1, 2,..., k+l. 127! BETA (output) REAL array, shape (:) with size(BETA) = n. 128! The real scalars beta(i) , i = 1, 2, ..., k+l. 129! Note: The generalized singular values of the pair (A,B) are 130! sigma(i) = ALPHA(i)/BETA(i), i = 1, 2, ..., k+l. 131! If k + l < n, then ALPHA(k+l+1:n) = BETA(k+l+1:n) = 0. 132! K, L Optional (output) INTEGER. 133! The dimension parameters k and l. 134! U Optional (output) REAL or COMPLEX square array, shape (:,:) with 135! size(U,1) = m. 136! The matrix U . 137! V Optional (output) REAL or COMPLEX square array, shape (:,:) with 138! size(V,1) = p. 139! The matrix V . 140! Q Optional (output) REAL or COMPLEX square array, shape (:,:) with 141! size(Q,1) = n. 142! The matrix Q. 143! IWORK Optional (output) INTEGER array, shape(:) with size(IWORK) = n. 144! IWORK contains sorting information. More precisely, the loop 145! for i = k + 1, min(m, k + l) 146! swap ALPHA(i) and ALPHA(IWORK(i)) 147! end 148! will sort ALPHA so that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(n). 149! INFO Optional (output) INTEGER. 150! = 0: successful exit. 151! < 0: if INFO = -i, the i-th argument had an illegal value. 152! > 0: if INFO = 1, the algorithm failed to converge. 153! If INFO is not present and an error occurs, then the program is 154! terminated with an error message. 155!---------------------------------------------------------------------- 156! .. LOCAL PARAMETERS .. 157 CHARACTER(LEN=8), PARAMETER :: SRNAME = 'LA_GGSVD' 158! .. LOCAL SCALARS .. 159 CHARACTER(LEN=1) :: LJOBU, LJOBV, LJOBQ 160 INTEGER :: M, N, P, LINFO, ISTAT, ISTAT1, S1U, S2U, S1V, S2V, & 161 S1Q, S2Q, LK, LL 162! .. LOCAL ARRAYS .. 163 REAL(WP), TARGET :: LLU(1,1), LLV(1,1), LLQ(1,1) 164 INTEGER, POINTER :: LIWORK(:) 165 REAL(WP), POINTER :: WORK(:) 166! .. INTRINSIC FUNCTIONS .. 167 INTRINSIC MAX, PRESENT, SIZE 168! .. EXECUTABLE STATEMENTS .. 169 LINFO = 0; ISTAT = 0; M = SIZE(A,1); N = SIZE(A,2); P = SIZE(B,1) 170 IF( PRESENT(U) )THEN; S1U = SIZE(U,1); S2U = SIZE(U,2); LJOBU = 'U' 171 ELSE; S1U = 1; S2U = 1; LJOBU = 'N'; END IF 172 IF( PRESENT(V) )THEN; S1V = SIZE(V,1); S2V = SIZE(V,2); LJOBV = 'V' 173 ELSE; S1V = 1; S2V = 1; LJOBV = 'N'; END IF 174 IF( PRESENT(Q) )THEN; S1Q = SIZE(Q,1); S2Q = SIZE(Q,2); LJOBQ = 'Q' 175 ELSE; S1Q = 1; S2Q = 1; LJOBQ = 'N'; END IF 176! .. TEST THE ARGUMENTS 177 IF( M < 0 .OR. N < N )THEN; LINFO = -1 178 ELSE IF( SIZE(B,2) /= N .OR. P < 0 ) THEN; LINFO = -2 179 ELSE IF( SIZE( ALPHA ) /= N )THEN; LINFO = -3 180 ELSE IF( SIZE( BETA ) /= N )THEN; LINFO = -4 181 ELSE IF( PRESENT(U) .AND. ( S1U /= MAX(1,M) .OR. S2U /= M ) )THEN; LINFO = -7 182 ELSE IF( PRESENT(V) .AND. ( S1V /= MAX(1,P) .OR. S2V /= P ) )THEN; LINFO = -8 183 ELSE IF( PRESENT(Q) .AND. ( S1Q /= MAX(1,N) .OR. S2Q /= N ) )THEN; LINFO = -9 184 ELSE 185 IF (PRESENT(IWORK)) THEN 186 LIWORK => IWORK 187 ELSE 188 ALLOCATE( LIWORK(MAX(1,N)), STAT=ISTAT ) 189 ENDIF 190 IF( ISTAT == 0 ) THEN 191 ALLOCATE( WORK( MAX(1,MAX(3*N,M,P)+N) ), STAT=ISTAT) 192 END IF 193 IF( ISTAT == 0 ) THEN 194 IF( PRESENT(U) )THEN 195 IF( PRESENT(V) )THEN 196 IF( PRESENT(Q) )THEN 197 CALL GGSVD_F77( LJOBU, LJOBV, LJOBQ, M, N, P, LK, LL, A, MAX(1,M), & 198 B, MAX(1,P), ALPHA, BETA, U, MAX(1,S1U), & 199 V, MAX(1,S1V), Q, MAX(1,S1Q), WORK, LIWORK, LINFO ) 200 ELSE 201 CALL GGSVD_F77( LJOBU, LJOBV, LJOBQ, M, N, P, LK, LL, A, MAX(1,M), & 202 B, MAX(1,P), ALPHA, BETA, U, MAX(1,S1U), & 203 V, MAX(1,S1V), LLQ, MAX(1,S1Q), WORK, LIWORK, LINFO ) 204 ENDIF 205 ELSE 206 IF( PRESENT(Q) )THEN 207 CALL GGSVD_F77( LJOBU, LJOBV, LJOBQ, M, N, P, LK, LL, A, MAX(1,M), & 208 B, MAX(1,P), ALPHA, BETA, U, MAX(1,S1U), & 209 LLV, MAX(1,S1V), Q, MAX(1,S1Q), WORK, LIWORK, LINFO ) 210 ELSE 211 CALL GGSVD_F77( LJOBU, LJOBV, LJOBQ, M, N, P, LK, LL, A, MAX(1,M), & 212 B, MAX(1,P), ALPHA, BETA, U, MAX(1,S1U), & 213 LLV, MAX(1,S1V), LLQ, MAX(1,S1Q), WORK, LIWORK, LINFO ) 214 ENDIF 215 ENDIF 216 ELSE 217 IF( PRESENT(V) )THEN 218 IF( PRESENT(Q) )THEN 219 CALL GGSVD_F77( LJOBU, LJOBV, LJOBQ, M, N, P, LK, LL, A, MAX(1,M), & 220 B, MAX(1,P), ALPHA, BETA, LLU, MAX(1,S1U), & 221 V, MAX(1,S1V), Q, MAX(1,S1Q), WORK, LIWORK, LINFO ) 222 ELSE 223 CALL GGSVD_F77( LJOBU, LJOBV, LJOBQ, M, N, P, LK, LL, A, MAX(1,M), & 224 B, MAX(1,P), ALPHA, BETA, LLU, MAX(1,S1U), & 225 V, MAX(1,S1V), LLQ, MAX(1,S1Q), WORK, LIWORK, LINFO ) 226 ENDIF 227 ELSE 228 IF( PRESENT(Q) )THEN 229 CALL GGSVD_F77( LJOBU, LJOBV, LJOBQ, M, N, P, LK, LL, A, MAX(1,M), & 230 B, MAX(1,P), ALPHA, BETA, LLU, MAX(1,S1U), & 231 LLV, MAX(1,S1V), Q, MAX(1,S1Q), WORK, LIWORK, LINFO ) 232 ELSE 233 CALL GGSVD_F77( LJOBU, LJOBV, LJOBQ, M, N, P, LK, LL, A, MAX(1,M), & 234 B, MAX(1,P), ALPHA, BETA, LLU, MAX(1,S1U), & 235 LLV, MAX(1,S1V), LLQ, MAX(1,S1Q), WORK, LIWORK, LINFO ) 236 ENDIF 237 ENDIF 238 ENDIF 239 IF( PRESENT(K) ) K = LK 240 IF( PRESENT(L) ) L = LL 241 ELSE; LINFO = -100; ENDIF 242 DEALLOCATE(WORK, STAT=ISTAT1) 243 IF (.NOT. PRESENT(IWORK)) DEALLOCATE(LIWORK, STAT=ISTAT1) 244 ENDIF 245 CALL ERINFO(LINFO,SRNAME,INFO,ISTAT) 246END SUBROUTINE DGGSVD_F95 247