|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, 00002 $ LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, 00003 $ INFO ) 00004 * 00005 * -- LAPACK auxiliary routine (version 3.2) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * November 2006 00009 * 00010 * .. Scalar Arguments .. 00011 CHARACTER TRANS 00012 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N 00013 REAL RDSCAL, RDSUM, SCALE 00014 * .. 00015 * .. Array Arguments .. 00016 COMPLEX A( LDA, * ), B( LDB, * ), C( LDC, * ), 00017 $ D( LDD, * ), E( LDE, * ), F( LDF, * ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * CTGSY2 solves the generalized Sylvester equation 00024 * 00025 * A * R - L * B = scale * C (1) 00026 * D * R - L * E = scale * F 00027 * 00028 * using Level 1 and 2 BLAS, where R and L are unknown M-by-N matrices, 00029 * (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M, 00030 * N-by-N and M-by-N, respectively. A, B, D and E are upper triangular 00031 * (i.e., (A,D) and (B,E) in generalized Schur form). 00032 * 00033 * The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output 00034 * scaling factor chosen to avoid overflow. 00035 * 00036 * In matrix notation solving equation (1) corresponds to solve 00037 * Zx = scale * b, where Z is defined as 00038 * 00039 * Z = [ kron(In, A) -kron(B', Im) ] (2) 00040 * [ kron(In, D) -kron(E', Im) ], 00041 * 00042 * Ik is the identity matrix of size k and X' is the transpose of X. 00043 * kron(X, Y) is the Kronecker product between the matrices X and Y. 00044 * 00045 * If TRANS = 'C', y in the conjugate transposed system Z'y = scale*b 00046 * is solved for, which is equivalent to solve for R and L in 00047 * 00048 * A' * R + D' * L = scale * C (3) 00049 * R * B' + L * E' = scale * -F 00050 * 00051 * This case is used to compute an estimate of Dif[(A, D), (B, E)] = 00052 * = sigma_min(Z) using reverse communicaton with CLACON. 00053 * 00054 * CTGSY2 also (IJOB >= 1) contributes to the computation in CTGSYL 00055 * of an upper bound on the separation between to matrix pairs. Then 00056 * the input (A, D), (B, E) are sub-pencils of two matrix pairs in 00057 * CTGSYL. 00058 * 00059 * Arguments 00060 * ========= 00061 * 00062 * TRANS (input) CHARACTER*1 00063 * = 'N', solve the generalized Sylvester equation (1). 00064 * = 'T': solve the 'transposed' system (3). 00065 * 00066 * IJOB (input) INTEGER 00067 * Specifies what kind of functionality to be performed. 00068 * =0: solve (1) only. 00069 * =1: A contribution from this subsystem to a Frobenius 00070 * norm-based estimate of the separation between two matrix 00071 * pairs is computed. (look ahead strategy is used). 00072 * =2: A contribution from this subsystem to a Frobenius 00073 * norm-based estimate of the separation between two matrix 00074 * pairs is computed. (SGECON on sub-systems is used.) 00075 * Not referenced if TRANS = 'T'. 00076 * 00077 * M (input) INTEGER 00078 * On entry, M specifies the order of A and D, and the row 00079 * dimension of C, F, R and L. 00080 * 00081 * N (input) INTEGER 00082 * On entry, N specifies the order of B and E, and the column 00083 * dimension of C, F, R and L. 00084 * 00085 * A (input) COMPLEX array, dimension (LDA, M) 00086 * On entry, A contains an upper triangular matrix. 00087 * 00088 * LDA (input) INTEGER 00089 * The leading dimension of the matrix A. LDA >= max(1, M). 00090 * 00091 * B (input) COMPLEX array, dimension (LDB, N) 00092 * On entry, B contains an upper triangular matrix. 00093 * 00094 * LDB (input) INTEGER 00095 * The leading dimension of the matrix B. LDB >= max(1, N). 00096 * 00097 * C (input/output) COMPLEX array, dimension (LDC, N) 00098 * On entry, C contains the right-hand-side of the first matrix 00099 * equation in (1). 00100 * On exit, if IJOB = 0, C has been overwritten by the solution 00101 * R. 00102 * 00103 * LDC (input) INTEGER 00104 * The leading dimension of the matrix C. LDC >= max(1, M). 00105 * 00106 * D (input) COMPLEX array, dimension (LDD, M) 00107 * On entry, D contains an upper triangular matrix. 00108 * 00109 * LDD (input) INTEGER 00110 * The leading dimension of the matrix D. LDD >= max(1, M). 00111 * 00112 * E (input) COMPLEX array, dimension (LDE, N) 00113 * On entry, E contains an upper triangular matrix. 00114 * 00115 * LDE (input) INTEGER 00116 * The leading dimension of the matrix E. LDE >= max(1, N). 00117 * 00118 * F (input/output) COMPLEX array, dimension (LDF, N) 00119 * On entry, F contains the right-hand-side of the second matrix 00120 * equation in (1). 00121 * On exit, if IJOB = 0, F has been overwritten by the solution 00122 * L. 00123 * 00124 * LDF (input) INTEGER 00125 * The leading dimension of the matrix F. LDF >= max(1, M). 00126 * 00127 * SCALE (output) REAL 00128 * On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions 00129 * R and L (C and F on entry) will hold the solutions to a 00130 * slightly perturbed system but the input matrices A, B, D and 00131 * E have not been changed. If SCALE = 0, R and L will hold the 00132 * solutions to the homogeneous system with C = F = 0. 00133 * Normally, SCALE = 1. 00134 * 00135 * RDSUM (input/output) REAL 00136 * On entry, the sum of squares of computed contributions to 00137 * the Dif-estimate under computation by CTGSYL, where the 00138 * scaling factor RDSCAL (see below) has been factored out. 00139 * On exit, the corresponding sum of squares updated with the 00140 * contributions from the current sub-system. 00141 * If TRANS = 'T' RDSUM is not touched. 00142 * NOTE: RDSUM only makes sense when CTGSY2 is called by 00143 * CTGSYL. 00144 * 00145 * RDSCAL (input/output) REAL 00146 * On entry, scaling factor used to prevent overflow in RDSUM. 00147 * On exit, RDSCAL is updated w.r.t. the current contributions 00148 * in RDSUM. 00149 * If TRANS = 'T', RDSCAL is not touched. 00150 * NOTE: RDSCAL only makes sense when CTGSY2 is called by 00151 * CTGSYL. 00152 * 00153 * INFO (output) INTEGER 00154 * On exit, if INFO is set to 00155 * =0: Successful exit 00156 * <0: If INFO = -i, input argument number i is illegal. 00157 * >0: The matrix pairs (A, D) and (B, E) have common or very 00158 * close eigenvalues. 00159 * 00160 * Further Details 00161 * =============== 00162 * 00163 * Based on contributions by 00164 * Bo Kagstrom and Peter Poromaa, Department of Computing Science, 00165 * Umea University, S-901 87 Umea, Sweden. 00166 * 00167 * ===================================================================== 00168 * 00169 * .. Parameters .. 00170 REAL ZERO, ONE 00171 INTEGER LDZ 00172 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, LDZ = 2 ) 00173 * .. 00174 * .. Local Scalars .. 00175 LOGICAL NOTRAN 00176 INTEGER I, IERR, J, K 00177 REAL SCALOC 00178 COMPLEX ALPHA 00179 * .. 00180 * .. Local Arrays .. 00181 INTEGER IPIV( LDZ ), JPIV( LDZ ) 00182 COMPLEX RHS( LDZ ), Z( LDZ, LDZ ) 00183 * .. 00184 * .. External Functions .. 00185 LOGICAL LSAME 00186 EXTERNAL LSAME 00187 * .. 00188 * .. External Subroutines .. 00189 EXTERNAL CAXPY, CGESC2, CGETC2, CSCAL, CLATDF, XERBLA 00190 * .. 00191 * .. Intrinsic Functions .. 00192 INTRINSIC CMPLX, CONJG, MAX 00193 * .. 00194 * .. Executable Statements .. 00195 * 00196 * Decode and test input parameters 00197 * 00198 INFO = 0 00199 IERR = 0 00200 NOTRAN = LSAME( TRANS, 'N' ) 00201 IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN 00202 INFO = -1 00203 ELSE IF( NOTRAN ) THEN 00204 IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.2 ) ) THEN 00205 INFO = -2 00206 END IF 00207 END IF 00208 IF( INFO.EQ.0 ) THEN 00209 IF( M.LE.0 ) THEN 00210 INFO = -3 00211 ELSE IF( N.LE.0 ) THEN 00212 INFO = -4 00213 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00214 INFO = -5 00215 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00216 INFO = -8 00217 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 00218 INFO = -10 00219 ELSE IF( LDD.LT.MAX( 1, M ) ) THEN 00220 INFO = -12 00221 ELSE IF( LDE.LT.MAX( 1, N ) ) THEN 00222 INFO = -14 00223 ELSE IF( LDF.LT.MAX( 1, M ) ) THEN 00224 INFO = -16 00225 END IF 00226 END IF 00227 IF( INFO.NE.0 ) THEN 00228 CALL XERBLA( 'CTGSY2', -INFO ) 00229 RETURN 00230 END IF 00231 * 00232 IF( NOTRAN ) THEN 00233 * 00234 * Solve (I, J) - system 00235 * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J) 00236 * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J) 00237 * for I = M, M - 1, ..., 1; J = 1, 2, ..., N 00238 * 00239 SCALE = ONE 00240 SCALOC = ONE 00241 DO 30 J = 1, N 00242 DO 20 I = M, 1, -1 00243 * 00244 * Build 2 by 2 system 00245 * 00246 Z( 1, 1 ) = A( I, I ) 00247 Z( 2, 1 ) = D( I, I ) 00248 Z( 1, 2 ) = -B( J, J ) 00249 Z( 2, 2 ) = -E( J, J ) 00250 * 00251 * Set up right hand side(s) 00252 * 00253 RHS( 1 ) = C( I, J ) 00254 RHS( 2 ) = F( I, J ) 00255 * 00256 * Solve Z * x = RHS 00257 * 00258 CALL CGETC2( LDZ, Z, LDZ, IPIV, JPIV, IERR ) 00259 IF( IERR.GT.0 ) 00260 $ INFO = IERR 00261 IF( IJOB.EQ.0 ) THEN 00262 CALL CGESC2( LDZ, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) 00263 IF( SCALOC.NE.ONE ) THEN 00264 DO 10 K = 1, N 00265 CALL CSCAL( M, CMPLX( SCALOC, ZERO ), C( 1, K ), 00266 $ 1 ) 00267 CALL CSCAL( M, CMPLX( SCALOC, ZERO ), F( 1, K ), 00268 $ 1 ) 00269 10 CONTINUE 00270 SCALE = SCALE*SCALOC 00271 END IF 00272 ELSE 00273 CALL CLATDF( IJOB, LDZ, Z, LDZ, RHS, RDSUM, RDSCAL, 00274 $ IPIV, JPIV ) 00275 END IF 00276 * 00277 * Unpack solution vector(s) 00278 * 00279 C( I, J ) = RHS( 1 ) 00280 F( I, J ) = RHS( 2 ) 00281 * 00282 * Substitute R(I, J) and L(I, J) into remaining equation. 00283 * 00284 IF( I.GT.1 ) THEN 00285 ALPHA = -RHS( 1 ) 00286 CALL CAXPY( I-1, ALPHA, A( 1, I ), 1, C( 1, J ), 1 ) 00287 CALL CAXPY( I-1, ALPHA, D( 1, I ), 1, F( 1, J ), 1 ) 00288 END IF 00289 IF( J.LT.N ) THEN 00290 CALL CAXPY( N-J, RHS( 2 ), B( J, J+1 ), LDB, 00291 $ C( I, J+1 ), LDC ) 00292 CALL CAXPY( N-J, RHS( 2 ), E( J, J+1 ), LDE, 00293 $ F( I, J+1 ), LDF ) 00294 END IF 00295 * 00296 20 CONTINUE 00297 30 CONTINUE 00298 ELSE 00299 * 00300 * Solve transposed (I, J) - system: 00301 * A(I, I)' * R(I, J) + D(I, I)' * L(J, J) = C(I, J) 00302 * R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J) 00303 * for I = 1, 2, ..., M, J = N, N - 1, ..., 1 00304 * 00305 SCALE = ONE 00306 SCALOC = ONE 00307 DO 80 I = 1, M 00308 DO 70 J = N, 1, -1 00309 * 00310 * Build 2 by 2 system Z' 00311 * 00312 Z( 1, 1 ) = CONJG( A( I, I ) ) 00313 Z( 2, 1 ) = -CONJG( B( J, J ) ) 00314 Z( 1, 2 ) = CONJG( D( I, I ) ) 00315 Z( 2, 2 ) = -CONJG( E( J, J ) ) 00316 * 00317 * 00318 * Set up right hand side(s) 00319 * 00320 RHS( 1 ) = C( I, J ) 00321 RHS( 2 ) = F( I, J ) 00322 * 00323 * Solve Z' * x = RHS 00324 * 00325 CALL CGETC2( LDZ, Z, LDZ, IPIV, JPIV, IERR ) 00326 IF( IERR.GT.0 ) 00327 $ INFO = IERR 00328 CALL CGESC2( LDZ, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) 00329 IF( SCALOC.NE.ONE ) THEN 00330 DO 40 K = 1, N 00331 CALL CSCAL( M, CMPLX( SCALOC, ZERO ), C( 1, K ), 00332 $ 1 ) 00333 CALL CSCAL( M, CMPLX( SCALOC, ZERO ), F( 1, K ), 00334 $ 1 ) 00335 40 CONTINUE 00336 SCALE = SCALE*SCALOC 00337 END IF 00338 * 00339 * Unpack solution vector(s) 00340 * 00341 C( I, J ) = RHS( 1 ) 00342 F( I, J ) = RHS( 2 ) 00343 * 00344 * Substitute R(I, J) and L(I, J) into remaining equation. 00345 * 00346 DO 50 K = 1, J - 1 00347 F( I, K ) = F( I, K ) + RHS( 1 )*CONJG( B( K, J ) ) + 00348 $ RHS( 2 )*CONJG( E( K, J ) ) 00349 50 CONTINUE 00350 DO 60 K = I + 1, M 00351 C( K, J ) = C( K, J ) - CONJG( A( I, K ) )*RHS( 1 ) - 00352 $ CONJG( D( I, K ) )*RHS( 2 ) 00353 60 CONTINUE 00354 * 00355 70 CONTINUE 00356 80 CONTINUE 00357 END IF 00358 RETURN 00359 * 00360 * End of CTGSY2 00361 * 00362 END
1.7.3 
