dtrsyl
C DTRSYL SOURCE FANDEUR 22/05/02 21:15:18 11359 *> \brief \b DTRSYL * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download DTRSYL + dependencies *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsyl.f"> *> [TGZ]</a> *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsyl.f"> *> [ZIP]</a> *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsyl.f"> *> [TXT]</a> *> \endhtmlonly * * Definition: * =========== * * SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, * LDC, SCALE, INFO ) * * .. Scalar Arguments .. * CHARACTER TRANA, TRANB * INTEGER INFO, ISGN, LDA, LDB, LDC, M, N * REAL*8 SCALE * .. * .. Array Arguments .. * REAL*8 A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> DTRSYL solves the real Sylvester matrix equation: *> *> op(A)*X + X*op(B) = scale*C or *> op(A)*X - X*op(B) = scale*C, *> *> where op(A) = A or A**T, and A and B are both upper quasi- *> triangular. A is M-by-M and B is N-by-N; the right hand side C and *> the solution X are M-by-N; and scale is an output scale factor, set *> <= 1 to avoid overflow in X. *> *> A and B must be in Schur canonical form (as returned by DHSEQR), that *> is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; *> each 2-by-2 diagonal block has its diagonal elements equal and its *> off-diagonal elements of opposite sign. *> \endverbatim * * Arguments: * ========== * *> \param[in] TRANA *> \verbatim *> TRANA is CHARACTER*1 *> Specifies the option op(A): *> = 'N': op(A) = A (No transpose) *> = 'T': op(A) = A**T (Transpose) *> = 'C': op(A) = A**H (Conjugate transpose = Transpose) *> \endverbatim *> *> \param[in] TRANB *> \verbatim *> TRANB is CHARACTER*1 *> Specifies the option op(B): *> = 'N': op(B) = B (No transpose) *> = 'T': op(B) = B**T (Transpose) *> = 'C': op(B) = B**H (Conjugate transpose = Transpose) *> \endverbatim *> *> \param[in] ISGN *> \verbatim *> ISGN is INTEGER *> Specifies the sign in the equation: *> = +1: solve op(A)*X + X*op(B) = scale*C *> = -1: solve op(A)*X - X*op(B) = scale*C *> \endverbatim *> *> \param[in] M *> \verbatim *> M is INTEGER *> The order of the matrix A, and the number of rows in the *> matrices X and C. M >= 0. *> \endverbatim *> *> \param[in] N *> \verbatim *> N is INTEGER *> The order of the matrix B, and the number of columns in the *> matrices X and C. N >= 0. *> \endverbatim *> *> \param[in] A *> \verbatim *> A is DOUBLE PRECISION array, dimension (LDA,M) *> The upper quasi-triangular matrix A, in Schur canonical form. *> \endverbatim *> *> \param[in] LDA *> \verbatim *> LDA is INTEGER *> The leading dimension of the array A. LDA >= max(1,M). *> \endverbatim *> *> \param[in] B *> \verbatim *> B is DOUBLE PRECISION array, dimension (LDB,N) *> The upper quasi-triangular matrix B, in Schur canonical form. *> \endverbatim *> *> \param[in] LDB *> \verbatim *> LDB is INTEGER *> The leading dimension of the array B. LDB >= max(1,N). *> \endverbatim *> *> \param[in,out] C *> \verbatim *> C is DOUBLE PRECISION array, dimension (LDC,N) *> On entry, the M-by-N right hand side matrix C. *> On exit, C is overwritten by the solution matrix X. *> \endverbatim *> *> \param[in] LDC *> \verbatim *> LDC is INTEGER *> The leading dimension of the array C. LDC >= max(1,M) *> \endverbatim *> *> \param[out] SCALE *> \verbatim *> SCALE is DOUBLE PRECISION *> The scale factor, scale, set <= 1 to avoid overflow in X. *> \endverbatim *> *> \param[out] INFO *> \verbatim *> INFO is INTEGER *> = 0: successful exit *> < 0: if INFO = -i, the i-th argument had an illegal value *> = 1: A and B have common or very close eigenvalues; perturbed *> values were used to solve the equation (but the matrices *> A and B are unchanged). *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date December 2016 * *> \ingroup doubleSYcomputational * * ===================================================================== $ LDC, SCALE, INFO ) * * -- LAPACK computational routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * * .. Scalar Arguments .. CHARACTER TRANA, TRANB INTEGER INFO, ISGN, LDA, LDB, LDC, M, N REAL*8 SCALE * .. * .. Array Arguments .. REAL*8 A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * ===================================================================== * * .. Parameters .. * .. * .. Local Scalars .. LOGICAL NOTRNA, NOTRNB INTEGER IERR, J, K, K1, K2, KNEXT, L, L1, L2, LNEXT REAL*8 A11, BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN, $ SMLNUM, SUML, SUMR, XNORM * .. * .. Local Arrays .. REAL*8 DUM( 1 ), VEC( 2, 2 ), X( 2, 2 ) * .. * .. External Functions .. LOGICAL LSAME * .. * .. External Subroutines .. * .. ** .. Intrinsic Functions .. * INTRINSIC ABS, DBLE, MAX, MIN ** .. ** .. Executable Statements .. * * Decode and Test input parameters * * INFO = 0 IF( .NOT.NOTRNA .AND. INFO = -1 ELSE IF( .NOT.NOTRNB .AND. INFO = -2 ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN INFO = -3 ELSE IF( M.LT.0 ) THEN INFO = -4 ELSE IF( N.LT.0 ) THEN INFO = -5 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -7 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN INFO = -9 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN INFO = -11 END IF IF( INFO.NE.0 ) THEN RETURN END IF * * Quick return if possible * SCALE = ONE IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * * Set constants to control overflow * BIGNUM = ONE / SMLNUM SMLNUM = SMLNUM*DBLE( M*N ) / EPS BIGNUM = ONE / SMLNUM * * SGN = ISGN * IF( NOTRNA .AND. NOTRNB ) THEN * * Solve A*X + ISGN*X*B = scale*C. * * The (K,L)th block of X is determined starting from * bottom-left corner column by column by * * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) * * Where * M L-1 * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)]. * I=K+1 J=1 * * Start column loop (index = L) * L1 (L2) : column index of the first (first) row of X(K,L). * LNEXT = 1 DO 60 L = 1, N IF( L.LT.LNEXT ) $ GO TO 60 IF( L.EQ.N ) THEN L1 = L L2 = L ELSE L1 = L L2 = L + 1 LNEXT = L + 2 ELSE L1 = L L2 = L LNEXT = L + 1 END IF END IF * * Start row loop (index = K) * K1 (K2): row index of the first (last) row of X(K,L). * KNEXT = M DO 50 K = M, 1, -1 IF( K.GT.KNEXT ) $ GO TO 50 IF( K.EQ.1 ) THEN K1 = K K2 = K ELSE K1 = K - 1 K2 = K KNEXT = K - 2 ELSE K1 = K K2 = K KNEXT = K - 1 END IF END IF * IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN $ LDA, C( MIN( K1+1, M ), L1 ), 1 ) $ LDC, B( 1, L1 ), 1 ) VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) SCALOC = ONE * A11 = A( K1, K1 ) + SGN*B( L1, L1 ) DA11 = ABS( A11 ) IF( DA11.LE.SMIN ) THEN A11 = SMIN DA11 = SMIN INFO = 1 END IF DB = ABS( VEC( 1, 1 ) ) IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN IF( DB.GT.BIGNUM*DA11 ) $ SCALOC = ONE / DB END IF X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 * IF( SCALOC.NE.ONE ) THEN DO 10 J = 1, N 10 CONTINUE SCALE = SCALE*SCALOC END IF C( K1, L1 ) = X( 1, 1 ) * ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN * $ LDA, C( MIN( K2+1, M ), L1 ), 1 ) $ LDC, B( 1, L1 ), 1 ) VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) * $ LDA, C( MIN( K2+1, M ), L1 ), 1 ) $ LDC, B( 1, L1 ), 1 ) VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) * $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), IF( IERR.NE.0 ) $ INFO = 1 * IF( SCALOC.NE.ONE ) THEN DO 20 J = 1, N 20 CONTINUE SCALE = SCALE*SCALOC END IF C( K1, L1 ) = X( 1, 1 ) C( K2, L1 ) = X( 2, 1 ) * ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN * $ LDA, C( MIN( K1+1, M ), L1 ), 1 ) $ B( 1, L1 ), 1 ) VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) * $ LDA, C( MIN( K1+1, M ), L2 ), 1 ) $ B( 1, L2 ), 1 ) VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) * $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), IF( IERR.NE.0 ) $ INFO = 1 * IF( SCALOC.NE.ONE ) THEN DO 30 J = 1, N 30 CONTINUE SCALE = SCALE*SCALOC END IF C( K1, L1 ) = X( 1, 1 ) C( K1, L2 ) = X( 2, 1 ) * ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN * $ LDA, C( MIN( K2+1, M ), L1 ), 1 ) $ B( 1, L1 ), 1 ) VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) * $ LDA, C( MIN( K2+1, M ), L2 ), 1 ) $ B( 1, L2 ), 1 ) VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) * $ LDA, C( MIN( K2+1, M ), L1 ), 1 ) $ B( 1, L1 ), 1 ) VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) * $ LDA, C( MIN( K2+1, M ), L2 ), 1 ) $ B( 1, L2 ), 1 ) VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) * $ A( K1, K1 ), LDA, B( L1, L1 ), LDB, VEC, $ 2, SCALOC, X, 2, XNORM, IERR ) IF( IERR.NE.0 ) $ INFO = 1 * IF( SCALOC.NE.ONE ) THEN DO 40 J = 1, N 40 CONTINUE SCALE = SCALE*SCALOC END IF C( K1, L1 ) = X( 1, 1 ) C( K1, L2 ) = X( 1, 2 ) C( K2, L1 ) = X( 2, 1 ) C( K2, L2 ) = X( 2, 2 ) END IF * 50 CONTINUE * 60 CONTINUE * ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN * * Solve A**T *X + ISGN*X*B = scale*C. * * The (K,L)th block of X is determined starting from * upper-left corner column by column by * * A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) * * Where * K-1 T L-1 * R(K,L) = SUM [A(I,K)**T*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)] * I=1 J=1 * * Start column loop (index = L) * L1 (L2): column index of the first (last) row of X(K,L) * LNEXT = 1 DO 120 L = 1, N IF( L.LT.LNEXT ) $ GO TO 120 IF( L.EQ.N ) THEN L1 = L L2 = L ELSE L1 = L L2 = L + 1 LNEXT = L + 2 ELSE L1 = L L2 = L LNEXT = L + 1 END IF END IF * * Start row loop (index = K) * K1 (K2): row index of the first (last) row of X(K,L) * KNEXT = 1 DO 110 K = 1, M IF( K.LT.KNEXT ) $ GO TO 110 IF( K.EQ.M ) THEN K1 = K K2 = K ELSE K1 = K K2 = K + 1 KNEXT = K + 2 ELSE K1 = K K2 = K KNEXT = K + 1 END IF END IF * IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN $ C( 1, L1 ), 1 ) $ B( 1, L1 ), 1 ) VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) SCALOC = ONE * A11 = A( K1, K1 ) + SGN*B( L1, L1 ) DA11 = ABS( A11 ) IF( DA11.LE.SMIN ) THEN A11 = SMIN DA11 = SMIN INFO = 1 END IF DB = ABS( VEC( 1, 1 ) ) IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN IF( DB.GT.BIGNUM*DA11 ) $ SCALOC = ONE / DB END IF X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 * IF( SCALOC.NE.ONE ) THEN DO 70 J = 1, N 70 CONTINUE SCALE = SCALE*SCALOC END IF C( K1, L1 ) = X( 1, 1 ) * ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN * $ C( 1, L1 ), 1 ) $ B( 1, L1 ), 1 ) VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) * $ C( 1, L1 ), 1 ) $ B( 1, L1 ), 1 ) VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) * $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), IF( IERR.NE.0 ) $ INFO = 1 * IF( SCALOC.NE.ONE ) THEN DO 80 J = 1, N 80 CONTINUE SCALE = SCALE*SCALOC END IF C( K1, L1 ) = X( 1, 1 ) C( K2, L1 ) = X( 2, 1 ) * ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN * $ C( 1, L1 ), 1 ) $ B( 1, L1 ), 1 ) VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) * $ C( 1, L2 ), 1 ) $ B( 1, L2 ), 1 ) VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) * $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), IF( IERR.NE.0 ) $ INFO = 1 * IF( SCALOC.NE.ONE ) THEN DO 90 J = 1, N 90 CONTINUE SCALE = SCALE*SCALOC END IF C( K1, L1 ) = X( 1, 1 ) C( K1, L2 ) = X( 2, 1 ) * ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN * $ C( 1, L1 ), 1 ) $ B( 1, L1 ), 1 ) VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) * $ C( 1, L2 ), 1 ) $ B( 1, L2 ), 1 ) VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) * $ C( 1, L1 ), 1 ) $ B( 1, L1 ), 1 ) VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) * $ C( 1, L2 ), 1 ) $ B( 1, L2 ), 1 ) VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) * $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X, $ 2, XNORM, IERR ) IF( IERR.NE.0 ) $ INFO = 1 * IF( SCALOC.NE.ONE ) THEN DO 100 J = 1, N 100 CONTINUE SCALE = SCALE*SCALOC END IF C( K1, L1 ) = X( 1, 1 ) C( K1, L2 ) = X( 1, 2 ) C( K2, L1 ) = X( 2, 1 ) C( K2, L2 ) = X( 2, 2 ) END IF * 110 CONTINUE 120 CONTINUE * ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN * * Solve A**T*X + ISGN*X*B**T = scale*C. * * The (K,L)th block of X is determined starting from * top-right corner column by column by * * A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L) * * Where * K-1 N * R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T]. * I=1 J=L+1 * * Start column loop (index = L) * L1 (L2): column index of the first (last) row of X(K,L) * LNEXT = N DO 180 L = N, 1, -1 IF( L.GT.LNEXT ) $ GO TO 180 IF( L.EQ.1 ) THEN L1 = L L2 = L ELSE L1 = L - 1 L2 = L LNEXT = L - 2 ELSE L1 = L L2 = L LNEXT = L - 1 END IF END IF * * Start row loop (index = K) * K1 (K2): row index of the first (last) row of X(K,L) * KNEXT = 1 DO 170 K = 1, M IF( K.LT.KNEXT ) $ GO TO 170 IF( K.EQ.M ) THEN K1 = K K2 = K ELSE K1 = K K2 = K + 1 KNEXT = K + 2 ELSE K1 = K K2 = K KNEXT = K + 1 END IF END IF * IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN $ C( 1, L1 ), 1 ) $ LDC, B( L1, MIN( L1+1, N ) ), LDB ) VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) SCALOC = ONE * A11 = A( K1, K1 ) + SGN*B( L1, L1 ) DA11 = ABS( A11 ) IF( DA11.LE.SMIN ) THEN A11 = SMIN DA11 = SMIN INFO = 1 END IF DB = ABS( VEC( 1, 1 ) ) IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN IF( DB.GT.BIGNUM*DA11 ) $ SCALOC = ONE / DB END IF X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 * IF( SCALOC.NE.ONE ) THEN DO 130 J = 1, N 130 CONTINUE SCALE = SCALE*SCALOC END IF C( K1, L1 ) = X( 1, 1 ) * ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN * $ C( 1, L1 ), 1 ) $ B( L1, MIN( L2+1, N ) ), LDB ) VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) * $ B( L1, MIN( L2+1, N ) ), LDB ) VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) * $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), IF( IERR.NE.0 ) $ INFO = 1 * IF( SCALOC.NE.ONE ) THEN DO 140 J = 1, N 140 CONTINUE SCALE = SCALE*SCALOC END IF C( K1, L1 ) = X( 1, 1 ) C( K2, L1 ) = X( 2, 1 ) * ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN * $ B( L1, MIN( L2+1, N ) ), LDB ) VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) * $ B( L2, MIN( L2+1, N ) ), LDB ) VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) * $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), IF( IERR.NE.0 ) $ INFO = 1 * IF( SCALOC.NE.ONE ) THEN DO 150 J = 1, N 150 CONTINUE SCALE = SCALE*SCALOC END IF C( K1, L1 ) = X( 1, 1 ) C( K1, L2 ) = X( 2, 1 ) * ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN * $ B( L1, MIN( L2+1, N ) ), LDB ) VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) * $ B( L2, MIN( L2+1, N ) ), LDB ) VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) * $ B( L1, MIN( L2+1, N ) ), LDB ) VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) * $ B( L2, MIN( L2+1, N ) ), LDB ) VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) * $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X, $ 2, XNORM, IERR ) IF( IERR.NE.0 ) $ INFO = 1 * IF( SCALOC.NE.ONE ) THEN DO 160 J = 1, N 160 CONTINUE SCALE = SCALE*SCALOC END IF C( K1, L1 ) = X( 1, 1 ) C( K1, L2 ) = X( 1, 2 ) C( K2, L1 ) = X( 2, 1 ) C( K2, L2 ) = X( 2, 2 ) END IF * 170 CONTINUE 180 CONTINUE * ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN * * Solve A*X + ISGN*X*B**T = scale*C. * * The (K,L)th block of X is determined starting from * bottom-right corner column by column by * * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L) * * Where * M N * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T]. * I=K+1 J=L+1 * * Start column loop (index = L) * L1 (L2): column index of the first (last) row of X(K,L) * LNEXT = N DO 240 L = N, 1, -1 IF( L.GT.LNEXT ) $ GO TO 240 IF( L.EQ.1 ) THEN L1 = L L2 = L ELSE L1 = L - 1 L2 = L LNEXT = L - 2 ELSE L1 = L L2 = L LNEXT = L - 1 END IF END IF * * Start row loop (index = K) * K1 (K2): row index of the first (last) row of X(K,L) * KNEXT = M DO 230 K = M, 1, -1 IF( K.GT.KNEXT ) $ GO TO 230 IF( K.EQ.1 ) THEN K1 = K K2 = K ELSE K1 = K - 1 K2 = K KNEXT = K - 2 ELSE K1 = K K2 = K KNEXT = K - 1 END IF END IF * IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN $ C( MIN( K1+1, M ), L1 ), 1 ) $ B( L1, MIN( L1+1, N ) ), LDB ) VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) SCALOC = ONE * A11 = A( K1, K1 ) + SGN*B( L1, L1 ) DA11 = ABS( A11 ) IF( DA11.LE.SMIN ) THEN A11 = SMIN DA11 = SMIN INFO = 1 END IF DB = ABS( VEC( 1, 1 ) ) IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN IF( DB.GT.BIGNUM*DA11 ) $ SCALOC = ONE / DB END IF X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 * IF( SCALOC.NE.ONE ) THEN DO 190 J = 1, N 190 CONTINUE SCALE = SCALE*SCALOC END IF C( K1, L1 ) = X( 1, 1 ) * ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN * $ C( MIN( K2+1, M ), L1 ), 1 ) $ B( L1, MIN( L2+1, N ) ), LDB ) VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) * $ C( MIN( K2+1, M ), L1 ), 1 ) $ B( L1, MIN( L2+1, N ) ), LDB ) VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) * $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), IF( IERR.NE.0 ) $ INFO = 1 * IF( SCALOC.NE.ONE ) THEN DO 200 J = 1, N 200 CONTINUE SCALE = SCALE*SCALOC END IF C( K1, L1 ) = X( 1, 1 ) C( K2, L1 ) = X( 2, 1 ) * ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN * $ C( MIN( K1+1, M ), L1 ), 1 ) $ B( L1, MIN( L2+1, N ) ), LDB ) VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) * $ C( MIN( K1+1, M ), L2 ), 1 ) $ B( L2, MIN( L2+1, N ) ), LDB ) VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) * $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), IF( IERR.NE.0 ) $ INFO = 1 * IF( SCALOC.NE.ONE ) THEN DO 210 J = 1, N 210 CONTINUE SCALE = SCALE*SCALOC END IF C( K1, L1 ) = X( 1, 1 ) C( K1, L2 ) = X( 2, 1 ) * ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN * $ C( MIN( K2+1, M ), L1 ), 1 ) $ B( L1, MIN( L2+1, N ) ), LDB ) VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) * $ C( MIN( K2+1, M ), L2 ), 1 ) $ B( L2, MIN( L2+1, N ) ), LDB ) VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) * $ C( MIN( K2+1, M ), L1 ), 1 ) $ B( L1, MIN( L2+1, N ) ), LDB ) VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) * $ C( MIN( K2+1, M ), L2 ), 1 ) $ B( L2, MIN( L2+1, N ) ), LDB ) VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) * $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X, $ 2, XNORM, IERR ) IF( IERR.NE.0 ) $ INFO = 1 * IF( SCALOC.NE.ONE ) THEN DO 220 J = 1, N 220 CONTINUE SCALE = SCALE*SCALOC END IF C( K1, L1 ) = X( 1, 1 ) C( K1, L2 ) = X( 1, 2 ) C( K2, L1 ) = X( 2, 1 ) C( K2, L2 ) = X( 2, 2 ) END IF * 230 CONTINUE 240 CONTINUE * END IF * RETURN * * End of DTRSYL * END
© Cast3M 2003 - Tous droits réservés.
Mentions légales