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: *> \date November 2011
106: *
107: *> \ingroup complex16OTHERcomputational
108: *
109: *> \par Further Details:
110: * =====================
111: *>
112: *> \verbatim
113: *>
114: *> The method used is described in Nicholas J. Higham, "Efficient
115: *> Algorithms for Computing the Condition Number of a Tridiagonal
116: *> Matrix", SIAM J. Sci. Stat. Comput., Vol. 7, No. 1, January 1986.
117: *> \endverbatim
118: *>
119: * =====================================================================
120: SUBROUTINE ZPTCON( N, D, E, ANORM, RCOND, RWORK, INFO )
121: *
122: * -- LAPACK computational routine (version 3.4.0) --
123: * -- LAPACK is a software package provided by Univ. of Tennessee, --
124: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
125: * November 2011
126: *
127: * .. Scalar Arguments ..
128: INTEGER INFO, N
129: DOUBLE PRECISION ANORM, RCOND
130: * ..
131: * .. Array Arguments ..
132: DOUBLE PRECISION D( * ), RWORK( * )
133: COMPLEX*16 E( * )
134: * ..
135: *
136: * =====================================================================
137: *
138: * .. Parameters ..
139: DOUBLE PRECISION ONE, ZERO
140: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
141: * ..
142: * .. Local Scalars ..
143: INTEGER I, IX
144: DOUBLE PRECISION AINVNM
145: * ..
146: * .. External Functions ..
147: INTEGER IDAMAX
148: EXTERNAL IDAMAX
149: * ..
150: * .. External Subroutines ..
151: EXTERNAL XERBLA
152: * ..
153: * .. Intrinsic Functions ..
154: INTRINSIC ABS
155: * ..
156: * .. Executable Statements ..
157: *
158: * Test the input arguments.
159: *
160: INFO = 0
161: IF( N.LT.0 ) THEN
162: INFO = -1
163: ELSE IF( ANORM.LT.ZERO ) THEN
164: INFO = -4
165: END IF
166: IF( INFO.NE.0 ) THEN
167: CALL XERBLA( 'ZPTCON', -INFO )
168: RETURN
169: END IF
170: *
171: * Quick return if possible
172: *
173: RCOND = ZERO
174: IF( N.EQ.0 ) THEN
175: RCOND = ONE
176: RETURN
177: ELSE IF( ANORM.EQ.ZERO ) THEN
178: RETURN
179: END IF
180: *
181: * Check that D(1:N) is positive.
182: *
183: DO 10 I = 1, N
184: IF( D( I ).LE.ZERO )
185: $ RETURN
186: 10 CONTINUE
187: *
188: * Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
189: *
190: * m(i,j) = abs(A(i,j)), i = j,
191: * m(i,j) = -abs(A(i,j)), i .ne. j,
192: *
193: * and e = [ 1, 1, ..., 1 ]**T. Note M(A) = M(L)*D*M(L)**H.
194: *
195: * Solve M(L) * x = e.
196: *
197: RWORK( 1 ) = ONE
198: DO 20 I = 2, N
199: RWORK( I ) = ONE + RWORK( I-1 )*ABS( E( I-1 ) )
200: 20 CONTINUE
201: *
202: * Solve D * M(L)**H * x = b.
203: *
204: RWORK( N ) = RWORK( N ) / D( N )
205: DO 30 I = N - 1, 1, -1
206: RWORK( I ) = RWORK( I ) / D( I ) + RWORK( I+1 )*ABS( E( I ) )
207: 30 CONTINUE
208: *
209: * Compute AINVNM = max(x(i)), 1<=i<=n.
210: *
211: IX = IDAMAX( N, RWORK, 1 )
212: AINVNM = ABS( RWORK( IX ) )
213: *
214: * Compute the reciprocal condition number.
215: *
216: IF( AINVNM.NE.ZERO )
217: $ RCOND = ( ONE / AINVNM ) / ANORM
218: *
219: RETURN
220: *
221: * End of ZPTCON
222: *
223: END
CVSweb interface <joel.bertrand@systella.fr>