LAPACK: SRC/ddisna.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DDISNA( JOB, M, N, D, SEP, INFO )
00002 *
00003 *  -- LAPACK 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       CHARACTER          JOB
00010       INTEGER            INFO, M, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   D( * ), SEP( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DDISNA computes the reciprocal condition numbers for the eigenvectors
00020 *  of a real symmetric or complex Hermitian matrix or for the left or
00021 *  right singular vectors of a general m-by-n matrix. The reciprocal
00022 *  condition number is the 'gap' between the corresponding eigenvalue or
00023 *  singular value and the nearest other one.
00024 *
00025 *  The bound on the error, measured by angle in radians, in the I-th
00026 *  computed vector is given by
00027 *
00028 *         DLAMCH( 'E' ) * ( ANORM / SEP( I ) )
00029 *
00030 *  where ANORM = 2-norm(A) = max( abs( D(j) ) ).  SEP(I) is not allowed
00031 *  to be smaller than DLAMCH( 'E' )*ANORM in order to limit the size of
00032 *  the error bound.
00033 *
00034 *  DDISNA may also be used to compute error bounds for eigenvectors of
00035 *  the generalized symmetric definite eigenproblem.
00036 *
00037 *  Arguments
00038 *  =========
00039 *
00040 *  JOB     (input) CHARACTER*1
00041 *          Specifies for which problem the reciprocal condition numbers
00042 *          should be computed:
00043 *          = 'E':  the eigenvectors of a symmetric/Hermitian matrix;
00044 *          = 'L':  the left singular vectors of a general matrix;
00045 *          = 'R':  the right singular vectors of a general matrix.
00046 *
00047 *  M       (input) INTEGER
00048 *          The number of rows of the matrix. M >= 0.
00049 *
00050 *  N       (input) INTEGER
00051 *          If JOB = 'L' or 'R', the number of columns of the matrix,
00052 *          in which case N >= 0. Ignored if JOB = 'E'.
00053 *
00054 *  D       (input) DOUBLE PRECISION array, dimension (M) if JOB = 'E'
00055 *                              dimension (min(M,N)) if JOB = 'L' or 'R'
00056 *          The eigenvalues (if JOB = 'E') or singular values (if JOB =
00057 *          'L' or 'R') of the matrix, in either increasing or decreasing
00058 *          order. If singular values, they must be non-negative.
00059 *
00060 *  SEP     (output) DOUBLE PRECISION array, dimension (M) if JOB = 'E'
00061 *                               dimension (min(M,N)) if JOB = 'L' or 'R'
00062 *          The reciprocal condition numbers of the vectors.
00063 *
00064 *  INFO    (output) INTEGER
00065 *          = 0:  successful exit.
00066 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00067 *
00068 *  =====================================================================
00069 *
00070 *     .. Parameters ..
00071       DOUBLE PRECISION   ZERO
00072       PARAMETER          ( ZERO = 0.0D+0 )
00073 *     ..
00074 *     .. Local Scalars ..
00075       LOGICAL            DECR, EIGEN, INCR, LEFT, RIGHT, SING
00076       INTEGER            I, K
00077       DOUBLE PRECISION   ANORM, EPS, NEWGAP, OLDGAP, SAFMIN, THRESH
00078 *     ..
00079 *     .. External Functions ..
00080       LOGICAL            LSAME
00081       DOUBLE PRECISION   DLAMCH
00082       EXTERNAL           LSAME, DLAMCH
00083 *     ..
00084 *     .. Intrinsic Functions ..
00085       INTRINSIC          ABS, MAX, MIN
00086 *     ..
00087 *     .. External Subroutines ..
00088       EXTERNAL           XERBLA
00089 *     ..
00090 *     .. Executable Statements ..
00091 *
00092 *     Test the input arguments
00093 *
00094       INFO = 0
00095       EIGEN = LSAME( JOB, 'E' )
00096       LEFT = LSAME( JOB, 'L' )
00097       RIGHT = LSAME( JOB, 'R' )
00098       SING = LEFT .OR. RIGHT
00099       IF( EIGEN ) THEN
00100          K = M
00101       ELSE IF( SING ) THEN
00102          K = MIN( M, N )
00103       END IF
00104       IF( .NOT.EIGEN .AND. .NOT.SING ) THEN
00105          INFO = -1
00106       ELSE IF( M.LT.0 ) THEN
00107          INFO = -2
00108       ELSE IF( K.LT.0 ) THEN
00109          INFO = -3
00110       ELSE
00111          INCR = .TRUE.
00112          DECR = .TRUE.
00113          DO 10 I = 1, K - 1
00114             IF( INCR )
00115      $         INCR = INCR .AND. D( I ).LE.D( I+1 )
00116             IF( DECR )
00117      $         DECR = DECR .AND. D( I ).GE.D( I+1 )
00118    10    CONTINUE
00119          IF( SING .AND. K.GT.0 ) THEN
00120             IF( INCR )
00121      $         INCR = INCR .AND. ZERO.LE.D( 1 )
00122             IF( DECR )
00123      $         DECR = DECR .AND. D( K ).GE.ZERO
00124          END IF
00125          IF( .NOT.( INCR .OR. DECR ) )
00126      $      INFO = -4
00127       END IF
00128       IF( INFO.NE.0 ) THEN
00129          CALL XERBLA( 'DDISNA', -INFO )
00130          RETURN
00131       END IF
00132 *
00133 *     Quick return if possible
00134 *
00135       IF( K.EQ.0 )
00136      $   RETURN
00137 *
00138 *     Compute reciprocal condition numbers
00139 *
00140       IF( K.EQ.1 ) THEN
00141          SEP( 1 ) = DLAMCH( 'O' )
00142       ELSE
00143          OLDGAP = ABS( D( 2 )-D( 1 ) )
00144          SEP( 1 ) = OLDGAP
00145          DO 20 I = 2, K - 1
00146             NEWGAP = ABS( D( I+1 )-D( I ) )
00147             SEP( I ) = MIN( OLDGAP, NEWGAP )
00148             OLDGAP = NEWGAP
00149    20    CONTINUE
00150          SEP( K ) = OLDGAP
00151       END IF
00152       IF( SING ) THEN
00153          IF( ( LEFT .AND. M.GT.N ) .OR. ( RIGHT .AND. M.LT.N ) ) THEN
00154             IF( INCR )
00155      $         SEP( 1 ) = MIN( SEP( 1 ), D( 1 ) )
00156             IF( DECR )
00157      $         SEP( K ) = MIN( SEP( K ), D( K ) )
00158          END IF
00159       END IF
00160 *
00161 *     Ensure that reciprocal condition numbers are not less than
00162 *     threshold, in order to limit the size of the error bound
00163 *
00164       EPS = DLAMCH( 'E' )
00165       SAFMIN = DLAMCH( 'S' )
00166       ANORM = MAX( ABS( D( 1 ) ), ABS( D( K ) ) )
00167       IF( ANORM.EQ.ZERO ) THEN
00168          THRESH = EPS
00169       ELSE
00170          THRESH = MAX( EPS*ANORM, SAFMIN )
00171       END IF
00172       DO 30 I = 1, K
00173          SEP( I ) = MAX( SEP( I ), THRESH )
00174    30 CONTINUE
00175 *
00176       RETURN
00177 *
00178 *     End of DDISNA
00179 *
00180       END