dgehrd.f.html

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

      EXTERNAL           <a name="ILAENV.131"></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 = MIN( NBMAX, <a name="ILAENV.138"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DGEHRD.138"></a><a href="dgehrd.f.html#DGEHRD.1">DGEHRD</a>'</span>, <span class="string">' '</span>, N, ILO, IHI, -1 ) )
      LWKOPT = N*NB
      WORK( 1 ) = LWKOPT
      LQUERY = ( LWORK.EQ.-1 )
      IF( N.LT.0 ) THEN
         INFO = -1
      ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
         INFO = -2
      ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -5
      ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
         INFO = -8
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.154"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DGEHRD.154"></a><a href="dgehrd.f.html#DGEHRD.1">DGEHRD</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">     Set elements 1:ILO-1 and IHI:N-1 of TAU to zero
</span><span class="comment">*</span><span class="comment">
</span>      DO 10 I = 1, ILO - 1
         TAU( I ) = ZERO
   10 CONTINUE
      DO 20 I = MAX( 1, IHI ), N - 1
         TAU( I ) = ZERO
   20 CONTINUE
<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>      NH = IHI - ILO + 1
      IF( NH.LE.1 ) THEN
         WORK( 1 ) = 1
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine the block size
</span><span class="comment">*</span><span class="comment">
</span>      NB = MIN( NBMAX, <a name="ILAENV.179"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DGEHRD.179"></a><a href="dgehrd.f.html#DGEHRD.1">DGEHRD</a>'</span>, <span class="string">' '</span>, N, ILO, IHI, -1 ) )
      NBMIN = 2
      IWS = 1
      IF( NB.GT.1 .AND. NB.LT.NH ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Determine when to cross over from blocked to unblocked code
</span><span class="comment">*</span><span class="comment">        (last block is always handled by unblocked code)
</span><span class="comment">*</span><span class="comment">
</span>         NX = MAX( NB, <a name="ILAENV.187"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 3, <span class="string">'<a name="DGEHRD.187"></a><a href="dgehrd.f.html#DGEHRD.1">DGEHRD</a>'</span>, <span class="string">' '</span>, N, ILO, IHI, -1 ) )
         IF( NX.LT.NH ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Determine if workspace is large enough for blocked code
</span><span class="comment">*</span><span class="comment">
</span>            IWS = N*NB
            IF( LWORK.LT.IWS ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Not enough workspace to use optimal NB:  determine the
</span><span class="comment">*</span><span class="comment">              minimum value of NB, and reduce NB or force use of
</span><span class="comment">*</span><span class="comment">              unblocked code
</span><span class="comment">*</span><span class="comment">
</span>               NBMIN = MAX( 2, <a name="ILAENV.199"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 2, <span class="string">'<a name="DGEHRD.199"></a><a href="dgehrd.f.html#DGEHRD.1">DGEHRD</a>'</span>, <span class="string">' '</span>, N, ILO, IHI,
     $                 -1 ) )
               IF( LWORK.GE.N*NBMIN ) THEN
                  NB = LWORK / N
               ELSE
                  NB = 1
               END IF
            END IF
         END IF
      END IF
      LDWORK = N
<span class="comment">*</span><span class="comment">
</span>      IF( NB.LT.NBMIN .OR. NB.GE.NH ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use unblocked code below
</span><span class="comment">*</span><span class="comment">
</span>         I = ILO
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use blocked code
</span><span class="comment">*</span><span class="comment">
</span>         DO 40 I = ILO, IHI - 1 - NX, NB
            IB = MIN( NB, IHI-I )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Reduce columns i:i+ib-1 to Hessenberg form, returning the
</span><span class="comment">*</span><span class="comment">           matrices V and T of the block reflector H = I - V*T*V'
</span><span class="comment">*</span><span class="comment">           which performs the reduction, and also the matrix Y = A*V*T
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DLAHR2.228"></a><a href="dlahr2.f.html#DLAHR2.1">DLAHR2</a>( IHI, I, IB, A( 1, I ), LDA, TAU( I ), T, LDT,
     $                   WORK, LDWORK )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Apply the block reflector H to A(1:ihi,i+ib:ihi) from the
</span><span class="comment">*</span><span class="comment">           right, computing  A := A - Y * V'. V(i+ib,ib-1) must be set
</span><span class="comment">*</span><span class="comment">           to 1
</span><span class="comment">*</span><span class="comment">
</span>            EI = A( I+IB, I+IB-1 )
            A( I+IB, I+IB-1 ) = ONE
            CALL DGEMM( <span class="string">'No transpose'</span>, <span class="string">'Transpose'</span>, 
     $                  IHI, IHI-I-IB+1,
     $                  IB, -ONE, WORK, LDWORK, A( I+IB, I ), LDA, ONE,
     $                  A( 1, I+IB ), LDA )
            A( I+IB, I+IB-1 ) = EI
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Apply the block reflector H to A(1:i,i+1:i+ib-1) from the
</span><span class="comment">*</span><span class="comment">           right
</span><span class="comment">*</span><span class="comment">
</span>            CALL DTRMM( <span class="string">'Right'</span>, <span class="string">'Lower'</span>, <span class="string">'Transpose'</span>,
     $                  <span class="string">'Unit'</span>, I, IB-1,
     $                  ONE, A( I+1, I ), LDA, WORK, LDWORK )
            DO 30 J = 0, IB-2
               CALL DAXPY( I, -ONE, WORK( LDWORK*J+1 ), 1,
     $                     A( 1, I+J+1 ), 1 )
   30       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Apply the block reflector H to A(i+1:ihi,i+ib:n) from the
</span><span class="comment">*</span><span class="comment">           left
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DLARFB.257"></a><a href="dlarfb.f.html#DLARFB.1">DLARFB</a>( <span class="string">'Left'</span>, <span class="string">'Transpose'</span>, <span class="string">'Forward'</span>,
     $                   <span class="string">'Columnwise'</span>,
     $                   IHI-I, N-I-IB+1, IB, A( I+1, I ), LDA, T, LDT,
     $                   A( I+1, I+IB ), LDA, WORK, LDWORK )
   40    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Use unblocked code to reduce the rest of the matrix
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="DGEHD2.266"></a><a href="dgehd2.f.html#DGEHD2.1">DGEHD2</a>( N, I, IHI, A, LDA, TAU, WORK, IINFO )
      WORK( 1 ) = IWS
<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="DGEHRD.271"></a><a href="dgehrd.f.html#DGEHRD.1">DGEHRD</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>