1: SUBROUTINE DPTCON( N, D, E, ANORM, RCOND, WORK, INFO )
2: *
3: * -- LAPACK routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: INTEGER INFO, N
10: DOUBLE PRECISION ANORM, RCOND
11: * ..
12: * .. Array Arguments ..
13: DOUBLE PRECISION D( * ), E( * ), WORK( * )
14: * ..
15: *
16: * Purpose
17: * =======
18: *
19: * DPTCON computes the reciprocal of the condition number (in the
20: * 1-norm) of a real symmetric positive definite tridiagonal matrix
21: * using the factorization A = L*D*L**T or A = U**T*D*U computed by
22: * DPTTRF.
23: *
24: * Norm(inv(A)) is computed by a direct method, and the reciprocal of
25: * the condition number is computed as
26: * RCOND = 1 / (ANORM * norm(inv(A))).
27: *
28: * Arguments
29: * =========
30: *
31: * N (input) INTEGER
32: * The order of the matrix A. N >= 0.
33: *
34: * D (input) DOUBLE PRECISION array, dimension (N)
35: * The n diagonal elements of the diagonal matrix D from the
36: * factorization of A, as computed by DPTTRF.
37: *
38: * E (input) DOUBLE PRECISION array, dimension (N-1)
39: * The (n-1) off-diagonal elements of the unit bidiagonal factor
40: * U or L from the factorization of A, as computed by DPTTRF.
41: *
42: * ANORM (input) DOUBLE PRECISION
43: * The 1-norm of the original matrix A.
44: *
45: * RCOND (output) DOUBLE PRECISION
46: * The reciprocal of the condition number of the matrix A,
47: * computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is the
48: * 1-norm of inv(A) computed in this routine.
49: *
50: * WORK (workspace) DOUBLE PRECISION array, dimension (N)
51: *
52: * INFO (output) INTEGER
53: * = 0: successful exit
54: * < 0: if INFO = -i, the i-th argument had an illegal value
55: *
56: * Further Details
57: * ===============
58: *
59: * The method used is described in Nicholas J. Higham, "Efficient
60: * Algorithms for Computing the Condition Number of a Tridiagonal
61: * Matrix", SIAM J. Sci. Stat. Comput., Vol. 7, No. 1, January 1986.
62: *
63: * =====================================================================
64: *
65: * .. Parameters ..
66: DOUBLE PRECISION ONE, ZERO
67: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
68: * ..
69: * .. Local Scalars ..
70: INTEGER I, IX
71: DOUBLE PRECISION AINVNM
72: * ..
73: * .. External Functions ..
74: INTEGER IDAMAX
75: EXTERNAL IDAMAX
76: * ..
77: * .. External Subroutines ..
78: EXTERNAL XERBLA
79: * ..
80: * .. Intrinsic Functions ..
81: INTRINSIC ABS
82: * ..
83: * .. Executable Statements ..
84: *
85: * Test the input arguments.
86: *
87: INFO = 0
88: IF( N.LT.0 ) THEN
89: INFO = -1
90: ELSE IF( ANORM.LT.ZERO ) THEN
91: INFO = -4
92: END IF
93: IF( INFO.NE.0 ) THEN
94: CALL XERBLA( 'DPTCON', -INFO )
95: RETURN
96: END IF
97: *
98: * Quick return if possible
99: *
100: RCOND = ZERO
101: IF( N.EQ.0 ) THEN
102: RCOND = ONE
103: RETURN
104: ELSE IF( ANORM.EQ.ZERO ) THEN
105: RETURN
106: END IF
107: *
108: * Check that D(1:N) is positive.
109: *
110: DO 10 I = 1, N
111: IF( D( I ).LE.ZERO )
112: $ RETURN
113: 10 CONTINUE
114: *
115: * Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
116: *
117: * m(i,j) = abs(A(i,j)), i = j,
118: * m(i,j) = -abs(A(i,j)), i .ne. j,
119: *
120: * and e = [ 1, 1, ..., 1 ]'. Note M(A) = M(L)*D*M(L)'.
121: *
122: * Solve M(L) * x = e.
123: *
124: WORK( 1 ) = ONE
125: DO 20 I = 2, N
126: WORK( I ) = ONE + WORK( I-1 )*ABS( E( I-1 ) )
127: 20 CONTINUE
128: *
129: * Solve D * M(L)' * x = b.
130: *
131: WORK( N ) = WORK( N ) / D( N )
132: DO 30 I = N - 1, 1, -1
133: WORK( I ) = WORK( I ) / D( I ) + WORK( I+1 )*ABS( E( I ) )
134: 30 CONTINUE
135: *
136: * Compute AINVNM = max(x(i)), 1<=i<=n.
137: *
138: IX = IDAMAX( N, WORK, 1 )
139: AINVNM = ABS( WORK( IX ) )
140: *
141: * Compute the reciprocal condition number.
142: *
143: IF( AINVNM.NE.ZERO )
144: $ RCOND = ( ONE / AINVNM ) / ANORM
145: *
146: RETURN
147: *
148: * End of DPTCON
149: *
150: END
CVSweb interface <joel.bertrand@systella.fr>