dsbtrd.f.html

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

</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">
</span>               DO 200 K = KDN + 1, 2, -1
                  J1 = J1 + KDN
                  J2 = J2 + KDN
<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="DLARGV.379"></a><a href="dlargv.f.html#DLARGV.1">DLARGV</a>( NR, AB( KD1, J1-KD1 ), INCA,
     $                            WORK( J1 ), KD1, D( J1 ), KD1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    apply plane rotations from one side
</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="DLARTV.386"></a><a href="dlartv.f.html#DLARTV.1">DLARTV</a> or DROT is used
</span><span class="comment">*</span><span class="comment">
</span>                     IF( NR.GT.2*KD-1 ) THEN
                        DO 130 L = 1, KD - 1
                           CALL <a name="DLARTV.390"></a><a href="dlartv.f.html#DLARTV.1">DLARTV</a>( NR, AB( KD1-L, J1-KD1+L ), INCA,
     $                                  AB( KD1-L+1, J1-KD1+L ), INCA,
     $                                  D( J1 ), WORK( J1 ), KD1 )
  130                   CONTINUE
                     ELSE
                        JEND = J1 + KD1*( NR-1 )
                        DO 140 JINC = J1, JEND, KD1
                           CALL DROT( KDM1, AB( KD, JINC-KD ), INCX,
     $                                AB( KD1, JINC-KD ), INCX,
     $                                D( JINC ), WORK( JINC ) )
  140                   CONTINUE
                     END IF
<span class="comment">*</span><span class="comment">
</span>                  END IF
<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+k-1,i)
</span><span class="comment">*</span><span class="comment">                       within the band
</span><span class="comment">*</span><span class="comment">
</span>                        CALL <a name="DLARTG.411"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>( AB( K-1, I ), AB( K, I ),
     $                               D( I+K-1 ), WORK( I+K-1 ), TEMP )
                        AB( K-1, I ) = TEMP
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       apply rotation from the left
</span><span class="comment">*</span><span class="comment">
</span>                        CALL DROT( K-3, AB( K-2, I+1 ), LDAB-1,
     $                             AB( K-1, I+1 ), LDAB-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="DLAR2V.429"></a><a href="dlar2v.f.html#DLAR2V.1">DLAR2V</a>( NR, AB( 1, J1-1 ), AB( 1, J1 ),
     $                            AB( 2, J1-1 ), 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 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="DLARTV.437"></a><a href="dlartv.f.html#DLARTV.1">DLARTV</a> or DROT is used
</span><span class="comment">*</span><span class="comment">
</span>                  IF( NR.GT.0 ) THEN
                     IF( NR.GT.2*KD-1 ) THEN
                        DO 150 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="DLARTV.448"></a><a href="dlartv.f.html#DLARTV.1">DLARTV</a>( NRT, AB( L+2, J1-1 ), INCA,
     $                                     AB( L+1, J1 ), INCA, D( J1 ),
     $                                     WORK( J1 ), KD1 )
  150                   CONTINUE
                     ELSE
                        J1END = J1 + KD1*( NR-2 )
                        IF( J1END.GE.J1 ) THEN
                           DO 160 J1INC = J1, J1END, KD1
                              CALL DROT( KDM1, AB( 3, J1INC-1 ), 1,
     $                                   AB( 2, J1INC ), 1, D( J1INC ),
     $                                   WORK( J1INC ) )
  160                      CONTINUE
                        END IF
                        LEND = MIN( KDM1, N-J2 )
                        LAST = J1END + KD1
                        IF( LEND.GT.0 )
     $                     CALL DROT( LEND, AB( 3, LAST-1 ), 1,
     $                                AB( 2, LAST ), 1, D( LAST ),
     $                                WORK( LAST ) )
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span><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 170 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 DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
     $                                1, D( J ), WORK( J ) )
  170                   CONTINUE
                     ELSE
<span class="comment">*</span><span class="comment">
</span>                        DO 180 J = J1, J2, KD1
                           CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
     $                                D( J ), WORK( J ) )
  180                   CONTINUE
                     END IF
                  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 190 J = J1, J2, KD1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    create nonzero element a(j+kd,j-1) outside the
</span><span class="comment">*</span><span class="comment">                    band and store it in WORK
</span><span class="comment">*</span><span class="comment">
</span>                     WORK( J+KD ) = WORK( J )*AB( KD1, J )
                     AB( KD1, J ) = D( J )*AB( KD1, J )
  190             CONTINUE
  200          CONTINUE
  210       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">           copy off-diagonal elements to E
</span><span class="comment">*</span><span class="comment">
</span>            DO 220 I = 1, N - 1
               E( I ) = AB( 2, I )
  220       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 230 I = 1, N - 1
               E( I ) = ZERO
  230       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 240 I = 1, N
            D( I ) = AB( 1, I )
  240    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DSBTRD.550"></a><a href="dsbtrd.f.html#DSBTRD.1">DSBTRD</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>