LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cpbcon.f
Go to the documentation of this file.
1*> \brief \b CPBCON
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CPBCON + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cpbcon.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cpbcon.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cpbcon.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CPBCON( UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK,
22* RWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER UPLO
26* INTEGER INFO, KD, LDAB, N
27* REAL ANORM, RCOND
28* ..
29* .. Array Arguments ..
30* REAL RWORK( * )
31* COMPLEX AB( LDAB, * ), WORK( * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> CPBCON estimates the reciprocal of the condition number (in the
41*> 1-norm) of a complex Hermitian positive definite band matrix using
42*> the Cholesky factorization A = U**H*U or A = L*L**H computed by
43*> CPBTRF.
44*>
45*> An estimate is obtained for norm(inv(A)), and the reciprocal of the
46*> condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
47*> \endverbatim
48*
49* Arguments:
50* ==========
51*
52*> \param[in] UPLO
53*> \verbatim
54*> UPLO is CHARACTER*1
55*> = 'U': Upper triangular factor stored in AB;
56*> = 'L': Lower triangular factor stored in AB.
57*> \endverbatim
58*>
59*> \param[in] N
60*> \verbatim
61*> N is INTEGER
62*> The order of the matrix A. N >= 0.
63*> \endverbatim
64*>
65*> \param[in] KD
66*> \verbatim
67*> KD is INTEGER
68*> The number of superdiagonals of the matrix A if UPLO = 'U',
69*> or the number of sub-diagonals if UPLO = 'L'. KD >= 0.
70*> \endverbatim
71*>
72*> \param[in] AB
73*> \verbatim
74*> AB is COMPLEX array, dimension (LDAB,N)
75*> The triangular factor U or L from the Cholesky factorization
76*> A = U**H*U or A = L*L**H of the band matrix A, stored in the
77*> first KD+1 rows of the array. The j-th column of U or L is
78*> stored in the j-th column of the array AB as follows:
79*> if UPLO ='U', AB(kd+1+i-j,j) = U(i,j) for max(1,j-kd)<=i<=j;
80*> if UPLO ='L', AB(1+i-j,j) = L(i,j) for j<=i<=min(n,j+kd).
81*> \endverbatim
82*>
83*> \param[in] LDAB
84*> \verbatim
85*> LDAB is INTEGER
86*> The leading dimension of the array AB. LDAB >= KD+1.
87*> \endverbatim
88*>
89*> \param[in] ANORM
90*> \verbatim
91*> ANORM is REAL
92*> The 1-norm (or infinity-norm) of the Hermitian band matrix A.
93*> \endverbatim
94*>
95*> \param[out] RCOND
96*> \verbatim
97*> RCOND is REAL
98*> The reciprocal of the condition number of the matrix A,
99*> computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
100*> estimate of the 1-norm of inv(A) computed in this routine.
101*> \endverbatim
102*>
103*> \param[out] WORK
104*> \verbatim
105*> WORK is COMPLEX array, dimension (2*N)
106*> \endverbatim
107*>
108*> \param[out] RWORK
109*> \verbatim
110*> RWORK is REAL array, dimension (N)
111*> \endverbatim
112*>
113*> \param[out] INFO
114*> \verbatim
115*> INFO is INTEGER
116*> = 0: successful exit
117*> < 0: if INFO = -i, the i-th argument had an illegal value
118*> \endverbatim
119*
120* Authors:
121* ========
122*
123*> \author Univ. of Tennessee
124*> \author Univ. of California Berkeley
125*> \author Univ. of Colorado Denver
126*> \author NAG Ltd.
127*
128*> \ingroup complexOTHERcomputational
129*
130* =====================================================================
131 SUBROUTINE cpbcon( UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK,
132 $ RWORK, INFO )
133*
134* -- LAPACK computational routine --
135* -- LAPACK is a software package provided by Univ. of Tennessee, --
136* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
137*
138* .. Scalar Arguments ..
139 CHARACTER UPLO
140 INTEGER INFO, KD, LDAB, N
141 REAL ANORM, RCOND
142* ..
143* .. Array Arguments ..
144 REAL RWORK( * )
145 COMPLEX AB( LDAB, * ), WORK( * )
146* ..
147*
148* =====================================================================
149*
150* .. Parameters ..
151 REAL ONE, ZERO
152 parameter( one = 1.0e+0, zero = 0.0e+0 )
153* ..
154* .. Local Scalars ..
155 LOGICAL UPPER
156 CHARACTER NORMIN
157 INTEGER IX, KASE
158 REAL AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
159 COMPLEX ZDUM
160* ..
161* .. Local Arrays ..
162 INTEGER ISAVE( 3 )
163* ..
164* .. External Functions ..
165 LOGICAL LSAME
166 INTEGER ICAMAX
167 REAL SLAMCH
168 EXTERNAL lsame, icamax, slamch
169* ..
170* .. External Subroutines ..
171 EXTERNAL clacn2, clatbs, csrscl, xerbla
172* ..
173* .. Intrinsic Functions ..
174 INTRINSIC abs, aimag, real
175* ..
176* .. Statement Functions ..
177 REAL CABS1
178* ..
179* .. Statement Function definitions ..
180 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
181* ..
182* .. Executable Statements ..
183*
184* Test the input parameters.
185*
186 info = 0
187 upper = lsame( uplo, 'U' )
188 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
189 info = -1
190 ELSE IF( n.LT.0 ) THEN
191 info = -2
192 ELSE IF( kd.LT.0 ) THEN
193 info = -3
194 ELSE IF( ldab.LT.kd+1 ) THEN
195 info = -5
196 ELSE IF( anorm.LT.zero ) THEN
197 info = -6
198 END IF
199 IF( info.NE.0 ) THEN
200 CALL xerbla( 'CPBCON', -info )
201 RETURN
202 END IF
203*
204* Quick return if possible
205*
206 rcond = zero
207 IF( n.EQ.0 ) THEN
208 rcond = one
209 RETURN
210 ELSE IF( anorm.EQ.zero ) THEN
211 RETURN
212 END IF
213*
214 smlnum = slamch( 'Safe minimum' )
215*
216* Estimate the 1-norm of the inverse.
217*
218 kase = 0
219 normin = 'N'
220 10 CONTINUE
221 CALL clacn2( n, work( n+1 ), work, ainvnm, kase, isave )
222 IF( kase.NE.0 ) THEN
223 IF( upper ) THEN
224*
225* Multiply by inv(U**H).
226*
227 CALL clatbs( 'Upper', 'Conjugate transpose', 'Non-unit',
228 $ normin, n, kd, ab, ldab, work, scalel, rwork,
229 $ info )
230 normin = 'Y'
231*
232* Multiply by inv(U).
233*
234 CALL clatbs( 'Upper', 'No transpose', 'Non-unit', normin, n,
235 $ kd, ab, ldab, work, scaleu, rwork, info )
236 ELSE
237*
238* Multiply by inv(L).
239*
240 CALL clatbs( 'Lower', 'No transpose', 'Non-unit', normin, n,
241 $ kd, ab, ldab, work, scalel, rwork, info )
242 normin = 'Y'
243*
244* Multiply by inv(L**H).
245*
246 CALL clatbs( 'Lower', 'Conjugate transpose', 'Non-unit',
247 $ normin, n, kd, ab, ldab, work, scaleu, rwork,
248 $ info )
249 END IF
250*
251* Multiply by 1/SCALE if doing so will not cause overflow.
252*
253 scale = scalel*scaleu
254 IF( scale.NE.one ) THEN
255 ix = icamax( n, work, 1 )
256 IF( scale.LT.cabs1( work( ix ) )*smlnum .OR. scale.EQ.zero )
257 $ GO TO 20
258 CALL csrscl( n, scale, work, 1 )
259 END IF
260 GO TO 10
261 END IF
262*
263* Compute the estimate of the reciprocal condition number.
264*
265 IF( ainvnm.NE.zero )
266 $ rcond = ( one / ainvnm ) / anorm
267*
268 20 CONTINUE
269*
270 RETURN
271*
272* End of CPBCON
273*
274 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine csrscl(N, SA, SX, INCX)
CSRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: csrscl.f:84
subroutine clatbs(UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, SCALE, CNORM, INFO)
CLATBS solves a triangular banded system of equations.
Definition: clatbs.f:243
subroutine clacn2(N, V, X, EST, KASE, ISAVE)
CLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: clacn2.f:133
subroutine cpbcon(UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK, RWORK, INFO)
CPBCON
Definition: cpbcon.f:133