stgsyl.f.html

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

</span>      Q = P + 1
      J = 1
   60 CONTINUE
      IF( J.GT.N )
     $   GO TO 70
      Q = Q + 1
      IWORK( Q ) = J
      J = J + NB
      IF( J.GE.N )
     $   GO TO 70
      IF( B( J, J-1 ).NE.ZERO )
     $   J = J + 1
      GO TO 60
   70 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      IWORK( Q+1 ) = N + 1
      IF( IWORK( Q ).EQ.IWORK( Q+1 ) )
     $   Q = Q - 1
<span class="comment">*</span><span class="comment">
</span>      IF( NOTRAN ) THEN
<span class="comment">*</span><span class="comment">
</span>         DO 150 IROUND = 1, ISOLVE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Solve (I, J)-subsystem
</span><span class="comment">*</span><span class="comment">               A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
</span><span class="comment">*</span><span class="comment">               D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
</span><span class="comment">*</span><span class="comment">           for I = P, P - 1,..., 1; J = 1, 2,..., Q
</span><span class="comment">*</span><span class="comment">
</span>            DSCALE = ZERO
            DSUM = ONE
            PQ = 0
            SCALE = ONE
            DO 130 J = P + 2, Q
               JS = IWORK( J )
               JE = IWORK( J+1 ) - 1
               NB = JE - JS + 1
               DO 120 I = P, 1, -1
                  IS = IWORK( I )
                  IE = IWORK( I+1 ) - 1
                  MB = IE - IS + 1
                  PPQQ = 0
                  CALL <a name="STGSY2.411"></a><a href="stgsy2.f.html#STGSY2.1">STGSY2</a>( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA,
     $                         B( JS, JS ), LDB, C( IS, JS ), LDC,
     $                         D( IS, IS ), LDD, E( JS, JS ), LDE,
     $                         F( IS, JS ), LDF, SCALOC, DSUM, DSCALE,
     $                         IWORK( Q+2 ), PPQQ, LINFO )
                  IF( LINFO.GT.0 )
     $               INFO = LINFO
<span class="comment">*</span><span class="comment">
</span>                  PQ = PQ + PPQQ
                  IF( SCALOC.NE.ONE ) THEN
                     DO 80 K = 1, JS - 1
                        CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
                        CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
   80                CONTINUE
                     DO 90 K = JS, JE
                        CALL SSCAL( IS-1, SCALOC, C( 1, K ), 1 )
                        CALL SSCAL( IS-1, SCALOC, F( 1, K ), 1 )
   90                CONTINUE
                     DO 100 K = JS, JE
                        CALL SSCAL( M-IE, SCALOC, C( IE+1, K ), 1 )
                        CALL SSCAL( M-IE, SCALOC, F( IE+1, K ), 1 )
  100                CONTINUE
                     DO 110 K = JE + 1, N
                        CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
                        CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
  110                CONTINUE
                     SCALE = SCALE*SCALOC
                  END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Substitute R(I, J) and L(I, J) into remaining
</span><span class="comment">*</span><span class="comment">                 equation.
</span><span class="comment">*</span><span class="comment">
</span>                  IF( I.GT.1 ) THEN
                     CALL SGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, IS-1, NB, MB, -ONE,
     $                           A( 1, IS ), LDA, C( IS, JS ), LDC, ONE,
     $                           C( 1, JS ), LDC )
                     CALL SGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, IS-1, NB, MB, -ONE,
     $                           D( 1, IS ), LDD, C( IS, JS ), LDC, ONE,
     $                           F( 1, JS ), LDF )
                  END IF
                  IF( J.LT.Q ) THEN
                     CALL SGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, MB, N-JE, NB, ONE,
     $                           F( IS, JS ), LDF, B( JS, JE+1 ), LDB,
     $                           ONE, C( IS, JE+1 ), LDC )
                     CALL SGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, MB, N-JE, NB, ONE,
     $                           F( IS, JS ), LDF, E( JS, JE+1 ), LDE,
     $                           ONE, F( IS, JE+1 ), LDF )
                  END IF
  120          CONTINUE
  130       CONTINUE
            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="SLACPY.473"></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.474"></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.475"></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.476"></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.478"></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.479"></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
  150    CONTINUE
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Solve transposed (I, J)-subsystem
</span><span class="comment">*</span><span class="comment">             A(I, I)' * R(I, J)  + D(I, I)' * L(I, J)  =  C(I, J)
</span><span class="comment">*</span><span class="comment">             R(I, J)  * B(J, J)' + L(I, J)  * E(J, J)' = -F(I, J)
</span><span class="comment">*</span><span class="comment">        for I = 1,2,..., P; J = Q, Q-1,..., 1
</span><span class="comment">*</span><span class="comment">
</span>         SCALE = ONE
         DO 210 I = 1, P
            IS = IWORK( I )
            IE = IWORK( I+1 ) - 1
            MB = IE - IS + 1
            DO 200 J = Q, P + 2, -1
               JS = IWORK( J )
               JE = IWORK( J+1 ) - 1
               NB = JE - JS + 1
               CALL <a name="STGSY2.500"></a><a href="stgsy2.f.html#STGSY2.1">STGSY2</a>( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA,
     $                      B( JS, JS ), LDB, C( IS, JS ), LDC,
     $                      D( IS, IS ), LDD, E( JS, JS ), LDE,
     $                      F( IS, JS ), LDF, SCALOC, DSUM, DSCALE,
     $                      IWORK( Q+2 ), PPQQ, LINFO )
               IF( LINFO.GT.0 )
     $            INFO = LINFO
               IF( SCALOC.NE.ONE ) THEN
                  DO 160 K = 1, JS - 1
                     CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
                     CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
  160             CONTINUE
                  DO 170 K = JS, JE
                     CALL SSCAL( IS-1, SCALOC, C( 1, K ), 1 )
                     CALL SSCAL( IS-1, SCALOC, F( 1, K ), 1 )
  170             CONTINUE
                  DO 180 K = JS, JE
                     CALL SSCAL( M-IE, SCALOC, C( IE+1, K ), 1 )
                     CALL SSCAL( M-IE, SCALOC, F( IE+1, K ), 1 )
  180             CONTINUE
                  DO 190 K = JE + 1, N
                     CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
                     CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
  190             CONTINUE
                  SCALE = SCALE*SCALOC
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Substitute R(I, J) and L(I, J) into remaining equation.
</span><span class="comment">*</span><span class="comment">
</span>               IF( J.GT.P+2 ) THEN
                  CALL SGEMM( <span class="string">'N'</span>, <span class="string">'T'</span>, MB, JS-1, NB, ONE, C( IS, JS ),
     $                        LDC, B( 1, JS ), LDB, ONE, F( IS, 1 ),
     $                        LDF )
                  CALL SGEMM( <span class="string">'N'</span>, <span class="string">'T'</span>, MB, JS-1, NB, ONE, F( IS, JS ),
     $                        LDF, E( 1, JS ), LDE, ONE, F( IS, 1 ),
     $                        LDF )
               END IF
               IF( I.LT.P ) THEN
                  CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, M-IE, NB, MB, -ONE,
     $                        A( IS, IE+1 ), LDA, C( IS, JS ), LDC, ONE,
     $                        C( IE+1, JS ), LDC )
                  CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, M-IE, NB, MB, -ONE,
     $                        D( IS, IE+1 ), LDD, F( IS, JS ), LDF, ONE,
     $                        C( IE+1, JS ), LDC )
               END IF
  200       CONTINUE
  210    CONTINUE
<span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = LWMIN
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="STGSYL.554"></a><a href="stgsyl.f.html#STGSYL.1">STGSYL</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>