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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ddisna.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack