dstemr.f.html

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

</span><span class="comment">*</span><span class="comment">          The support of the eigenvectors in Z, i.e., the indices
</span><span class="comment">*</span><span class="comment">          indicating the nonzero elements in Z. The i-th computed eigenvector
</span><span class="comment">*</span><span class="comment">          is nonzero only in elements ISUPPZ( 2*i-1 ) through
</span><span class="comment">*</span><span class="comment">          ISUPPZ( 2*i ). This is relevant in the case when the matrix
</span><span class="comment">*</span><span class="comment">          is split. ISUPPZ is only accessed when JOBZ is 'V' and N &gt; 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  TRYRAC  (input/output) LOGICAL
</span><span class="comment">*</span><span class="comment">          If TRYRAC.EQ..TRUE., indicates that the code should check whether
</span><span class="comment">*</span><span class="comment">          the tridiagonal matrix defines its eigenvalues to high relative
</span><span class="comment">*</span><span class="comment">          accuracy.  If so, the code uses relative-accuracy preserving
</span><span class="comment">*</span><span class="comment">          algorithms that might be (a bit) slower depending on the matrix.
</span><span class="comment">*</span><span class="comment">          If the matrix does not define its eigenvalues to high relative
</span><span class="comment">*</span><span class="comment">          accuracy, the code can uses possibly faster algorithms.
</span><span class="comment">*</span><span class="comment">          If TRYRAC.EQ..FALSE., the code is not required to guarantee
</span><span class="comment">*</span><span class="comment">          relatively accurate eigenvalues and can use the fastest possible
</span><span class="comment">*</span><span class="comment">          techniques.
</span><span class="comment">*</span><span class="comment">          On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix
</span><span class="comment">*</span><span class="comment">          does not define its eigenvalues to high relative accuracy.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
</span><span class="comment">*</span><span class="comment">          On exit, if INFO = 0, WORK(1) returns the optimal
</span><span class="comment">*</span><span class="comment">          (and minimal) LWORK.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LWORK   (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The dimension of the array WORK. LWORK &gt;= max(1,18*N)
</span><span class="comment">*</span><span class="comment">          if JOBZ = 'V', and LWORK &gt;= max(1,12*N) if JOBZ = 'N'.
</span><span class="comment">*</span><span class="comment">          If LWORK = -1, then a workspace query is assumed; the routine
</span><span class="comment">*</span><span class="comment">          only calculates the optimal size of the WORK array, returns
</span><span class="comment">*</span><span class="comment">          this value as the first entry of the WORK array, and no error
</span><span class="comment">*</span><span class="comment">          message related to LWORK is issued by <a name="XERBLA.179"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IWORK   (workspace/output) INTEGER array, dimension (LIWORK)
</span><span class="comment">*</span><span class="comment">          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LIWORK  (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The dimension of the array IWORK.  LIWORK &gt;= max(1,10*N)
</span><span class="comment">*</span><span class="comment">          if the eigenvectors are desired, and LIWORK &gt;= max(1,8*N)
</span><span class="comment">*</span><span class="comment">          if only the eigenvalues are to be computed.
</span><span class="comment">*</span><span class="comment">          If LIWORK = -1, then a workspace query is assumed; the
</span><span class="comment">*</span><span class="comment">          routine only calculates the optimal size of the IWORK array,
</span><span class="comment">*</span><span class="comment">          returns this value as the first entry of the IWORK array, and
</span><span class="comment">*</span><span class="comment">          no error message related to LIWORK is issued by <a name="XERBLA.191"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          On exit, INFO
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit
</span><span class="comment">*</span><span class="comment">          &lt; 0:  if INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment">          &gt; 0:  if INFO = 1X, internal error in <a name="DLARRE.197"></a><a href="dlarre.f.html#DLARRE.1">DLARRE</a>,
</span><span class="comment">*</span><span class="comment">                if INFO = 2X, internal error in <a name="DLARRV.198"></a><a href="dlarrv.f.html#DLARRV.1">DLARRV</a>.
</span><span class="comment">*</span><span class="comment">                Here, the digit X = ABS( IINFO ) &lt; 10, where IINFO is
</span><span class="comment">*</span><span class="comment">                the nonzero error code returned by <a name="DLARRE.200"></a><a href="dlarre.f.html#DLARRE.1">DLARRE</a> or
</span><span class="comment">*</span><span class="comment">                <a name="DLARRV.201"></a><a href="dlarrv.f.html#DLARRV.1">DLARRV</a>, respectively.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Based on contributions by
</span><span class="comment">*</span><span class="comment">     Beresford Parlett, University of California, Berkeley, USA
</span><span class="comment">*</span><span class="comment">     Jim Demmel, University of California, Berkeley, USA
</span><span class="comment">*</span><span class="comment">     Inderjit Dhillon, University of Texas, Austin, USA
</span><span class="comment">*</span><span class="comment">     Osni Marques, LBNL/NERSC, USA
</span><span class="comment">*</span><span class="comment">     Christof Voemel, University of California, Berkeley, USA
</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>      DOUBLE PRECISION   ZERO, ONE, FOUR, MINRGP
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0,
     $                     FOUR = 4.0D0,
     $                     MINRGP = 1.0D-3 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY
      INTEGER            I, IBEGIN, IEND, IFIRST, IIL, IINDBL, IINDW,
     $                   IINDWK, IINFO, IINSPL, IIU, ILAST, IN, INDD,
     $                   INDE2, INDERR, INDGP, INDGRS, INDWRK, ITMP,
     $                   ITMP2, J, JBLK, JJ, LIWMIN, LWMIN, NSPLIT,
     $                   NZCMIN, OFFSET, WBEGIN, WEND
      DOUBLE PRECISION   BIGNUM, CS, EPS, PIVMIN, R1, R2, RMAX, RMIN,
     $                   RTOL1, RTOL2, SAFMIN, SCALE, SMLNUM, SN,
     $                   THRESH, TMP, TNRM, WL, WU
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.235"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      DOUBLE PRECISION   <a name="DLAMCH.236"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANST.236"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>
      EXTERNAL           <a name="LSAME.237"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="DLAMCH.237"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANST.237"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           DCOPY, <a name="DLAE2.240"></a><a href="dlae2.f.html#DLAE2.1">DLAE2</a>, <a name="DLAEV2.240"></a><a href="dlaev2.f.html#DLAEV2.1">DLAEV2</a>, <a name="DLARRC.240"></a><a href="dlarrc.f.html#DLARRC.1">DLARRC</a>, <a name="DLARRE.240"></a><a href="dlarre.f.html#DLARRE.1">DLARRE</a>, <a name="DLARRJ.240"></a><a href="dlarrj.f.html#DLARRJ.1">DLARRJ</a>,
     $                   <a name="DLARRR.241"></a><a href="dlarrr.f.html#DLARRR.1">DLARRR</a>, <a name="DLARRV.241"></a><a href="dlarrv.f.html#DLARRV.1">DLARRV</a>, <a name="DLASRT.241"></a><a href="dlasrt.f.html#DLASRT.1">DLASRT</a>, DSCAL, DSWAP, <a name="XERBLA.241"></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, SQRT


<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>      WANTZ = <a name="LSAME.252"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBZ, <span class="string">'V'</span> )
      ALLEIG = <a name="LSAME.253"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'A'</span> )
      VALEIG = <a name="LSAME.254"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'V'</span> )
      INDEIG = <a name="LSAME.255"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'I'</span> )
