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