1: *> \brief \b DPBCON
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DPBCON + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dpbcon.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dpbcon.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dpbcon.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DPBCON( UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK,
22: * IWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER UPLO
26: * INTEGER INFO, KD, LDAB, N
27: * DOUBLE PRECISION ANORM, RCOND
28: * ..
29: * .. Array Arguments ..
30: * INTEGER IWORK( * )
31: * DOUBLE PRECISION AB( LDAB, * ), WORK( * )
32: * ..
33: *
34: *
35: *> \par Purpose:
36: * =============
37: *>
38: *> \verbatim
39: *>
40: *> DPBCON estimates the reciprocal of the condition number (in the
41: *> 1-norm) of a real symmetric positive definite band matrix using the
42: *> Cholesky factorization A = U**T*U or A = L*L**T computed by DPBTRF.
43: *>
44: *> An estimate is obtained for norm(inv(A)), and the reciprocal of the
45: *> condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
46: *> \endverbatim
47: *
48: * Arguments:
49: * ==========
50: *
51: *> \param[in] UPLO
52: *> \verbatim
53: *> UPLO is CHARACTER*1
54: *> = 'U': Upper triangular factor stored in AB;
55: *> = 'L': Lower triangular factor stored in AB.
56: *> \endverbatim
57: *>
58: *> \param[in] N
59: *> \verbatim
60: *> N is INTEGER
61: *> The order of the matrix A. N >= 0.
62: *> \endverbatim
63: *>
64: *> \param[in] KD
65: *> \verbatim
66: *> KD is INTEGER
67: *> The number of superdiagonals of the matrix A if UPLO = 'U',
68: *> or the number of subdiagonals if UPLO = 'L'. KD >= 0.
69: *> \endverbatim
70: *>
71: *> \param[in] AB
72: *> \verbatim
73: *> AB is DOUBLE PRECISION array, dimension (LDAB,N)
74: *> The triangular factor U or L from the Cholesky factorization
75: *> A = U**T*U or A = L*L**T of the band matrix A, stored in the
76: *> first KD+1 rows of the array. The j-th column of U or L is
77: *> stored in the j-th column of the array AB as follows:
78: *> if UPLO ='U', AB(kd+1+i-j,j) = U(i,j) for max(1,j-kd)<=i<=j;
79: *> if UPLO ='L', AB(1+i-j,j) = L(i,j) for j<=i<=min(n,j+kd).
80: *> \endverbatim
81: *>
82: *> \param[in] LDAB
83: *> \verbatim
84: *> LDAB is INTEGER
85: *> The leading dimension of the array AB. LDAB >= KD+1.
86: *> \endverbatim
87: *>
88: *> \param[in] ANORM
89: *> \verbatim
90: *> ANORM is DOUBLE PRECISION
91: *> The 1-norm (or infinity-norm) of the symmetric band matrix A.
92: *> \endverbatim
93: *>
94: *> \param[out] RCOND
95: *> \verbatim
96: *> RCOND is DOUBLE PRECISION
97: *> The reciprocal of the condition number of the matrix A,
98: *> computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
99: *> estimate of the 1-norm of inv(A) computed in this routine.
100: *> \endverbatim
101: *>
102: *> \param[out] WORK
103: *> \verbatim
104: *> WORK is DOUBLE PRECISION array, dimension (3*N)
105: *> \endverbatim
106: *>
107: *> \param[out] IWORK
108: *> \verbatim
109: *> IWORK is INTEGER array, dimension (N)
110: *> \endverbatim
111: *>
112: *> \param[out] INFO
113: *> \verbatim
114: *> INFO is INTEGER
115: *> = 0: successful exit
116: *> < 0: if INFO = -i, the i-th argument had an illegal value
117: *> \endverbatim
118: *
119: * Authors:
120: * ========
121: *
122: *> \author Univ. of Tennessee
123: *> \author Univ. of California Berkeley
124: *> \author Univ. of Colorado Denver
125: *> \author NAG Ltd.
126: *
127: *> \ingroup doubleOTHERcomputational
128: *
129: * =====================================================================
130: SUBROUTINE DPBCON( UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK,
131: $ IWORK, INFO )
132: *
133: * -- LAPACK computational routine --
134: * -- LAPACK is a software package provided by Univ. of Tennessee, --
135: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
136: *
137: * .. Scalar Arguments ..
138: CHARACTER UPLO
139: INTEGER INFO, KD, LDAB, N
140: DOUBLE PRECISION ANORM, RCOND
141: * ..
142: * .. Array Arguments ..
143: INTEGER IWORK( * )
144: DOUBLE PRECISION AB( LDAB, * ), WORK( * )
145: * ..
146: *
147: * =====================================================================
148: *
149: * .. Parameters ..
150: DOUBLE PRECISION ONE, ZERO
151: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
152: * ..
153: * .. Local Scalars ..
154: LOGICAL UPPER
155: CHARACTER NORMIN
156: INTEGER IX, KASE
157: DOUBLE PRECISION AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
158: * ..
159: * .. Local Arrays ..
160: INTEGER ISAVE( 3 )
161: * ..
162: * .. External Functions ..
163: LOGICAL LSAME
164: INTEGER IDAMAX
165: DOUBLE PRECISION DLAMCH
166: EXTERNAL LSAME, IDAMAX, DLAMCH
167: * ..
168: * .. External Subroutines ..
169: EXTERNAL DLACN2, DLATBS, DRSCL, XERBLA
170: * ..
171: * .. Intrinsic Functions ..
172: INTRINSIC ABS
173: * ..
174: * .. Executable Statements ..
175: *
176: * Test the input parameters.
177: *
178: INFO = 0
179: UPPER = LSAME( UPLO, 'U' )
180: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
181: INFO = -1
182: ELSE IF( N.LT.0 ) THEN
183: INFO = -2
184: ELSE IF( KD.LT.0 ) THEN
185: INFO = -3
186: ELSE IF( LDAB.LT.KD+1 ) THEN
187: INFO = -5
188: ELSE IF( ANORM.LT.ZERO ) THEN
189: INFO = -6
190: END IF
191: IF( INFO.NE.0 ) THEN
192: CALL XERBLA( 'DPBCON', -INFO )
193: RETURN
194: END IF
195: *
196: * Quick return if possible
197: *
198: RCOND = ZERO
199: IF( N.EQ.0 ) THEN
200: RCOND = ONE
201: RETURN
202: ELSE IF( ANORM.EQ.ZERO ) THEN
203: RETURN
204: END IF
205: *
206: SMLNUM = DLAMCH( 'Safe minimum' )
207: *
208: * Estimate the 1-norm of the inverse.
209: *
210: KASE = 0
211: NORMIN = 'N'
212: 10 CONTINUE
213: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
214: IF( KASE.NE.0 ) THEN
215: IF( UPPER ) THEN
216: *
217: * Multiply by inv(U**T).
218: *
219: CALL DLATBS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N,
220: $ KD, AB, LDAB, WORK, SCALEL, WORK( 2*N+1 ),
221: $ INFO )
222: NORMIN = 'Y'
223: *
224: * Multiply by inv(U).
225: *
226: CALL DLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
227: $ KD, AB, LDAB, WORK, SCALEU, WORK( 2*N+1 ),
228: $ INFO )
229: ELSE
230: *
231: * Multiply by inv(L).
232: *
233: CALL DLATBS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
234: $ KD, AB, LDAB, WORK, SCALEL, WORK( 2*N+1 ),
235: $ INFO )
236: NORMIN = 'Y'
237: *
238: * Multiply by inv(L**T).
239: *
240: CALL DLATBS( 'Lower', 'Transpose', 'Non-unit', NORMIN, N,
241: $ KD, AB, LDAB, WORK, SCALEU, WORK( 2*N+1 ),
242: $ INFO )
243: END IF
244: *
245: * Multiply by 1/SCALE if doing so will not cause overflow.
246: *
247: SCALE = SCALEL*SCALEU
248: IF( SCALE.NE.ONE ) THEN
249: IX = IDAMAX( N, WORK, 1 )
250: IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
251: $ GO TO 20
252: CALL DRSCL( N, SCALE, WORK, 1 )
253: END IF
254: GO TO 10
255: END IF
256: *
257: * Compute the estimate of the reciprocal condition number.
258: *
259: IF( AINVNM.NE.ZERO )
260: $ RCOND = ( ONE / AINVNM ) / ANORM
261: *
262: 20 CONTINUE
263: *
264: RETURN
265: *
266: * End of DPBCON
267: *
268: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dpbcon.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack