dsbgst.f.html

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

HTML
1,028 行
字号
12345
</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="DLARGV.604"></a><a href="dlargv.f.html#DLARGV.1">DLARGV</a>( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ),
     $                      KA1, WORK( N+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="DLARTV.611"></a><a href="dlartv.f.html#DLARTV.1">DLARTV</a>( NR, AB( L+1, J2-L ), INCA,
     $                         AB( L+2, J2-L ), INCA, WORK( N+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="DLAR2V.619"></a><a href="dlar2v.f.html#DLAR2V.1">DLAR2V</a>( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
     $                      INCA, WORK( N+J2-M ), WORK( J2-M ), KA1 )
<span class="comment">*</span><span class="comment">
</span>            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="DLARTV.629"></a><a href="dlartv.f.html#DLARTV.1">DLARTV</a>( NRT, AB( KA1-L+1, J2 ), INCA,
     $                         AB( KA1-L, J2+1 ), INCA, WORK( N+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 DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
     $                       WORK( N+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="DLARTV.667"></a><a href="dlartv.f.html#DLARTV.1">DLARTV</a>( NRT, AB( KA1-L+1, J2-KA ), INCA,
     $                         AB( KA1-L, J2-KA+1 ), INCA,
     $                         WORK( N+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 )
               WORK( N+J ) = WORK( N+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 ) = WORK( N+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="DLARGV.700"></a><a href="dlargv.f.html#DLARGV.1">DLARGV</a>( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1,
     $                      WORK( N+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="DLARTV.706"></a><a href="dlartv.f.html#DLARTV.1">DLARTV</a>( NR, AB( L+1, J2-L ), INCA,
     $                         AB( L+2, J2-L ), INCA, WORK( N+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="DLAR2V.714"></a><a href="dlar2v.f.html#DLAR2V.1">DLAR2V</a>( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
     $                      INCA, WORK( N+J2 ), WORK( J2 ), KA1 )
<span class="comment">*</span><span class="comment">
</span>            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="DLARTV.724"></a><a href="dlartv.f.html#DLARTV.1">DLARTV</a>( NRT, AB( KA1-L+1, J2 ), INCA,
     $                         AB( KA1-L, J2+1 ), INCA, WORK( N+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 DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
     $                       WORK( N+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="DLARTV.748"></a><a href="dlartv.f.html#DLARTV.1">DLARTV</a>( NRT, AB( KA1-L+1, J2 ), INCA,
     $                         AB( KA1-L, J2+1 ), INCA, WORK( N+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, I - KB + 2*KA + 1, -1
               WORK( N+J-M ) = WORK( N+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">
</span><span class="comment">*</span><span class="comment">     To avoid duplicating code, the two loops are merged.
</span><span class="comment">*</span><span class="comment">
</span>      UPDATE = .TRUE.
      I = 0
  490 CONTINUE
      IF( UPDATE ) THEN
         I = I + 1
         KBT = MIN( KB, M-I )
         I0 = I + 1
         I1 = MAX( 1, I-KA )
         I2 = I + KBT - KA1
         IF( I.GT.M ) THEN
            UPDATE = .FALSE.
            I = I - 1
            I0 = M + 1
            IF( KA.EQ.0 )
     $         RETURN
            GO TO 490
         END IF
      ELSE
         I = I - KA
         IF( I.LT.2 )
     $      RETURN
      END IF
<span class="comment">*</span><span class="comment">
12345

dsbgst.f.html - 源码说明

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

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

⌨️ 快捷键说明

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