chbgst.f.html

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

HTML
1,017 行
字号
12345
                  RA1 = RA
               END IF
            END IF
            J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
            NR = ( N-J2+KA ) / KA1
            J1 = J2 + ( NR-1 )*KA1
            IF( UPDATE ) THEN
               J2T = MAX( J2, I+2*KA-K+1 )
            ELSE
               J2T = J2
            END IF
            NRT = ( N-J2T+KA ) / KA1
            DO 320 J = J2T, J1, KA1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              create nonzero element a(j+1,j-ka) outside the band
</span><span class="comment">*</span><span class="comment">              and store it in WORK(j-m)
</span><span class="comment">*</span><span class="comment">
</span>               WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 )
               AB( KA1, J-KA+1 ) = RWORK( J-M )*AB( KA1, J-KA+1 )
  320       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           generate rotations in 1st set to annihilate elements which
</span><span class="comment">*</span><span class="comment">           have been created outside the band
</span><span class="comment">*</span><span class="comment">
</span>            IF( NRT.GT.0 )
     $         CALL <a name="CLARGV.619"></a><a href="clargv.f.html#CLARGV.1">CLARGV</a>( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ),
     $                      KA1, RWORK( J2T-M ), KA1 )
            IF( NR.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              apply rotations in 1st set from the left
</span><span class="comment">*</span><span class="comment">
</span>               DO 330 L = 1, KA - 1
                  CALL <a name="CLARTV.626"></a><a href="clartv.f.html#CLARTV.1">CLARTV</a>( NR, AB( L+1, J2-L ), INCA,
     $                         AB( L+2, J2-L ), INCA, RWORK( J2-M ),
     $                         WORK( J2-M ), KA1 )
  330          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              apply rotations in 1st set from both sides to diagonal
</span><span class="comment">*</span><span class="comment">              blocks
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="CLAR2V.634"></a><a href="clar2v.f.html#CLAR2V.1">CLAR2V</a>( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
     $                      INCA, RWORK( J2-M ), WORK( J2-M ), KA1 )
<span class="comment">*</span><span class="comment">
</span>               CALL <a name="CLACGV.637"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( NR, WORK( J2-M ), KA1 )
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           start applying rotations in 1st set from the right
</span><span class="comment">*</span><span class="comment">
</span>            DO 340 L = KA - 1, KB - K + 1, -1
               NRT = ( N-J2+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL <a name="CLARTV.645"></a><a href="clartv.f.html#CLARTV.1">CLARTV</a>( NRT, AB( KA1-L+1, J2 ), INCA,
     $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ),
     $                         WORK( J2-M ), KA1 )
  340       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            IF( WANTX ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              post-multiply X by product of rotations in 1st set
</span><span class="comment">*</span><span class="comment">
</span>               DO 350 J = J2, J1, KA1
                  CALL <a name="CROT.655"></a><a href="crot.f.html#CROT.1">CROT</a>( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
     $                       RWORK( J-M ), WORK( J-M ) )
  350          CONTINUE
            END IF
  360    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         IF( UPDATE ) THEN
            IF( I2.LE.N .AND. KBT.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              create nonzero element a(i-kbt+ka+1,i-kbt) outside the
</span><span class="comment">*</span><span class="comment">              band and store it in WORK(i-kbt)
</span><span class="comment">*</span><span class="comment">
</span>               WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span>         DO 400 K = KB, 1, -1
            IF( UPDATE ) THEN
               J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
            ELSE
               J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           finish applying rotations in 2nd set from the right
</span><span class="comment">*</span><span class="comment">
</span>            DO 370 L = KB - K, 1, -1
               NRT = ( N-J2+KA+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL <a name="CLARTV.683"></a><a href="clartv.f.html#CLARTV.1">CLARTV</a>( NRT, AB( KA1-L+1, J2-KA ), INCA,
     $                         AB( KA1-L, J2-KA+1 ), INCA,
     $                         RWORK( J2-KA ), WORK( J2-KA ), KA1 )
  370       CONTINUE
            NR = ( N-J2+KA ) / KA1
            J1 = J2 + ( NR-1 )*KA1
            DO 380 J = J1, J2, -KA1
               WORK( J ) = WORK( J-KA )
               RWORK( J ) = RWORK( J-KA )
  380       CONTINUE
            DO 390 J = J2, J1, KA1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              create nonzero element a(j+1,j-ka) outside the band
</span><span class="comment">*</span><span class="comment">              and store it in WORK(j)
</span><span class="comment">*</span><span class="comment">
</span>               WORK( J ) = WORK( J )*AB( KA1, J-KA+1 )
               AB( KA1, J-KA+1 ) = RWORK( J )*AB( KA1, J-KA+1 )
  390       CONTINUE
            IF( UPDATE ) THEN
               IF( I-K.LT.N-KA .AND. K.LE.KBT )
     $            WORK( I-K+KA ) = WORK( I-K )
            END IF
  400    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         DO 440 K = KB, 1, -1
            J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
            NR = ( N-J2+KA ) / KA1
            J1 = J2 + ( NR-1 )*KA1
            IF( NR.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              generate rotations in 2nd set to annihilate elements
</span><span class="comment">*</span><span class="comment">              which have been created outside the band
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="CLARGV.716"></a><a href="clargv.f.html#CLARGV.1">CLARGV</a>( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1,
     $                      RWORK( J2 ), KA1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              apply rotations in 2nd set from the left
</span><span class="comment">*</span><span class="comment">
</span>               DO 410 L = 1, KA - 1
                  CALL <a name="CLARTV.722"></a><a href="clartv.f.html#CLARTV.1">CLARTV</a>( NR, AB( L+1, J2-L ), INCA,
     $                         AB( L+2, J2-L ), INCA, RWORK( J2 ),
     $                         WORK( J2 ), KA1 )
  410          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              apply rotations in 2nd set from both sides to diagonal
</span><span class="comment">*</span><span class="comment">              blocks
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="CLAR2V.730"></a><a href="clar2v.f.html#CLAR2V.1">CLAR2V</a>( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
     $                      INCA, RWORK( J2 ), WORK( J2 ), KA1 )
<span class="comment">*</span><span class="comment">
</span>               CALL <a name="CLACGV.733"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( NR, WORK( J2 ), KA1 )
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           start applying rotations in 2nd set from the right
</span><span class="comment">*</span><span class="comment">
</span>            DO 420 L = KA - 1, KB - K + 1, -1
               NRT = ( N-J2+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL <a name="CLARTV.741"></a><a href="clartv.f.html#CLARTV.1">CLARTV</a>( NRT, AB( KA1-L+1, J2 ), INCA,
     $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2 ),
     $                         WORK( J2 ), KA1 )
  420       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            IF( WANTX ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              post-multiply X by product of rotations in 2nd set
</span><span class="comment">*</span><span class="comment">
</span>               DO 430 J = J2, J1, KA1
                  CALL <a name="CROT.751"></a><a href="crot.f.html#CROT.1">CROT</a>( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
     $                       RWORK( J ), WORK( J ) )
  430          CONTINUE
            END IF
  440    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         DO 460 K = 1, KB - 1
            J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           finish applying rotations in 1st set from the right
</span><span class="comment">*</span><span class="comment">
</span>            DO 450 L = KB - K, 1, -1
               NRT = ( N-J2+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL <a name="CLARTV.765"></a><a href="clartv.f.html#CLARTV.1">CLARTV</a>( NRT, AB( KA1-L+1, J2 ), INCA,
     $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ),
     $                         WORK( J2-M ), KA1 )
  450       CONTINUE
  460    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         IF( KB.GT.1 ) THEN
            DO 470 J = N - 1, I2 + KA, -1
               RWORK( J-M ) = RWORK( J-KA-M )
               WORK( J-M ) = WORK( J-KA-M )
  470       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span>      GO TO 10
<span class="comment">*</span><span class="comment">
</span>  480 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     **************************** Phase 2 *****************************
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     The logical structure of this phase is:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     UPDATE = .TRUE.
</span><span class="comment">*</span><span class="comment">     DO I = 1, M
</span><span class="comment">*</span><span class="comment">        use S(i) to update A and create a new bulge
</span><span class="comment">*</span><span class="comment">        apply rotations to push all bulges KA positions upward
</span><span class="comment">*</span><span class="comment">     END DO
</span><span class="comment">*</span><span class="comment">     UPDATE = .FALSE.
</span><span class="comment">*</span><span class="comment">     DO I = M - KA - 1, 2, -1
</span><span class="comment">*</span><span class="comment">        apply rotations to push all bulges KA positions upward
</span><span class="comment">*</span><span class="comment">     END DO
</span><span class="comment">*</span><span class="comment">
12345

chbgst.f.html - 源码说明

本页面展示了「famous linear algebra library (LAPACK) ports to windows」中的 chbgst.f.html 源码文件，采用 HTML 编程语言编写，共 1,017 行代码。您可以在线阅读完整代码内容，也可以返回资源详情页下载完整源码包进行本地学习和开发。

虫虫下载站收录了大量与algebra相关的技术资源，包括源代码、技术文档、电路图等，是电子工程师和嵌入式开发者的专业学习平台。

⌨️ 快捷键说明

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