LAPACK 3.3.1
Linear Algebra PACKage

dptcon.f

Go to the documentation of this file.
00001       SUBROUTINE DPTCON( N, D, E, ANORM, RCOND, WORK, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.3.1) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *  -- April 2011                                                      --
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            INFO, N
00010       DOUBLE PRECISION   ANORM, RCOND
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   D( * ), E( * ), WORK( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DPTCON computes the reciprocal of the condition number (in the
00020 *  1-norm) of a real symmetric positive definite tridiagonal matrix
00021 *  using the factorization A = L*D*L**T or A = U**T*D*U computed by
00022 *  DPTTRF.
00023 *
00024 *  Norm(inv(A)) is computed by a direct method, and the reciprocal of
00025 *  the condition number is computed as
00026 *               RCOND = 1 / (ANORM * norm(inv(A))).
00027 *
00028 *  Arguments
00029 *  =========
00030 *
00031 *  N       (input) INTEGER
00032 *          The order of the matrix A.  N >= 0.
00033 *
00034 *  D       (input) DOUBLE PRECISION array, dimension (N)
00035 *          The n diagonal elements of the diagonal matrix D from the
00036 *          factorization of A, as computed by DPTTRF.
00037 *
00038 *  E       (input) DOUBLE PRECISION array, dimension (N-1)
00039 *          The (n-1) off-diagonal elements of the unit bidiagonal factor
00040 *          U or L from the factorization of A,  as computed by DPTTRF.
00041 *
00042 *  ANORM   (input) DOUBLE PRECISION
00043 *          The 1-norm of the original matrix A.
00044 *
00045 *  RCOND   (output) DOUBLE PRECISION
00046 *          The reciprocal of the condition number of the matrix A,
00047 *          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is the
00048 *          1-norm of inv(A) computed in this routine.
00049 *
00050 *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
00051 *
00052 *  INFO    (output) INTEGER
00053 *          = 0:  successful exit
00054 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00055 *
00056 *  Further Details
00057 *  ===============
00058 *
00059 *  The method used is described in Nicholas J. Higham, "Efficient
00060 *  Algorithms for Computing the Condition Number of a Tridiagonal
00061 *  Matrix", SIAM J. Sci. Stat. Comput., Vol. 7, No. 1, January 1986.
00062 *
00063 *  =====================================================================
00064 *
00065 *     .. Parameters ..
00066       DOUBLE PRECISION   ONE, ZERO
00067       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00068 *     ..
00069 *     .. Local Scalars ..
00070       INTEGER            I, IX
00071       DOUBLE PRECISION   AINVNM
00072 *     ..
00073 *     .. External Functions ..
00074       INTEGER            IDAMAX
00075       EXTERNAL           IDAMAX
00076 *     ..
00077 *     .. External Subroutines ..
00078       EXTERNAL           XERBLA
00079 *     ..
00080 *     .. Intrinsic Functions ..
00081       INTRINSIC          ABS
00082 *     ..
00083 *     .. Executable Statements ..
00084 *
00085 *     Test the input arguments.
00086 *
00087       INFO = 0
00088       IF( N.LT.0 ) THEN
00089          INFO = -1
00090       ELSE IF( ANORM.LT.ZERO ) THEN
00091          INFO = -4
00092       END IF
00093       IF( INFO.NE.0 ) THEN
00094          CALL XERBLA( 'DPTCON', -INFO )
00095          RETURN
00096       END IF
00097 *
00098 *     Quick return if possible
00099 *
00100       RCOND = ZERO
00101       IF( N.EQ.0 ) THEN
00102          RCOND = ONE
00103          RETURN
00104       ELSE IF( ANORM.EQ.ZERO ) THEN
00105          RETURN
00106       END IF
00107 *
00108 *     Check that D(1:N) is positive.
00109 *
00110       DO 10 I = 1, N
00111          IF( D( I ).LE.ZERO )
00112      $      RETURN
00113    10 CONTINUE
00114 *
00115 *     Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
00116 *
00117 *        m(i,j) =  abs(A(i,j)), i = j,
00118 *        m(i,j) = -abs(A(i,j)), i .ne. j,
00119 *
00120 *     and e = [ 1, 1, ..., 1 ]**T.  Note M(A) = M(L)*D*M(L)**T.
00121 *
00122 *     Solve M(L) * x = e.
00123 *
00124       WORK( 1 ) = ONE
00125       DO 20 I = 2, N
00126          WORK( I ) = ONE + WORK( I-1 )*ABS( E( I-1 ) )
00127    20 CONTINUE
00128 *
00129 *     Solve D * M(L)**T * x = b.
00130 *
00131       WORK( N ) = WORK( N ) / D( N )
00132       DO 30 I = N - 1, 1, -1
00133          WORK( I ) = WORK( I ) / D( I ) + WORK( I+1 )*ABS( E( I ) )
00134    30 CONTINUE
00135 *
00136 *     Compute AINVNM = max(x(i)), 1<=i<=n.
00137 *
00138       IX = IDAMAX( N, WORK, 1 )
00139       AINVNM = ABS( WORK( IX ) )
00140 *
00141 *     Compute the reciprocal condition number.
00142 *
00143       IF( AINVNM.NE.ZERO )
00144      $   RCOND = ( ONE / AINVNM ) / ANORM
00145 *
00146       RETURN
00147 *
00148 *     End of DPTCON
00149 *
00150       END
 All Files Functions