ctgsyl.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 597 行 · 第 1/3 页

</span><span class="comment">*</span><span class="comment">      Note 75.  To appear in ACM Trans. on Math. Software, Vol 22,
</span><span class="comment">*</span><span class="comment">      No 1, 1996.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  [2] B. Kagstrom, A Perturbation Analysis of the Generalized Sylvester
</span><span class="comment">*</span><span class="comment">      Equation (AR - LB, DR - LE ) = (C, F), SIAM J. Matrix Anal.
</span><span class="comment">*</span><span class="comment">      Appl., 15(4):1045-1060, 1994.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  [3] B. Kagstrom and L. Westin, Generalized Schur Methods with
</span><span class="comment">*</span><span class="comment">      Condition Estimators for Solving the Generalized Sylvester
</span><span class="comment">*</span><span class="comment">      Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7,
</span><span class="comment">*</span><span class="comment">      July 1989, pp 745-751.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">  Replaced various illegal calls to CCOPY by calls to <a name="CLASET.194"></a><a href="claset.f.html#CLASET.1">CLASET</a>.
</span><span class="comment">*</span><span class="comment">  Sven Hammarling, 1/5/02.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               ZERO, ONE
      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
      COMPLEX            CZERO
      PARAMETER          ( CZERO = (0.0E+0, 0.0E+0) )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            LQUERY, NOTRAN
      INTEGER            I, IE, IFUNC, IROUND, IS, ISOLVE, J, JE, JS, K,
     $                   LINFO, LWMIN, MB, NB, P, PQ, Q
      REAL               DSCALE, DSUM, SCALE2, SCALOC
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.210"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            <a name="ILAENV.211"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
      EXTERNAL           <a name="LSAME.212"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="ILAENV.212"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           CGEMM, <a name="CLACPY.215"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>, <a name="CLASET.215"></a><a href="claset.f.html#CLASET.1">CLASET</a>, CSCAL, <a name="CTGSY2.215"></a><a href="ctgsy2.f.html#CTGSY2.1">CTGSY2</a>, <a name="XERBLA.215"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          CMPLX, MAX, REAL, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Decode and test input parameters
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      NOTRAN = <a name="LSAME.225"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'N'</span> )
      LQUERY = ( LWORK.EQ.-1 )
<span class="comment">*</span><span class="comment">
</span>      IF( .NOT.NOTRAN .AND. .NOT.<a name="LSAME.228"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'C'</span> ) ) THEN
         INFO = -1
      ELSE IF( NOTRAN ) THEN
         IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.4 ) ) THEN
            INFO = -2
         END IF
      END IF
      IF( INFO.EQ.0 ) THEN
         IF( M.LE.0 ) THEN
            INFO = -3
         ELSE IF( N.LE.0 ) THEN
            INFO = -4
         ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
            INFO = -6
         ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
            INFO = -8
         ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
            INFO = -10
         ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
            INFO = -12
         ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
            INFO = -14
         ELSE IF( LDF.LT.MAX( 1, M ) ) THEN
            INFO = -16
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.EQ.0 ) THEN
         IF( NOTRAN ) THEN
            IF( IJOB.EQ.1 .OR. IJOB.EQ.2 ) THEN
               LWMIN = MAX( 1, 2*M*N )
            ELSE
               LWMIN = 1
            END IF
         ELSE
            LWMIN = 1
         END IF
         WORK( 1 ) = LWMIN
<span class="comment">*</span><span class="comment">
</span>         IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
            INFO = -20
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.273"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CTGSYL.273"></a><a href="ctgsyl.f.html#CTGSYL.1">CTGSYL</a>'</span>, -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( M.EQ.0 .OR. N.EQ.0 ) THEN
         SCALE = 1
         IF( NOTRAN ) THEN
            IF( IJOB.NE.0 ) THEN
               DIF = 0
            END IF
         END IF
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine  optimal block sizes MB and NB
</span><span class="comment">*</span><span class="comment">
</span>      MB = <a name="ILAENV.293"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 2, <span class="string">'<a name="CTGSYL.293"></a><a href="ctgsyl.f.html#CTGSYL.1">CTGSYL</a>'</span>, TRANS, M, N, -1, -1 )
      NB = <a name="ILAENV.294"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 5, <span class="string">'<a name="CTGSYL.294"></a><a href="ctgsyl.f.html#CTGSYL.1">CTGSYL</a>'</span>, TRANS, M, N, -1, -1 )
