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

Diff for /rpl/lapack/lapack/dlarre.f between versions 1.7 and 1.24

version 1.7, 2010/08/13 21:03:51 version 1.24, 2023/08/07 08:38:57
Line 1 Line 1
   *> \brief \b DLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduced block Ti, finds base representations and eigenvalues.
   *
   *  =========== DOCUMENTATION ===========
   *
   * Online html documentation available at
   *            http://www.netlib.org/lapack/explore-html/
   *
   *> \htmlonly
   *> Download DLARRE + dependencies
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarre.f">
   *> [TGZ]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarre.f">
   *> [ZIP]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarre.f">
   *> [TXT]</a>
   *> \endhtmlonly
   *
   *  Definition:
   *  ===========
   *
   *       SUBROUTINE DLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2,
   *                           RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,
   *                           W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
   *                           WORK, IWORK, INFO )
   *
   *       .. Scalar Arguments ..
   *       CHARACTER          RANGE
   *       INTEGER            IL, INFO, IU, M, N, NSPLIT
   *       DOUBLE PRECISION  PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
   *       ..
   *       .. Array Arguments ..
   *       INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * ),
   *      $                   INDEXW( * )
   *       DOUBLE PRECISION   D( * ), E( * ), E2( * ), GERS( * ),
   *      $                   W( * ),WERR( * ), WGAP( * ), WORK( * )
   *       ..
   *
   *
   *> \par Purpose:
   *  =============
   *>
   *> \verbatim
   *>
   *> To find the desired eigenvalues of a given real symmetric
   *> tridiagonal matrix T, DLARRE sets any "small" off-diagonal
   *> elements to zero, and for each unreduced block T_i, it finds
   *> (a) a suitable shift at one end of the block's spectrum,
   *> (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
   *> (c) eigenvalues of each L_i D_i L_i^T.
   *> The representations and eigenvalues found are then used by
   *> DSTEMR to compute the eigenvectors of T.
   *> The accuracy varies depending on whether bisection is used to
   *> find a few eigenvalues or the dqds algorithm (subroutine DLASQ2) to
   *> conpute all and then discard any unwanted one.
   *> As an added benefit, DLARRE also outputs the n
   *> Gerschgorin intervals for the matrices L_i D_i L_i^T.
   *> \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] N
   *> \verbatim
   *>          N is INTEGER
   *>          The order of the matrix. N > 0.
   *> \endverbatim
   *>
   *> \param[in,out] VL
   *> \verbatim
   *>          VL is DOUBLE PRECISION
   *>          If RANGE='V', the lower bound for the eigenvalues.
   *>          Eigenvalues less than or equal to VL, or greater than VU,
   *>          will not be returned.  VL < VU.
   *>          If RANGE='I' or ='A', DLARRE computes bounds on the desired
   *>          part of the spectrum.
   *> \endverbatim
   *>
   *> \param[in,out] VU
   *> \verbatim
   *>          VU is DOUBLE PRECISION
   *>          If RANGE='V', the upper bound for the eigenvalues.
   *>          Eigenvalues less than or equal to VL, or greater than VU,
   *>          will not be returned.  VL < VU.
   *>          If RANGE='I' or ='A', DLARRE computes bounds on the desired
   *>          part of the spectrum.
   *> \endverbatim
   *>
   *> \param[in] IL
   *> \verbatim
   *>          IL is INTEGER
   *>          If RANGE='I', the index of the
   *>          smallest eigenvalue to be returned.
   *>          1 <= IL <= IU <= N.
   *> \endverbatim
   *>
   *> \param[in] IU
   *> \verbatim
   *>          IU is INTEGER
   *>          If RANGE='I', the index of the
   *>          largest eigenvalue to be returned.
   *>          1 <= IL <= IU <= N.
   *> \endverbatim
   *>
   *> \param[in,out] D
   *> \verbatim
   *>          D is DOUBLE PRECISION array, dimension (N)
   *>          On entry, the N diagonal elements of the tridiagonal
   *>          matrix T.
   *>          On exit, the N diagonal elements of the diagonal
   *>          matrices D_i.
   *> \endverbatim
   *>
   *> \param[in,out] E
   *> \verbatim
   *>          E is DOUBLE PRECISION array, dimension (N)
   *>          On entry, the first (N-1) entries contain the subdiagonal
   *>          elements of the tridiagonal matrix T; E(N) need not be set.
   *>          On exit, E contains the subdiagonal elements of the unit
   *>          bidiagonal matrices L_i. The entries E( ISPLIT( I ) ),
   *>          1 <= I <= NSPLIT, contain the base points sigma_i on output.
   *> \endverbatim
   *>
   *> \param[in,out] E2
   *> \verbatim
   *>          E2 is DOUBLE PRECISION array, dimension (N)
   *>          On entry, the first (N-1) entries contain the SQUARES of the
   *>          subdiagonal elements of the tridiagonal matrix T;
   *>          E2(N) need not be set.
   *>          On exit, the entries E2( ISPLIT( I ) ),
   *>          1 <= I <= NSPLIT, have been set to zero
   *> \endverbatim
   *>
   *> \param[in] RTOL1
   *> \verbatim
   *>          RTOL1 is DOUBLE PRECISION
   *> \endverbatim
   *>
   *> \param[in] RTOL2
   *> \verbatim
   *>          RTOL2 is DOUBLE PRECISION
   *>           Parameters for bisection.
   *>           An interval [LEFT,RIGHT] has converged if
   *>           RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
   *> \endverbatim
   *>
   *> \param[in] SPLTOL
   *> \verbatim
   *>          SPLTOL is DOUBLE PRECISION
   *>          The threshold for splitting.
   *> \endverbatim
   *>
   *> \param[out] NSPLIT
   *> \verbatim
   *>          NSPLIT is INTEGER
   *>          The number of blocks T splits into. 1 <= NSPLIT <= N.
   *> \endverbatim
   *>
   *> \param[out] ISPLIT
   *> \verbatim
   *>          ISPLIT is INTEGER array, dimension (N)
   *>          The splitting points, at which T breaks up into blocks.
   *>          The first block 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.
   *> \endverbatim
   *>
   *> \param[out] M
   *> \verbatim
   *>          M is INTEGER
   *>          The total number of eigenvalues (of all L_i D_i L_i^T)
   *>          found.
   *> \endverbatim
   *>
   *> \param[out] W
   *> \verbatim
   *>          W is DOUBLE PRECISION array, dimension (N)
   *>          The first M elements contain the eigenvalues. The
   *>          eigenvalues of each of the blocks, L_i D_i L_i^T, are
   *>          sorted in ascending order ( DLARRE may use the
   *>          remaining N-M elements as workspace).
   *> \endverbatim
   *>
   *> \param[out] WERR
   *> \verbatim
   *>          WERR is DOUBLE PRECISION array, dimension (N)
   *>          The error bound on the corresponding eigenvalue in W.
   *> \endverbatim
   *>
   *> \param[out] WGAP
   *> \verbatim
   *>          WGAP is DOUBLE PRECISION array, dimension (N)
   *>          The separation from the right neighbor eigenvalue in W.
   *>          The gap is only with respect to the eigenvalues of the same block
   *>          as each block has its own representation tree.
   *>          Exception: at the right end of a block we store the left gap
   *> \endverbatim
   *>
   *> \param[out] IBLOCK
   *> \verbatim
   *>          IBLOCK is INTEGER array, dimension (N)
   *>          The indices of the blocks (submatrices) associated with the
   *>          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
   *>          W(i) belongs to the first block from the top, =2 if W(i)
   *>          belongs to the second block, etc.
   *> \endverbatim
   *>
   *> \param[out] INDEXW
   *> \verbatim
   *>          INDEXW is INTEGER array, dimension (N)
   *>          The indices of the eigenvalues within each block (submatrix);
   *>          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
   *>          i-th eigenvalue W(i) is the 10-th eigenvalue in block 2
   *> \endverbatim
   *>
   *> \param[out] GERS
   *> \verbatim
   *>          GERS is DOUBLE PRECISION array, dimension (2*N)
   *>          The N Gerschgorin intervals (the i-th Gerschgorin interval
   *>          is (GERS(2*i-1), GERS(2*i)).
   *> \endverbatim
   *>
   *> \param[out] PIVMIN
   *> \verbatim
   *>          PIVMIN is DOUBLE PRECISION
   *>          The minimum pivot in the Sturm sequence for T.
   *> \endverbatim
   *>
   *> \param[out] WORK
   *> \verbatim
   *>          WORK is DOUBLE PRECISION array, dimension (6*N)
   *>          Workspace.
   *> \endverbatim
   *>
   *> \param[out] IWORK
   *> \verbatim
   *>          IWORK is INTEGER array, dimension (5*N)
   *>          Workspace.
   *> \endverbatim
   *>
   *> \param[out] INFO
   *> \verbatim
   *>          INFO is INTEGER
   *>          = 0:  successful exit
   *>          > 0:  A problem occurred in DLARRE.
   *>          < 0:  One of the called subroutines signaled an internal problem.
   *>                Needs inspection of the corresponding parameter IINFO
   *>                for further information.
   *>
   *>          =-1:  Problem in DLARRD.
   *>          = 2:  No base representation could be found in MAXTRY iterations.
   *>                Increasing MAXTRY and recompilation might be a remedy.
   *>          =-3:  Problem in DLARRB when computing the refined root
   *>                representation for DLASQ2.
   *>          =-4:  Problem in DLARRB when preforming bisection on the
   *>                desired part of the spectrum.
   *>          =-5:  Problem in DLASQ2.
   *>          =-6:  Problem in DLASQ2.
   *> \endverbatim
   *
   *  Authors:
   *  ========
   *
   *> \author Univ. of Tennessee
   *> \author Univ. of California Berkeley
   *> \author Univ. of Colorado Denver
   *> \author NAG Ltd.
   *
   *> \ingroup OTHERauxiliary
   *
   *> \par Further Details:
   *  =====================
   *>
   *> \verbatim
   *>
   *>  The base representations are required to suffer very little
   *>  element growth and consequently define all their eigenvalues to
   *>  high relative accuracy.
   *> \endverbatim
   *
   *> \par Contributors:
   *  ==================
   *>
   *>     Beresford Parlett, University of California, Berkeley, USA \n
   *>     Jim Demmel, University of California, Berkeley, USA \n
   *>     Inderjit Dhillon, University of Texas, Austin, USA \n
   *>     Osni Marques, LBNL/NERSC, USA \n
   *>     Christof Voemel, University of California, Berkeley, USA \n
   *>
   *  =====================================================================
       SUBROUTINE DLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2,        SUBROUTINE DLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2,
      $                    RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,       $                    RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,
      $                    W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,       $                    W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
      $                    WORK, IWORK, INFO )       $                    WORK, IWORK, INFO )
       IMPLICIT NONE  
 *  *
 *  -- LAPACK auxiliary routine (version 3.2.2) --  *  -- LAPACK auxiliary routine --
 *  -- 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..--
 *     June 2010  
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          RANGE        CHARACTER          RANGE
Line 21 Line 319
      $                   W( * ),WERR( * ), WGAP( * ), WORK( * )       $                   W( * ),WERR( * ), WGAP( * ), WORK( * )
 *     ..  *     ..
 *  *
 *  Purpose  
 *  =======  
 *  
 *  To find the desired eigenvalues of a given real symmetric  
 *  tridiagonal matrix T, DLARRE sets any "small" off-diagonal  
 *  elements to zero, and for each unreduced block T_i, it finds  
 *  (a) a suitable shift at one end of the block's spectrum,  
 *  (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and  
 *  (c) eigenvalues of each L_i D_i L_i^T.  
 *  The representations and eigenvalues found are then used by  
 *  DSTEMR to compute the eigenvectors of T.  
 *  The accuracy varies depending on whether bisection is used to  
 *  find a few eigenvalues or the dqds algorithm (subroutine DLASQ2) to  
 *  conpute all and then discard any unwanted one.  
 *  As an added benefit, DLARRE also outputs the n  
 *  Gerschgorin intervals for the matrices L_i D_i L_i^T.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  RANGE   (input) CHARACTER  
 *          = '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.  
 *  
 *  N       (input) INTEGER  
 *          The order of the matrix. N > 0.  
 *  
 *  VL      (input/output) DOUBLE PRECISION  
 *  VU      (input/output) DOUBLE PRECISION  
 *          If RANGE='V', the lower and upper bounds for the eigenvalues.  
 *          Eigenvalues less than or equal to VL, or greater than VU,  
 *          will not be returned.  VL < VU.  
 *          If RANGE='I' or ='A', DLARRE computes bounds on the desired  
 *          part of the spectrum.  
 *  
 *  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.  
 *  
 *  D       (input/output) DOUBLE PRECISION array, dimension (N)  
 *          On entry, the N diagonal elements of the tridiagonal  
 *          matrix T.  
 *          On exit, the N diagonal elements of the diagonal  
 *          matrices D_i.  
 *  
 *  E       (input/output) DOUBLE PRECISION array, dimension (N)  
 *          On entry, the first (N-1) entries contain the subdiagonal  
 *          elements of the tridiagonal matrix T; E(N) need not be set.  
 *          On exit, E contains the subdiagonal elements of the unit  
 *          bidiagonal matrices L_i. The entries E( ISPLIT( I ) ),  
 *          1 <= I <= NSPLIT, contain the base points sigma_i on output.  
 *  
 *  E2      (input/output) DOUBLE PRECISION array, dimension (N)  
 *          On entry, the first (N-1) entries contain the SQUARES of the  
 *          subdiagonal elements of the tridiagonal matrix T;  
 *          E2(N) need not be set.  
 *          On exit, the entries E2( ISPLIT( I ) ),  
 *          1 <= I <= NSPLIT, have been set to zero  
 *  
 *  RTOL1   (input) DOUBLE PRECISION  
 *  RTOL2   (input) DOUBLE PRECISION  
 *           Parameters for bisection.  
 *           An interval [LEFT,RIGHT] has converged if  
 *           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )  
 *  
 *  SPLTOL  (input) DOUBLE PRECISION  
 *          The threshold for splitting.  
 *  
 *  NSPLIT  (output) INTEGER  
 *          The number of blocks T splits into. 1 <= NSPLIT <= N.  
 *  
 *  ISPLIT  (output) INTEGER array, dimension (N)  
 *          The splitting points, at which T breaks up into blocks.  
 *          The first block 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.  
 *  
 *  M       (output) INTEGER  
 *          The total number of eigenvalues (of all L_i D_i L_i^T)  
 *          found.  
 *  
 *  W       (output) DOUBLE PRECISION array, dimension (N)  
 *          The first M elements contain the eigenvalues. The  
 *          eigenvalues of each of the blocks, L_i D_i L_i^T, are  
 *          sorted in ascending order ( DLARRE may use the  
 *          remaining N-M elements as workspace).  
 *  
 *  WERR    (output) DOUBLE PRECISION array, dimension (N)  
 *          The error bound on the corresponding eigenvalue in W.  
 *  
 *  WGAP    (output) DOUBLE PRECISION array, dimension (N)  
 *          The separation from the right neighbor eigenvalue in W.  
 *          The gap is only with respect to the eigenvalues of the same block  
 *          as each block has its own representation tree.  
 *          Exception: at the right end of a block we store the left gap  
 *  
 *  IBLOCK  (output) INTEGER array, dimension (N)  
 *          The indices of the blocks (submatrices) associated with the  
 *          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue  
 *          W(i) belongs to the first block from the top, =2 if W(i)  
 *          belongs to the second block, etc.  
 *  
 *  INDEXW  (output) INTEGER array, dimension (N)  
 *          The indices of the eigenvalues within each block (submatrix);  
 *          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the  
 *          i-th eigenvalue W(i) is the 10-th eigenvalue in block 2  
 *  
 *  GERS    (output) DOUBLE PRECISION array, dimension (2*N)  
 *          The N Gerschgorin intervals (the i-th Gerschgorin interval  
 *          is (GERS(2*i-1), GERS(2*i)).  
 *  
 *  PIVMIN  (output) DOUBLE PRECISION  
 *          The minimum pivot in the Sturm sequence for T.  
 *  
 *  WORK    (workspace) DOUBLE PRECISION array, dimension (6*N)  
 *          Workspace.  
 *  
 *  IWORK   (workspace) INTEGER array, dimension (5*N)  
 *          Workspace.  
 *  
 *  INFO    (output) INTEGER  
 *          = 0:  successful exit  
 *          > 0:  A problem occured in DLARRE.  
 *          < 0:  One of the called subroutines signaled an internal problem.  
 *                Needs inspection of the corresponding parameter IINFO  
 *                for further information.  
 *  
 *          =-1:  Problem in DLARRD.  
 *          = 2:  No base representation could be found in MAXTRY iterations.  
 *                Increasing MAXTRY and recompilation might be a remedy.  
 *          =-3:  Problem in DLARRB when computing the refined root  
 *                representation for DLASQ2.  
 *          =-4:  Problem in DLARRB when preforming bisection on the  
 *                desired part of the spectrum.  
 *          =-5:  Problem in DLASQ2.  
 *          =-6:  Problem in DLASQ2.  
 *  
 *  Further Details  
 *  The base representations are required to suffer very little  
 *  element growth and consequently define all their eigenvalues to  
 *  high relative accuracy.  
 *  ===============  
 *  
 *  Based on contributions by  
 *     Beresford Parlett, University of California, Berkeley, USA  
 *     Jim Demmel, University of California, Berkeley, USA  
 *     Inderjit Dhillon, University of Texas, Austin, USA  
 *     Osni Marques, LBNL/NERSC, USA  
 *     Christof Voemel, University of California, Berkeley, USA  
 *  
 *  =====================================================================  *  =====================================================================
 *  *
 *     .. Parameters ..  *     .. Parameters ..
Line 215 Line 357
 *     ..  *     ..
 *     .. External Subroutines ..  *     .. External Subroutines ..
       EXTERNAL           DCOPY, DLARNV, DLARRA, DLARRB, DLARRC, DLARRD,        EXTERNAL           DCOPY, DLARNV, DLARRA, DLARRB, DLARRC, DLARRD,
      $                   DLASQ2       $                   DLASQ2, DLARRK
 *     ..  *     ..
 *     .. Intrinsic Functions ..  *     .. Intrinsic Functions ..
       INTRINSIC          ABS, MAX, MIN        INTRINSIC          ABS, MAX, MIN
Line 225 Line 367
 *  *
   
       INFO = 0        INFO = 0
   *
   *     Quick return if possible
   *
         IF( N.LE.0 ) THEN
            RETURN
         END IF
 *  *
 *     Decode RANGE  *     Decode RANGE
 *  *
Line 532 Line 679
 *           The initial SIGMA was to the outer end of the spectrum  *           The initial SIGMA was to the outer end of the spectrum
 *           the matrix is definite and we need not retreat.  *           the matrix is definite and we need not retreat.
             TAU = SPDIAM*EPS*N + TWO*PIVMIN              TAU = SPDIAM*EPS*N + TWO*PIVMIN
               TAU = MAX( TAU,TWO*EPS*ABS(SIGMA) )
          ELSE           ELSE
             IF(MB.GT.1) THEN              IF(MB.GT.1) THEN
                CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN)                 CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN)
Line 748 Line 896
   
       RETURN        RETURN
 *  *
 *     end of DLARRE  *     End of DLARRE
 *  *
       END        END
Removed from v.1.7  
changed lines
  Added in v.1.24

CVSweb interface <joel.bertrand@systella.fr>