sgbbrd.f.html

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

</span>               IF( NR.GT.0 )
     $            CALL <a name="SLARGV.295"></a><a href="slargv.f.html#SLARGV.1">SLARGV</a>( NR, AB( 1, J1+KUN-1 ), INCA,
     $                         WORK( J1+KUN ), KB1, WORK( MN+J1+KUN ),
     $                         KB1 )
<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>               DO 50 L = 1, KB
                  IF( J2+L-1.GT.M ) THEN
                     NRT = NR - 1
                  ELSE
                     NRT = NR
                  END IF
                  IF( NRT.GT.0 )
     $               CALL <a name="SLARTV.308"></a><a href="slartv.f.html#SLARTV.1">SLARTV</a>( NRT, AB( L+1, J1+KUN-1 ), INCA,
     $                            AB( L, J1+KUN ), INCA,
     $                            WORK( MN+J1+KUN ), WORK( J1+KUN ),
     $                            KB1 )
   50          CONTINUE
<span class="comment">*</span><span class="comment">
</span>               IF( ML.EQ.ML0 .AND. MU.GT.MU0 ) THEN
                  IF( MU.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+mu-1)
</span><span class="comment">*</span><span class="comment">                    within the band, and apply rotation from the right
</span><span class="comment">*</span><span class="comment">
</span>                     CALL <a name="SLARTG.320"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( AB( KU-MU+3, I+MU-2 ),
     $                            AB( KU-MU+2, I+MU-1 ),
     $                            WORK( MN+I+MU-1 ), WORK( I+MU-1 ),
     $                            RA )
                     AB( KU-MU+3, I+MU-2 ) = RA
                     CALL SROT( MIN( KL+MU-2, M-I ),
     $                          AB( KU-MU+4, I+MU-2 ), 1,
     $                          AB( KU-MU+3, I+MU-1 ), 1,
     $                          WORK( MN+I+MU-1 ), WORK( I+MU-1 ) )
                  END IF
                  NR = NR + 1
                  J1 = J1 - KB1
               END IF
<span class="comment">*</span><span class="comment">
</span>               IF( WANTPT ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 accumulate product of plane rotations in P'
</span><span class="comment">*</span><span class="comment">
</span>                  DO 60 J = J1, J2, KB1
                     CALL SROT( N, PT( J+KUN-1, 1 ), LDPT,
     $                          PT( J+KUN, 1 ), LDPT, WORK( MN+J+KUN ),
     $                          WORK( J+KUN ) )
   60             CONTINUE
               END IF
<span class="comment">*</span><span class="comment">
</span>               IF( J2+KB.GT.M ) 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 - KB1
               END IF
<span class="comment">*</span><span class="comment">
</span>               DO 70 J = J1, J2, KB1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 create nonzero element a(j+kl+ku,j+ku-1) below the
</span><span class="comment">*</span><span class="comment">                 band and store it in WORK(1:n)
</span><span class="comment">*</span><span class="comment">
</span>                  WORK( J+KB ) = WORK( J+KUN )*AB( KLU1, J+KUN )
                  AB( KLU1, J+KUN ) = WORK( MN+J+KUN )*AB( KLU1, J+KUN )
   70          CONTINUE
<span class="comment">*</span><span class="comment">
</span>               IF( ML.GT.ML0 ) THEN
                  ML = ML - 1
               ELSE
                  MU = MU - 1
               END IF
   80       CONTINUE
   90    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( KU.EQ.0 .AND. KL.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        A has been reduced to lower bidiagonal form
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Transform lower bidiagonal form to upper bidiagonal by applying
</span><span class="comment">*</span><span class="comment">        plane rotations from the left, storing diagonal elements in D
</span><span class="comment">*</span><span class="comment">        and off-diagonal elements in E
</span><span class="comment">*</span><span class="comment">
</span>         DO 100 I = 1, MIN( M-1, N )
            CALL <a name="SLARTG.380"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( AB( 1, I ), AB( 2, I ), RC, RS, RA )
            D( I ) = RA
            IF( I.LT.N ) THEN
               E( I ) = RS*AB( 1, I+1 )
               AB( 1, I+1 ) = RC*AB( 1, I+1 )
            END IF
            IF( WANTQ )
     $         CALL SROT( M, Q( 1, I ), 1, Q( 1, I+1 ), 1, RC, RS )
            IF( WANTC )
     $         CALL SROT( NCC, C( I, 1 ), LDC, C( I+1, 1 ), LDC, RC,
     $                    RS )
  100    CONTINUE
         IF( M.LE.N )
     $      D( M ) = AB( 1, M )
      ELSE IF( KU.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        A has been reduced to upper bidiagonal form
</span><span class="comment">*</span><span class="comment">
</span>         IF( M.LT.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Annihilate a(m,m+1) by applying plane rotations from the
</span><span class="comment">*</span><span class="comment">           right, storing diagonal elements in D and off-diagonal
</span><span class="comment">*</span><span class="comment">           elements in E
</span><span class="comment">*</span><span class="comment">
</span>            RB = AB( KU, M+1 )
            DO 110 I = M, 1, -1
               CALL <a name="SLARTG.406"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( AB( KU+1, I ), RB, RC, RS, RA )
               D( I ) = RA
               IF( I.GT.1 ) THEN
                  RB = -RS*AB( KU, I )
                  E( I-1 ) = RC*AB( KU, I )
               END IF
               IF( WANTPT )
     $            CALL SROT( N, PT( I, 1 ), LDPT, PT( M+1, 1 ), LDPT,
     $                       RC, RS )
  110       CONTINUE
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Copy off-diagonal elements to E and diagonal elements to D
</span><span class="comment">*</span><span class="comment">
</span>            DO 120 I = 1, MINMN - 1
               E( I ) = AB( KU, I+1 )
  120       CONTINUE
            DO 130 I = 1, MINMN
               D( I ) = AB( KU+1, I )
  130       CONTINUE
         END IF
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        A is diagonal. Set elements of E to zero and copy diagonal
</span><span class="comment">*</span><span class="comment">        elements to D.
</span><span class="comment">*</span><span class="comment">
</span>         DO 140 I = 1, MINMN - 1
            E( I ) = ZERO
  140    CONTINUE
         DO 150 I = 1, MINMN
            D( I ) = AB( 1, I )
  150    CONTINUE
      END IF
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SGBBRD.441"></a><a href="sgbbrd.f.html#SGBBRD.1">SGBBRD</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>