<span class="comment">*</span><span class="comment">
</span>      LQUERY = ( ( LWORK.EQ.-1 ).OR.( LIWORK.EQ.-1 ) )
      ZQUERY = ( NZC.EQ.-1 )
      TRYRAC = ( INFO.NE.0 )

<span class="comment">*</span><span class="comment">     <a name="DSTEMR.261"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a> needs WORK of size 6*N, IWORK of size 3*N.
</span><span class="comment">*</span><span class="comment">     In addition, <a name="DLARRE.262"></a><a href="dlarre.f.html#DLARRE.1">DLARRE</a> needs WORK of size 6*N, IWORK of size 5*N.
</span><span class="comment">*</span><span class="comment">     Furthermore, <a name="DLARRV.263"></a><a href="dlarrv.f.html#DLARRV.1">DLARRV</a> needs WORK of size 12*N, IWORK of size 7*N.
</span>      IF( WANTZ ) THEN
         LWMIN = 18*N
         LIWMIN = 10*N
      ELSE
<span class="comment">*</span><span class="comment">        need less workspace if only the eigenvalues are wanted
</span>         LWMIN = 12*N
         LIWMIN = 8*N
      ENDIF

      WL = ZERO
      WU = ZERO
      IIL = 0
      IIU = 0

      IF( VALEIG ) THEN
<span class="comment">*</span><span class="comment">        We do not reference VL, VU in the cases RANGE = 'I','A'
</span><span class="comment">*</span><span class="comment">        The interval (WL, WU] contains all the wanted eigenvalues.
</span><span class="comment">*</span><span class="comment">        It is either given by the user or computed in <a name="DLARRE.281"></a><a href="dlarre.f.html#DLARRE.1">DLARRE</a>.
</span>         WL = VL
         WU = VU
      ELSEIF( INDEIG ) THEN
<span class="comment">*</span><span class="comment">        We do not reference IL, IU in the cases RANGE = 'V','A'
</span>         IIL = IL
         IIU = IU
      ENDIF
<span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( .NOT.( WANTZ .OR. <a name="LSAME.291"></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( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( VALEIG .AND. N.GT.0 .AND. WU.LE.WL ) THEN
         INFO = -7
      ELSE IF( INDEIG .AND. ( IIL.LT.1 .OR. IIL.GT.N ) ) THEN
         INFO = -8
      ELSE IF( INDEIG .AND. ( IIU.LT.IIL .OR. IIU.GT.N ) ) THEN
         INFO = -9
      ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
         INFO = -13
      ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
         INFO = -17
      ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
         INFO = -19
      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.313"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Safe minimum'</span> )
      EPS = <a name="DLAMCH.314"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Precision'</span> )
      SMLNUM = SAFMIN / EPS
      BIGNUM = ONE / SMLNUM
      RMIN = SQRT( SMLNUM )