205 lines
		
	
	
		
			5.8 KiB
		
	
	
	
		
			Fortran
		
	
	
	
			
		
		
	
	
			205 lines
		
	
	
		
			5.8 KiB
		
	
	
	
		
			Fortran
		
	
	
	
| *> \brief \b SLARRR performs tests to decide whether the symmetric tridiagonal matrix T warrants expensive computations which guarantee high relative accuracy in the eigenvalues.
 | |
| *
 | |
| *  =========== DOCUMENTATION ===========
 | |
| *
 | |
| * Online html documentation available at 
 | |
| *            http://www.netlib.org/lapack/explore-html/ 
 | |
| *
 | |
| *> \htmlonly
 | |
| *> Download SLARRR + dependencies 
 | |
| *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarrr.f"> 
 | |
| *> [TGZ]</a> 
 | |
| *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarrr.f"> 
 | |
| *> [ZIP]</a> 
 | |
| *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarrr.f"> 
 | |
| *> [TXT]</a>
 | |
| *> \endhtmlonly 
 | |
| *
 | |
| *  Definition:
 | |
| *  ===========
 | |
| *
 | |
| *       SUBROUTINE SLARRR( N, D, E, INFO )
 | |
| * 
 | |
| *       .. Scalar Arguments ..
 | |
| *       INTEGER            N, INFO
 | |
| *       ..
 | |
| *       .. Array Arguments ..
 | |
| *       REAL               D( * ), E( * )
 | |
| *       ..
 | |
| *  
 | |
| *  
 | |
| *
 | |
| *> \par Purpose:
 | |
| *  =============
 | |
| *>
 | |
| *> \verbatim
 | |
| *>
 | |
| *> Perform tests to decide whether the symmetric tridiagonal matrix T
 | |
| *> warrants expensive computations which guarantee high relative accuracy
 | |
| *> in the eigenvalues.
 | |
| *> \endverbatim
 | |
| *
 | |
| *  Arguments:
 | |
| *  ==========
 | |
| *
 | |
| *> \param[in] N
 | |
| *> \verbatim
 | |
| *>          N is INTEGER
 | |
| *>          The order of the matrix. N > 0.
 | |
| *> \endverbatim
 | |
| *>
 | |
| *> \param[in] D
 | |
| *> \verbatim
 | |
| *>          D is REAL array, dimension (N)
 | |
| *>          The N diagonal elements of the tridiagonal matrix T.
 | |
| *> \endverbatim
 | |
| *>
 | |
| *> \param[in,out] E
 | |
| *> \verbatim
 | |
| *>          E is REAL array, dimension (N)
 | |
| *>          On entry, the first (N-1) entries contain the subdiagonal
 | |
| *>          elements of the tridiagonal matrix T; E(N) is set to ZERO.
 | |
| *> \endverbatim
 | |
| *>
 | |
| *> \param[out] INFO
 | |
| *> \verbatim
 | |
| *>          INFO is INTEGER
 | |
| *>          INFO = 0(default) : the matrix warrants computations preserving
 | |
| *>                              relative accuracy.
 | |
| *>          INFO = 1          : the matrix warrants computations guaranteeing
 | |
| *>                              only absolute accuracy.
 | |
| *> \endverbatim
 | |
| *
 | |
| *  Authors:
 | |
| *  ========
 | |
| *
 | |
| *> \author Univ. of Tennessee 
 | |
| *> \author Univ. of California Berkeley 
 | |
| *> \author Univ. of Colorado Denver 
 | |
| *> \author NAG Ltd. 
 | |
| *
 | |
| *> \date September 2012
 | |
| *
 | |
| *> \ingroup auxOTHERauxiliary
 | |
| *
 | |
| *> \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
 | |
| *
 | |
| *  =====================================================================
 | |
|       SUBROUTINE SLARRR( N, D, E, INFO )
 | |
| *
 | |
| *  -- LAPACK auxiliary routine (version 3.4.2) --
 | |
| *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 | |
| *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 | |
| *     September 2012
 | |
| *
 | |
| *     .. Scalar Arguments ..
 | |
|       INTEGER            N, INFO
 | |
| *     ..
 | |
| *     .. Array Arguments ..
 | |
|       REAL               D( * ), E( * )
 | |
| *     ..
 | |
| *
 | |
| *
 | |
| *  =====================================================================
 | |
| *
 | |
| *     .. Parameters ..
 | |
|       REAL               ZERO, RELCOND
 | |
