SRC/cpocon.f

   1 *> \brief \b CPOCON
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CPOCON + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cpocon.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cpocon.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cpocon.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, RWORK,
  22 *                          INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          UPLO
  26 *       INTEGER            INFO, LDA, N
  27 *       REAL               ANORM, RCOND
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       REAL               RWORK( * )
  31 *       COMPLEX            A( LDA, * ), WORK( * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> CPOCON estimates the reciprocal of the condition number (in the
  41 *> 1-norm) of a complex Hermitian positive definite matrix using the
  42 *> Cholesky factorization A = U**H*U or A = L*L**H computed by CPOTRF.
  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 COMPLEX array, dimension (LDA,N)
  67 *>          The triangular factor U or L from the Cholesky factorization
  68 *>          A = U**H*U or A = L*L**H, as computed by CPOTRF.
  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 REAL
  80 *>          The 1-norm (or infinity-norm) of the Hermitian matrix A.
  81 *> \endverbatim
  82 *>
  83 *> \param[out] RCOND
  84 *> \verbatim
  85 *>          RCOND is REAL
  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 COMPLEX array, dimension (2*N)
  94 *> \endverbatim
  95 *>
  96 *> \param[out] RWORK
  97 *> \verbatim
  98 *>          RWORK is REAL 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 complexPOcomputational
 119 *
 120 *  =====================================================================
 121       SUBROUTINE CPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, RWORK,
 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       REAL               ANORM, RCOND
 133 *     ..
 134 *     .. Array Arguments ..
 135       REAL               RWORK( * )
 136       COMPLEX            A( LDA, * ), WORK( * )
 137 *     ..
 138 *
 139 *  =====================================================================
 140 *
 141 *     .. Parameters ..
 142       REAL               ONE, ZERO
 143       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
 144 *     ..
 145 *     .. Local Scalars ..
 146       LOGICAL            UPPER
 147       CHARACTER          NORMIN
 148       INTEGER            IX, KASE
 149       REAL               AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
 150       COMPLEX            ZDUM
 151 *     ..
 152 *     .. Local Arrays ..
 153       INTEGER            ISAVE( 3 )
 154 *     ..
 155 *     .. External Functions ..
 156       LOGICAL            LSAME
 157       INTEGER            ICAMAX
 158       REAL               SLAMCH
 159       EXTERNAL           LSAME, ICAMAX, SLAMCH
 160 *     ..
 161 *     .. External Subroutines ..
 162       EXTERNAL           CLACN2, CLATRS, CSRSCL, XERBLA
 163 *     ..
 164 *     .. Intrinsic Functions ..
 165       INTRINSIC          ABS, AIMAG, MAX, REAL
 166 *     ..
 167 *     .. Statement Functions ..
 168       REAL               CABS1
 169 *     ..
 170 *     .. Statement Function definitions ..
 171       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
 172 *     ..
 173 *     .. Executable Statements ..
 174 *
 175 *     Test the input parameters.
 176 *
 177       INFO = 0
 178       UPPER = LSAME( UPLO, 'U' )
 179       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 180          INFO = -1
 181       ELSE IF( N.LT.0 ) THEN
 182          INFO = -2
 183       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 184          INFO = -4
 185       ELSE IF( ANORM.LT.ZERO ) THEN
 186          INFO = -5
 187       END IF
 188       IF( INFO.NE.0 ) THEN
 189          CALL XERBLA( 'CPOCON', -INFO )
 190          RETURN
 191       END IF
 192 *
 193 *     Quick return if possible
 194 *
 195       RCOND = ZERO
 196       IF( N.EQ.0 ) THEN
 197          RCOND = ONE
 198          RETURN
 199       ELSE IF( ANORM.EQ.ZERO ) THEN
 200          RETURN
 201       END IF
 202 *
 203       SMLNUM = SLAMCH( 'Safe minimum' )
 204 *
 205 *     Estimate the 1-norm of inv(A).
 206 *
 207       KASE = 0
 208       NORMIN = 'N'
 209    10 CONTINUE
 210       CALL CLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
 211       IF( KASE.NE.0 ) THEN
 212          IF( UPPER ) THEN
 213 *
 214 *           Multiply by inv(U**H).
 215 *
 216             CALL CLATRS( 'Upper', 'Conjugate transpose', 'Non-unit',
 217      $                   NORMIN, N, A, LDA, WORK, SCALEL, RWORK, INFO )
 218             NORMIN = 'Y'
 219 *
 220 *           Multiply by inv(U).
 221 *
 222             CALL CLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
 223      $                   A, LDA, WORK, SCALEU, RWORK, INFO )
 224          ELSE
 225 *
 226 *           Multiply by inv(L).
 227 *
 228             CALL CLATRS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
 229      $                   A, LDA, WORK, SCALEL, RWORK, INFO )
 230             NORMIN = 'Y'
 231 *
 232 *           Multiply by inv(L**H).
 233 *
 234             CALL CLATRS( 'Lower', 'Conjugate transpose', 'Non-unit',
 235      $                   NORMIN, N, A, LDA, WORK, SCALEU, RWORK, INFO )
 236          END IF
 237 *
 238 *        Multiply by 1/SCALE if doing so will not cause overflow.
 239 *
 240          SCALE = SCALEL*SCALEU
 241          IF( SCALE.NE.ONE ) THEN
 242             IX = ICAMAX( N, WORK, 1 )
 243             IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
 244      $         GO TO 20
 245             CALL CSRSCL( N, SCALE, WORK, 1 )
 246          END IF
 247          GO TO 10
 248       END IF
 249 *
 250 *     Compute the estimate of the reciprocal condition number.
 251 *
 252       IF( AINVNM.NE.ZERO )
 253      $   RCOND = ( ONE / AINVNM ) / ANORM
 254 *
 255    20 CONTINUE
 256       RETURN
 257 *
 258 *     End of CPOCON
 259 *
 260       END