|
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 |
SUBROUTINE SDISNA( JOB, M, N, D, SEP, INFO )
* * -- LAPACK routine (version 3.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2006 * * .. Scalar Arguments .. CHARACTER JOB INTEGER INFO, M, N * .. * .. Array Arguments .. REAL D( * ), SEP( * ) * .. * * Purpose * ======= * * SDISNA computes the reciprocal condition numbers for the eigenvectors * of a real symmetric or complex Hermitian matrix or for the left or * right singular vectors of a general m-by-n matrix. The reciprocal * condition number is the 'gap' between the corresponding eigenvalue or * singular value and the nearest other one. * * The bound on the error, measured by angle in radians, in the I-th * computed vector is given by * * SLAMCH( 'E' ) * ( ANORM / SEP( I ) ) * * where ANORM = 2-norm(A) = max( abs( D(j) ) ). SEP(I) is not allowed * to be smaller than SLAMCH( 'E' )*ANORM in order to limit the size of * the error bound. * * SDISNA may also be used to compute error bounds for eigenvectors of * the generalized symmetric definite eigenproblem. * * Arguments * ========= * * JOB (input) CHARACTER*1 * Specifies for which problem the reciprocal condition numbers * should be computed: * = 'E': the eigenvectors of a symmetric/Hermitian matrix; * = 'L': the left singular vectors of a general matrix; * = 'R': the right singular vectors of a general matrix. * * M (input) INTEGER * The number of rows of the matrix. M >= 0. * * N (input) INTEGER * If JOB = 'L' or 'R', the number of columns of the matrix, * in which case N >= 0. Ignored if JOB = 'E'. * * D (input) REAL array, dimension (M) if JOB = 'E' * dimension (min(M,N)) if JOB = 'L' or 'R' * The eigenvalues (if JOB = 'E') or singular values (if JOB = * 'L' or 'R') of the matrix, in either increasing or decreasing * order. If singular values, they must be non-negative. * * SEP (output) REAL array, dimension (M) if JOB = 'E' * dimension (min(M,N)) if JOB = 'L' or 'R' * The reciprocal condition numbers of the vectors. * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * * ===================================================================== * * .. Parameters .. REAL ZERO PARAMETER ( ZERO = 0.0E+0 ) * .. * .. Local Scalars .. LOGICAL DECR, EIGEN, INCR, LEFT, RIGHT, SING INTEGER I, K REAL ANORM, EPS, NEWGAP, OLDGAP, SAFMIN, THRESH * .. * .. External Functions .. LOGICAL LSAME REAL SLAMCH EXTERNAL LSAME, SLAMCH * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN * .. * .. External Subroutines .. EXTERNAL XERBLA * .. * .. Executable Statements .. * * Test the input arguments * INFO = 0 EIGEN = LSAME( JOB, 'E' ) LEFT = LSAME( JOB, 'L' ) RIGHT = LSAME( JOB, 'R' ) SING = LEFT .OR. RIGHT IF( EIGEN ) THEN K = M ELSE IF( SING ) THEN K = MIN( M, N ) END IF IF( .NOT.EIGEN .AND. .NOT.SING ) THEN INFO = -1 ELSE IF( M.LT.0 ) THEN INFO = -2 ELSE IF( K.LT.0 ) THEN INFO = -3 ELSE INCR = .TRUE. DECR = .TRUE. DO 10 I = 1, K - 1 IF( INCR ) $ INCR = INCR .AND. D( I ).LE.D( I+1 ) IF( DECR ) $ DECR = DECR .AND. D( I ).GE.D( I+1 ) 10 CONTINUE IF( SING .AND. K.GT.0 ) THEN IF( INCR ) $ INCR = INCR .AND. ZERO.LE.D( 1 ) IF( DECR ) $ DECR = DECR .AND. D( K ).GE.ZERO END IF IF( .NOT.( INCR .OR. DECR ) ) $ INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'SDISNA', -INFO ) RETURN END IF * * Quick return if possible * IF( K.EQ.0 ) $ RETURN * * Compute reciprocal condition numbers * IF( K.EQ.1 ) THEN SEP( 1 ) = SLAMCH( 'O' ) ELSE OLDGAP = ABS( D( 2 )-D( 1 ) ) SEP( 1 ) = OLDGAP DO 20 I = 2, K - 1 NEWGAP = ABS( D( I+1 )-D( I ) ) SEP( I ) = MIN( OLDGAP, NEWGAP ) OLDGAP = NEWGAP 20 CONTINUE SEP( K ) = OLDGAP END IF IF( SING ) THEN IF( ( LEFT .AND. M.GT.N ) .OR. ( RIGHT .AND. M.LT.N ) ) THEN IF( INCR ) $ SEP( 1 ) = MIN( SEP( 1 ), D( 1 ) ) IF( DECR ) $ SEP( K ) = MIN( SEP( K ), D( K ) ) END IF END IF * * Ensure that reciprocal condition numbers are not less than * threshold, in order to limit the size of the error bound * EPS = SLAMCH( 'E' ) SAFMIN = SLAMCH( 'S' ) ANORM = MAX( ABS( D( 1 ) ), ABS( D( K ) ) ) IF( ANORM.EQ.ZERO ) THEN THRESH = EPS ELSE THRESH = MAX( EPS*ANORM, SAFMIN ) END IF DO 30 I = 1, K SEP( I ) = MAX( SEP( I ), THRESH ) 30 CONTINUE * RETURN * * End of SDISNA * END |