LAPACK 3.3.0

strsyl.f

Go to the documentation of this file.
00001       SUBROUTINE STRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
00002      $                   LDC, SCALE, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
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' *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)'*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)'*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'*X + ISGN*X*B' = 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)'*X(K,L) + ISGN*X(K,L)*B(L,L)' = C(K,L) - R(K,L)
00539 *
00540 *        Where
00541 *                     K-1                          N
00542 *            R(K,L) = SUM [A(I,K)'*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)'].
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' = 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)' = 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)'].
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
 All Files Functions