stgsyl.f.html

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

</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  [1] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
</span><span class="comment">*</span><span class="comment">      for Solving the Generalized Sylvester Equation and Estimating the
</span><span class="comment">*</span><span class="comment">      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
</span><span class="comment">*</span><span class="comment">      Department of Computing Science, Umea University, S-901 87 Umea,
</span><span class="comment">*</span><span class="comment">      Sweden, December 1993, Revised April 1994, Also as LAPACK Working
</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 SCOPY by calls to <a name="SLASET.195"></a><a href="slaset.f.html#SLASET.1">SLASET</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 )
<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, PPQQ, 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.209"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            <a name="ILAENV.210"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
      EXTERNAL           <a name="LSAME.211"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="ILAENV.211"></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           SGEMM, <a name="SLACPY.214"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>, <a name="SLASET.214"></a><a href="slaset.f.html#SLASET.1">SLASET</a>, SSCAL, <a name="STGSY2.214"></a><a href="stgsy2.f.html#STGSY2.1">STGSY2</a>, <a name="XERBLA.214"></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          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.224"></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.227"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'T'</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.272"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="STGSYL.272"></a><a href="stgsyl.f.html#STGSYL.1">STGSYL</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.292"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 2, <span class="string">'<a name="STGSYL.292"></a><a href="stgsyl.f.html#STGSYL.1">STGSYL</a>'</span>, TRANS, M, N, -1, -1 )
      NB = <a name="ILAENV.293"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 5, <span class="string">'<a name="STGSYL.293"></a><a href="stgsyl.f.html#STGSYL.1">STGSYL</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="SLASET.300"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'F'</span>, M, N, ZERO, ZERO, C, LDC )
            CALL <a name="SLASET.301"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'F'</span>, M, N, ZERO, ZERO, 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>         DO 30 IROUND = 1, ISOLVE
<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>            DSCALE = ZERO
            DSUM = ONE
            PQ = 0
            CALL <a name="STGSY2.317"></a><a href="stgsy2.f.html#STGSY2.1">STGSY2</a>( TRANS, IFUNC, M, N, A, LDA, B, LDB, C, LDC, D,
     $                   LDD, E, LDE, F, LDF, SCALE, DSUM, DSCALE,
     $                   IWORK, PQ, 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
<span class="comment">*</span><span class="comment">
</span>            IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN
               IF( NOTRAN ) THEN
                  IFUNC = IJOB
               END IF
               SCALE2 = SCALE
               CALL <a name="SLACPY.333"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'F'</span>, M, N, C, LDC, WORK, M )
               CALL <a name="SLACPY.334"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'F'</span>, M, N, F, LDF, WORK( M*N+1 ), M )
               CALL <a name="SLASET.335"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'F'</span>, M, N, ZERO, ZERO, C, LDC )
               CALL <a name="SLASET.336"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'F'</span>, M, N, ZERO, ZERO, F, LDF )
            ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
               CALL <a name="SLACPY.338"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'F'</span>, M, N, WORK, M, C, LDC )
               CALL <a name="SLACPY.339"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</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
      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
      IF( A( I, I-1 ).NE.ZERO )
     $   I = I + 1
      GO TO 40
   50 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      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">