stgsy2.f.html

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

HTML
981
字号
     $                          RHS( 1 ), 1, C( IE+1, JS ), LDC )
                     CALL SGER( M-IE, NB, -ONE, D( IS, IE+1 ), LDD,
     $                          RHS( 3 ), 1, C( IE+1, JS ), LDC )
                  END IF
<span class="comment">*</span><span class="comment">
</span>               ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) 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 ) = A( IS, ISP1 )
                  Z( 3, 1 ) = -B( JS, JS )
                  Z( 4, 1 ) = ZERO
<span class="comment">*</span><span class="comment">
</span>                  Z( 1, 2 ) = A( ISP1, IS )
                  Z( 2, 2 ) = A( ISP1, ISP1 )
                  Z( 3, 2 ) = ZERO
                  Z( 4, 2 ) = -B( JS, JS )
<span class="comment">*</span><span class="comment">
</span>                  Z( 1, 3 ) = D( IS, IS )
                  Z( 2, 3 ) = D( IS, ISP1 )
                  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( ISP1, ISP1 )
                  Z( 3, 4 ) = ZERO
                  Z( 4, 4 ) = -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 ) = C( ISP1, JS )
                  RHS( 3 ) = F( IS, JS )
                  RHS( 4 ) = F( ISP1, 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.808"></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.812"></a><a href="sgesc2.f.html#SGESC2.1">SGESC2</a>( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
                  IF( SCALOC.NE.ONE ) THEN
                     DO 150 K = 1, N
                        CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
                        CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
  150                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( ISP1, JS ) = RHS( 2 )
                  F( IS, JS ) = RHS( 3 )
                  F( ISP1, JS ) = 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 SGER( MB, JS-1, ONE, RHS( 1 ), 1, B( 1, JS ),
     $                          1, F( IS, 1 ), LDF )
                     CALL SGER( MB, JS-1, ONE, RHS( 3 ), 1, E( 1, JS ),
     $                          1, F( IS, 1 ), LDF )
                  END IF
                  IF( I.LT.P ) THEN
                     CALL SGEMV( <span class="string">'T'</span>, MB, M-IE, -ONE, A( IS, IE+1 ),
     $                           LDA, RHS( 1 ), 1, ONE, C( IE+1, JS ),
     $                           1 )
                     CALL SGEMV( <span class="string">'T'</span>, MB, M-IE, -ONE, D( IS, IE+1 ),
     $                           LDD, RHS( 3 ), 1, ONE, C( IE+1, JS ),
     $                           1 )
                  END IF
<span class="comment">*</span><span class="comment">
</span>               ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Build an 8-by-8 system Z' * x = RHS
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="SLASET.850"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'F'</span>, LDZ, LDZ, ZERO, ZERO, Z, LDZ )
<span class="comment">*</span><span class="comment">
</span>                  Z( 1, 1 ) = A( IS, IS )
                  Z( 2, 1 ) = A( IS, ISP1 )
                  Z( 5, 1 ) = -B( JS, JS )
                  Z( 7, 1 ) = -B( JSP1, JS )
<span class="comment">*</span><span class="comment">
</span>                  Z( 1, 2 ) = A( ISP1, IS )
                  Z( 2, 2 ) = A( ISP1, ISP1 )
                  Z( 6, 2 ) = -B( JS, JS )
                  Z( 8, 2 ) = -B( JSP1, JS )
<span class="comment">*</span><span class="comment">
</span>                  Z( 3, 3 ) = A( IS, IS )
                  Z( 4, 3 ) = A( IS, ISP1 )
                  Z( 5, 3 ) = -B( JS, JSP1 )
                  Z( 7, 3 ) = -B( JSP1, JSP1 )
<span class="comment">*</span><span class="comment">
</span>                  Z( 3, 4 ) = A( ISP1, IS )
                  Z( 4, 4 ) = A( ISP1, ISP1 )
                  Z( 6, 4 ) = -B( JS, JSP1 )
                  Z( 8, 4 ) = -B( JSP1, JSP1 )
<span class="comment">*</span><span class="comment">
</span>                  Z( 1, 5 ) = D( IS, IS )
                  Z( 2, 5 ) = D( IS, ISP1 )
                  Z( 5, 5 ) = -E( JS, JS )
<span class="comment">*</span><span class="comment">
</span>                  Z( 2, 6 ) = D( ISP1, ISP1 )
                  Z( 6, 6 ) = -E( JS, JS )
<span class="comment">*</span><span class="comment">
</span>                  Z( 3, 7 ) = D( IS, IS )
                  Z( 4, 7 ) = D( IS, ISP1 )
                  Z( 5, 7 ) = -E( JS, JSP1 )
                  Z( 7, 7 ) = -E( JSP1, JSP1 )
<span class="comment">*</span><span class="comment">
</span>                  Z( 4, 8 ) = D( ISP1, ISP1 )
                  Z( 6, 8 ) = -E( JS, JSP1 )
                  Z( 8, 8 ) = -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>                  K = 1
                  II = MB*NB + 1
                  DO 160 JJ = 0, NB - 1
                     CALL SCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 )
                     CALL SCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 )
                     K = K + MB
                     II = II + MB
  160             CONTINUE
<span class="comment">*</span><span class="comment">
</span><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.902"></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.906"></a><a href="sgesc2.f.html#SGESC2.1">SGESC2</a>( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
                  IF( SCALOC.NE.ONE ) THEN
                     DO 170 K = 1, N
                        CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
                        CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
  170                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>                  K = 1
                  II = MB*NB + 1
                  DO 180 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
  180             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( 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
<span class="comment">*</span><span class="comment">
</span>               END IF
<span class="comment">*</span><span class="comment">
</span>  190       CONTINUE
  200    CONTINUE
<span class="comment">*</span><span class="comment">
</span>      END IF
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="STGSY2.954"></a><a href="stgsy2.f.html#STGSY2.1">STGSY2</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>

⌨️ 快捷键说明

复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?