 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, 00002 $ LDC, SCALE, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.3.1) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * -- April 2011 -- 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER TRANA, TRANB 00011 INTEGER INFO, ISGN, LDA, LDB, LDC, M, N 00012 DOUBLE PRECISION SCALE 00013 * .. 00014 * .. Array Arguments .. 00015 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DTRSYL solves the real Sylvester matrix equation: 00022 * 00023 * op(A)*X + X*op(B) = scale*C or 00024 * op(A)*X - X*op(B) = scale*C, 00025 * 00026 * where op(A) = A or A**T, and A and B are both upper quasi- 00027 * triangular. A is M-by-M and B is N-by-N; the right hand side C and 00028 * the solution X are M-by-N; and scale is an output scale factor, set 00029 * <= 1 to avoid overflow in X. 00030 * 00031 * A and B must be in Schur canonical form (as returned by DHSEQR), that 00032 * is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; 00033 * each 2-by-2 diagonal block has its diagonal elements equal and its 00034 * off-diagonal elements of opposite sign. 00035 * 00036 * Arguments 00037 * ========= 00038 * 00039 * TRANA (input) CHARACTER*1 00040 * Specifies the option op(A): 00041 * = 'N': op(A) = A (No transpose) 00042 * = 'T': op(A) = A**T (Transpose) 00043 * = 'C': op(A) = A**H (Conjugate transpose = Transpose) 00044 * 00045 * TRANB (input) CHARACTER*1 00046 * Specifies the option op(B): 00047 * = 'N': op(B) = B (No transpose) 00048 * = 'T': op(B) = B**T (Transpose) 00049 * = 'C': op(B) = B**H (Conjugate transpose = Transpose) 00050 * 00051 * ISGN (input) INTEGER 00052 * Specifies the sign in the equation: 00053 * = +1: solve op(A)*X + X*op(B) = scale*C 00054 * = -1: solve op(A)*X - X*op(B) = scale*C 00055 * 00056 * M (input) INTEGER 00057 * The order of the matrix A, and the number of rows in the 00058 * matrices X and C. M >= 0. 00059 * 00060 * N (input) INTEGER 00061 * The order of the matrix B, and the number of columns in the 00062 * matrices X and C. N >= 0. 00063 * 00064 * A (input) DOUBLE PRECISION array, dimension (LDA,M) 00065 * The upper quasi-triangular matrix A, in Schur canonical form. 00066 * 00067 * LDA (input) INTEGER 00068 * The leading dimension of the array A. LDA >= max(1,M). 00069 * 00070 * B (input) DOUBLE PRECISION array, dimension (LDB,N) 00071 * The upper quasi-triangular matrix B, in Schur canonical form. 00072 * 00073 * LDB (input) INTEGER 00074 * The leading dimension of the array B. LDB >= max(1,N). 00075 * 00076 * C (input/output) DOUBLE PRECISION array, dimension (LDC,N) 00077 * On entry, the M-by-N right hand side matrix C. 00078 * On exit, C is overwritten by the solution matrix X. 00079 * 00080 * LDC (input) INTEGER 00081 * The leading dimension of the array C. LDC >= max(1,M) 00082 * 00083 * SCALE (output) DOUBLE PRECISION 00084 * The scale factor, scale, set <= 1 to avoid overflow in X. 00085 * 00086 * INFO (output) INTEGER 00087 * = 0: successful exit 00088 * < 0: if INFO = -i, the i-th argument had an illegal value 00089 * = 1: A and B have common or very close eigenvalues; perturbed 00090 * values were used to solve the equation (but the matrices 00091 * A and B are unchanged). 00092 * 00093 * ===================================================================== 00094 * 00095 * .. Parameters .. 00096 DOUBLE PRECISION ZERO, ONE 00097 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00098 * .. 00099 * .. Local Scalars .. 00100 LOGICAL NOTRNA, NOTRNB 00101 INTEGER IERR, J, K, K1, K2, KNEXT, L, L1, L2, LNEXT 00102 DOUBLE PRECISION A11, BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN, 00103 $ SMLNUM, SUML, SUMR, XNORM 00104 * .. 00105 * .. Local Arrays .. 00106 DOUBLE PRECISION DUM( 1 ), VEC( 2, 2 ), X( 2, 2 ) 00107 * .. 00108 * .. External Functions .. 00109 LOGICAL LSAME 00110 DOUBLE PRECISION DDOT, DLAMCH, DLANGE 00111 EXTERNAL LSAME, DDOT, DLAMCH, DLANGE 00112 * .. 00113 * .. External Subroutines .. 00114 EXTERNAL DLABAD, DLALN2, DLASY2, DSCAL, XERBLA 00115 * .. 00116 * .. Intrinsic Functions .. 00117 INTRINSIC ABS, DBLE, MAX, MIN 00118 * .. 00119 * .. Executable Statements .. 00120 * 00121 * Decode and Test input parameters 00122 * 00123 NOTRNA = LSAME( TRANA, 'N' ) 00124 NOTRNB = LSAME( TRANB, 'N' ) 00125 * 00126 INFO = 0 00127 IF( .NOT.NOTRNA .AND. .NOT.LSAME( TRANA, 'T' ) .AND. .NOT. 00128 $ LSAME( TRANA, 'C' ) ) THEN 00129 INFO = -1 00130 ELSE IF( .NOT.NOTRNB .AND. .NOT.LSAME( TRANB, 'T' ) .AND. .NOT. 00131 $ LSAME( TRANB, 'C' ) ) THEN 00132 INFO = -2 00133 ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN 00134 INFO = -3 00135 ELSE IF( M.LT.0 ) THEN 00136 INFO = -4 00137 ELSE IF( N.LT.0 ) THEN 00138 INFO = -5 00139 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00140 INFO = -7 00141 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00142 INFO = -9 00143 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 00144 INFO = -11 00145 END IF 00146 IF( INFO.NE.0 ) THEN 00147 CALL XERBLA( 'DTRSYL', -INFO ) 00148 RETURN 00149 END IF 00150 * 00151 * Quick return if possible 00152 * 00153 SCALE = ONE 00154 IF( M.EQ.0 .OR. N.EQ.0 ) 00155 $ RETURN 00156 * 00157 * Set constants to control overflow 00158 * 00159 EPS = DLAMCH( 'P' ) 00160 SMLNUM = DLAMCH( 'S' ) 00161 BIGNUM = ONE / SMLNUM 00162 CALL DLABAD( SMLNUM, BIGNUM ) 00163 SMLNUM = SMLNUM*DBLE( M*N ) / EPS 00164 BIGNUM = ONE / SMLNUM 00165 * 00166 SMIN = MAX( SMLNUM, EPS*DLANGE( 'M', M, M, A, LDA, DUM ), 00167 $ EPS*DLANGE( 'M', N, N, B, LDB, DUM ) ) 00168 * 00169 SGN = ISGN 00170 * 00171 IF( NOTRNA .AND. NOTRNB ) THEN 00172 * 00173 * Solve A*X + ISGN*X*B = scale*C. 00174 * 00175 * The (K,L)th block of X is determined starting from 00176 * bottom-left corner column by column by 00177 * 00178 * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) 00179 * 00180 * Where 00181 * M L-1 00182 * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)]. 00183 * I=K+1 J=1 00184 * 00185 * Start column loop (index = L) 00186 * L1 (L2) : column index of the first (first) row of X(K,L). 00187 * 00188 LNEXT = 1 00189 DO 60 L = 1, N 00190 IF( L.LT.LNEXT ) 00191 $ GO TO 60 00192 IF( L.EQ.N ) THEN 00193 L1 = L 00194 L2 = L 00195 ELSE 00196 IF( B( L+1, L ).NE.ZERO ) THEN 00197 L1 = L 00198 L2 = L + 1 00199 LNEXT = L + 2 00200 ELSE 00201 L1 = L 00202 L2 = L 00203 LNEXT = L + 1 00204 END IF 00205 END IF 00206 * 00207 * Start row loop (index = K) 00208 * K1 (K2): row index of the first (last) row of X(K,L). 00209 * 00210 KNEXT = M 00211 DO 50 K = M, 1, -1 00212 IF( K.GT.KNEXT ) 00213 $ GO TO 50 00214 IF( K.EQ.1 ) THEN 00215 K1 = K 00216 K2 = K 00217 ELSE 00218 IF( A( K, K-1 ).NE.ZERO ) THEN 00219 K1 = K - 1 00220 K2 = K 00221 KNEXT = K - 2 00222 ELSE 00223 K1 = K 00224 K2 = K 00225 KNEXT = K - 1 00226 END IF 00227 END IF 00228 * 00229 IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN 00230 SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 00231 $ C( MIN( K1+1, M ), L1 ), 1 ) 00232 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 00233 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 00234 SCALOC = ONE 00235 * 00236 A11 = A( K1, K1 ) + SGN*B( L1, L1 ) 00237 DA11 = ABS( A11 ) 00238 IF( DA11.LE.SMIN ) THEN 00239 A11 = SMIN 00240 DA11 = SMIN 00241 INFO = 1 00242 END IF 00243 DB = ABS( VEC( 1, 1 ) ) 00244 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 00245 IF( DB.GT.BIGNUM*DA11 ) 00246 $ SCALOC = ONE / DB 00247 END IF 00248 X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 00249 * 00250 IF( SCALOC.NE.ONE ) THEN 00251 DO 10 J = 1, N 00252 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 00253 10 CONTINUE 00254 SCALE = SCALE*SCALOC 00255 END IF 00256 C( K1, L1 ) = X( 1, 1 ) 00257 * 00258 ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN 00259 * 00260 SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 00261 $ C( MIN( K2+1, M ), L1 ), 1 ) 00262 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 00263 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 00264 * 00265 SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 00266 $ C( MIN( K2+1, M ), L1 ), 1 ) 00267 SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 00268 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 00269 * 00270 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ), 00271 $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), 00272 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 00273 IF( IERR.NE.0 ) 00274 $ INFO = 1 00275 * 00276 IF( SCALOC.NE.ONE ) THEN 00277 DO 20 J = 1, N 00278 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 00279 20 CONTINUE 00280 SCALE = SCALE*SCALOC 00281 END IF 00282 C( K1, L1 ) = X( 1, 1 ) 00283 C( K2, L1 ) = X( 2, 1 ) 00284 * 00285 ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN 00286 * 00287 SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 00288 $ C( MIN( K1+1, M ), L1 ), 1 ) 00289 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 00290 VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) 00291 * 00292 SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 00293 $ C( MIN( K1+1, M ), L2 ), 1 ) 00294 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 00295 VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) 00296 * 00297 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, B( L1, L1 ), 00298 $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), 00299 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 00300 IF( IERR.NE.0 ) 00301 $ INFO = 1 00302 * 00303 IF( SCALOC.NE.ONE ) THEN 00304 DO 30 J = 1, N 00305 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 00306 30 CONTINUE 00307 SCALE = SCALE*SCALOC 00308 END IF 00309 C( K1, L1 ) = X( 1, 1 ) 00310 C( K1, L2 ) = X( 2, 1 ) 00311 * 00312 ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN 00313 * 00314 SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 00315 $ C( MIN( K2+1, M ), L1 ), 1 ) 00316 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 00317 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 00318 * 00319 SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 00320 $ C( MIN( K2+1, M ), L2 ), 1 ) 00321 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 00322 VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) 00323 * 00324 SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 00325 $ C( MIN( K2+1, M ), L1 ), 1 ) 00326 SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 00327 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 00328 * 00329 SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 00330 $ C( MIN( K2+1, M ), L2 ), 1 ) 00331 SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 ) 00332 VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) 00333 * 00334 CALL DLASY2( .FALSE., .FALSE., ISGN, 2, 2, 00335 $ A( K1, K1 ), LDA, B( L1, L1 ), LDB, VEC, 00336 $ 2, SCALOC, X, 2, XNORM, IERR ) 00337 IF( IERR.NE.0 ) 00338 $ INFO = 1 00339 * 00340 IF( SCALOC.NE.ONE ) THEN 00341 DO 40 J = 1, N 00342 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 00343 40 CONTINUE 00344 SCALE = SCALE*SCALOC 00345 END IF 00346 C( K1, L1 ) = X( 1, 1 ) 00347 C( K1, L2 ) = X( 1, 2 ) 00348 C( K2, L1 ) = X( 2, 1 ) 00349 C( K2, L2 ) = X( 2, 2 ) 00350 END IF 00351 * 00352 50 CONTINUE 00353 * 00354 60 CONTINUE 00355 * 00356 ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN 00357 * 00358 * Solve A**T *X + ISGN*X*B = scale*C. 00359 * 00360 * The (K,L)th block of X is determined starting from 00361 * upper-left corner column by column by 00362 * 00363 * A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) 00364 * 00365 * Where 00366 * K-1 L-1 00367 * R(K,L) = SUM [A(I,K)**T*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)] 00368 * I=1 J=1 00369 * 00370 * Start column loop (index = L) 00371 * L1 (L2): column index of the first (last) row of X(K,L) 00372 * 00373 LNEXT = 1 00374 DO 120 L = 1, N 00375 IF( L.LT.LNEXT ) 00376 $ GO TO 120 00377 IF( L.EQ.N ) THEN 00378 L1 = L 00379 L2 = L 00380 ELSE 00381 IF( B( L+1, L ).NE.ZERO ) THEN 00382 L1 = L 00383 L2 = L + 1 00384 LNEXT = L + 2 00385 ELSE 00386 L1 = L 00387 L2 = L 00388 LNEXT = L + 1 00389 END IF 00390 END IF 00391 * 00392 * Start row loop (index = K) 00393 * K1 (K2): row index of the first (last) row of X(K,L) 00394 * 00395 KNEXT = 1 00396 DO 110 K = 1, M 00397 IF( K.LT.KNEXT ) 00398 $ GO TO 110 00399 IF( K.EQ.M ) THEN 00400 K1 = K 00401 K2 = K 00402 ELSE 00403 IF( A( K+1, K ).NE.ZERO ) THEN 00404 K1 = K 00405 K2 = K + 1 00406 KNEXT = K + 2 00407 ELSE 00408 K1 = K 00409 K2 = K 00410 KNEXT = K + 1 00411 END IF 00412 END IF 00413 * 00414 IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN 00415 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 00416 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 00417 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 00418 SCALOC = ONE 00419 * 00420 A11 = A( K1, K1 ) + SGN*B( L1, L1 ) 00421 DA11 = ABS( A11 ) 00422 IF( DA11.LE.SMIN ) THEN 00423 A11 = SMIN 00424 DA11 = SMIN 00425 INFO = 1 00426 END IF 00427 DB = ABS( VEC( 1, 1 ) ) 00428 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 00429 IF( DB.GT.BIGNUM*DA11 ) 00430 $ SCALOC = ONE / DB 00431 END IF 00432 X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 00433 * 00434 IF( SCALOC.NE.ONE ) THEN 00435 DO 70 J = 1, N 00436 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 00437 70 CONTINUE 00438 SCALE = SCALE*SCALOC 00439 END IF 00440 C( K1, L1 ) = X( 1, 1 ) 00441 * 00442 ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN 00443 * 00444 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 00445 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 00446 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 00447 * 00448 SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 00449 SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 00450 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 00451 * 00452 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ), 00453 $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), 00454 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 00455 IF( IERR.NE.0 ) 00456 $ INFO = 1 00457 * 00458 IF( SCALOC.NE.ONE ) THEN 00459 DO 80 J = 1, N 00460 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 00461 80 CONTINUE 00462 SCALE = SCALE*SCALOC 00463 END IF 00464 C( K1, L1 ) = X( 1, 1 ) 00465 C( K2, L1 ) = X( 2, 1 ) 00466 * 00467 ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN 00468 * 00469 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 00470 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 00471 VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) 00472 * 00473 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 00474 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 00475 VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) 00476 * 00477 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, B( L1, L1 ), 00478 $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), 00479 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 00480 IF( IERR.NE.0 ) 00481 $ INFO = 1 00482 * 00483 IF( SCALOC.NE.ONE ) THEN 00484 DO 90 J = 1, N 00485 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 00486 90 CONTINUE 00487 SCALE = SCALE*SCALOC 00488 END IF 00489 C( K1, L1 ) = X( 1, 1 ) 00490 C( K1, L2 ) = X( 2, 1 ) 00491 * 00492 ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN 00493 * 00494 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 00495 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 00496 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 00497 * 00498 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 00499 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 00500 VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) 00501 * 00502 SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 00503 SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 00504 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 00505 * 00506 SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 ) 00507 SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 ) 00508 VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) 00509 * 00510 CALL DLASY2( .TRUE., .FALSE., ISGN, 2, 2, A( K1, K1 ), 00511 $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X, 00512 $ 2, XNORM, IERR ) 00513 IF( IERR.NE.0 ) 00514 $ INFO = 1 00515 * 00516 IF( SCALOC.NE.ONE ) THEN 00517 DO 100 J = 1, N 00518 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 00519 100 CONTINUE 00520 SCALE = SCALE*SCALOC 00521 END IF 00522 C( K1, L1 ) = X( 1, 1 ) 00523 C( K1, L2 ) = X( 1, 2 ) 00524 C( K2, L1 ) = X( 2, 1 ) 00525 C( K2, L2 ) = X( 2, 2 ) 00526 END IF 00527 * 00528 110 CONTINUE 00529 120 CONTINUE 00530 * 00531 ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN 00532 * 00533 * Solve A**T*X + ISGN*X*B**T = scale*C. 00534 * 00535 * The (K,L)th block of X is determined starting from 00536 * top-right corner column by column by 00537 * 00538 * A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L) 00539 * 00540 * Where 00541 * K-1 N 00542 * R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T]. 00543 * I=1 J=L+1 00544 * 00545 * Start column loop (index = L) 00546 * L1 (L2): column index of the first (last) row of X(K,L) 00547 * 00548 LNEXT = N 00549 DO 180 L = N, 1, -1 00550 IF( L.GT.LNEXT ) 00551 $ GO TO 180 00552 IF( L.EQ.1 ) THEN 00553 L1 = L 00554 L2 = L 00555 ELSE 00556 IF( B( L, L-1 ).NE.ZERO ) THEN 00557 L1 = L - 1 00558 L2 = L 00559 LNEXT = L - 2 00560 ELSE 00561 L1 = L 00562 L2 = L 00563 LNEXT = L - 1 00564 END IF 00565 END IF 00566 * 00567 * Start row loop (index = K) 00568 * K1 (K2): row index of the first (last) row of X(K,L) 00569 * 00570 KNEXT = 1 00571 DO 170 K = 1, M 00572 IF( K.LT.KNEXT ) 00573 $ GO TO 170 00574 IF( K.EQ.M ) THEN 00575 K1 = K 00576 K2 = K 00577 ELSE 00578 IF( A( K+1, K ).NE.ZERO ) THEN 00579 K1 = K 00580 K2 = K + 1 00581 KNEXT = K + 2 00582 ELSE 00583 K1 = K 00584 K2 = K 00585 KNEXT = K + 1 00586 END IF 00587 END IF 00588 * 00589 IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN 00590 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 00591 SUMR = DDOT( N-L1, C( K1, MIN( L1+1, N ) ), LDC, 00592 $ B( L1, MIN( L1+1, N ) ), LDB ) 00593 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 00594 SCALOC = ONE 00595 * 00596 A11 = A( K1, K1 ) + SGN*B( L1, L1 ) 00597 DA11 = ABS( A11 ) 00598 IF( DA11.LE.SMIN ) THEN 00599 A11 = SMIN 00600 DA11 = SMIN 00601 INFO = 1 00602 END IF 00603 DB = ABS( VEC( 1, 1 ) ) 00604 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 00605 IF( DB.GT.BIGNUM*DA11 ) 00606 $ SCALOC = ONE / DB 00607 END IF 00608 X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 00609 * 00610 IF( SCALOC.NE.ONE ) THEN 00611 DO 130 J = 1, N 00612 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 00613 130 CONTINUE 00614 SCALE = SCALE*SCALOC 00615 END IF 00616 C( K1, L1 ) = X( 1, 1 ) 00617 * 00618 ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN 00619 * 00620 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 00621 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 00622 $ B( L1, MIN( L2+1, N ) ), LDB ) 00623 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 00624 * 00625 SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 00626 SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 00627 $ B( L1, MIN( L2+1, N ) ), LDB ) 00628 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 00629 * 00630 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ), 00631 $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), 00632 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 00633 IF( IERR.NE.0 ) 00634 $ INFO = 1 00635 * 00636 IF( SCALOC.NE.ONE ) THEN 00637 DO 140 J = 1, N 00638 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 00639 140 CONTINUE 00640 SCALE = SCALE*SCALOC 00641 END IF 00642 C( K1, L1 ) = X( 1, 1 ) 00643 C( K2, L1 ) = X( 2, 1 ) 00644 * 00645 ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN 00646 * 00647 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 00648 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 00649 $ B( L1, MIN( L2+1, N ) ), LDB ) 00650 VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) 00651 * 00652 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 00653 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 00654 $ B( L2, MIN( L2+1, N ) ), LDB ) 00655 VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) 00656 * 00657 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, B( L1, L1 ), 00658 $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), 00659 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 00660 IF( IERR.NE.0 ) 00661 $ INFO = 1 00662 * 00663 IF( SCALOC.NE.ONE ) THEN 00664 DO 150 J = 1, N 00665 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 00666 150 CONTINUE 00667 SCALE = SCALE*SCALOC 00668 END IF 00669 C( K1, L1 ) = X( 1, 1 ) 00670 C( K1, L2 ) = X( 2, 1 ) 00671 * 00672 ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN 00673 * 00674 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 00675 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 00676 $ B( L1, MIN( L2+1, N ) ), LDB ) 00677 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 00678 * 00679 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 00680 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 00681 $ B( L2, MIN( L2+1, N ) ), LDB ) 00682 VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) 00683 * 00684 SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 00685 SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 00686 $ B( L1, MIN( L2+1, N ) ), LDB ) 00687 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 00688 * 00689 SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 ) 00690 SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 00691 $ B( L2, MIN( L2+1, N ) ), LDB ) 00692 VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) 00693 * 00694 CALL DLASY2( .TRUE., .TRUE., ISGN, 2, 2, A( K1, K1 ), 00695 $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X, 00696 $ 2, XNORM, IERR ) 00697 IF( IERR.NE.0 ) 00698 $ INFO = 1 00699 * 00700 IF( SCALOC.NE.ONE ) THEN 00701 DO 160 J = 1, N 00702 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 00703 160 CONTINUE 00704 SCALE = SCALE*SCALOC 00705 END IF 00706 C( K1, L1 ) = X( 1, 1 ) 00707 C( K1, L2 ) = X( 1, 2 ) 00708 C( K2, L1 ) = X( 2, 1 ) 00709 C( K2, L2 ) = X( 2, 2 ) 00710 END IF 00711 * 00712 170 CONTINUE 00713 180 CONTINUE 00714 * 00715 ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN 00716 * 00717 * Solve A*X + ISGN*X*B**T = scale*C. 00718 * 00719 * The (K,L)th block of X is determined starting from 00720 * bottom-right corner column by column by 00721 * 00722 * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L) 00723 * 00724 * Where 00725 * M N 00726 * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T]. 00727 * I=K+1 J=L+1 00728 * 00729 * Start column loop (index = L) 00730 * L1 (L2): column index of the first (last) row of X(K,L) 00731 * 00732 LNEXT = N 00733 DO 240 L = N, 1, -1 00734 IF( L.GT.LNEXT ) 00735 $ GO TO 240 00736 IF( L.EQ.1 ) THEN 00737 L1 = L 00738 L2 = L 00739 ELSE 00740 IF( B( L, L-1 ).NE.ZERO ) THEN 00741 L1 = L - 1 00742 L2 = L 00743 LNEXT = L - 2 00744 ELSE 00745 L1 = L 00746 L2 = L 00747 LNEXT = L - 1 00748 END IF 00749 END IF 00750 * 00751 * Start row loop (index = K) 00752 * K1 (K2): row index of the first (last) row of X(K,L) 00753 * 00754 KNEXT = M 00755 DO 230 K = M, 1, -1 00756 IF( K.GT.KNEXT ) 00757 $ GO TO 230 00758 IF( K.EQ.1 ) THEN 00759 K1 = K 00760 K2 = K 00761 ELSE 00762 IF( A( K, K-1 ).NE.ZERO ) THEN 00763 K1 = K - 1 00764 K2 = K 00765 KNEXT = K - 2 00766 ELSE 00767 K1 = K 00768 K2 = K 00769 KNEXT = K - 1 00770 END IF 00771 END IF 00772 * 00773 IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN 00774 SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 00775 $ C( MIN( K1+1, M ), L1 ), 1 ) 00776 SUMR = DDOT( N-L1, C( K1, MIN( L1+1, N ) ), LDC, 00777 $ B( L1, MIN( L1+1, N ) ), LDB ) 00778 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 00779 SCALOC = ONE 00780 * 00781 A11 = A( K1, K1 ) + SGN*B( L1, L1 ) 00782 DA11 = ABS( A11 ) 00783 IF( DA11.LE.SMIN ) THEN 00784 A11 = SMIN 00785 DA11 = SMIN 00786 INFO = 1 00787 END IF 00788 DB = ABS( VEC( 1, 1 ) ) 00789 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 00790 IF( DB.GT.BIGNUM*DA11 ) 00791 $ SCALOC = ONE / DB 00792 END IF 00793 X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 00794 * 00795 IF( SCALOC.NE.ONE ) THEN 00796 DO 190 J = 1, N 00797 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 00798 190 CONTINUE 00799 SCALE = SCALE*SCALOC 00800 END IF 00801 C( K1, L1 ) = X( 1, 1 ) 00802 * 00803 ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN 00804 * 00805 SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 00806 $ C( MIN( K2+1, M ), L1 ), 1 ) 00807 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 00808 $ B( L1, MIN( L2+1, N ) ), LDB ) 00809 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 00810 * 00811 SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 00812 $ C( MIN( K2+1, M ), L1 ), 1 ) 00813 SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 00814 $ B( L1, MIN( L2+1, N ) ), LDB ) 00815 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 00816 * 00817 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ), 00818 $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), 00819 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 00820 IF( IERR.NE.0 ) 00821 $ INFO = 1 00822 * 00823 IF( SCALOC.NE.ONE ) THEN 00824 DO 200 J = 1, N 00825 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 00826 200 CONTINUE 00827 SCALE = SCALE*SCALOC 00828 END IF 00829 C( K1, L1 ) = X( 1, 1 ) 00830 C( K2, L1 ) = X( 2, 1 ) 00831 * 00832 ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN 00833 * 00834 SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 00835 $ C( MIN( K1+1, M ), L1 ), 1 ) 00836 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 00837 $ B( L1, MIN( L2+1, N ) ), LDB ) 00838 VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) 00839 * 00840 SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 00841 $ C( MIN( K1+1, M ), L2 ), 1 ) 00842 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 00843 $ B( L2, MIN( L2+1, N ) ), LDB ) 00844 VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) 00845 * 00846 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, B( L1, L1 ), 00847 $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), 00848 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 00849 IF( IERR.NE.0 ) 00850 $ INFO = 1 00851 * 00852 IF( SCALOC.NE.ONE ) THEN 00853 DO 210 J = 1, N 00854 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 00855 210 CONTINUE 00856 SCALE = SCALE*SCALOC 00857 END IF 00858 C( K1, L1 ) = X( 1, 1 ) 00859 C( K1, L2 ) = X( 2, 1 ) 00860 * 00861 ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN 00862 * 00863 SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 00864 $ C( MIN( K2+1, M ), L1 ), 1 ) 00865 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 00866 $ B( L1, MIN( L2+1, N ) ), LDB ) 00867 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 00868 * 00869 SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 00870 $ C( MIN( K2+1, M ), L2 ), 1 ) 00871 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 00872 $ B( L2, MIN( L2+1, N ) ), LDB ) 00873 VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) 00874 * 00875 SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 00876 $ C( MIN( K2+1, M ), L1 ), 1 ) 00877 SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 00878 $ B( L1, MIN( L2+1, N ) ), LDB ) 00879 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 00880 * 00881 SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 00882 $ C( MIN( K2+1, M ), L2 ), 1 ) 00883 SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 00884 $ B( L2, MIN( L2+1, N ) ), LDB ) 00885 VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) 00886 * 00887 CALL DLASY2( .FALSE., .TRUE., ISGN, 2, 2, A( K1, K1 ), 00888 $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X, 00889 $ 2, XNORM, IERR ) 00890 IF( IERR.NE.0 ) 00891 $ INFO = 1 00892 * 00893 IF( SCALOC.NE.ONE ) THEN 00894 DO 220 J = 1, N 00895 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 00896 220 CONTINUE 00897 SCALE = SCALE*SCALOC 00898 END IF 00899 C( K1, L1 ) = X( 1, 1 ) 00900 C( K1, L2 ) = X( 1, 2 ) 00901 C( K2, L1 ) = X( 2, 1 ) 00902 C( K2, L2 ) = X( 2, 2 ) 00903 END IF 00904 * 00905 230 CONTINUE 00906 240 CONTINUE 00907 * 00908 END IF 00909 * 00910 RETURN 00911 * 00912 * End of DTRSYL 00913 * 00914 END
1.7.3 
