SRC/ddisna.f

   1 *> \brief \b DDISNA
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DDISNA + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ddisna.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ddisna.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ddisna.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DDISNA( JOB, M, N, D, SEP, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          JOB
  25 *       INTEGER            INFO, M, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       DOUBLE PRECISION   D( * ), SEP( * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> DDISNA computes the reciprocal condition numbers for the eigenvectors
  38 *> of a real symmetric or complex Hermitian matrix or for the left or
  39 *> right singular vectors of a general m-by-n matrix. The reciprocal
  40 *> condition number is the 'gap' between the corresponding eigenvalue or
  41 *> singular value and the nearest other one.
  42 *>
  43 *> The bound on the error, measured by angle in radians, in the I-th
  44 *> computed vector is given by
  45 *>
  46 *>        DLAMCH( 'E' ) * ( ANORM / SEP( I ) )
  47 *>
  48 *> where ANORM = 2-norm(A) = max( abs( D(j) ) ).  SEP(I) is not allowed
  49 *> to be smaller than DLAMCH( 'E' )*ANORM in order to limit the size of
  50 *> the error bound.
  51 *>
  52 *> DDISNA may also be used to compute error bounds for eigenvectors of
  53 *> the generalized symmetric definite eigenproblem.
  54 *> \endverbatim
  55 *
  56 *  Arguments:
  57 *  ==========
  58 *
  59 *> \param[in] JOB
  60 *> \verbatim
  61 *>          JOB is CHARACTER*1
  62 *>          Specifies for which problem the reciprocal condition numbers
  63 *>          should be computed:
  64 *>          = 'E':  the eigenvectors of a symmetric/Hermitian matrix;
  65 *>          = 'L':  the left singular vectors of a general matrix;
  66 *>          = 'R':  the right singular vectors of a general matrix.
  67 *> \endverbatim
  68 *>
  69 *> \param[in] M
  70 *> \verbatim
  71 *>          M is INTEGER
  72 *>          The number of rows of the matrix. M >= 0.
  73 *> \endverbatim
  74 *>
  75 *> \param[in] N
  76 *> \verbatim
  77 *>          N is INTEGER
  78 *>          If JOB = 'L' or 'R', the number of columns of the matrix,
  79 *>          in which case N >= 0. Ignored if JOB = 'E'.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] D
  83 *> \verbatim
  84 *>          D is DOUBLE PRECISION array, dimension (M) if JOB = 'E'
  85 *>                              dimension (min(M,N)) if JOB = 'L' or 'R'
  86 *>          The eigenvalues (if JOB = 'E') or singular values (if JOB =
  87 *>          'L' or 'R') of the matrix, in either increasing or decreasing
  88 *>          order. If singular values, they must be non-negative.
  89 *> \endverbatim
  90 *>
  91 *> \param[out] SEP
  92 *> \verbatim
  93 *>          SEP is DOUBLE PRECISION array, dimension (M) if JOB = 'E'
  94 *>                               dimension (min(M,N)) if JOB = 'L' or 'R'
  95 *>          The reciprocal condition numbers of the vectors.
  96 *> \endverbatim
  97 *>
  98 *> \param[out] INFO
  99 *> \verbatim
 100 *>          INFO is INTEGER
 101 *>          = 0:  successful exit.
 102 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 103 *> \endverbatim
 104 *
 105 *  Authors:
 106 *  ========
 107 *
 108 *> \author Univ. of Tennessee
 109 *> \author Univ. of California Berkeley
 110 *> \author Univ. of Colorado Denver
 111 *> \author NAG Ltd.
 112 *
 113 *> \date November 2011
 114 *
 115 *> \ingroup auxOTHERcomputational
 116 *
 117 *  =====================================================================
 118       SUBROUTINE DDISNA( JOB, M, N, D, SEP, INFO )
 119 *
 120 *  -- LAPACK computational routine (version 3.4.0) --
 121 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 122 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 123 *     November 2011
 124 *
 125 *     .. Scalar Arguments ..
 126       CHARACTER          JOB
 127       INTEGER            INFO, M, N
 128 *     ..
 129 *     .. Array Arguments ..
 130       DOUBLE PRECISION   D( * ), SEP( * )
 131 *     ..
 132 *
 133 *  =====================================================================
 134 *
 135 *     .. Parameters ..
 136       DOUBLE PRECISION   ZERO
 137       PARAMETER          ( ZERO = 0.0D+0 )
 138 *     ..
 139 *     .. Local Scalars ..
 140       LOGICAL            DECR, EIGEN, INCR, LEFT, RIGHT, SING
 141       INTEGER            I, K
 142       DOUBLE PRECISION   ANORM, EPS, NEWGAP, OLDGAP, SAFMIN, THRESH
 143 *     ..
 144 *     .. External Functions ..
 145       LOGICAL            LSAME
 146       DOUBLE PRECISION   DLAMCH
 147       EXTERNAL           LSAME, DLAMCH
 148 *     ..
 149 *     .. Intrinsic Functions ..
 150       INTRINSIC          ABS, MAX, MIN
 151 *     ..
 152 *     .. External Subroutines ..
 153       EXTERNAL           XERBLA
 154 *     ..
 155 *     .. Executable Statements ..
 156 *
 157 *     Test the input arguments
 158 *
 159       INFO = 0
 160       EIGEN = LSAME( JOB, 'E' )
 161       LEFT = LSAME( JOB, 'L' )
 162       RIGHT = LSAME( JOB, 'R' )
 163       SING = LEFT .OR. RIGHT
 164       IF( EIGEN ) THEN
 165          K = M
 166       ELSE IF( SING ) THEN
 167          K = MIN( M, N )
 168       END IF
 169       IF( .NOT.EIGEN .AND. .NOT.SING ) THEN
 170          INFO = -1
 171       ELSE IF( M.LT.0 ) THEN
 172          INFO = -2
 173       ELSE IF( K.LT.0 ) THEN
 174          INFO = -3
 175       ELSE
 176          INCR = .TRUE.
 177          DECR = .TRUE.
 178          DO 10 I = 1, K - 1
 179             IF( INCR )
 180      $         INCR = INCR .AND. D( I ).LE.D( I+1 )
 181             IF( DECR )
 182      $         DECR = DECR .AND. D( I ).GE.D( I+1 )
 183    10    CONTINUE
 184          IF( SING .AND. K.GT.0 ) THEN
 185             IF( INCR )
 186      $         INCR = INCR .AND. ZERO.LE.D( 1 )
 187             IF( DECR )
 188      $         DECR = DECR .AND. D( K ).GE.ZERO
 189          END IF
 190          IF( .NOT.( INCR .OR. DECR ) )
 191      $      INFO = -4
 192       END IF
 193       IF( INFO.NE.0 ) THEN
 194          CALL XERBLA( 'DDISNA', -INFO )
 195          RETURN
 196       END IF
 197 *
 198 *     Quick return if possible
 199 *
 200       IF( K.EQ.0 )
 201      $   RETURN
 202 *
 203 *     Compute reciprocal condition numbers
 204 *
 205       IF( K.EQ.1 ) THEN
 206          SEP( 1 ) = DLAMCH( 'O' )
 207       ELSE
 208          OLDGAP = ABS( D( 2 )-D( 1 ) )
 209          SEP( 1 ) = OLDGAP
 210          DO 20 I = 2, K - 1
 211             NEWGAP = ABS( D( I+1 )-D( I ) )
 212             SEP( I ) = MIN( OLDGAP, NEWGAP )
 213             OLDGAP = NEWGAP
 214    20    CONTINUE
 215          SEP( K ) = OLDGAP
 216       END IF
 217       IF( SING ) THEN
 218          IF( ( LEFT .AND. M.GT.N ) .OR. ( RIGHT .AND. M.LT.N ) ) THEN
 219             IF( INCR )
 220      $         SEP( 1 ) = MIN( SEP( 1 ), D( 1 ) )
 221             IF( DECR )
 222      $         SEP( K ) = MIN( SEP( K ), D( K ) )
 223          END IF
 224       END IF
 225 *
 226 *     Ensure that reciprocal condition numbers are not less than
 227 *     threshold, in order to limit the size of the error bound
 228 *
 229       EPS = DLAMCH( 'E' )
 230       SAFMIN = DLAMCH( 'S' )
 231       ANORM = MAX( ABS( D( 1 ) ), ABS( D( K ) ) )
 232       IF( ANORM.EQ.ZERO ) THEN
 233          THRESH = EPS
 234       ELSE
 235          THRESH = MAX( EPS*ANORM, SAFMIN )
 236       END IF
 237       DO 30 I = 1, K
 238          SEP( I ) = MAX( SEP( I ), THRESH )
 239    30 CONTINUE
 240 *
 241       RETURN
 242 *
 243 *     End of DDISNA
 244 *
 245       END