dsbgst.f.html

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

HTML
1,028 行
字号
12345
</span>      IF( I.LT.M-KBT ) THEN
         NX = M
      ELSE
         NX = N
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( UPPER ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Transform A, working with the upper triangle
</span><span class="comment">*</span><span class="comment">
</span>         IF( UPDATE ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Form  inv(S(i))**T * A * inv(S(i))
</span><span class="comment">*</span><span class="comment">
</span>            BII = BB( KB1, I )
            DO 500 J = I1, I
               AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
  500       CONTINUE
            DO 510 J = I, MIN( N, I+KA )
               AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
  510       CONTINUE
            DO 540 K = I + 1, I + KBT
               DO 520 J = K, I + KBT
                  AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
     $                               BB( I-J+KB1, J )*AB( I-K+KA1, K ) -
     $                               BB( I-K+KB1, K )*AB( I-J+KA1, J ) +
     $                               AB( KA1, I )*BB( I-J+KB1, J )*
     $                               BB( I-K+KB1, K )
  520          CONTINUE
               DO 530 J = I + KBT + 1, MIN( N, I+KA )
                  AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
     $                               BB( I-K+KB1, K )*AB( I-J+KA1, J )
  530          CONTINUE
  540       CONTINUE
            DO 560 J = I1, I
               DO 550 K = I + 1, MIN( J+KA, I+KBT )
                  AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
     $                               BB( I-K+KB1, K )*AB( J-I+KA1, I )
  550          CONTINUE
  560       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 inv(S(i))
</span><span class="comment">*</span><span class="comment">
</span>               CALL DSCAL( NX, ONE / BII, X( 1, I ), 1 )
               IF( KBT.GT.0 )
     $            CALL DGER( NX, KBT, -ONE, X( 1, I ), 1, BB( KB, I+1 ),
     $                       LDBB-1, X( 1, I+1 ), LDX )
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           store a(i1,i) in RA1 for use in next loop over K
</span><span class="comment">*</span><span class="comment">
</span>            RA1 = AB( I1-I+KA1, I )
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Generate and apply vectors of rotations to chase all the
</span><span class="comment">*</span><span class="comment">        existing bulges KA positions up toward the top of the band
</span><span class="comment">*</span><span class="comment">
</span>         DO 610 K = 1, KB - 1
            IF( UPDATE ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Determine the rotations which would annihilate the bulge
</span><span class="comment">*</span><span class="comment">              which has in theory just been created
</span><span class="comment">*</span><span class="comment">
</span>               IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 generate rotation to annihilate a(i+k-ka-1,i)
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="DLARTG.875"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>( AB( K+1, I ), RA1, WORK( N+I+K-KA ),
     $                         WORK( I+K-KA ), RA )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 create nonzero element a(i+k-ka-1,i+k) outside the
</span><span class="comment">*</span><span class="comment">                 band and store it in WORK(m-kb+i+k)
</span><span class="comment">*</span><span class="comment">
</span>                  T = -BB( KB1-K, I+K )*RA1
                  WORK( M-KB+I+K ) = WORK( N+I+K-KA )*T -
     $                               WORK( I+K-KA )*AB( 1, I+K )
                  AB( 1, I+K ) = WORK( I+K-KA )*T +
     $                           WORK( N+I+K-KA )*AB( 1, I+K )
                  RA1 = RA
               END IF
            END IF
            J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
            NR = ( J2+KA-1 ) / KA1
            J1 = J2 - ( NR-1 )*KA1
            IF( UPDATE ) THEN
               J2T = MIN( J2, I-2*KA+K-1 )
            ELSE
               J2T = J2
            END IF
            NRT = ( J2T+KA-1 ) / KA1
            DO 570 J = J1, J2T, 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( 1, J+KA-1 )
               AB( 1, J+KA-1 ) = WORK( N+J )*AB( 1, J+KA-1 )
  570       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="DLARGV.911"></a><a href="dlargv.f.html#DLARGV.1">DLARGV</a>( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1,
     $                      WORK( N+J1 ), 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 580 L = 1, KA - 1
                  CALL <a name="DLARTV.918"></a><a href="dlartv.f.html#DLARTV.1">DLARTV</a>( NR, AB( KA1-L, J1+L ), INCA,
     $                         AB( KA-L, J1+L ), INCA, WORK( N+J1 ),
     $                         WORK( J1 ), KA1 )
  580          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.926"></a><a href="dlar2v.f.html#DLAR2V.1">DLAR2V</a>( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
     $                      AB( KA, J1 ), INCA, WORK( N+J1 ),
     $                      WORK( J1 ), 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 590 L = KA - 1, KB - K + 1, -1
               NRT = ( J2+L-1 ) / KA1
               J1T = J2 - ( NRT-1 )*KA1
               IF( NRT.GT.0 )
     $            CALL <a name="DLARTV.938"></a><a href="dlartv.f.html#DLARTV.1">DLARTV</a>( NRT, AB( L, J1T ), INCA,
     $                         AB( L+1, J1T-1 ), INCA, WORK( N+J1T ),
     $                         WORK( J1T ), KA1 )
  590       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 600 J = J1, J2, KA1
                  CALL DROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
     $                       WORK( N+J ), WORK( J ) )
  600          CONTINUE
            END IF
  610    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         IF( UPDATE ) THEN
            IF( I2.GT.0 .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(m-kb+i+kbt)
</span><span class="comment">*</span><span class="comment">
</span>               WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span>         DO 650 K = KB, 1, -1
            IF( UPDATE ) THEN
               J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
            ELSE
               J2 = I + K + 1 - MAX( 1, K+I0-M )*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 620 L = KB - K, 1, -1
               NRT = ( J2+KA+L-1 ) / KA1
               J1T = J2 - ( NRT-1 )*KA1
               IF( NRT.GT.0 )
     $            CALL <a name="DLARTV.977"></a><a href="dlartv.f.html#DLARTV.1">DLARTV</a>( NRT, AB( L, J1T+KA ), INCA,
     $                         AB( L+1, J1T+KA-1 ), INCA,
     $                         WORK( N+M-KB+J1T+KA ),
     $                         WORK( M-KB+J1T+KA ), KA1 )
  620       CONTINUE
            NR = ( J2+KA-1 ) / KA1
            J1 = J2 - ( NR-1 )*KA1
            DO 630 J = J1, J2, KA1
               WORK( M-KB+J ) = WORK( M-KB+J+KA )
               WORK( N+M-KB+J ) = WORK( N+M-KB+J+KA )
  630       CONTINUE
            DO 640 J = J1, J2, 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(m-kb+j)
</span><span class="comment">*</span><span class="comment">
</span>               WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 )
               AB( 1, J+KA-1 ) = WORK( N+M-KB+J )*AB( 1, J+KA-1 )
  640       CONTINUE
            IF( UPDATE ) THEN
               IF( I+K.GT.KA1 .AND. K.LE.KBT )
     $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
            END IF
  650    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         DO 690 K = KB, 1, -1
            J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
            NR = ( J2+KA-1 ) / 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">
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 + -
显示快捷键?