SRC/dptcon.f

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