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