LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dptcon.f
Go to the documentation of this file.
00001 *> \brief \b DPTCON
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DPTCON + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dptcon.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dptcon.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dptcon.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DPTCON( N, D, E, ANORM, RCOND, WORK, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            INFO, N
00025 *       DOUBLE PRECISION   ANORM, RCOND
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       DOUBLE PRECISION   D( * ), E( * ), WORK( * )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> DPTCON computes the reciprocal of the condition number (in the
00038 *> 1-norm) of a real symmetric positive definite tridiagonal matrix
00039 *> using the factorization A = L*D*L**T or A = U**T*D*U computed by
00040 *> DPTTRF.
00041 *>
00042 *> Norm(inv(A)) is computed by a direct method, and the reciprocal of
00043 *> the condition number is computed as
00044 *>              RCOND = 1 / (ANORM * norm(inv(A))).
00045 *> \endverbatim
00046 *
00047 *  Arguments:
00048 *  ==========
00049 *
00050 *> \param[in] N
00051 *> \verbatim
00052 *>          N is INTEGER
00053 *>          The order of the matrix A.  N >= 0.
00054 *> \endverbatim
00055 *>
00056 *> \param[in] D
00057 *> \verbatim
00058 *>          D is DOUBLE PRECISION array, dimension (N)
00059 *>          The n diagonal elements of the diagonal matrix D from the
00060 *>          factorization of A, as computed by DPTTRF.
00061 *> \endverbatim
00062 *>
00063 *> \param[in] E
00064 *> \verbatim
00065 *>          E is DOUBLE PRECISION array, dimension (N-1)
00066 *>          The (n-1) off-diagonal elements of the unit bidiagonal factor
00067 *>          U or L from the factorization of A,  as computed by DPTTRF.
00068 *> \endverbatim
00069 *>
00070 *> \param[in] ANORM
00071 *> \verbatim
00072 *>          ANORM is DOUBLE PRECISION
00073 *>          The 1-norm of the original matrix A.
00074 *> \endverbatim
00075 *>
00076 *> \param[out] RCOND
00077 *> \verbatim
00078 *>          RCOND is DOUBLE PRECISION
00079 *>          The reciprocal of the condition number of the matrix A,
00080 *>          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is the
00081 *>          1-norm of inv(A) computed in this routine.
00082 *> \endverbatim
00083 *>
00084 *> \param[out] WORK
00085 *> \verbatim
00086 *>          WORK is DOUBLE PRECISION array, dimension (N)
00087 *> \endverbatim
00088 *>
00089 *> \param[out] INFO
00090 *> \verbatim
00091 *>          INFO is INTEGER
00092 *>          = 0:  successful exit
00093 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00094 *> \endverbatim
00095 *
00096 *  Authors:
00097 *  ========
00098 *
00099 *> \author Univ. of Tennessee 
00100 *> \author Univ. of California Berkeley 
00101 *> \author Univ. of Colorado Denver 
00102 *> \author NAG Ltd. 
00103 *
00104 *> \date November 2011
00105 *
00106 *> \ingroup doubleOTHERcomputational
00107 *
00108 *> \par Further Details:
00109 *  =====================
00110 *>
00111 *> \verbatim
00112 *>
00113 *>  The method used is described in Nicholas J. Higham, "Efficient
00114 *>  Algorithms for Computing the Condition Number of a Tridiagonal
00115 *>  Matrix", SIAM J. Sci. Stat. Comput., Vol. 7, No. 1, January 1986.
00116 *> \endverbatim
00117 *>
00118 *  =====================================================================
00119       SUBROUTINE DPTCON( N, D, E, ANORM, RCOND, WORK, INFO )
00120 *
00121 *  -- LAPACK computational routine (version 3.4.0) --
00122 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00123 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00124 *     November 2011
00125 *
00126 *     .. Scalar Arguments ..
00127       INTEGER            INFO, N
00128       DOUBLE PRECISION   ANORM, RCOND
00129 *     ..
00130 *     .. Array Arguments ..
00131       DOUBLE PRECISION   D( * ), E( * ), WORK( * )
00132 *     ..
00133 *
00134 *  =====================================================================
00135 *
00136 *     .. Parameters ..
00137       DOUBLE PRECISION   ONE, ZERO
00138       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00139 *     ..
00140 *     .. Local Scalars ..
00141       INTEGER            I, IX
00142       DOUBLE PRECISION   AINVNM
00143 *     ..
00144 *     .. External Functions ..
00145       INTEGER            IDAMAX
00146       EXTERNAL           IDAMAX
00147 *     ..
00148 *     .. External Subroutines ..
00149       EXTERNAL           XERBLA
00150 *     ..
00151 *     .. Intrinsic Functions ..
00152       INTRINSIC          ABS
00153 *     ..
00154 *     .. Executable Statements ..
00155 *
00156 *     Test the input arguments.
00157 *
00158       INFO = 0
00159       IF( N.LT.0 ) THEN
00160          INFO = -1
00161       ELSE IF( ANORM.LT.ZERO ) THEN
00162          INFO = -4
00163       END IF
00164       IF( INFO.NE.0 ) THEN
00165          CALL XERBLA( 'DPTCON', -INFO )
00166          RETURN
00167       END IF
00168 *
00169 *     Quick return if possible
00170 *
00171       RCOND = ZERO
00172       IF( N.EQ.0 ) THEN
00173          RCOND = ONE
00174          RETURN
00175       ELSE IF( ANORM.EQ.ZERO ) THEN
00176          RETURN
00177       END IF
00178 *
00179 *     Check that D(1:N) is positive.
00180 *
00181       DO 10 I = 1, N
00182          IF( D( I ).LE.ZERO )
00183      $      RETURN
00184    10 CONTINUE
00185 *
00186 *     Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
00187 *
00188 *        m(i,j) =  abs(A(i,j)), i = j,
00189 *        m(i,j) = -abs(A(i,j)), i .ne. j,
00190 *
00191 *     and e = [ 1, 1, ..., 1 ]**T.  Note M(A) = M(L)*D*M(L)**T.
00192 *
00193 *     Solve M(L) * x = e.
00194 *
00195       WORK( 1 ) = ONE
00196       DO 20 I = 2, N
00197          WORK( I ) = ONE + WORK( I-1 )*ABS( E( I-1 ) )
00198    20 CONTINUE
00199 *
00200 *     Solve D * M(L)**T * x = b.
00201 *
00202       WORK( N ) = WORK( N ) / D( N )
00203       DO 30 I = N - 1, 1, -1
00204          WORK( I ) = WORK( I ) / D( I ) + WORK( I+1 )*ABS( E( I ) )
00205    30 CONTINUE
00206 *
00207 *     Compute AINVNM = max(x(i)), 1<=i<=n.
00208 *
00209       IX = IDAMAX( N, WORK, 1 )
00210       AINVNM = ABS( WORK( IX ) )
00211 *
00212 *     Compute the reciprocal condition number.
00213 *
00214       IF( AINVNM.NE.ZERO )
00215      $   RCOND = ( ONE / AINVNM ) / ANORM
00216 *
00217       RETURN
00218 *
00219 *     End of DPTCON
00220 *
00221       END
 All Files Functions