LAPACK 3.3.0

cptcon.f

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