zgbbrd.f.html

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

         INFO = -3
      ELSE IF( NCC.LT.0 ) THEN
         INFO = -4
      ELSE IF( KL.LT.0 ) THEN
         INFO = -5
      ELSE IF( KU.LT.0 ) THEN
         INFO = -6
      ELSE IF( LDAB.LT.KLU1 ) THEN
         INFO = -8
      ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. LDQ.LT.MAX( 1, M ) ) THEN
         INFO = -12
      ELSE IF( LDPT.LT.1 .OR. WANTPT .AND. LDPT.LT.MAX( 1, N ) ) THEN
         INFO = -14
      ELSE IF( LDC.LT.1 .OR. WANTC .AND. LDC.LT.MAX( 1, M ) ) THEN
         INFO = -16
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.163"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZGBBRD.163"></a><a href="zgbbrd.f.html#ZGBBRD.1">ZGBBRD</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize Q and P' to the unit matrix, if needed
</span><span class="comment">*</span><span class="comment">
</span>      IF( WANTQ )
     $   CALL <a name="ZLASET.170"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'Full'</span>, M, M, CZERO, CONE, Q, LDQ )
      IF( WANTPT )
     $   CALL <a name="ZLASET.172"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'Full'</span>, N, N, CZERO, CONE, PT, LDPT )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible.
</span><span class="comment">*</span><span class="comment">
</span>      IF( M.EQ.0 .OR. N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span>      MINMN = MIN( M, N )
<span class="comment">*</span><span class="comment">
</span>      IF( KL+KU.GT.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Reduce to upper bidiagonal form if KU &gt; 0; if KU = 0, reduce
</span><span class="comment">*</span><span class="comment">        first to lower bidiagonal form and then transform to upper
</span><span class="comment">*</span><span class="comment">        bidiagonal
</span><span class="comment">*</span><span class="comment">
</span>         IF( KU.GT.0 ) THEN
            ML0 = 1
            MU0 = 2
         ELSE
            ML0 = 2
            MU0 = 1
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Wherever possible, plane rotations are generated and applied in
</span><span class="comment">*</span><span class="comment">        vector operations of length NR over the index set J1:J2:KLU1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        The complex sines of the plane rotations are stored in WORK,
</span><span class="comment">*</span><span class="comment">        and the real cosines in RWORK.
</span><span class="comment">*</span><span class="comment">
</span>         KLM = MIN( M-1, KL )
         KUN = MIN( N-1, KU )
         KB = KLM + KUN
         KB1 = KB + 1
         INCA = KB1*LDAB
         NR = 0
         J1 = KLM + 2
         J2 = 1 - KUN
<span class="comment">*</span><span class="comment">
</span>         DO 90 I = 1, MINMN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Reduce i-th column and i-th row of matrix to bidiagonal form
</span><span class="comment">*</span><span class="comment">
</span>            ML = KLM + 1
            MU = KUN + 1
            DO 80 KK = 1, KB
               J1 = J1 + KB
               J2 = J2 + KB
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              generate plane rotations to annihilate nonzero elements
</span><span class="comment">*</span><span class="comment">              which have been created below the band
</span><span class="comment">*</span><span class="comment">
</span>               IF( NR.GT.0 )
     $            CALL <a name="ZLARGV.224"></a><a href="zlargv.f.html#ZLARGV.1">ZLARGV</a>( NR, AB( KLU1, J1-KLM-1 ), INCA,
     $                         WORK( J1 ), KB1, RWORK( J1 ), KB1 )
<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>               DO 10 L = 1, KB
                  IF( J2-KLM+L-1.GT.N ) THEN
                     NRT = NR - 1
                  ELSE
                     NRT = NR
                  END IF
                  IF( NRT.GT.0 )
     $               CALL <a name="ZLARTV.236"></a><a href="zlartv.f.html#ZLARTV.1">ZLARTV</a>( NRT, AB( KLU1-L, J1-KLM+L-1 ), INCA,
     $                            AB( KLU1-L+1, J1-KLM+L-1 ), INCA,
     $                            RWORK( J1 ), WORK( J1 ), KB1 )
   10          CONTINUE
<span class="comment">*</span><span class="comment">
</span>               IF( ML.GT.ML0 ) THEN
                  IF( ML.LE.M-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+ml-1,i)
</span><span class="comment">*</span><span class="comment">                    within the band, and apply rotation from the left
</span><span class="comment">*</span><span class="comment">
</span>                     CALL <a name="ZLARTG.247"></a><a href="zlartg.f.html#ZLARTG.1">ZLARTG</a>( AB( KU+ML-1, I ), AB( KU+ML, I ),
     $                            RWORK( I+ML-1 ), WORK( I+ML-1 ), RA )
                     AB( KU+ML-1, I ) = RA
                     IF( I.LT.N )
     $                  CALL <a name="ZROT.251"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( MIN( KU+ML-2, N-I ),
     $                             AB( KU+ML-2, I+1 ), LDAB-1,
     $                             AB( KU+ML-1, I+1 ), LDAB-1,
     $                             RWORK( I+ML-1 ), WORK( I+ML-1 ) )
                  END IF
                  NR = NR + 1
                  J1 = J1 - KB1
               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>                  DO 20 J = J1, J2, KB1
                     CALL <a name="ZROT.265"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( M, Q( 1, J-1 ), 1, Q( 1, J ), 1,
     $                          RWORK( J ), DCONJG( WORK( J ) ) )
   20             CONTINUE
               END IF
<span class="comment">*</span><span class="comment">
</span>               IF( WANTC ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 apply plane rotations to C
</span><span class="comment">*</span><span class="comment">
</span>                  DO 30 J = J1, J2, KB1
                     CALL <a name="ZROT.275"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( NCC, C( J-1, 1 ), LDC, C( J, 1 ), LDC,
     $                          RWORK( J ), WORK( J ) )
   30             CONTINUE
               END IF
<span class="comment">*</span><span class="comment">
</span>               IF( J2+KUN.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 - KB1
               END IF
<span class="comment">*</span><span class="comment">
</span>               DO 40 J = J1, J2, KB1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 create nonzero element a(j-1,j+ku) above the band
</span><span class="comment">*</span><span class="comment">                 and store it in WORK(n+1:2*n)
</span><span class="comment">*</span><span class="comment">
</span>                  WORK( J+KUN ) = WORK( J )*AB( 1, J+KUN )
                  AB( 1, J+KUN ) = RWORK( J )*AB( 1, J+KUN )
   40          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              generate plane rotations to annihilate nonzero elements
</span><span class="comment">*</span><span class="comment">              which have been generated above the band
</span><span class="comment">*</span><span class="comment">
</span>               IF( NR.GT.0 )
     $            CALL <a name="ZLARGV.301"></a><a href="zlargv.f.html#ZLARGV.1">ZLARGV</a>( NR, AB( 1, J1+KUN-1 ), INCA,
     $                         WORK( J1+KUN ), KB1, RWORK( 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