sgebrd.f.html

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

</span><span class="comment">*</span><span class="comment">    (  v1  v2  v3  v4  d  )           (  v1  v2  v3  e   d   u5 )
</span><span class="comment">*</span><span class="comment">    (  v1  v2  v3  v4  v5 )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where d and e denote diagonal and off-diagonal elements of B, vi
</span><span class="comment">*</span><span class="comment">  denotes an element of the vector defining H(i), and ui an element of
</span><span class="comment">*</span><span class="comment">  the vector defining G(i).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               ONE
      PARAMETER          ( ONE = 1.0E+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            LQUERY
      INTEGER            I, IINFO, J, LDWRKX, LDWRKY, LWKOPT, MINMN, NB,
     $                   NBMIN, NX
      REAL               WS
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="SGEBD2.149"></a><a href="sgebd2.f.html#SGEBD2.1">SGEBD2</a>, SGEMM, <a name="SLABRD.149"></a><a href="slabrd.f.html#SLABRD.1">SLABRD</a>, <a name="XERBLA.149"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          MAX, MIN, REAL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      INTEGER            <a name="ILAENV.155"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
      EXTERNAL           <a name="ILAENV.156"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input parameters
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      NB = MAX( 1, <a name="ILAENV.163"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="SGEBRD.163"></a><a href="sgebrd.f.html#SGEBRD.1">SGEBRD</a>'</span>, <span class="string">' '</span>, M, N, -1, -1 ) )
      LWKOPT = ( M+N )*NB
      WORK( 1 ) = REAL( LWKOPT )
      LQUERY = ( LWORK.EQ.-1 )
      IF( M.LT.0 ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -4
      ELSE IF( LWORK.LT.MAX( 1, M, N ) .AND. .NOT.LQUERY ) THEN
         INFO = -10
      END IF
      IF( INFO.LT.0 ) THEN
         CALL <a name="XERBLA.177"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SGEBRD.177"></a><a href="sgebrd.f.html#SGEBRD.1">SGEBRD</a>'</span>, -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
<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>      MINMN = MIN( M, N )
      IF( MINMN.EQ.0 ) THEN
         WORK( 1 ) = 1
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      WS = MAX( M, N )
      LDWRKX = M
      LDWRKY = N
<span class="comment">*</span><span class="comment">
</span>      IF( NB.GT.1 .AND. NB.LT.MINMN ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Set the crossover point NX.
</span><span class="comment">*</span><span class="comment">
</span>         NX = MAX( NB, <a name="ILAENV.199"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 3, <span class="string">'<a name="SGEBRD.199"></a><a href="sgebrd.f.html#SGEBRD.1">SGEBRD</a>'</span>, <span class="string">' '</span>, M, N, -1, -1 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Determine when to switch from blocked to unblocked code.
</span><span class="comment">*</span><span class="comment">
</span>         IF( NX.LT.MINMN ) THEN
            WS = ( M+N )*NB
            IF( LWORK.LT.WS ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Not enough work space for the optimal NB, consider using
</span><span class="comment">*</span><span class="comment">              a smaller block size.
</span><span class="comment">*</span><span class="comment">
</span>               NBMIN = <a name="ILAENV.210"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 2, <span class="string">'<a name="SGEBRD.210"></a><a href="sgebrd.f.html#SGEBRD.1">SGEBRD</a>'</span>, <span class="string">' '</span>, M, N, -1, -1 )
               IF( LWORK.GE.( M+N )*NBMIN ) THEN
                  NB = LWORK / ( M+N )
               ELSE
                  NB = 1
                  NX = MINMN
               END IF
            END IF
         END IF
      ELSE
         NX = MINMN
      END IF
<span class="comment">*</span><span class="comment">
</span>      DO 30 I = 1, MINMN - NX, NB
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Reduce rows and columns i:i+nb-1 to bidiagonal form and return
</span><span class="comment">*</span><span class="comment">        the matrices X and Y which are needed to update the unreduced
</span><span class="comment">*</span><span class="comment">        part of the matrix
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLABRD.229"></a><a href="slabrd.f.html#SLABRD.1">SLABRD</a>( M-I+1, N-I+1, NB, A( I, I ), LDA, D( I ), E( I ),
     $                TAUQ( I ), TAUP( I ), WORK, LDWRKX,
     $                WORK( LDWRKX*NB+1 ), LDWRKY )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Update the trailing submatrix A(i+nb:m,i+nb:n), using an update
</span><span class="comment">*</span><span class="comment">        of the form  A := A - V*Y' - X*U'
</span><span class="comment">*</span><span class="comment">
</span>         CALL SGEMM( <span class="string">'No transpose'</span>, <span class="string">'Transpose'</span>, M-I-NB+1, N-I-NB+1,
     $               NB, -ONE, A( I+NB, I ), LDA,
     $               WORK( LDWRKX*NB+NB+1 ), LDWRKY, ONE,
     $               A( I+NB, I+NB ), LDA )
         CALL SGEMM( <span class="string">'No transpose'</span>, <span class="string">'No transpose'</span>, M-I-NB+1, N-I-NB+1,
     $               NB, -ONE, WORK( NB+1 ), LDWRKX, A( I, I+NB ), LDA,
     $               ONE, A( I+NB, I+NB ), LDA )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Copy diagonal and off-diagonal elements of B back into A
</span><span class="comment">*</span><span class="comment">
</span>         IF( M.GE.N ) THEN
            DO 10 J = I, I + NB - 1
               A( J, J ) = D( J )
               A( J, J+1 ) = E( J )
   10       CONTINUE
         ELSE
            DO 20 J = I, I + NB - 1
               A( J, J ) = D( J )
               A( J+1, J ) = E( J )
   20       CONTINUE
         END IF
   30 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Use unblocked code to reduce the remainder of the matrix
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="SGEBD2.261"></a><a href="sgebd2.f.html#SGEBD2.1">SGEBD2</a>( M-I+1, N-I+1, A( I, I ), LDA, D( I ), E( I ),
     $             TAUQ( I ), TAUP( I ), WORK, IINFO )
      WORK( 1 ) = WS
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SGEBRD.266"></a><a href="sgebrd.f.html#SGEBRD.1">SGEBRD</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>