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