00001 SUBROUTINE SPTCON( N, D, E, ANORM, RCOND, WORK, 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( * ), E( * ), WORK( * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * SPTCON 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 * SPTTRF. 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) REAL array, dimension (N) 00035 * The n diagonal elements of the diagonal matrix D from the 00036 * factorization of A, as computed by SPTTRF. 00037 * 00038 * E (input) REAL 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 SPTTRF. 00041 * 00042 * ANORM (input) REAL 00043 * The 1-norm of the original matrix A. 00044 * 00045 * RCOND (output) REAL 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) REAL 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 REAL ONE, ZERO 00067 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00068 * .. 00069 * .. Local Scalars .. 00070 INTEGER I, IX 00071 REAL AINVNM 00072 * .. 00073 * .. External Functions .. 00074 INTEGER ISAMAX 00075 EXTERNAL ISAMAX 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( 'SPTCON', -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 ]'. Note M(A) = M(L)*D*M(L)'. 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)' * 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 = ISAMAX( 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 SPTCON 00149 * 00150 END
1.7.2