chbtrd.f.html

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

HTML
613 行
字号
123
<span class="comment">*</span><span class="comment">
</span>                  IF( NR.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    generate plane rotations to annihilate nonzero
</span><span class="comment">*</span><span class="comment">                    elements which have been created outside the band
</span><span class="comment">*</span><span class="comment">
</span>                     CALL <a name="CLARGV.193"></a><a href="clargv.f.html#CLARGV.1">CLARGV</a>( NR, AB( 1, J1-1 ), INCA, WORK( J1 ),
     $                            KD1, D( J1 ), KD1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    apply rotations from the right
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Dependent on the the number of diagonals either
</span><span class="comment">*</span><span class="comment">                    <a name="CLARTV.200"></a><a href="clartv.f.html#CLARTV.1">CLARTV</a> or <a name="CROT.200"></a><a href="crot.f.html#CROT.1">CROT</a> is used
</span><span class="comment">*</span><span class="comment">
</span>                     IF( NR.GE.2*KD-1 ) THEN
                        DO 10 L = 1, KD - 1
                           CALL <a name="CLARTV.204"></a><a href="clartv.f.html#CLARTV.1">CLARTV</a>( NR, AB( L+1, J1-1 ), INCA,
     $                                  AB( L, J1 ), INCA, D( J1 ),
     $                                  WORK( J1 ), KD1 )
   10                   CONTINUE
<span class="comment">*</span><span class="comment">
</span>                     ELSE
                        JEND = J1 + ( NR-1 )*KD1
                        DO 20 JINC = J1, JEND, KD1
                           CALL <a name="CROT.212"></a><a href="crot.f.html#CROT.1">CROT</a>( KDM1, AB( 2, JINC-1 ), 1,
     $                                AB( 1, JINC ), 1, D( JINC ),
     $                                WORK( JINC ) )
   20                   CONTINUE
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span>                  IF( K.GT.2 ) THEN
                     IF( K.LE.N-I+1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       generate plane rotation to annihilate a(i,i+k-1)
</span><span class="comment">*</span><span class="comment">                       within the band
</span><span class="comment">*</span><span class="comment">
</span>                        CALL <a name="CLARTG.226"></a><a href="clartg.f.html#CLARTG.1">CLARTG</a>( AB( KD-K+3, I+K-2 ),
     $                               AB( KD-K+2, I+K-1 ), D( I+K-1 ),
     $                               WORK( I+K-1 ), TEMP )
                        AB( KD-K+3, I+K-2 ) = TEMP
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       apply rotation from the right
</span><span class="comment">*</span><span class="comment">
</span>                        CALL <a name="CROT.233"></a><a href="crot.f.html#CROT.1">CROT</a>( K-3, AB( KD-K+4, I+K-2 ), 1,
     $                             AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ),
     $                             WORK( I+K-1 ) )
                     END IF
                     NR = NR + 1
                     J1 = J1 - KDN - 1
                  END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 apply plane rotations from both sides to diagonal
</span><span class="comment">*</span><span class="comment">                 blocks
</span><span class="comment">*</span><span class="comment">
</span>                  IF( NR.GT.0 )
     $               CALL <a name="CLAR2V.245"></a><a href="clar2v.f.html#CLAR2V.1">CLAR2V</a>( NR, AB( KD1, J1-1 ), AB( KD1, J1 ),
     $                            AB( KD, J1 ), INCA, D( J1 ),
     $                            WORK( J1 ), KD1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 apply plane rotations from the left
</span><span class="comment">*</span><span class="comment">
</span>                  IF( NR.GT.0 ) THEN
                     CALL <a name="CLACGV.252"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( NR, WORK( J1 ), KD1 )
                     IF( 2*KD-1.LT.NR ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Dependent on the the number of diagonals either
</span><span class="comment">*</span><span class="comment">                    <a name="CLARTV.256"></a><a href="clartv.f.html#CLARTV.1">CLARTV</a> or <a name="CROT.256"></a><a href="crot.f.html#CROT.1">CROT</a> is used
</span><span class="comment">*</span><span class="comment">
</span>                        DO 30 L = 1, KD - 1
                           IF( J2+L.GT.N ) THEN
                              NRT = NR - 1
                           ELSE
                              NRT = NR
                           END IF
                           IF( NRT.GT.0 )
     $                        CALL <a name="CLARTV.265"></a><a href="clartv.f.html#CLARTV.1">CLARTV</a>( NRT, AB( KD-L, J1+L ), INCA,
     $                                     AB( KD-L+1, J1+L ), INCA,
     $                                     D( J1 ), WORK( J1 ), KD1 )
   30                   CONTINUE
                     ELSE
                        J1END = J1 + KD1*( NR-2 )
                        IF( J1END.GE.J1 ) THEN
                           DO 40 JIN = J1, J1END, KD1
                              CALL <a name="CROT.273"></a><a href="crot.f.html#CROT.1">CROT</a>( KD-1, AB( KD-1, JIN+1 ), INCX,
     $                                   AB( KD, JIN+1 ), INCX,
     $                                   D( JIN ), WORK( JIN ) )
   40                      CONTINUE
                        END IF
                        LEND = MIN( KDM1, N-J2 )
                        LAST = J1END + KD1
                        IF( LEND.GT.0 )
     $                     CALL <a name="CROT.281"></a><a href="crot.f.html#CROT.1">CROT</a>( LEND, AB( KD-1, LAST+1 ), INCX,
     $                                AB( KD, LAST+1 ), INCX, D( LAST ),
     $                                WORK( LAST ) )
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  IF( WANTQ ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    accumulate product of plane rotations in Q
</span><span class="comment">*</span><span class="comment">
</span>                     IF( INITQ ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 take advantage of the fact that Q was
</span><span class="comment">*</span><span class="comment">                 initially the Identity matrix
</span><span class="comment">*</span><span class="comment">
</span>                        IQEND = MAX( IQEND, J2 )
                        I2 = MAX( 0, K-3 )
                        IQAEND = 1 + I*KD
                        IF( K.EQ.2 )
     $                     IQAEND = IQAEND + KD
                        IQAEND = MIN( IQAEND, IQEND )
                        DO 50 J = J1, J2, KD1
                           IBL = I - I2 / KDM1
                           I2 = I2 + 1
                           IQB = MAX( 1, J-IBL )
                           NQ = 1 + IQAEND - IQB
                           IQAEND = MIN( IQAEND+KD, IQEND )
                           CALL <a name="CROT.308"></a><a href="crot.f.html#CROT.1">CROT</a>( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
     $                                1, D( J ), CONJG( WORK( J ) ) )
   50                   CONTINUE
                     ELSE
<span class="comment">*</span><span class="comment">
</span>                        DO 60 J = J1, J2, KD1
                           CALL <a name="CROT.314"></a><a href="crot.f.html#CROT.1">CROT</a>( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
     $                                D( J ), CONJG( WORK( J ) ) )
   60                   CONTINUE
                     END IF
<span class="comment">*</span><span class="comment">
</span>                  END IF
<span class="comment">*</span><span class="comment">
</span>                  IF( J2+KDN.GT.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    adjust J2 to keep within the bounds of the matrix
</span><span class="comment">*</span><span class="comment">
</span>                     NR = NR - 1
                     J2 = J2 - KDN - 1
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  DO 70 J = J1, J2, KD1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    create nonzero element a(j-1,j+kd) outside the band
</span><span class="comment">*</span><span class="comment">                    and store it in WORK
</span><span class="comment">*</span><span class="comment">
</span>                     WORK( J+KD ) = WORK( J )*AB( 1, J+KD )
                     AB( 1, J+KD ) = D( J )*AB( 1, J+KD )
   70             CONTINUE
   80          CONTINUE
   90       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( KD.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           make off-diagonal elements real and copy them to E
</span><span class="comment">*</span><span class="comment">
</span>            DO 100 I = 1, N - 1
               T = AB( KD, I+1 )
               ABST = ABS( T )
               AB( KD, I+1 ) = ABST
               E( I ) = ABST
               IF( ABST.NE.ZERO ) THEN
                  T = T / ABST
               ELSE
                  T = CONE
               END IF
               IF( I.LT.N-1 )
     $            AB( KD, I+2 ) = AB( KD, I+2 )*T
               IF( WANTQ ) THEN
                  CALL CSCAL( N, CONJG( T ), Q( 1, I+1 ), 1 )
               END IF
  100       CONTINUE
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           set E to zero if original matrix was diagonal
</span><span class="comment">*</span><span class="comment">
</span>            DO 110 I = 1, N - 1
               E( I ) = ZERO
  110       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        copy diagonal elements to D
</span><span class="comment">*</span><span class="comment">
</span>         DO 120 I = 1, N
            D( I ) = AB( KD1, I )
  120    CONTINUE
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span>         IF( KD.GT.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Reduce to complex Hermitian tridiagonal form, working with
</span><span class="comment">*</span><span class="comment">           the lower triangle
</span><span class="comment">*</span><span class="comment">
</span>            NR = 0
            J1 = KDN + 2
            J2 = 1
<span class="comment">*</span><span class="comment">
</span>            AB( 1, 1 ) = REAL( AB( 1, 1 ) )
            DO 210 I = 1, N - 2
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Reduce i-th column of matrix to tridiagonal form
</span><span class="comment">*</span><span class="comment">
123

chbtrd.f.html - 源码说明

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

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

⌨️ 快捷键说明

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