SRC/cptcon.f

   1 *> \brief \b CPTCON
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CPTCON + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cptcon.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cptcon.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cptcon.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CPTCON( N, D, E, ANORM, RCOND, RWORK, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            INFO, N
  25 *       REAL               ANORM, RCOND
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       REAL               D( * ), RWORK( * )
  29 *       COMPLEX            E( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> CPTCON computes the reciprocal of the condition number (in the
  39 *> 1-norm) of a complex Hermitian positive definite tridiagonal matrix
  40 *> using the factorization A = L*D*L**H or A = U**H*D*U computed by
  41 *> CPTTRF.
  42 *>
  43 *> Norm(inv(A)) is computed by a direct method, and the reciprocal of
  44 *> the condition number is computed as
  45 *>                  RCOND = 1 / (ANORM * norm(inv(A))).
  46 *> \endverbatim
  47 *
  48 *  Arguments:
  49 *  ==========
  50 *
  51 *> \param[in] N
  52 *> \verbatim
  53 *>          N is INTEGER
  54 *>          The order of the matrix A.  N >= 0.
  55 *> \endverbatim
  56 *>
  57 *> \param[in] D
  58 *> \verbatim
  59 *>          D is REAL array, dimension (N)
  60 *>          The n diagonal elements of the diagonal matrix D from the
  61 *>          factorization of A, as computed by CPTTRF.
  62 *> \endverbatim
  63 *>
  64 *> \param[in] E
  65 *> \verbatim
  66 *>          E is COMPLEX array, dimension (N-1)
  67 *>          The (n-1) off-diagonal elements of the unit bidiagonal factor
  68 *>          U or L from the factorization of A, as computed by CPTTRF.
  69 *> \endverbatim
  70 *>
  71 *> \param[in] ANORM
  72 *> \verbatim
  73 *>          ANORM is REAL
  74 *>          The 1-norm of the original matrix A.
  75 *> \endverbatim
  76 *>
  77 *> \param[out] RCOND
  78 *> \verbatim
  79 *>          RCOND is REAL
  80 *>          The reciprocal of the condition number of the matrix A,
  81 *>          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is the
  82 *>          1-norm of inv(A) computed in this routine.
  83 *> \endverbatim
  84 *>
  85 *> \param[out] RWORK
  86 *> \verbatim
  87 *>          RWORK is REAL array, dimension (N)
  88 *> \endverbatim
  89 *>
  90 *> \param[out] INFO
  91 *> \verbatim
  92 *>          INFO is INTEGER
  93 *>          = 0:  successful exit
  94 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  95 *> \endverbatim
  96 *
  97 *  Authors:
  98 *  ========
  99 *
 100 *> \author Univ. of Tennessee
 101 *> \author Univ. of California Berkeley
 102 *> \author Univ. of Colorado Denver
 103 *> \author NAG Ltd.
 104 *
 105 *> \date September 2012
 106 *
 107 *> \ingroup complexPTcomputational
 108 *
 109 *> \par Further Details:
 110 *  =====================
 111 *>
 112 *> \verbatim
 113 *>
 114 *>  The method used is described in Nicholas J. Higham, "Efficient
 115 *>  Algorithms for Computing the Condition Number of a Tridiagonal
 116 *>  Matrix", SIAM J. Sci. Stat. Comput., Vol. 7, No. 1, January 1986.
 117 *> \endverbatim
 118 *>
 119 *  =====================================================================
 120       SUBROUTINE CPTCON( N, D, E, ANORM, RCOND, RWORK, INFO )
 121 *
 122 *  -- LAPACK computational routine (version 3.4.2) --
 123 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 124 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 125 *     September 2012
 126 *
 127 *     .. Scalar Arguments ..
 128       INTEGER            INFO, N
 129       REAL               ANORM, RCOND
 130 *     ..
 131 *     .. Array Arguments ..
 132       REAL               D( * ), RWORK( * )
 133       COMPLEX            E( * )
 134 *     ..
 135 *
 136 *  =====================================================================
 137 *
 138 *     .. Parameters ..
 139       REAL               ONE, ZERO
 140       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
 141 *     ..
 142 *     .. Local Scalars ..
 143       INTEGER            I, IX
 144       REAL               AINVNM
 145 *     ..
 146 *     .. External Functions ..
 147       INTEGER            ISAMAX
 148       EXTERNAL           ISAMAX
 149 *     ..
 150 *     .. External Subroutines ..
 151       EXTERNAL           XERBLA
 152 *     ..
 153 *     .. Intrinsic Functions ..
 154       INTRINSIC          ABS
 155 *     ..
 156 *     .. Executable Statements ..
 157 *
 158 *     Test the input arguments.
 159 *
 160       INFO = 0
 161       IF( N.LT.0 ) THEN
 162          INFO = -1
 163       ELSE IF( ANORM.LT.ZERO ) THEN
 164          INFO = -4
 165       END IF
 166       IF( INFO.NE.0 ) THEN
 167          CALL XERBLA( 'CPTCON', -INFO )
 168          RETURN
 169       END IF
 170 *
 171 *     Quick return if possible
 172 *
 173       RCOND = ZERO
 174       IF( N.EQ.0 ) THEN
 175          RCOND = ONE
 176          RETURN
 177       ELSE IF( ANORM.EQ.ZERO ) THEN
 178          RETURN
 179       END IF
 180 *
 181 *     Check that D(1:N) is positive.
 182 *
 183       DO 10 I = 1, N
 184          IF( D( I ).LE.ZERO )
 185      $      RETURN
 186    10 CONTINUE
 187 *
 188 *     Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
 189 *
 190 *        m(i,j) =  abs(A(i,j)), i = j,
 191 *        m(i,j) = -abs(A(i,j)), i .ne. j,
 192 *
 193 *     and e = [ 1, 1, ..., 1 ]**T.  Note M(A) = M(L)*D*M(L)**H.
 194 *
 195 *     Solve M(L) * x = e.
 196 *
 197       RWORK( 1 ) = ONE
 198       DO 20 I = 2, N
 199          RWORK( I ) = ONE + RWORK( I-1 )*ABS( E( I-1 ) )
 200    20 CONTINUE
 201 *
 202 *     Solve D * M(L)**H * x = b.
 203 *
 204       RWORK( N ) = RWORK( N ) / D( N )
 205       DO 30 I = N - 1, 1, -1
 206          RWORK( I ) = RWORK( I ) / D( I ) + RWORK( I+1 )*ABS( E( I ) )
 207    30 CONTINUE
 208 *
 209 *     Compute AINVNM = max(x(i)), 1<=i<=n.
 210 *
 211       IX = ISAMAX( N, RWORK, 1 )
 212       AINVNM = ABS( RWORK( IX ) )
 213 *
 214 *     Compute the reciprocal condition number.
 215 *
 216       IF( AINVNM.NE.ZERO )
 217      $   RCOND = ( ONE / AINVNM ) / ANORM
 218 *
 219       RETURN
 220 *
 221 *     End of CPTCON
 222 *
 223       END