[BACK]Return to dstebz.f CVS log [TXT] [DIR] Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/dstebz.f between versions 1.8 and 1.9

version 1.8, 2011/07/22 07:38:11 version 1.9, 2011/11/21 20:43:04
Line 1 Line 1
   *> \brief \b DSTEBZ
   *
   *  =========== DOCUMENTATION ===========
   *
   * Online html documentation available at 
   *            http://www.netlib.org/lapack/explore-html/ 
   *
   *> \htmlonly
   *> Download DSTEBZ + dependencies 
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstebz.f"> 
   *> [TGZ]</a> 
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstebz.f"> 
   *> [ZIP]</a> 
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstebz.f"> 
   *> [TXT]</a>
   *> \endhtmlonly 
   *
   *  Definition:
   *  ===========
   *
   *       SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E,
   *                          M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK,
   *                          INFO )
   * 
   *       .. Scalar Arguments ..
   *       CHARACTER          ORDER, RANGE
   *       INTEGER            IL, INFO, IU, M, N, NSPLIT
   *       DOUBLE PRECISION   ABSTOL, VL, VU
   *       ..
   *       .. Array Arguments ..
   *       INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * )
   *       DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )
   *       ..
   *  
   *
   *> \par Purpose:
   *  =============
   *>
   *> \verbatim
   *>
   *> DSTEBZ computes the eigenvalues of a symmetric tridiagonal
   *> matrix T.  The user may ask for all eigenvalues, all eigenvalues
   *> in the half-open interval (VL, VU], or the IL-th through IU-th
   *> eigenvalues.
   *>
   *> To avoid overflow, the matrix must be scaled so that its
   *> largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
   *> accuracy, it should not be much smaller than that.
   *>
   *> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
   *> Matrix", Report CS41, Computer Science Dept., Stanford
   *> University, July 21, 1966.
   *> \endverbatim
   *
   *  Arguments:
   *  ==========
   *
   *> \param[in] RANGE
   *> \verbatim
   *>          RANGE is CHARACTER*1
   *>          = 'A': ("All")   all eigenvalues will be found.
   *>          = 'V': ("Value") all eigenvalues in the half-open interval
   *>                           (VL, VU] will be found.
   *>          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
   *>                           entire matrix) will be found.
   *> \endverbatim
   *>
   *> \param[in] ORDER
   *> \verbatim
   *>          ORDER is CHARACTER*1
   *>          = 'B': ("By Block") the eigenvalues will be grouped by
   *>                              split-off block (see IBLOCK, ISPLIT) and
   *>                              ordered from smallest to largest within
   *>                              the block.
   *>          = 'E': ("Entire matrix")
   *>                              the eigenvalues for the entire matrix
   *>                              will be ordered from smallest to
   *>                              largest.
   *> \endverbatim
   *>
   *> \param[in] N
   *> \verbatim
   *>          N is INTEGER
   *>          The order of the tridiagonal matrix T.  N >= 0.
   *> \endverbatim
   *>
   *> \param[in] VL
   *> \verbatim
   *>          VL is DOUBLE PRECISION
   *> \endverbatim
   *>
   *> \param[in] VU
   *> \verbatim
   *>          VU is DOUBLE PRECISION
   *>
   *>          If RANGE='V', the lower and upper bounds of the interval to
   *>          be searched for eigenvalues.  Eigenvalues less than or equal
   *>          to VL, or greater than VU, will not be returned.  VL < VU.
   *>          Not referenced if RANGE = 'A' or 'I'.
   *> \endverbatim
   *>
   *> \param[in] IL
   *> \verbatim
   *>          IL is INTEGER
   *> \endverbatim
   *>
   *> \param[in] IU
   *> \verbatim
   *>          IU is INTEGER
   *>
   *>          If RANGE='I', the indices (in ascending order) of the
   *>          smallest and largest eigenvalues to be returned.
   *>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
   *>          Not referenced if RANGE = 'A' or 'V'.
   *> \endverbatim
   *>
   *> \param[in] ABSTOL
   *> \verbatim
   *>          ABSTOL is DOUBLE PRECISION
   *>          The absolute tolerance for the eigenvalues.  An eigenvalue
   *>          (or cluster) is considered to be located if it has been
   *>          determined to lie in an interval whose width is ABSTOL or
   *>          less.  If ABSTOL is less than or equal to zero, then ULP*|T|
   *>          will be used, where |T| means the 1-norm of T.
   *>
   *>          Eigenvalues will be computed most accurately when ABSTOL is
   *>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
   *> \endverbatim
   *>
   *> \param[in] D
   *> \verbatim
   *>          D is DOUBLE PRECISION array, dimension (N)
   *>          The n diagonal elements of the tridiagonal matrix T.
   *> \endverbatim
   *>
   *> \param[in] E
   *> \verbatim
   *>          E is DOUBLE PRECISION array, dimension (N-1)
   *>          The (n-1) off-diagonal elements of the tridiagonal matrix T.
   *> \endverbatim
   *>
   *> \param[out] M
   *> \verbatim
   *>          M is INTEGER
   *>          The actual number of eigenvalues found. 0 <= M <= N.
   *>          (See also the description of INFO=2,3.)
   *> \endverbatim
   *>
   *> \param[out] NSPLIT
   *> \verbatim
   *>          NSPLIT is INTEGER
   *>          The number of diagonal blocks in the matrix T.
   *>          1 <= NSPLIT <= N.
   *> \endverbatim
   *>
   *> \param[out] W
   *> \verbatim
   *>          W is DOUBLE PRECISION array, dimension (N)
   *>          On exit, the first M elements of W will contain the
   *>          eigenvalues.  (DSTEBZ may use the remaining N-M elements as
   *>          workspace.)
   *> \endverbatim
   *>
   *> \param[out] IBLOCK
   *> \verbatim
   *>          IBLOCK is INTEGER array, dimension (N)
   *>          At each row/column j where E(j) is zero or small, the
   *>          matrix T is considered to split into a block diagonal
   *>          matrix.  On exit, if INFO = 0, IBLOCK(i) specifies to which
   *>          block (from 1 to the number of blocks) the eigenvalue W(i)
   *>          belongs.  (DSTEBZ may use the remaining N-M elements as
   *>          workspace.)
   *> \endverbatim
   *>
   *> \param[out] ISPLIT
   *> \verbatim
   *>          ISPLIT is INTEGER array, dimension (N)
   *>          The splitting points, at which T breaks up into submatrices.
   *>          The first submatrix consists of rows/columns 1 to ISPLIT(1),
   *>          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
   *>          etc., and the NSPLIT-th consists of rows/columns
   *>          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
   *>          (Only the first NSPLIT elements will actually be used, but
   *>          since the user cannot know a priori what value NSPLIT will
   *>          have, N words must be reserved for ISPLIT.)
   *> \endverbatim
   *>
   *> \param[out] WORK
   *> \verbatim
   *>          WORK is DOUBLE PRECISION array, dimension (4*N)
   *> \endverbatim
   *>
   *> \param[out] IWORK
   *> \verbatim
   *>          IWORK is INTEGER array, dimension (3*N)
   *> \endverbatim
   *>
   *> \param[out] INFO
   *> \verbatim
   *>          INFO is INTEGER
   *>          = 0:  successful exit
   *>          < 0:  if INFO = -i, the i-th argument had an illegal value
   *>          > 0:  some or all of the eigenvalues failed to converge or
   *>                were not computed:
   *>                =1 or 3: Bisection failed to converge for some
   *>                        eigenvalues; these eigenvalues are flagged by a
   *>                        negative block number.  The effect is that the
   *>                        eigenvalues may not be as accurate as the
   *>                        absolute and relative tolerances.  This is
   *>                        generally caused by unexpectedly inaccurate
   *>                        arithmetic.
   *>                =2 or 3: RANGE='I' only: Not all of the eigenvalues
   *>                        IL:IU were found.
   *>                        Effect: M < IU+1-IL
   *>                        Cause:  non-monotonic arithmetic, causing the
   *>                                Sturm sequence to be non-monotonic.
   *>                        Cure:   recalculate, using RANGE='A', and pick
   *>                                out eigenvalues IL:IU.  In some cases,
   *>                                increasing the PARAMETER "FUDGE" may
   *>                                make things work.
   *>                = 4:    RANGE='I', and the Gershgorin interval
   *>                        initially used was too small.  No eigenvalues
   *>                        were computed.
   *>                        Probable cause: your machine has sloppy
   *>                                        floating-point arithmetic.
   *>                        Cure: Increase the PARAMETER "FUDGE",
   *>                              recompile, and try again.
   *> \endverbatim
   *
   *> \par Internal Parameters:
   *  =========================
   *>
   *> \verbatim
   *>  RELFAC  DOUBLE PRECISION, default = 2.0e0
   *>          The relative tolerance.  An interval (a,b] lies within
   *>          "relative tolerance" if  b-a < RELFAC*ulp*max(|a|,|b|),
   *>          where "ulp" is the machine precision (distance from 1 to
   *>          the next larger floating point number.)
   *>
   *>  FUDGE   DOUBLE PRECISION, default = 2
   *>          A "fudge factor" to widen the Gershgorin intervals.  Ideally,
   *>          a value of 1 should work, but on machines with sloppy
   *>          arithmetic, this needs to be larger.  The default for
   *>          publicly released versions should be large enough to handle
   *>          the worst machine around.  Note that this has no effect
   *>          on accuracy of the solution.
   *> \endverbatim
   *
   *  Authors:
   *  ========
   *
   *> \author Univ. of Tennessee 
   *> \author Univ. of California Berkeley 
   *> \author Univ. of Colorado Denver 
   *> \author NAG Ltd. 
   *
   *> \date November 2011
   *
   *> \ingroup auxOTHERcomputational
   *
   *  =====================================================================
       SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E,        SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E,
      $                   M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK,       $                   M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK,
      $                   INFO )       $                   INFO )
 *  *
 *  -- LAPACK routine (version 3.3.1) --  *  -- LAPACK computational routine (version 3.4.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 *  -- April 2011                                                      --  *     November 2011
 *     8-18-00:  Increase FUDGE factor for T3E (eca)  
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          ORDER, RANGE        CHARACTER          ORDER, RANGE
Line 18 Line 278
       DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )        DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )
 *     ..  *     ..
 *  *
 *  Purpose  
 *  =======  
 *  
 *  DSTEBZ computes the eigenvalues of a symmetric tridiagonal  
 *  matrix T.  The user may ask for all eigenvalues, all eigenvalues  
 *  in the half-open interval (VL, VU], or the IL-th through IU-th  
 *  eigenvalues.  
 *  
 *  To avoid overflow, the matrix must be scaled so that its  
 *  largest element is no greater than overflow**(1/2) *  
 *  underflow**(1/4) in absolute value, and for greatest  
 *  accuracy, it should not be much smaller than that.  
 *  
 *  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal  
 *  Matrix", Report CS41, Computer Science Dept., Stanford  
 *  University, July 21, 1966.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  RANGE   (input) CHARACTER*1  
 *          = 'A': ("All")   all eigenvalues will be found.  
 *          = 'V': ("Value") all eigenvalues in the half-open interval  
 *                           (VL, VU] will be found.  
 *          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the  
 *                           entire matrix) will be found.  
 *  
 *  ORDER   (input) CHARACTER*1  
 *          = 'B': ("By Block") the eigenvalues will be grouped by  
 *                              split-off block (see IBLOCK, ISPLIT) and  
 *                              ordered from smallest to largest within  
 *                              the block.  
 *          = 'E': ("Entire matrix")  
 *                              the eigenvalues for the entire matrix  
 *                              will be ordered from smallest to  
 *                              largest.  
 *  
 *  N       (input) INTEGER  
 *          The order of the tridiagonal matrix T.  N >= 0.  
 *  
 *  VL      (input) DOUBLE PRECISION  
 *  VU      (input) DOUBLE PRECISION  
 *          If RANGE='V', the lower and upper bounds of the interval to  
 *          be searched for eigenvalues.  Eigenvalues less than or equal  
 *          to VL, or greater than VU, will not be returned.  VL < VU.  
 *          Not referenced if RANGE = 'A' or 'I'.  
 *  
 *  IL      (input) INTEGER  
 *  IU      (input) INTEGER  
 *          If RANGE='I', the indices (in ascending order) of the  
 *          smallest and largest eigenvalues to be returned.  
 *          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.  
 *          Not referenced if RANGE = 'A' or 'V'.  
 *  
 *  ABSTOL  (input) DOUBLE PRECISION  
 *          The absolute tolerance for the eigenvalues.  An eigenvalue  
 *          (or cluster) is considered to be located if it has been  
 *          determined to lie in an interval whose width is ABSTOL or  
 *          less.  If ABSTOL is less than or equal to zero, then ULP*|T|  
 *          will be used, where |T| means the 1-norm of T.  
 *  
 *          Eigenvalues will be computed most accurately when ABSTOL is  
 *          set to twice the underflow threshold 2*DLAMCH('S'), not zero.  
 *  
 *  D       (input) DOUBLE PRECISION array, dimension (N)  
 *          The n diagonal elements of the tridiagonal matrix T.  
 *  
 *  E       (input) DOUBLE PRECISION array, dimension (N-1)  
 *          The (n-1) off-diagonal elements of the tridiagonal matrix T.  
 *  
 *  M       (output) INTEGER  
 *          The actual number of eigenvalues found. 0 <= M <= N.  
 *          (See also the description of INFO=2,3.)  
 *  
 *  NSPLIT  (output) INTEGER  
 *          The number of diagonal blocks in the matrix T.  
 *          1 <= NSPLIT <= N.  
 *  
 *  W       (output) DOUBLE PRECISION array, dimension (N)  
 *          On exit, the first M elements of W will contain the  
 *          eigenvalues.  (DSTEBZ may use the remaining N-M elements as  
 *          workspace.)  
 *  
 *  IBLOCK  (output) INTEGER array, dimension (N)  
 *          At each row/column j where E(j) is zero or small, the  
 *          matrix T is considered to split into a block diagonal  
 *          matrix.  On exit, if INFO = 0, IBLOCK(i) specifies to which  
 *          block (from 1 to the number of blocks) the eigenvalue W(i)  
 *          belongs.  (DSTEBZ may use the remaining N-M elements as  
 *          workspace.)  
 *  
 *  ISPLIT  (output) INTEGER array, dimension (N)  
 *          The splitting points, at which T breaks up into submatrices.  
 *          The first submatrix consists of rows/columns 1 to ISPLIT(1),  
 *          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),  
 *          etc., and the NSPLIT-th consists of rows/columns  
 *          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.  
 *          (Only the first NSPLIT elements will actually be used, but  
 *          since the user cannot know a priori what value NSPLIT will  
 *          have, N words must be reserved for ISPLIT.)  
 *  
 *  WORK    (workspace) DOUBLE PRECISION array, dimension (4*N)  
 *  
 *  IWORK   (workspace) INTEGER array, dimension (3*N)  
 *  
 *  INFO    (output) INTEGER  
 *          = 0:  successful exit  
 *          < 0:  if INFO = -i, the i-th argument had an illegal value  
 *          > 0:  some or all of the eigenvalues failed to converge or  
 *                were not computed:  
 *                =1 or 3: Bisection failed to converge for some  
 *                        eigenvalues; these eigenvalues are flagged by a  
 *                        negative block number.  The effect is that the  
 *                        eigenvalues may not be as accurate as the  
 *                        absolute and relative tolerances.  This is  
 *                        generally caused by unexpectedly inaccurate  
 *                        arithmetic.  
 *                =2 or 3: RANGE='I' only: Not all of the eigenvalues  
 *                        IL:IU were found.  
 *                        Effect: M < IU+1-IL  
 *                        Cause:  non-monotonic arithmetic, causing the  
 *                                Sturm sequence to be non-monotonic.  
 *                        Cure:   recalculate, using RANGE='A', and pick  
 *                                out eigenvalues IL:IU.  In some cases,  
 *                                increasing the PARAMETER "FUDGE" may  
 *                                make things work.  
 *                = 4:    RANGE='I', and the Gershgorin interval  
 *                        initially used was too small.  No eigenvalues  
 *                        were computed.  
 *                        Probable cause: your machine has sloppy  
 *                                        floating-point arithmetic.  
 *                        Cure: Increase the PARAMETER "FUDGE",  
 *                              recompile, and try again.  
 *  
 *  Internal Parameters  
 *  ===================  
 *  
 *  RELFAC  DOUBLE PRECISION, default = 2.0e0  
 *          The relative tolerance.  An interval (a,b] lies within  
 *          "relative tolerance" if  b-a < RELFAC*ulp*max(|a|,|b|),  
 *          where "ulp" is the machine precision (distance from 1 to  
 *          the next larger floating point number.)  
 *  
 *  FUDGE   DOUBLE PRECISION, default = 2  
 *          A "fudge factor" to widen the Gershgorin intervals.  Ideally,  
 *          a value of 1 should work, but on machines with sloppy  
 *          arithmetic, this needs to be larger.  The default for  
 *          publicly released versions should be large enough to handle  
 *          the worst machine around.  Note that this has no effect  
 *          on accuracy of the solution.  
 *  
 *  =====================================================================  *  =====================================================================
 *  *
 *     .. Parameters ..  *     .. Parameters ..
Removed from v.1.8  
changed lines
  Added in v.1.9

CVSweb interface <joel.bertrand@systella.fr>