LAPACK 3.3.1
Linear Algebra PACKage

dlarrr.f

Go to the documentation of this file.
00001       SUBROUTINE DLARRR( N, D, E, INFO )
00002 *
00003 *  -- LAPACK auxiliary routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            N, INFO
00010 *     ..
00011 *     .. Array Arguments ..
00012       DOUBLE PRECISION   D( * ), E( * )
00013 *     ..
00014 *
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  Perform tests to decide whether the symmetric tridiagonal matrix T
00020 *  warrants expensive computations which guarantee high relative accuracy
00021 *  in the eigenvalues.
00022 *
00023 *  Arguments
00024 *  =========
00025 *
00026 *  N       (input) INTEGER
00027 *          The order of the matrix. N > 0.
00028 *
00029 *  D       (input) DOUBLE PRECISION array, dimension (N)
00030 *          The N diagonal elements of the tridiagonal matrix T.
00031 *
00032 *  E       (input/output) DOUBLE PRECISION array, dimension (N)
00033 *          On entry, the first (N-1) entries contain the subdiagonal
00034 *          elements of the tridiagonal matrix T; E(N) is set to ZERO.
00035 *
00036 *  INFO    (output) INTEGER
00037 *          INFO = 0(default) : the matrix warrants computations preserving
00038 *                              relative accuracy.
00039 *          INFO = 1          : the matrix warrants computations guaranteeing
00040 *                              only absolute accuracy.
00041 *
00042 *  Further Details
00043 *  ===============
00044 *
00045 *  Based on contributions by
00046 *     Beresford Parlett, University of California, Berkeley, USA
00047 *     Jim Demmel, University of California, Berkeley, USA
00048 *     Inderjit Dhillon, University of Texas, Austin, USA
00049 *     Osni Marques, LBNL/NERSC, USA
00050 *     Christof Voemel, University of California, Berkeley, USA
00051 *
00052 *  =====================================================================
00053 *
00054 *     .. Parameters ..
00055       DOUBLE PRECISION   ZERO, RELCOND
00056       PARAMETER          ( ZERO = 0.0D0,
00057      $                     RELCOND = 0.999D0 )
00058 *     ..
00059 *     .. Local Scalars ..
00060       INTEGER            I
00061       LOGICAL            YESREL
00062       DOUBLE PRECISION   EPS, SAFMIN, SMLNUM, RMIN, TMP, TMP2,
00063      $          OFFDIG, OFFDIG2
00064 
00065 *     ..
00066 *     .. External Functions ..
00067       DOUBLE PRECISION   DLAMCH
00068       EXTERNAL           DLAMCH
00069 *     ..
00070 *     .. Intrinsic Functions ..
00071       INTRINSIC          ABS
00072 *     ..
00073 *     .. Executable Statements ..
00074 *
00075 *     As a default, do NOT go for relative-accuracy preserving computations.
00076       INFO = 1
00077 
00078       SAFMIN = DLAMCH( 'Safe minimum' )
00079       EPS = DLAMCH( 'Precision' )
00080       SMLNUM = SAFMIN / EPS
00081       RMIN = SQRT( SMLNUM )
00082 
00083 *     Tests for relative accuracy
00084 *
00085 *     Test for scaled diagonal dominance
00086 *     Scale the diagonal entries to one and check whether the sum of the
00087 *     off-diagonals is less than one
00088 *
00089 *     The sdd relative error bounds have a 1/(1- 2*x) factor in them,
00090 *     x = max(OFFDIG + OFFDIG2), so when x is close to 1/2, no relative
00091 *     accuracy is promised.  In the notation of the code fragment below,
00092 *     1/(1 - (OFFDIG + OFFDIG2)) is the condition number.
00093 *     We don't think it is worth going into "sdd mode" unless the relative
00094 *     condition number is reasonable, not 1/macheps.
00095 *     The threshold should be compatible with other thresholds used in the
00096 *     code. We set  OFFDIG + OFFDIG2 <= .999 =: RELCOND, it corresponds
00097 *     to losing at most 3 decimal digits: 1 / (1 - (OFFDIG + OFFDIG2)) <= 1000
00098 *     instead of the current OFFDIG + OFFDIG2 < 1
00099 *
00100       YESREL = .TRUE.
00101       OFFDIG = ZERO
00102       TMP = SQRT(ABS(D(1)))
00103       IF (TMP.LT.RMIN) YESREL = .FALSE.
00104       IF(.NOT.YESREL) GOTO 11
00105       DO 10 I = 2, N
00106          TMP2 = SQRT(ABS(D(I)))
00107          IF (TMP2.LT.RMIN) YESREL = .FALSE.
00108          IF(.NOT.YESREL) GOTO 11
00109          OFFDIG2 = ABS(E(I-1))/(TMP*TMP2)
00110          IF(OFFDIG+OFFDIG2.GE.RELCOND) YESREL = .FALSE.
00111          IF(.NOT.YESREL) GOTO 11
00112          TMP = TMP2
00113          OFFDIG = OFFDIG2
00114  10   CONTINUE
00115  11   CONTINUE
00116 
00117       IF( YESREL ) THEN
00118          INFO = 0
00119          RETURN
00120       ELSE
00121       ENDIF
00122 *
00123 
00124 *
00125 *     *** MORE TO BE IMPLEMENTED ***
00126 *
00127 
00128 *
00129 *     Test if the lower bidiagonal matrix L from T = L D L^T
00130 *     (zero shift facto) is well conditioned
00131 *
00132 
00133 *
00134 *     Test if the upper bidiagonal matrix U from T = U D U^T
00135 *     (zero shift facto) is well conditioned.
00136 *     In this case, the matrix needs to be flipped and, at the end
00137 *     of the eigenvector computation, the flip needs to be applied
00138 *     to the computed eigenvectors (and the support)
00139 *
00140 
00141 *
00142       RETURN
00143 *
00144 *     END OF DLARRR
00145 *
00146       END
 All Files Functions