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: *> \date December 2016
128: *
129: *> \ingroup doubleOTHERcomputational
130: *
131: * =====================================================================
132: SUBROUTINE DPBCON( UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK,
133: $ IWORK, INFO )
134: *
135: * -- LAPACK computational routine (version 3.7.0) --
136: * -- LAPACK is a software package provided by Univ. of Tennessee, --
137: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138: * December 2016
139: *
140: * .. Scalar Arguments ..
141: CHARACTER UPLO
142: INTEGER INFO, KD, LDAB, N
143: DOUBLE PRECISION ANORM, RCOND
144: * ..
145: * .. Array Arguments ..
146: INTEGER IWORK( * )
147: DOUBLE PRECISION AB( LDAB, * ), WORK( * )
148: * ..
149: *
150: * =====================================================================
151: *
152: * .. Parameters ..
153: DOUBLE PRECISION ONE, ZERO
154: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
155: * ..
156: * .. Local Scalars ..
157: LOGICAL UPPER
158: CHARACTER NORMIN
159: INTEGER IX, KASE
160: DOUBLE PRECISION AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
161: * ..
162: * .. Local Arrays ..
163: INTEGER ISAVE( 3 )
164: * ..
165: * .. External Functions ..
166: LOGICAL LSAME
167: INTEGER IDAMAX
168: DOUBLE PRECISION DLAMCH
169: EXTERNAL LSAME, IDAMAX, DLAMCH
170: * ..
171: * .. External Subroutines ..
172: EXTERNAL DLACN2, DLATBS, DRSCL, XERBLA
173: * ..
174: * .. Intrinsic Functions ..
175: INTRINSIC ABS
176: * ..
177: * .. Executable Statements ..
178: *
179: * Test the input parameters.
180: *
181: INFO = 0
182: UPPER = LSAME( UPLO, 'U' )
183: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
184: INFO = -1
185: ELSE IF( N.LT.0 ) THEN
186: INFO = -2
187: ELSE IF( KD.LT.0 ) THEN
188: INFO = -3
189: ELSE IF( LDAB.LT.KD+1 ) THEN
190: INFO = -5
191: ELSE IF( ANORM.LT.ZERO ) THEN
192: INFO = -6
193: END IF
194: IF( INFO.NE.0 ) THEN
195: CALL XERBLA( 'DPBCON', -INFO )
196: RETURN
197: END IF
198: *
199: * Quick return if possible
200: *
201: RCOND = ZERO
202: IF( N.EQ.0 ) THEN
203: RCOND = ONE
204: RETURN
205: ELSE IF( ANORM.EQ.ZERO ) THEN
206: RETURN
207: END IF
208: *
209: SMLNUM = DLAMCH( 'Safe minimum' )
210: *
211: * Estimate the 1-norm of the inverse.
212: *
213: KASE = 0
214: NORMIN = 'N'
215: 10 CONTINUE
216: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
217: IF( KASE.NE.0 ) THEN
218: IF( UPPER ) THEN
219: *
220: * Multiply by inv(U**T).
221: *
222: CALL DLATBS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N,
223: $ KD, AB, LDAB, WORK, SCALEL, WORK( 2*N+1 ),
224: $ INFO )
225: NORMIN = 'Y'
226: *
227: * Multiply by inv(U).
228: *
229: CALL DLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
230: $ KD, AB, LDAB, WORK, SCALEU, WORK( 2*N+1 ),
231: $ INFO )
232: ELSE
233: *
234: * Multiply by inv(L).
235: *
236: CALL DLATBS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
237: $ KD, AB, LDAB, WORK, SCALEL, WORK( 2*N+1 ),
238: $ INFO )
239: NORMIN = 'Y'
240: *
241: * Multiply by inv(L**T).
242: *
243: CALL DLATBS( 'Lower', 'Transpose', 'Non-unit', NORMIN, N,
244: $ KD, AB, LDAB, WORK, SCALEU, WORK( 2*N+1 ),
245: $ INFO )
246: END IF
247: *
248: * Multiply by 1/SCALE if doing so will not cause overflow.
249: *
250: SCALE = SCALEL*SCALEU
251: IF( SCALE.NE.ONE ) THEN
252: IX = IDAMAX( N, WORK, 1 )
253: IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
254: $ GO TO 20
255: CALL DRSCL( N, SCALE, WORK, 1 )
256: END IF
257: GO TO 10
258: END IF
259: *
260: * Compute the estimate of the reciprocal condition number.
261: *
262: IF( AINVNM.NE.ZERO )
263: $ RCOND = ( ONE / AINVNM ) / ANORM
264: *
265: 20 CONTINUE
266: *
267: RETURN
268: *
269: * End of DPBCON
270: *
271: 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