1 SUBROUTINE DDISNA( JOB, M, N, D, SEP, INFO )
2 *
3 * -- LAPACK routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER JOB
10 INTEGER INFO, M, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION D( * ), SEP( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DDISNA computes the reciprocal condition numbers for the eigenvectors
20 * of a real symmetric or complex Hermitian matrix or for the left or
21 * right singular vectors of a general m-by-n matrix. The reciprocal
22 * condition number is the 'gap' between the corresponding eigenvalue or
23 * singular value and the nearest other one.
24 *
25 * The bound on the error, measured by angle in radians, in the I-th
26 * computed vector is given by
27 *
28 * DLAMCH( 'E' ) * ( ANORM / SEP( I ) )
29 *
30 * where ANORM = 2-norm(A) = max( abs( D(j) ) ). SEP(I) is not allowed
31 * to be smaller than DLAMCH( 'E' )*ANORM in order to limit the size of
32 * the error bound.
33 *
34 * DDISNA may also be used to compute error bounds for eigenvectors of
35 * the generalized symmetric definite eigenproblem.
36 *
37 * Arguments
38 * =========
39 *
40 * JOB (input) CHARACTER*1
41 * Specifies for which problem the reciprocal condition numbers
42 * should be computed:
43 * = 'E': the eigenvectors of a symmetric/Hermitian matrix;
44 * = 'L': the left singular vectors of a general matrix;
45 * = 'R': the right singular vectors of a general matrix.
46 *
47 * M (input) INTEGER
48 * The number of rows of the matrix. M >= 0.
49 *
50 * N (input) INTEGER
51 * If JOB = 'L' or 'R', the number of columns of the matrix,
52 * in which case N >= 0. Ignored if JOB = 'E'.
53 *
54 * D (input) DOUBLE PRECISION array, dimension (M) if JOB = 'E'
55 * dimension (min(M,N)) if JOB = 'L' or 'R'
56 * The eigenvalues (if JOB = 'E') or singular values (if JOB =
57 * 'L' or 'R') of the matrix, in either increasing or decreasing
58 * order. If singular values, they must be non-negative.
59 *
60 * SEP (output) DOUBLE PRECISION array, dimension (M) if JOB = 'E'
61 * dimension (min(M,N)) if JOB = 'L' or 'R'
62 * The reciprocal condition numbers of the vectors.
63 *
64 * INFO (output) INTEGER
65 * = 0: successful exit.
66 * < 0: if INFO = -i, the i-th argument had an illegal value.
67 *
68 * =====================================================================
69 *
70 * .. Parameters ..
71 DOUBLE PRECISION ZERO
72 PARAMETER ( ZERO = 0.0D+0 )
73 * ..
74 * .. Local Scalars ..
75 LOGICAL DECR, EIGEN, INCR, LEFT, RIGHT, SING
76 INTEGER I, K
77 DOUBLE PRECISION ANORM, EPS, NEWGAP, OLDGAP, SAFMIN, THRESH
78 * ..
79 * .. External Functions ..
80 LOGICAL LSAME
81 DOUBLE PRECISION DLAMCH
82 EXTERNAL LSAME, DLAMCH
83 * ..
84 * .. Intrinsic Functions ..
85 INTRINSIC ABS, MAX, MIN
86 * ..
87 * .. External Subroutines ..
88 EXTERNAL XERBLA
89 * ..
90 * .. Executable Statements ..
91 *
92 * Test the input arguments
93 *
94 INFO = 0
95 EIGEN = LSAME( JOB, 'E' )
96 LEFT = LSAME( JOB, 'L' )
97 RIGHT = LSAME( JOB, 'R' )
98 SING = LEFT .OR. RIGHT
99 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
2 *
3 * -- LAPACK routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER JOB
10 INTEGER INFO, M, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION D( * ), SEP( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DDISNA computes the reciprocal condition numbers for the eigenvectors
20 * of a real symmetric or complex Hermitian matrix or for the left or
21 * right singular vectors of a general m-by-n matrix. The reciprocal
22 * condition number is the 'gap' between the corresponding eigenvalue or
23 * singular value and the nearest other one.
24 *
25 * The bound on the error, measured by angle in radians, in the I-th
26 * computed vector is given by
27 *
28 * DLAMCH( 'E' ) * ( ANORM / SEP( I ) )
29 *
30 * where ANORM = 2-norm(A) = max( abs( D(j) ) ). SEP(I) is not allowed
31 * to be smaller than DLAMCH( 'E' )*ANORM in order to limit the size of
32 * the error bound.
33 *
34 * DDISNA may also be used to compute error bounds for eigenvectors of
35 * the generalized symmetric definite eigenproblem.
36 *
37 * Arguments
38 * =========
39 *
40 * JOB (input) CHARACTER*1
41 * Specifies for which problem the reciprocal condition numbers
42 * should be computed:
43 * = 'E': the eigenvectors of a symmetric/Hermitian matrix;
44 * = 'L': the left singular vectors of a general matrix;
45 * = 'R': the right singular vectors of a general matrix.
46 *
47 * M (input) INTEGER
48 * The number of rows of the matrix. M >= 0.
49 *
50 * N (input) INTEGER
51 * If JOB = 'L' or 'R', the number of columns of the matrix,
52 * in which case N >= 0. Ignored if JOB = 'E'.
53 *
54 * D (input) DOUBLE PRECISION array, dimension (M) if JOB = 'E'
55 * dimension (min(M,N)) if JOB = 'L' or 'R'
56 * The eigenvalues (if JOB = 'E') or singular values (if JOB =
57 * 'L' or 'R') of the matrix, in either increasing or decreasing
58 * order. If singular values, they must be non-negative.
59 *
60 * SEP (output) DOUBLE PRECISION array, dimension (M) if JOB = 'E'
61 * dimension (min(M,N)) if JOB = 'L' or 'R'
62 * The reciprocal condition numbers of the vectors.
63 *
64 * INFO (output) INTEGER
65 * = 0: successful exit.
66 * < 0: if INFO = -i, the i-th argument had an illegal value.
67 *
68 * =====================================================================
69 *
70 * .. Parameters ..
71 DOUBLE PRECISION ZERO
72 PARAMETER ( ZERO = 0.0D+0 )
73 * ..
74 * .. Local Scalars ..
75 LOGICAL DECR, EIGEN, INCR, LEFT, RIGHT, SING
76 INTEGER I, K
77 DOUBLE PRECISION ANORM, EPS, NEWGAP, OLDGAP, SAFMIN, THRESH
78 * ..
79 * .. External Functions ..
80 LOGICAL LSAME
81 DOUBLE PRECISION DLAMCH
82 EXTERNAL LSAME, DLAMCH
83 * ..
84 * .. Intrinsic Functions ..
85 INTRINSIC ABS, MAX, MIN
86 * ..
87 * .. External Subroutines ..
88 EXTERNAL XERBLA
89 * ..
90 * .. Executable Statements ..
91 *
92 * Test the input arguments
93 *
94 INFO = 0
95 EIGEN = LSAME( JOB, 'E' )
96 LEFT = LSAME( JOB, 'L' )
97 RIGHT = LSAME( JOB, 'R' )
98 SING = LEFT .OR. RIGHT
99 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