stgsy2.f.html

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

</span>                  CALL <a name="SGETC2.573"></a><a href="sgetc2.f.html#SGETC2.1">SGETC2</a>( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
                  IF( IERR.GT.0 )
     $               INFO = IERR
                  IF( IJOB.EQ.0 ) THEN
                     CALL <a name="SGESC2.577"></a><a href="sgesc2.f.html#SGESC2.1">SGESC2</a>( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
     $                            SCALOC )
                     IF( SCALOC.NE.ONE ) THEN
                        DO 90 K = 1, N
                           CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
                           CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
   90                   CONTINUE
                        SCALE = SCALE*SCALOC
                     END IF
                  ELSE
                     CALL <a name="SLATDF.587"></a><a href="slatdf.f.html#SLATDF.1">SLATDF</a>( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
     $                            RDSCAL, IPIV, JPIV )
                  END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Unpack solution vector(s)
</span><span class="comment">*</span><span class="comment">
</span>                  K = 1
                  II = MB*NB + 1
                  DO 100 JJ = 0, NB - 1
                     CALL SCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 )
                     CALL SCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 )
                     K = K + MB
                     II = II + MB
  100             CONTINUE
<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, RHS( 1 ), MB, 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, RHS( 1 ), MB, ONE,
     $                           F( 1, JS ), LDF )
                  END IF
                  IF( J.LT.Q ) THEN
                     K = MB*NB + 1
                     CALL SGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, MB, N-JE, NB, ONE, RHS( K ),
     $                           MB, 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, RHS( K ),
     $                           MB, E( JS, JE+1 ), LDE, ONE,
     $                           F( IS, JE+1 ), LDF )
                  END IF
<span class="comment">*</span><span class="comment">
</span>               END IF
<span class="comment">*</span><span class="comment">
</span>  110       CONTINUE
  120    CONTINUE
      ELSE
<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) + D(I, I)' * L(J, J)  =  C(I, J)
</span><span class="comment">*</span><span class="comment">             R(I, I)  * 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
         SCALOC = ONE
         DO 200 I = 1, P
<span class="comment">*</span><span class="comment">
</span>            IS = IWORK( I )
            ISP1 = IS + 1
            IE = IWORK( I+1 ) - 1
            MB = IE - IS + 1
            DO 190 J = Q, P + 2, -1
<span class="comment">*</span><span class="comment">
</span>               JS = IWORK( J )
               JSP1 = JS + 1
               JE = IWORK( J+1 ) - 1
               NB = JE - JS + 1
               ZDIM = MB*NB*2
               IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Build a 2-by-2 system Z' * x = RHS
</span><span class="comment">*</span><span class="comment">
</span>                  Z( 1, 1 ) = A( IS, IS )
                  Z( 2, 1 ) = -B( JS, JS )
                  Z( 1, 2 ) = D( IS, IS )
                  Z( 2, 2 ) = -E( JS, JS )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Set up right hand side(s)
