zheevr.f.html

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

</span>      LRWMIN = MAX( 1, 24*N )
      LIWMIN = MAX( 1, 10*N )
      LWMIN = MAX( 1, 2*N )
<span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( .NOT.( WANTZ .OR. <a name="LSAME.295"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBZ, <span class="string">'N'</span> ) ) ) THEN
         INFO = -1
      ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
         INFO = -2
      ELSE IF( .NOT.( LOWER .OR. <a name="LSAME.299"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> ) ) ) THEN
         INFO = -3
      ELSE IF( N.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -6
      ELSE
         IF( VALEIG ) THEN
            IF( N.GT.0 .AND. VU.LE.VL )
     $         INFO = -8
         ELSE IF( INDEIG ) THEN
            IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
               INFO = -9
            ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
               INFO = -10
            END IF
         END IF
      END IF
      IF( INFO.EQ.0 ) THEN
         IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
            INFO = -15
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.EQ.0 ) THEN
         NB = <a name="ILAENV.324"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZHETRD.324"></a><a href="zhetrd.f.html#ZHETRD.1">ZHETRD</a>'</span>, UPLO, N, -1, -1, -1 )
         NB = MAX( NB, <a name="ILAENV.325"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNMTR.325"></a><a href="zunmtr.f.html#ZUNMTR.1">ZUNMTR</a>'</span>, UPLO, N, -1, -1, -1 ) )
         LWKOPT = MAX( ( NB+1 )*N, LWMIN )
         WORK( 1 ) = LWKOPT
         RWORK( 1 ) = LRWMIN
         IWORK( 1 ) = LIWMIN
<span class="comment">*</span><span class="comment">
</span>         IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
            INFO = -18
         ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
            INFO = -20
         ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
            INFO = -22
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.341"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZHEEVR.341"></a><a href="zheevr.f.html#ZHEEVR.1">ZHEEVR</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>      M = 0
      IF( N.EQ.0 ) THEN
         WORK( 1 ) = 1
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.1 ) THEN
         WORK( 1 ) = 2
         IF( ALLEIG .OR. INDEIG ) THEN
            M = 1
            W( 1 ) = DBLE( A( 1, 1 ) )
         ELSE
            IF( VL.LT.DBLE( A( 1, 1 ) ) .AND. VU.GE.DBLE( A( 1, 1 ) ) )
     $           THEN
               M = 1
               W( 1 ) = DBLE( A( 1, 1 ) )
            END IF
         END IF
         IF( WANTZ )
     $      Z( 1, 1 ) = ONE
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Get machine constants.
</span><span class="comment">*</span><span class="comment">
</span>      SAFMIN = <a name="DLAMCH.374"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Safe minimum'</span> )
      EPS = <a name="DLAMCH.375"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Precision'</span> )
      SMLNUM = SAFMIN / EPS
      BIGNUM = ONE / SMLNUM
      RMIN = SQRT( SMLNUM )
      RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale matrix to allowable range, if necessary.
</span><span class="comment">*</span><span class="comment">
</span>      ISCALE = 0
      ABSTLL = ABSTOL
      IF (VALEIG) THEN
         VLL = VL
         VUU = VU
      END IF
      ANRM = <a name="ZLANSY.389"></a><a href="zlansy.f.html#ZLANSY.1">ZLANSY</a>( <span class="string">'M'</span>, UPLO, N, A, LDA, RWORK )
      IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
         ISCALE = 1
         SIGMA = RMIN / ANRM
      ELSE IF( ANRM.GT.RMAX ) THEN
         ISCALE = 1
         SIGMA = RMAX / ANRM
      END IF
      IF( ISCALE.EQ.1 ) THEN
         IF( LOWER ) THEN
            DO 10 J = 1, N
               CALL ZDSCAL( N-J+1, SIGMA, A( J, J ), 1 )
   10       CONTINUE
         ELSE
            DO 20 J = 1, N
               CALL ZDSCAL( J, SIGMA, A( 1, J ), 1 )
   20       CONTINUE
         END IF
         IF( ABSTOL.GT.0 )
     $      ABSTLL = ABSTOL*SIGMA
         IF( VALEIG ) THEN
            VLL = VL*SIGMA
            VUU = VU*SIGMA
         END IF
      END IF

<span class="comment">*</span><span class="comment">     Initialize indices into workspaces.  Note: The IWORK indices are
</span><span class="comment">*</span><span class="comment">     used only if <a name="DSTERF.416"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a> or <a name="ZSTEMR.416"></a><a href="zstemr.f.html#ZSTEMR.1">ZSTEMR</a> fail.
</span>
<span class="comment">*</span><span class="comment">     WORK(INDTAU:INDTAU+N-1) stores the complex scalar factors of the
</span><span class="comment">*</span><span class="comment">     elementary reflectors used in <a name="ZHETRD.419"></a><a href="zhetrd.f.html#ZHETRD.1">ZHETRD</a>.
</span>      INDTAU = 1
<span class="comment">*</span><span class="comment">     INDWK is the starting offset of the remaining complex workspace,
</span><span class="comment">*</span><span class="comment">     and LLWORK is the remaining complex workspace size.
</span>      INDWK = INDTAU + N
      LLWORK = LWORK - INDWK + 1

<span class="comment">*</span><span class="comment">     RWORK(INDRD:INDRD+N-1) stores the real tridiagonal's diagonal
</span><span class="comment">*</span><span class="comment">     entries.
</span>      INDRD = 1
<span class="comment">*</span><span class="comment">     RWORK(INDRE:INDRE+N-1) stores the off-diagonal entries of the
</span><span class="comment">*</span><span class="comment">     tridiagonal matrix from <a name="ZHETRD.430"></a><a href="zhetrd.f.html#ZHETRD.1">ZHETRD</a>.
</span>      INDRE = INDRD + N
<span class="comment">*</span><span class="comment">     RWORK(INDRDD:INDRDD+N-1) is a copy of the diagonal entries over
</span><span class="comment">*</span><span class="comment">     -written by <a name="ZSTEMR.433"></a><a href="zstemr.f.html#ZSTEMR.1">ZSTEMR</a> (the <a name="DSTERF.433"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a> path copies the diagonal to W).
</span>      INDRDD = INDRE + N
<span class="comment">*</span><span class="comment">     RWORK(INDREE:INDREE+N-1) is a copy of the off-diagonal entries over
</span><span class="comment">*</span><span class="comment">     -written while computing the eigenvalues in <a name="DSTERF.436"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a> and <a name="ZSTEMR.436"></a><a href="zstemr.f.html#ZSTEMR.1">ZSTEMR</a>.
</span>      INDREE = INDRDD + N
<span class="comment">*</span><span class="comment">     INDRWK is the starting offset of the left-over real workspace, and
</span><span class="comment">*</span><span class="comment">     LLRWORK is the remaining workspace size.
</span>      INDRWK = INDREE + N
      LLRWORK = LRWORK - INDRWK + 1

<span class="comment">*</span><span class="comment">     IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in <a name="DSTEBZ.443"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a> and