|       PARAMETER          ( ZERO = 0.0E0,
 | |
|      $                     RELCOND = 0.999E0 )
 | |
| *     ..
 | |
| *     .. Local Scalars ..
 | |
|       INTEGER            I
 | |
|       LOGICAL            YESREL
 | |
|       REAL               EPS, SAFMIN, SMLNUM, RMIN, TMP, TMP2,
 | |
|      $          OFFDIG, OFFDIG2
 | |
| 
 | |
| *     ..
 | |
| *     .. External Functions ..
 | |
|       REAL               SLAMCH
 | |
|       EXTERNAL           SLAMCH
 | |
| *     ..
 | |
| *     .. Intrinsic Functions ..
 | |
|       INTRINSIC          ABS
 | |
| *     ..
 | |
| *     .. Executable Statements ..
 | |
| *
 | |
| *     As a default, do NOT go for relative-accuracy preserving computations.
 | |
|       INFO = 1
 | |
| 
 | |
|       SAFMIN = SLAMCH( 'Safe minimum' )
 | |
|       EPS = SLAMCH( 'Precision' )
 | |
|       SMLNUM = SAFMIN / EPS
 | |
|       RMIN = SQRT( SMLNUM )
 | |
| 
 | |
| *     Tests for relative accuracy
 | |
| *
 | |
| *     Test for scaled diagonal dominance
 | |
| *     Scale the diagonal entries to one and check whether the sum of the
 | |
| *     off-diagonals is less than one
 | |
| *
 | |
| *     The sdd relative error bounds have a 1/(1- 2*x) factor in them,
 | |
| *     x = max(OFFDIG + OFFDIG2), so when x is close to 1/2, no relative
 | |
| *     accuracy is promised.  In the notation of the code fragment below,
 | |
| *     1/(1 - (OFFDIG + OFFDIG2)) is the condition number.
 | |
| *     We don't think it is worth going into "sdd mode" unless the relative
 | |
| *     condition number is reasonable, not 1/macheps.
 | |
| *     The threshold should be compatible with other thresholds used in the
 | |
| *     code. We set  OFFDIG + OFFDIG2 <= .999 =: RELCOND, it corresponds
 | |
| *     to losing at most 3 decimal digits: 1 / (1 - (OFFDIG + OFFDIG2)) <= 1000
 | |
| *     instead of the current OFFDIG + OFFDIG2 < 1
 | |
| *
 | |
|       YESREL = .TRUE.
 | |
|       OFFDIG = ZERO
 | |
|       TMP = SQRT(ABS(D(1)))
 | |
|       IF (TMP.LT.RMIN) YESREL = .FALSE.
 | |
|       IF(.NOT.YESREL) GOTO 11
 | |
|       DO 10 I = 2, N
 | |
|          TMP2 = SQRT(ABS(D(I)))
 | |
|          IF (TMP2.LT.RMIN) YESREL = .FALSE.
 | |
|          IF(.NOT.YESREL) GOTO 11
 | |
|          OFFDIG2 = ABS(E(I-1))/(TMP*TMP2)
 | |
|          IF(OFFDIG+OFFDIG2.GE.RELCOND) YESREL = .FALSE.
 | |
|          IF(.NOT.YESREL) GOTO 11
 | |
|          TMP = TMP2
 | |
|          OFFDIG = OFFDIG2
 | |
|  10   CONTINUE
 | |
|  11   CONTINUE
 | |
| 
 | |
|       IF( YESREL ) THEN
 | |
|          INFO = 0
 | |
|          RETURN
 | |
|       ELSE
 | |
|       ENDIF
 | |
| *
 | |
| 
 | |
| *
 | |
| *     *** MORE TO BE IMPLEMENTED ***
 | |
| *
 | |
| 
 | |
| *
 | |
| *     Test if the lower bidiagonal matrix L from T = L D L^T
 | |
| *     (zero shift facto) is well conditioned
 | |
| *
 | |
| 
 | |
| *
 | |
| *     Test if the upper bidiagonal matrix U from T = U D U^T
 | |
| *     (zero shift facto) is well conditioned.
 | |
| *     In this case, the matrix needs to be flipped and, at the end
 | |
| *     of the eigenvector computation, the flip needs to be applied
 | |
| *     to the computed eigenvectors (and the support)
 | |
| *
 | |
| 
 | |
| *
 | |
|       RETURN
 | |
| *
 | |
| *     END OF SLARRR
 | |
| *
 | |
|       END
 |