*> \brief \b DPTCON * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download DPTCON + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * SUBROUTINE DPTCON( N, D, E, ANORM, RCOND, WORK, INFO ) * * .. Scalar Arguments .. * INTEGER INFO, N * DOUBLE PRECISION ANORM, RCOND * .. * .. Array Arguments .. * DOUBLE PRECISION D( * ), E( * ), WORK( * ) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> DPTCON computes the reciprocal of the condition number (in the *> 1-norm) of a real symmetric positive definite tridiagonal matrix *> using the factorization A = L*D*L**T or A = U**T*D*U computed by *> DPTTRF. *> *> Norm(inv(A)) is computed by a direct method, and the reciprocal of *> the condition number is computed as *> RCOND = 1 / (ANORM * norm(inv(A))). *> \endverbatim * * Arguments: * ========== * *> \param[in] N *> \verbatim *> N is INTEGER *> The order of the matrix A. N >= 0. *> \endverbatim *> *> \param[in] D *> \verbatim *> D is DOUBLE PRECISION array, dimension (N) *> The n diagonal elements of the diagonal matrix D from the *> factorization of A, as computed by DPTTRF. *> \endverbatim *> *> \param[in] E *> \verbatim *> E is DOUBLE PRECISION array, dimension (N-1) *> The (n-1) off-diagonal elements of the unit bidiagonal factor *> U or L from the factorization of A, as computed by DPTTRF. *> \endverbatim *> *> \param[in] ANORM *> \verbatim *> ANORM is DOUBLE PRECISION *> The 1-norm of the original matrix A. *> \endverbatim *> *> \param[out] RCOND *> \verbatim *> RCOND is DOUBLE PRECISION *> The reciprocal of the condition number of the matrix A, *> computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is the *> 1-norm of inv(A) computed in this routine. *> \endverbatim *> *> \param[out] WORK *> \verbatim *> WORK is DOUBLE PRECISION array, dimension (N) *> \endverbatim *> *> \param[out] INFO *> \verbatim *> INFO is INTEGER *> = 0: successful exit *> < 0: if INFO = -i, the i-th argument had an illegal value *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \ingroup doublePTcomputational * *> \par Further Details: * ===================== *> *> \verbatim *> *> The method used is described in Nicholas J. Higham, "Efficient *> Algorithms for Computing the Condition Number of a Tridiagonal *> Matrix", SIAM J. Sci. Stat. Comput., Vol. 7, No. 1, January 1986. *> \endverbatim *> * ===================================================================== SUBROUTINE DPTCON( N, D, E, ANORM, RCOND, WORK, INFO ) * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * * .. Scalar Arguments .. INTEGER INFO, N DOUBLE PRECISION ANORM, RCOND * .. * .. Array Arguments .. DOUBLE PRECISION D( * ), E( * ), WORK( * ) * .. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. INTEGER I, IX DOUBLE PRECISION AINVNM * .. * .. External Functions .. INTEGER IDAMAX EXTERNAL IDAMAX * .. * .. External Subroutines .. EXTERNAL XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS * .. * .. Executable Statements .. * * Test the input arguments. * INFO = 0 IF( N.LT.0 ) THEN INFO = -1 ELSE IF( ANORM.LT.ZERO ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DPTCON', -INFO ) RETURN END IF * * Quick return if possible * RCOND = ZERO IF( N.EQ.0 ) THEN RCOND = ONE RETURN ELSE IF( ANORM.EQ.ZERO ) THEN RETURN END IF * * Check that D(1:N) is positive. * DO 10 I = 1, N IF( D( I ).LE.ZERO ) \$ RETURN 10 CONTINUE * * Solve M(A) * x = e, where M(A) = (m(i,j)) is given by * * m(i,j) = abs(A(i,j)), i = j, * m(i,j) = -abs(A(i,j)), i .ne. j, * * and e = [ 1, 1, ..., 1 ]**T. Note M(A) = M(L)*D*M(L)**T. * * Solve M(L) * x = e. * WORK( 1 ) = ONE DO 20 I = 2, N WORK( I ) = ONE + WORK( I-1 )*ABS( E( I-1 ) ) 20 CONTINUE * * Solve D * M(L)**T * x = b. * WORK( N ) = WORK( N ) / D( N ) DO 30 I = N - 1, 1, -1 WORK( I ) = WORK( I ) / D( I ) + WORK( I+1 )*ABS( E( I ) ) 30 CONTINUE * * Compute AINVNM = max(x(i)), 1<=i<=n. * IX = IDAMAX( N, WORK, 1 ) AINVNM = ABS( WORK( IX ) ) * * Compute the reciprocal condition number. * IF( AINVNM.NE.ZERO ) \$ RCOND = ( ONE / AINVNM ) / ANORM * RETURN * * End of DPTCON * END