001:       SUBROUTINE DDISNA( JOB, M, N, D, SEP, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          JOB
010:       INTEGER            INFO, M, N
011: *     ..
012: *     .. Array Arguments ..
013:       DOUBLE PRECISION   D( * ), SEP( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  DDISNA computes the reciprocal condition numbers for the eigenvectors
020: *  of a real symmetric or complex Hermitian matrix or for the left or
021: *  right singular vectors of a general m-by-n matrix. The reciprocal
022: *  condition number is the 'gap' between the corresponding eigenvalue or
023: *  singular value and the nearest other one.
024: *
025: *  The bound on the error, measured by angle in radians, in the I-th
026: *  computed vector is given by
027: *
028: *         DLAMCH( 'E' ) * ( ANORM / SEP( I ) )
029: *
030: *  where ANORM = 2-norm(A) = max( abs( D(j) ) ).  SEP(I) is not allowed
031: *  to be smaller than DLAMCH( 'E' )*ANORM in order to limit the size of
032: *  the error bound.
033: *
034: *  DDISNA may also be used to compute error bounds for eigenvectors of
035: *  the generalized symmetric definite eigenproblem.
036: *
037: *  Arguments
038: *  =========
039: *
040: *  JOB     (input) CHARACTER*1
041: *          Specifies for which problem the reciprocal condition numbers
042: *          should be computed:
043: *          = 'E':  the eigenvectors of a symmetric/Hermitian matrix;
044: *          = 'L':  the left singular vectors of a general matrix;
045: *          = 'R':  the right singular vectors of a general matrix.
046: *
047: *  M       (input) INTEGER
048: *          The number of rows of the matrix. M >= 0.
049: *
050: *  N       (input) INTEGER
051: *          If JOB = 'L' or 'R', the number of columns of the matrix,
052: *          in which case N >= 0. Ignored if JOB = 'E'.
053: *
054: *  D       (input) DOUBLE PRECISION array, dimension (M) if JOB = 'E'
055: *                              dimension (min(M,N)) if JOB = 'L' or 'R'
056: *          The eigenvalues (if JOB = 'E') or singular values (if JOB =
057: *          'L' or 'R') of the matrix, in either increasing or decreasing
058: *          order. If singular values, they must be non-negative.
059: *
060: *  SEP     (output) DOUBLE PRECISION array, dimension (M) if JOB = 'E'
061: *                               dimension (min(M,N)) if JOB = 'L' or 'R'
062: *          The reciprocal condition numbers of the vectors.
063: *
064: *  INFO    (output) INTEGER
065: *          = 0:  successful exit.
066: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
067: *
068: *  =====================================================================
069: *
070: *     .. Parameters ..
071:       DOUBLE PRECISION   ZERO
072:       PARAMETER          ( ZERO = 0.0D+0 )
073: *     ..
074: *     .. Local Scalars ..
075:       LOGICAL            DECR, EIGEN, INCR, LEFT, RIGHT, SING
076:       INTEGER            I, K
077:       DOUBLE PRECISION   ANORM, EPS, NEWGAP, OLDGAP, SAFMIN, THRESH
078: *     ..
079: *     .. External Functions ..
080:       LOGICAL            LSAME
081:       DOUBLE PRECISION   DLAMCH
082:       EXTERNAL           LSAME, DLAMCH
083: *     ..
084: *     .. Intrinsic Functions ..
085:       INTRINSIC          ABS, MAX, MIN
086: *     ..
087: *     .. External Subroutines ..
088:       EXTERNAL           XERBLA
089: *     ..
090: *     .. Executable Statements ..
091: *
092: *     Test the input arguments
093: *
094:       INFO = 0
095:       EIGEN = LSAME( JOB, 'E' )
096:       LEFT = LSAME( JOB, 'L' )
097:       RIGHT = LSAME( JOB, 'R' )
098:       SING = LEFT .OR. RIGHT
099:       IF( EIGEN ) THEN
100:          K = M
101:       ELSE IF( SING ) THEN
102:          K = MIN( M, N )
103:       END IF
104:       IF( .NOT.EIGEN .AND. .NOT.SING ) THEN
105:          INFO = -1
106:       ELSE IF( M.LT.0 ) THEN
107:          INFO = -2
108:       ELSE IF( K.LT.0 ) THEN
109:          INFO = -3
110:       ELSE
111:          INCR = .TRUE.
112:          DECR = .TRUE.
113:          DO 10 I = 1, K - 1
114:             IF( INCR )
115:      $         INCR = INCR .AND. D( I ).LE.D( I+1 )
116:             IF( DECR )
117:      $         DECR = DECR .AND. D( I ).GE.D( I+1 )
118:    10    CONTINUE
119:          IF( SING .AND. K.GT.0 ) THEN
120:             IF( INCR )
121:      $         INCR = INCR .AND. ZERO.LE.D( 1 )
122:             IF( DECR )
123:      $         DECR = DECR .AND. D( K ).GE.ZERO
124:          END IF
125:          IF( .NOT.( INCR .OR. DECR ) )
126:      $      INFO = -4
127:       END IF
128:       IF( INFO.NE.0 ) THEN
129:          CALL XERBLA( 'DDISNA', -INFO )
130:          RETURN
131:       END IF
132: *
133: *     Quick return if possible
134: *
135:       IF( K.EQ.0 )
136:      $   RETURN
137: *
138: *     Compute reciprocal condition numbers
139: *
140:       IF( K.EQ.1 ) THEN
141:          SEP( 1 ) = DLAMCH( 'O' )
142:       ELSE
143:          OLDGAP = ABS( D( 2 )-D( 1 ) )
144:          SEP( 1 ) = OLDGAP
145:          DO 20 I = 2, K - 1
146:             NEWGAP = ABS( D( I+1 )-D( I ) )
147:             SEP( I ) = MIN( OLDGAP, NEWGAP )
148:             OLDGAP = NEWGAP
149:    20    CONTINUE
150:          SEP( K ) = OLDGAP
151:       END IF
152:       IF( SING ) THEN
153:          IF( ( LEFT .AND. M.GT.N ) .OR. ( RIGHT .AND. M.LT.N ) ) THEN
154:             IF( INCR )
155:      $         SEP( 1 ) = MIN( SEP( 1 ), D( 1 ) )
156:             IF( DECR )
157:      $         SEP( K ) = MIN( SEP( K ), D( K ) )
158:          END IF
159:       END IF
160: *
161: *     Ensure that reciprocal condition numbers are not less than
162: *     threshold, in order to limit the size of the error bound
163: *
164:       EPS = DLAMCH( 'E' )
165:       SAFMIN = DLAMCH( 'S' )
166:       ANORM = MAX( ABS( D( 1 ) ), ABS( D( K ) ) )
167:       IF( ANORM.EQ.ZERO ) THEN
168:          THRESH = EPS
169:       ELSE
170:          THRESH = MAX( EPS*ANORM, SAFMIN )
171:       END IF
172:       DO 30 I = 1, K
173:          SEP( I ) = MAX( SEP( I ), THRESH )
174:    30 CONTINUE
175: *
176:       RETURN
177: *
178: *     End of DDISNA
179: *
180:       END
181: