chetrd.f.html

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

</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          MAX
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.147"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            <a name="ILAENV.148"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
      EXTERNAL           <a name="LSAME.149"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="ILAENV.149"></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
      UPPER = <a name="LSAME.156"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> )
      LQUERY = ( LWORK.EQ.-1 )
      IF( .NOT.UPPER .AND. .NOT.<a name="LSAME.158"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'L'</span> ) ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -4
      ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
         INFO = -9
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.EQ.0 ) THEN
<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 = <a name="ILAENV.172"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="CHETRD.172"></a><a href="chetrd.f.html#CHETRD.1">CHETRD</a>'</span>, UPLO, N, -1, -1, -1 )
         LWKOPT = N*NB
         WORK( 1 ) = LWKOPT
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.178"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CHETRD.178"></a><a href="chetrd.f.html#CHETRD.1">CHETRD</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>      IF( N.EQ.0 ) THEN
         WORK( 1 ) = 1
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      NX = N
      IWS = 1
      IF( NB.GT.1 .AND. NB.LT.N ) 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.198"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 3, <span class="string">'<a name="CHETRD.198"></a><a href="chetrd.f.html#CHETRD.1">CHETRD</a>'</span>, UPLO, N, -1, -1, -1 ) )
         IF( NX.LT.N ) 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>            LDWORK = N
            IWS = LDWORK*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 by setting NX = N.
</span><span class="comment">*</span><span class="comment">
</span>               NB = MAX( LWORK / LDWORK, 1 )
               NBMIN = <a name="ILAENV.212"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 2, <span class="string">'<a name="CHETRD.212"></a><a href="chetrd.f.html#CHETRD.1">CHETRD</a>'</span>, UPLO, N, -1, -1, -1 )
               IF( NB.LT.NBMIN )
     $            NX = N
            END IF
         ELSE
            NX = N
         END IF
      ELSE
         NB = 1
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( UPPER ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Reduce the upper triangle of A.
</span><span class="comment">*</span><span class="comment">        Columns 1:kk are handled by the unblocked method.
</span><span class="comment">*</span><span class="comment">
</span>         KK = N - ( ( N-NX+NB-1 ) / NB )*NB
         DO 20 I = N - NB + 1, KK + 1, -NB
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Reduce columns i:i+nb-1 to tridiagonal form and form the
</span><span class="comment">*</span><span class="comment">           matrix W which is needed to update the unreduced part of
</span><span class="comment">*</span><span class="comment">           the matrix
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="CLATRD.235"></a><a href="clatrd.f.html#CLATRD.1">CLATRD</a>( UPLO, I+NB-1, NB, A, LDA, E, TAU, WORK,
     $                   LDWORK )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update the unreduced submatrix A(1:i-1,1:i-1), using an
</span><span class="comment">*</span><span class="comment">           update of the form:  A := A - V*W' - W*V'
</span><span class="comment">*</span><span class="comment">
</span>            CALL CHER2K( UPLO, <span class="string">'No transpose'</span>, I-1, NB, -CONE,
     $                   A( 1, I ), LDA, WORK, LDWORK, ONE, A, LDA )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Copy superdiagonal elements back into A, and diagonal
</span><span class="comment">*</span><span class="comment">           elements into D
</span><span class="comment">*</span><span class="comment">
</span>            DO 10 J = I, I + NB - 1
               A( J-1, J ) = E( J-1 )
               D( J ) = A( J, J )
   10       CONTINUE
   20    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use unblocked code to reduce the last or only block
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CHETD2.255"></a><a href="chetd2.f.html#CHETD2.1">CHETD2</a>( UPLO, KK, A, LDA, D, E, TAU, IINFO )
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Reduce the lower triangle of A
</span><span class="comment">*</span><span class="comment">
</span>         DO 40 I = 1, N - NX, NB
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Reduce columns i:i+nb-1 to tridiagonal form and form the
</span><span class="comment">*</span><span class="comment">           matrix W which is needed to update the unreduced part of
</span><span class="comment">*</span><span class="comment">           the matrix
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="CLATRD.266"></a><a href="clatrd.f.html#CLATRD.1">CLATRD</a>( UPLO, N-I+1, NB, A( I, I ), LDA, E( I ),
     $                   TAU( I ), WORK, LDWORK )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update the unreduced submatrix A(i+nb:n,i+nb:n), using
</span><span class="comment">*</span><span class="comment">           an update of the form:  A := A - V*W' - W*V'
</span><span class="comment">*</span><span class="comment">
</span>            CALL CHER2K( UPLO, <span class="string">'No transpose'</span>, N-I-NB+1, NB, -CONE,
     $                   A( I+NB, I ), LDA, WORK( NB+1 ), LDWORK, ONE,
     $                   A( I+NB, I+NB ), LDA )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Copy subdiagonal elements back into A, and diagonal
</span><span class="comment">*</span><span class="comment">           elements into D
</span><span class="comment">*</span><span class="comment">
</span>            DO 30 J = I, I + NB - 1
               A( J+1, J ) = E( J )
               D( J ) = A( J, J )
   30       CONTINUE
   40    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use unblocked code to reduce the last or only block
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CHETD2.287"></a><a href="chetd2.f.html#CHETD2.1">CHETD2</a>( UPLO, N-I+1, A( I, I ), LDA, D( I ), E( I ),
     $                TAU( I ), IINFO )
      END IF
<span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = LWKOPT
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="CHETRD.294"></a><a href="chetrd.f.html#CHETRD.1">CHETRD</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>