1: *> \brief \b ZPTCON
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZPTCON + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zptcon.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zptcon.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zptcon.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZPTCON( N, D, E, ANORM, RCOND, RWORK, INFO )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER INFO, N
25: * DOUBLE PRECISION ANORM, RCOND
26: * ..
27: * .. Array Arguments ..
28: * DOUBLE PRECISION D( * ), RWORK( * )
29: * COMPLEX*16 E( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> ZPTCON computes the reciprocal of the condition number (in the
39: *> 1-norm) of a complex Hermitian positive definite tridiagonal matrix
40: *> using the factorization A = L*D*L**H or A = U**H*D*U computed by
41: *> ZPTTRF.
42: *>
43: *> Norm(inv(A)) is computed by a direct method, and the reciprocal of
44: *> the condition number is computed as
45: *> RCOND = 1 / (ANORM * norm(inv(A))).
46: *> \endverbatim
47: *
48: * Arguments:
49: * ==========
50: *
51: *> \param[in] N
52: *> \verbatim
53: *> N is INTEGER
54: *> The order of the matrix A. N >= 0.
55: *> \endverbatim
56: *>
57: *> \param[in] D
58: *> \verbatim
59: *> D is DOUBLE PRECISION array, dimension (N)
60: *> The n diagonal elements of the diagonal matrix D from the
61: *> factorization of A, as computed by ZPTTRF.
62: *> \endverbatim
63: *>
64: *> \param[in] E
65: *> \verbatim
66: *> E is COMPLEX*16 array, dimension (N-1)
67: *> The (n-1) off-diagonal elements of the unit bidiagonal factor
68: *> U or L from the factorization of A, as computed by ZPTTRF.
69: *> \endverbatim
70: *>
71: *> \param[in] ANORM
72: *> \verbatim
73: *> ANORM is DOUBLE PRECISION
74: *> The 1-norm of the original matrix A.
75: *> \endverbatim
76: *>
77: *> \param[out] RCOND
78: *> \verbatim
79: *> RCOND is DOUBLE PRECISION
80: *> The reciprocal of the condition number of the matrix A,
81: *> computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is the
82: *> 1-norm of inv(A) computed in this routine.
83: *> \endverbatim
84: *>
85: *> \param[out] RWORK
86: *> \verbatim
87: *> RWORK is DOUBLE PRECISION array, dimension (N)
88: *> \endverbatim
89: *>
90: *> \param[out] INFO
91: *> \verbatim
92: *> INFO is INTEGER
93: *> = 0: successful exit
94: *> < 0: if INFO = -i, the i-th argument had an illegal value
95: *> \endverbatim
96: *
97: * Authors:
98: * ========
99: *
100: *> \author Univ. of Tennessee
101: *> \author Univ. of California Berkeley
102: *> \author Univ. of Colorado Denver
103: *> \author NAG Ltd.
104: *
105: *> \ingroup complex16PTcomputational
106: *
107: *> \par Further Details:
108: * =====================
109: *>
110: *> \verbatim
111: *>
112: *> The method used is described in Nicholas J. Higham, "Efficient
113: *> Algorithms for Computing the Condition Number of a Tridiagonal
114: *> Matrix", SIAM J. Sci. Stat. Comput., Vol. 7, No. 1, January 1986.
115: *> \endverbatim
116: *>
117: * =====================================================================
118: SUBROUTINE ZPTCON( N, D, E, ANORM, RCOND, RWORK, INFO )
119: *
120: * -- LAPACK computational routine --
121: * -- LAPACK is a software package provided by Univ. of Tennessee, --
122: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
123: *
124: * .. Scalar Arguments ..
125: INTEGER INFO, N
126: DOUBLE PRECISION ANORM, RCOND
127: * ..
128: * .. Array Arguments ..
129: DOUBLE PRECISION D( * ), RWORK( * )
130: COMPLEX*16 E( * )
131: * ..
132: *
133: * =====================================================================
134: *
135: * .. Parameters ..
136: DOUBLE PRECISION ONE, ZERO
137: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
138: * ..
139: * .. Local Scalars ..
140: INTEGER I, IX
141: DOUBLE PRECISION AINVNM
142: * ..
143: * .. External Functions ..
144: INTEGER IDAMAX
145: EXTERNAL IDAMAX
146: * ..
147: * .. External Subroutines ..
148: EXTERNAL XERBLA
149: * ..
150: * .. Intrinsic Functions ..
151: INTRINSIC ABS
152: * ..
153: * .. Executable Statements ..
154: *
155: * Test the input arguments.
156: *
157: INFO = 0
158: IF( N.LT.0 ) THEN
159: INFO = -1
160: ELSE IF( ANORM.LT.ZERO ) THEN
161: INFO = -4
162: END IF
163: IF( INFO.NE.0 ) THEN
164: CALL XERBLA( 'ZPTCON', -INFO )
165: RETURN
166: END IF
167: *
168: * Quick return if possible
169: *
170: RCOND = ZERO
171: IF( N.EQ.0 ) THEN
172: RCOND = ONE
173: RETURN
174: ELSE IF( ANORM.EQ.ZERO ) THEN
175: RETURN
176: END IF
177: *
178: * Check that D(1:N) is positive.
179: *
180: DO 10 I = 1, N
181: IF( D( I ).LE.ZERO )
182: $ RETURN
183: 10 CONTINUE
184: *
185: * Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
186: *
187: * m(i,j) = abs(A(i,j)), i = j,
188: * m(i,j) = -abs(A(i,j)), i .ne. j,
189: *
190: * and e = [ 1, 1, ..., 1 ]**T. Note M(A) = M(L)*D*M(L)**H.
191: *
192: * Solve M(L) * x = e.
193: *
194: RWORK( 1 ) = ONE
195: DO 20 I = 2, N
196: RWORK( I ) = ONE + RWORK( I-1 )*ABS( E( I-1 ) )
197: 20 CONTINUE
198: *
199: * Solve D * M(L)**H * x = b.
200: *
201: RWORK( N ) = RWORK( N ) / D( N )
202: DO 30 I = N - 1, 1, -1
203: RWORK( I ) = RWORK( I ) / D( I ) + RWORK( I+1 )*ABS( E( I ) )
204: 30 CONTINUE
205: *
206: * Compute AINVNM = max(x(i)), 1<=i<=n.
207: *
208: IX = IDAMAX( N, RWORK, 1 )
209: AINVNM = ABS( RWORK( IX ) )
210: *
211: * Compute the reciprocal condition number.
212: *
213: IF( AINVNM.NE.ZERO )
214: $ RCOND = ( ONE / AINVNM ) / ANORM
215: *
216: RETURN
217: *
218: * End of ZPTCON
219: *
220: END
CVSweb interface <joel.bertrand@systella.fr>