</span><span class="comment">*</span><span class="comment">
</span>                  RHS( 1 ) = C( IS, JS )
                  RHS( 2 ) = F( IS, JS )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Solve Z' * x = RHS
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="SGETC2.665"></a><a href="sgetc2.f.html#SGETC2.1">SGETC2</a>( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
                  IF( IERR.GT.0 )
     $               INFO = IERR
<span class="comment">*</span><span class="comment">
</span>                  CALL <a name="SGESC2.669"></a><a href="sgesc2.f.html#SGESC2.1">SGESC2</a>( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
                  IF( SCALOC.NE.ONE ) THEN
                     DO 130 K = 1, N
                        CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
                        CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
  130                CONTINUE
                     SCALE = SCALE*SCALOC
                  END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Unpack solution vector(s)
</span><span class="comment">*</span><span class="comment">
</span>                  C( IS, JS ) = RHS( 1 )
                  F( IS, JS ) = RHS( 2 )
<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( J.GT.P+2 ) THEN
                     ALPHA = RHS( 1 )
                     CALL SAXPY( JS-1, ALPHA, B( 1, JS ), 1, F( IS, 1 ),
     $                           LDF )
                     ALPHA = RHS( 2 )
                     CALL SAXPY( JS-1, ALPHA, E( 1, JS ), 1, F( IS, 1 ),
     $                           LDF )
                  END IF
                  IF( I.LT.P ) THEN
                     ALPHA = -RHS( 1 )
                     CALL SAXPY( M-IE, ALPHA, A( IS, IE+1 ), LDA,
     $                           C( IE+1, JS ), 1 )
                     ALPHA = -RHS( 2 )
                     CALL SAXPY( M-IE, ALPHA, D( IS, IE+1 ), LDD,
     $                           C( IE+1, JS ), 1 )
                  END IF
<span class="comment">*</span><span class="comment">
</span>               ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Build a 4-by-4 system Z' * x = RHS
</span><span class="comment">*</span><span class="comment">
</span>                  Z( 1, 1 ) = A( IS, IS )
                  Z( 2, 1 ) = ZERO
                  Z( 3, 1 ) = -B( JS, JS )
                  Z( 4, 1 ) = -B( JSP1, JS )
<span class="comment">*</span><span class="comment">
</span>                  Z( 1, 2 ) = ZERO
                  Z( 2, 2 ) = A( IS, IS )
                  Z( 3, 2 ) = -B( JS, JSP1 )
                  Z( 4, 2 ) = -B( JSP1, JSP1 )
<span class="comment">*</span><span class="comment">
</span>                  Z( 1, 3 ) = D( IS, IS )
                  Z( 2, 3 ) = ZERO
                  Z( 3, 3 ) = -E( JS, JS )
                  Z( 4, 3 ) = ZERO
<span class="comment">*</span><span class="comment">
</span>                  Z( 1, 4 ) = ZERO
                  Z( 2, 4 ) = D( IS, IS )
                  Z( 3, 4 ) = -E( JS, JSP1 )
                  Z( 4, 4 ) = -E( JSP1, JSP1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Set up right hand side(s)
</span><span class="comment">*</span><span class="comment">
</span>                  RHS( 1 ) = C( IS, JS )
                  RHS( 2 ) = C( IS, JSP1 )
                  RHS( 3 ) = F( IS, JS )
                  RHS( 4 ) = F( IS, JSP1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Solve Z' * x = RHS
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="SGETC2.736"></a><a href="sgetc2.f.html#SGETC2.1">SGETC2</a>( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
                  IF( IERR.GT.0 )
     $               INFO = IERR
                  CALL <a name="SGESC2.739"></a><a href="sgesc2.f.html#SGESC2.1">SGESC2</a>( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
                  IF( SCALOC.NE.ONE ) THEN
                     DO 140 K = 1, N
                        CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
                        CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
  140                CONTINUE
                     SCALE = SCALE*SCALOC
                  END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Unpack solution vector(s)
</span><span class="comment">*</span><span class="comment">
</span>                  C( IS, JS ) = RHS( 1 )
                  C( IS, JSP1 ) = RHS( 2 )
                  F( IS, JS ) = RHS( 3 )
                  F( IS, JSP1 ) = RHS( 4 )
<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( J.GT.P+2 ) THEN
                     CALL SAXPY( JS-1, RHS( 1 ), B( 1, JS ), 1,
     $                           F( IS, 1 ), LDF )
                     CALL SAXPY( JS-1, RHS( 2 ), B( 1, JSP1 ), 1,
     $                           F( IS, 1 ), LDF )
                     CALL SAXPY( JS-1, RHS( 3 ), E( 1, JS ), 1,
     $                           F( IS, 1 ), LDF )
                     CALL SAXPY( JS-1, RHS( 4 ), E( 1, JSP1 ), 1,
     $                           F( IS, 1 ), LDF )
                  END IF
                  IF( I.LT.P ) THEN
                     CALL SGER( M-IE, NB, -ONE, A( IS, IE+1 ), LDA,