<span class="comment">*</span><span class="comment">
</span>      ISOLVE = 1
      IFUNC = 0
      IF( NOTRAN ) THEN
         IF( IJOB.GE.3 ) THEN
            IFUNC = IJOB - 2
            CALL <a name="CLASET.301"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'F'</span>, M, N, CZERO, CZERO, C, LDC )
            CALL <a name="CLASET.302"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'F'</span>, M, N, CZERO, CZERO, F, LDF )
         ELSE IF( IJOB.GE.1 .AND. NOTRAN ) THEN
            ISOLVE = 2
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( ( MB.LE.1 .AND. NB.LE.1 ) .OR. ( MB.GE.M .AND. NB.GE.N ) )
     $     THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use unblocked Level 2 solver
</span><span class="comment">*</span><span class="comment">
</span>         DO 30 IROUND = 1, ISOLVE
<span class="comment">*</span><span class="comment">
</span>            SCALE = ONE
            DSCALE = ZERO
            DSUM = ONE
            PQ = M*N
            CALL <a name="CTGSY2.319"></a><a href="ctgsy2.f.html#CTGSY2.1">CTGSY2</a>( TRANS, IFUNC, M, N, A, LDA, B, LDB, C, LDC, D,
     $                   LDD, E, LDE, F, LDF, SCALE, DSUM, DSCALE,
     $                   INFO )
            IF( DSCALE.NE.ZERO ) THEN
               IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN
                  DIF = SQRT( REAL( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
               ELSE
                  DIF = SQRT( REAL( PQ ) ) / ( DSCALE*SQRT( DSUM ) )
               END IF
            END IF
            IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN
               IF( NOTRAN ) THEN
                  IFUNC = IJOB
               END IF
               SCALE2 = SCALE
               CALL <a name="CLACPY.334"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'F'</span>, M, N, C, LDC, WORK, M )
               CALL <a name="CLACPY.335"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'F'</span>, M, N, F, LDF, WORK( M*N+1 ), M )
               CALL <a name="CLASET.336"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'F'</span>, M, N, CZERO, CZERO, C, LDC )
               CALL <a name="CLASET.337"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'F'</span>, M, N, CZERO, CZERO, F, LDF )
            ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
               CALL <a name="CLACPY.339"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'F'</span>, M, N, WORK, M, C, LDC )
               CALL <a name="CLACPY.340"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'F'</span>, M, N, WORK( M*N+1 ), M, F, LDF )
               SCALE = SCALE2
            END IF
   30    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         RETURN
<span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine block structure of A
</span><span class="comment">*</span><span class="comment">
</span>      P = 0
      I = 1
   40 CONTINUE
      IF( I.GT.M )
     $   GO TO 50
      P = P + 1
      IWORK( P ) = I
      I = I + MB
      IF( I.GE.M )
     $   GO TO 50
      GO TO 40
   50 CONTINUE
      IWORK( P+1 ) = M + 1
      IF( IWORK( P ).EQ.IWORK( P+1 ) )
     $   P = P - 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine block structure of B
</span><span class="comment">*</span><span class="comment">
</span>      Q = P + 1
      J = 1
   60 CONTINUE
      IF( J.GT.N )
     $   GO TO 70
<span class="comment">*</span><span class="comment">
</span>      Q = Q + 1
      IWORK( Q ) = J
      J = J + NB
      IF( J.GE.N )
     $   GO TO 70