 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE STRSYL( 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 REAL SCALE 00013 * .. 00014 * .. Array Arguments .. 00015 REAL A( LDA, * ), B( LDB, * ), C( LDC, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * STRSYL 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 SHSEQR), 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) REAL 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) REAL 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) REAL 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) REAL 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 REAL ZERO, ONE 00097 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00098 * .. 00099 * .. Local Scalars .. 00100 LOGICAL NOTRNA, NOTRNB 00101 INTEGER IERR, J, K, K1, K2, KNEXT, L, L1, L2, LNEXT 00102 REAL A11, BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN, 00103 $ SMLNUM, SUML, SUMR, XNORM 00104 * .. 00105 * .. Local Arrays .. 00106 REAL DUM( 1 ), VEC( 2, 2 ), X( 2, 2 ) 00107 * .. 00108 * .. External Functions .. 00109 LOGICAL LSAME 00110 REAL SDOT, SLAMCH, SLANGE 00111 EXTERNAL LSAME, SDOT, SLAMCH, SLANGE 00112 * .. 00113 * .. External Subroutines .. 00114 EXTERNAL SLABAD, SLALN2, SLASY2, SSCAL, XERBLA 00115 * .. 00116 * .. Intrinsic Functions .. 00117 INTRINSIC ABS, MAX, MIN, REAL 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( 'STRSYL', -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 = SLAMCH( 'P' ) 00160 SMLNUM = SLAMCH( 'S' ) 00161 BIGNUM = ONE / SMLNUM 00162 CALL SLABAD( SMLNUM, BIGNUM ) 00163 SMLNUM = SMLNUM*REAL( M*N ) / EPS 00164 BIGNUM = ONE / SMLNUM 00165 * 00166 SMIN = MAX( SMLNUM, EPS*SLANGE( 'M', M, M, A, LDA, DUM ), 00167 $ EPS*SLANGE( '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 70 L = 1, N 00190 IF( L.LT.LNEXT ) 00191 $ GO TO 70 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 60 K = M, 1, -1 00212 IF( K.GT.KNEXT ) 00213 $ GO TO 60 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 = SDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 00231 $ C( MIN( K1+1, M ), L1 ), 1 ) 00232 SUMR = SDOT( 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 SSCAL( 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 = SDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 00261 $ C( MIN( K2+1, M ), L1 ), 1 ) 00262 SUMR = SDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 00263 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 00264 * 00265 SUML = SDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 00266 $ C( MIN( K2+1, M ), L1 ), 1 ) 00267 SUMR = SDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 00268 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 00269 * 00270 CALL SLALN2( .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 SSCAL( 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 = SDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 00288 $ C( MIN( K1+1, M ), L1 ), 1 ) 00289 SUMR = SDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 00290 VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) 00291 * 00292 SUML = SDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 00293 $ C( MIN( K1+1, M ), L2 ), 1 ) 00294 SUMR = SDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 00295 VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) 00296 * 00297 CALL SLALN2( .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 40 J = 1, N 00305 CALL SSCAL( M, SCALOC, C( 1, J ), 1 ) 00306 40 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 = SDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 00315 $ C( MIN( K2+1, M ), L1 ), 1 ) 00316 SUMR = SDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 00317 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 00318 * 00319 SUML = SDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 00320 $ C( MIN( K2+1, M ), L2 ), 1 ) 00321 SUMR = SDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 00322 VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) 00323 * 00324 SUML = SDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 00325 $ C( MIN( K2+1, M ), L1 ), 1 ) 00326 SUMR = SDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 00327 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 00328 * 00329 SUML = SDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 00330 $ C( MIN( K2+1, M ), L2 ), 1 ) 00331 SUMR = SDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 ) 00332 VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) 00333 * 00334 CALL SLASY2( .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 50 J = 1, N 00342 CALL SSCAL( M, SCALOC, C( 1, J ), 1 ) 00343 50 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 60 CONTINUE 00353 * 00354 70 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 130 L = 1, N 00375 IF( L.LT.LNEXT ) 00376 $ GO TO 130 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 120 K = 1, M 00397 IF( K.LT.KNEXT ) 00398 $ GO TO 120 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 = SDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 00416 SUMR = SDOT( 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 80 J = 1, N 00436 CALL SSCAL( M, SCALOC, C( 1, J ), 1 ) 00437 80 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 = SDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 00445 SUMR = SDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 00446 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 00447 * 00448 SUML = SDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 00449 SUMR = SDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 00450 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 00451 * 00452 CALL SLALN2( .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 90 J = 1, N 00460 CALL SSCAL( M, SCALOC, C( 1, J ), 1 ) 00461 90 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 = SDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 00470 SUMR = SDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 00471 VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) 00472 * 00473 SUML = SDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 00474 SUMR = SDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 00475 VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) 00476 * 00477 CALL SLALN2( .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 100 J = 1, N 00485 CALL SSCAL( M, SCALOC, C( 1, J ), 1 ) 00486 100 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 = SDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 00495 SUMR = SDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 00496 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 00497 * 00498 SUML = SDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 00499 SUMR = SDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 00500 VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) 00501 * 00502 SUML = SDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 00503 SUMR = SDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 00504 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 00505 * 00506 SUML = SDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 ) 00507 SUMR = SDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 ) 00508 VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) 00509 * 00510 CALL SLASY2( .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 110 J = 1, N 00518 CALL SSCAL( M, SCALOC, C( 1, J ), 1 ) 00519 110 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 120 CONTINUE 00529 130 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 190 L = N, 1, -1 00550 IF( L.GT.LNEXT ) 00551 $ GO TO 190 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 180 K = 1, M 00572 IF( K.LT.KNEXT ) 00573 $ GO TO 180 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 = SDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 00591 SUMR = SDOT( 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 140 J = 1, N 00612 CALL SSCAL( M, SCALOC, C( 1, J ), 1 ) 00613 140 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 = SDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 00621 SUMR = SDOT( 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 = SDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 00626 SUMR = SDOT( 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 SLALN2( .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 150 J = 1, N 00638 CALL SSCAL( M, SCALOC, C( 1, J ), 1 ) 00639 150 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 = SDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 00648 SUMR = SDOT( 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 = SDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 00653 SUMR = SDOT( 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 SLALN2( .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 160 J = 1, N 00665 CALL SSCAL( M, SCALOC, C( 1, J ), 1 ) 00666 160 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 = SDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 00675 SUMR = SDOT( 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 = SDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 00680 SUMR = SDOT( 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 = SDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 00685 SUMR = SDOT( 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 = SDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 ) 00690 SUMR = SDOT( 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 SLASY2( .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 170 J = 1, N 00702 CALL SSCAL( M, SCALOC, C( 1, J ), 1 ) 00703 170 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 180 CONTINUE 00713 190 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 250 L = N, 1, -1 00734 IF( L.GT.LNEXT ) 00735 $ GO TO 250 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 240 K = M, 1, -1 00756 IF( K.GT.KNEXT ) 00757 $ GO TO 240 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 = SDOT( M-K1, A( K1, MIN(K1+1, M ) ), LDA, 00775 $ C( MIN( K1+1, M ), L1 ), 1 ) 00776 SUMR = SDOT( 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 200 J = 1, N 00797 CALL SSCAL( M, SCALOC, C( 1, J ), 1 ) 00798 200 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 = SDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 00806 $ C( MIN( K2+1, M ), L1 ), 1 ) 00807 SUMR = SDOT( 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 = SDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 00812 $ C( MIN( K2+1, M ), L1 ), 1 ) 00813 SUMR = SDOT( 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 SLALN2( .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 210 J = 1, N 00825 CALL SSCAL( M, SCALOC, C( 1, J ), 1 ) 00826 210 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 = SDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 00835 $ C( MIN( K1+1, M ), L1 ), 1 ) 00836 SUMR = SDOT( 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 = SDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 00841 $ C( MIN( K1+1, M ), L2 ), 1 ) 00842 SUMR = SDOT( 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 SLALN2( .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 220 J = 1, N 00854 CALL SSCAL( M, SCALOC, C( 1, J ), 1 ) 00855 220 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 = SDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 00864 $ C( MIN( K2+1, M ), L1 ), 1 ) 00865 SUMR = SDOT( 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 = SDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 00870 $ C( MIN( K2+1, M ), L2 ), 1 ) 00871 SUMR = SDOT( 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 = SDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 00876 $ C( MIN( K2+1, M ), L1 ), 1 ) 00877 SUMR = SDOT( 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 = SDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 00882 $ C( MIN( K2+1, M ), L2 ), 1 ) 00883 SUMR = SDOT( 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 SLASY2( .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 230 J = 1, N 00895 CALL SSCAL( M, SCALOC, C( 1, J ), 1 ) 00896 230 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 240 CONTINUE 00906 250 CONTINUE 00907 * 00908 END IF 00909 * 00910 RETURN 00911 * 00912 * End of STRSYL 00913 * 00914 END
1.7.3 
