SRC/dpocon.f

   1 *> \brief \b DPOCON
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DPOCON + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dpocon.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dpocon.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dpocon.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK,
  22 *                          INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          UPLO
  26 *       INTEGER            INFO, LDA, N
  27 *       DOUBLE PRECISION   ANORM, RCOND
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       INTEGER            IWORK( * )
  31 *       DOUBLE PRECISION   A( LDA, * ), WORK( * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> DPOCON estimates the reciprocal of the condition number (in the
  41 *> 1-norm) of a real symmetric positive definite matrix using the
  42 *> Cholesky factorization A = U**T*U or A = L*L**T computed by DPOTRF.
  43 *>
  44 *> An estimate is obtained for norm(inv(A)), and the reciprocal of the
  45 *> condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
  46 *> \endverbatim
  47 *
  48 *  Arguments:
  49 *  ==========
  50 *
  51 *> \param[in] UPLO
  52 *> \verbatim
  53 *>          UPLO is CHARACTER*1
  54 *>          = 'U':  Upper triangle of A is stored;
  55 *>          = 'L':  Lower triangle of A is stored.
  56 *> \endverbatim
  57 *>
  58 *> \param[in] N
  59 *> \verbatim
  60 *>          N is INTEGER
  61 *>          The order of the matrix A.  N >= 0.
  62 *> \endverbatim
  63 *>
  64 *> \param[in] A
  65 *> \verbatim
  66 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
  67 *>          The triangular factor U or L from the Cholesky factorization
  68 *>          A = U**T*U or A = L*L**T, as computed by DPOTRF.
  69 *> \endverbatim
  70 *>
  71 *> \param[in] LDA
  72 *> \verbatim
  73 *>          LDA is INTEGER
  74 *>          The leading dimension of the array A.  LDA >= max(1,N).
  75 *> \endverbatim
  76 *>
  77 *> \param[in] ANORM
  78 *> \verbatim
  79 *>          ANORM is DOUBLE PRECISION
  80 *>          The 1-norm (or infinity-norm) of the symmetric matrix A.
  81 *> \endverbatim
  82 *>
  83 *> \param[out] RCOND
  84 *> \verbatim
  85 *>          RCOND is DOUBLE PRECISION
  86 *>          The reciprocal of the condition number of the matrix A,
  87 *>          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
  88 *>          estimate of the 1-norm of inv(A) computed in this routine.
  89 *> \endverbatim
  90 *>
  91 *> \param[out] WORK
  92 *> \verbatim
  93 *>          WORK is DOUBLE PRECISION array, dimension (3*N)
  94 *> \endverbatim
  95 *>
  96 *> \param[out] IWORK
  97 *> \verbatim
  98 *>          IWORK is INTEGER array, dimension (N)
  99 *> \endverbatim
 100 *>
 101 *> \param[out] INFO
 102 *> \verbatim
 103 *>          INFO is INTEGER
 104 *>          = 0:  successful exit
 105 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 106 *> \endverbatim
 107 *
 108 *  Authors:
 109 *  ========
 110 *
 111 *> \author Univ. of Tennessee
 112 *> \author Univ. of California Berkeley
 113 *> \author Univ. of Colorado Denver
 114 *> \author NAG Ltd.
 115 *
 116 *> \date November 2011
 117 *
 118 *> \ingroup doublePOcomputational
 119 *
 120 *  =====================================================================
 121       SUBROUTINE DPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK,
 122      $                   INFO )
 123 *
 124 *  -- LAPACK computational routine (version 3.4.0) --
 125 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 126 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 127 *     November 2011
 128 *
 129 *     .. Scalar Arguments ..
 130       CHARACTER          UPLO
 131       INTEGER            INFO, LDA, N
 132       DOUBLE PRECISION   ANORM, RCOND
 133 *     ..
 134 *     .. Array Arguments ..
 135       INTEGER            IWORK( * )
 136       DOUBLE PRECISION   A( LDA, * ), WORK( * )
 137 *     ..
 138 *
 139 *  =====================================================================
 140 *
 141 *     .. Parameters ..
 142       DOUBLE PRECISION   ONE, ZERO
 143       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 144 *     ..
 145 *     .. Local Scalars ..
 146       LOGICAL            UPPER
 147       CHARACTER          NORMIN
 148       INTEGER            IX, KASE
 149       DOUBLE PRECISION   AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
 150 *     ..
 151 *     .. Local Arrays ..
 152       INTEGER            ISAVE( 3 )
 153 *     ..
 154 *     .. External Functions ..
 155       LOGICAL            LSAME
 156       INTEGER            IDAMAX
 157       DOUBLE PRECISION   DLAMCH
 158       EXTERNAL           LSAME, IDAMAX, DLAMCH
 159 *     ..
 160 *     .. External Subroutines ..
 161       EXTERNAL           DLACN2, DLATRS, DRSCL, XERBLA
 162 *     ..
 163 *     .. Intrinsic Functions ..
 164       INTRINSIC          ABS, MAX
 165 *     ..
 166 *     .. Executable Statements ..
 167 *
 168 *     Test the input parameters.
 169 *
 170       INFO = 0
 171       UPPER = LSAME( UPLO, 'U' )
 172       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 173          INFO = -1
 174       ELSE IF( N.LT.0 ) THEN
 175          INFO = -2
 176       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 177          INFO = -4
 178       ELSE IF( ANORM.LT.ZERO ) THEN
 179          INFO = -5
 180       END IF
 181       IF( INFO.NE.0 ) THEN
 182          CALL XERBLA( 'DPOCON', -INFO )
 183          RETURN
 184       END IF
 185 *
 186 *     Quick return if possible
 187 *
 188       RCOND = ZERO
 189       IF( N.EQ.0 ) THEN
 190          RCOND = ONE
 191          RETURN
 192       ELSE IF( ANORM.EQ.ZERO ) THEN
 193          RETURN
 194       END IF
 195 *
 196       SMLNUM = DLAMCH( 'Safe minimum' )
 197 *
 198 *     Estimate the 1-norm of inv(A).
 199 *
 200       KASE = 0
 201       NORMIN = 'N'
 202    10 CONTINUE
 203       CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
 204       IF( KASE.NE.0 ) THEN
 205          IF( UPPER ) THEN
 206 *
 207 *           Multiply by inv(U**T).
 208 *
 209             CALL DLATRS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, A,
 210      $                   LDA, WORK, SCALEL, WORK( 2*N+1 ), INFO )
 211             NORMIN = 'Y'
 212 *
 213 *           Multiply by inv(U).
 214 *
 215             CALL DLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
 216      $                   A, LDA, WORK, SCALEU, WORK( 2*N+1 ), INFO )
 217          ELSE
 218 *
 219 *           Multiply by inv(L).
 220 *
 221             CALL DLATRS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
 222      $                   A, LDA, WORK, SCALEL, WORK( 2*N+1 ), INFO )
 223             NORMIN = 'Y'
 224 *
 225 *           Multiply by inv(L**T).
 226 *
 227             CALL DLATRS( 'Lower', 'Transpose', 'Non-unit', NORMIN, N, A,
 228      $                   LDA, WORK, SCALEU, WORK( 2*N+1 ), INFO )
 229          END IF
 230 *
 231 *        Multiply by 1/SCALE if doing so will not cause overflow.
 232 *
 233          SCALE = SCALEL*SCALEU
 234          IF( SCALE.NE.ONE ) THEN
 235             IX = IDAMAX( N, WORK, 1 )
 236             IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
 237      $         GO TO 20
 238             CALL DRSCL( N, SCALE, WORK, 1 )
 239          END IF
 240          GO TO 10
 241       END IF
 242 *
 243 *     Compute the estimate of the reciprocal condition number.
 244 *
 245       IF( AINVNM.NE.ZERO )
 246      $   RCOND = ( ONE / AINVNM ) / ANORM
 247 *
 248    20 CONTINUE
 249       RETURN
 250 *
 251 *     End of DPOCON
 252 *
 253       END