1: SUBROUTINE ZPBCON( UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK,
2: $ RWORK, INFO )
3: *
4: * -- LAPACK routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH.
10: *
11: * .. Scalar Arguments ..
12: CHARACTER UPLO
13: INTEGER INFO, KD, LDAB, N
14: DOUBLE PRECISION ANORM, RCOND
15: * ..
16: * .. Array Arguments ..
17: DOUBLE PRECISION RWORK( * )
18: COMPLEX*16 AB( LDAB, * ), WORK( * )
19: * ..
20: *
21: * Purpose
22: * =======
23: *
24: * ZPBCON estimates the reciprocal of the condition number (in the
25: * 1-norm) of a complex Hermitian positive definite band matrix using
26: * the Cholesky factorization A = U**H*U or A = L*L**H computed by
27: * ZPBTRF.
28: *
29: * An estimate is obtained for norm(inv(A)), and the reciprocal of the
30: * condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
31: *
32: * Arguments
33: * =========
34: *
35: * UPLO (input) CHARACTER*1
36: * = 'U': Upper triangular factor stored in AB;
37: * = 'L': Lower triangular factor stored in AB.
38: *
39: * N (input) INTEGER
40: * The order of the matrix A. N >= 0.
41: *
42: * KD (input) INTEGER
43: * The number of superdiagonals of the matrix A if UPLO = 'U',
44: * or the number of sub-diagonals if UPLO = 'L'. KD >= 0.
45: *
46: * AB (input) COMPLEX*16 array, dimension (LDAB,N)
47: * The triangular factor U or L from the Cholesky factorization
48: * A = U**H*U or A = L*L**H of the band matrix A, stored in the
49: * first KD+1 rows of the array. The j-th column of U or L is
50: * stored in the j-th column of the array AB as follows:
51: * if UPLO ='U', AB(kd+1+i-j,j) = U(i,j) for max(1,j-kd)<=i<=j;
52: * if UPLO ='L', AB(1+i-j,j) = L(i,j) for j<=i<=min(n,j+kd).
53: *
54: * LDAB (input) INTEGER
55: * The leading dimension of the array AB. LDAB >= KD+1.
56: *
57: * ANORM (input) DOUBLE PRECISION
58: * The 1-norm (or infinity-norm) of the Hermitian band matrix A.
59: *
60: * RCOND (output) DOUBLE PRECISION
61: * The reciprocal of the condition number of the matrix A,
62: * computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
63: * estimate of the 1-norm of inv(A) computed in this routine.
64: *
65: * WORK (workspace) COMPLEX*16 array, dimension (2*N)
66: *
67: * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
68: *
69: * INFO (output) INTEGER
70: * = 0: successful exit
71: * < 0: if INFO = -i, the i-th argument had an illegal value
72: *
73: * =====================================================================
74: *
75: * .. Parameters ..
76: DOUBLE PRECISION ONE, ZERO
77: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
78: * ..
79: * .. Local Scalars ..
80: LOGICAL UPPER
81: CHARACTER NORMIN
82: INTEGER IX, KASE
83: DOUBLE PRECISION AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
84: COMPLEX*16 ZDUM
85: * ..
86: * .. Local Arrays ..
87: INTEGER ISAVE( 3 )
88: * ..
89: * .. External Functions ..
90: LOGICAL LSAME
91: INTEGER IZAMAX
92: DOUBLE PRECISION DLAMCH
93: EXTERNAL LSAME, IZAMAX, DLAMCH
94: * ..
95: * .. External Subroutines ..
96: EXTERNAL XERBLA, ZDRSCL, ZLACN2, ZLATBS
97: * ..
98: * .. Intrinsic Functions ..
99: INTRINSIC ABS, DBLE, DIMAG
100: * ..
101: * .. Statement Functions ..
102: DOUBLE PRECISION CABS1
103: * ..
104: * .. Statement Function definitions ..
105: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
106: * ..
107: * .. Executable Statements ..
108: *
109: * Test the input parameters.
110: *
111: INFO = 0
112: UPPER = LSAME( UPLO, 'U' )
113: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
114: INFO = -1
115: ELSE IF( N.LT.0 ) THEN
116: INFO = -2
117: ELSE IF( KD.LT.0 ) THEN
118: INFO = -3
119: ELSE IF( LDAB.LT.KD+1 ) THEN
120: INFO = -5
121: ELSE IF( ANORM.LT.ZERO ) THEN
122: INFO = -6
123: END IF
124: IF( INFO.NE.0 ) THEN
125: CALL XERBLA( 'ZPBCON', -INFO )
126: RETURN
127: END IF
128: *
129: * Quick return if possible
130: *
131: RCOND = ZERO
132: IF( N.EQ.0 ) THEN
133: RCOND = ONE
134: RETURN
135: ELSE IF( ANORM.EQ.ZERO ) THEN
136: RETURN
137: END IF
138: *
139: SMLNUM = DLAMCH( 'Safe minimum' )
140: *
141: * Estimate the 1-norm of the inverse.
142: *
143: KASE = 0
144: NORMIN = 'N'
145: 10 CONTINUE
146: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
147: IF( KASE.NE.0 ) THEN
148: IF( UPPER ) THEN
149: *
150: * Multiply by inv(U').
151: *
152: CALL ZLATBS( 'Upper', 'Conjugate transpose', 'Non-unit',
153: $ NORMIN, N, KD, AB, LDAB, WORK, SCALEL, RWORK,
154: $ INFO )
155: NORMIN = 'Y'
156: *
157: * Multiply by inv(U).
158: *
159: CALL ZLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
160: $ KD, AB, LDAB, WORK, SCALEU, RWORK, INFO )
161: ELSE
162: *
163: * Multiply by inv(L).
164: *
165: CALL ZLATBS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
166: $ KD, AB, LDAB, WORK, SCALEL, RWORK, INFO )
167: NORMIN = 'Y'
168: *
169: * Multiply by inv(L').
170: *
171: CALL ZLATBS( 'Lower', 'Conjugate transpose', 'Non-unit',
172: $ NORMIN, N, KD, AB, LDAB, WORK, SCALEU, RWORK,
173: $ INFO )
174: END IF
175: *
176: * Multiply by 1/SCALE if doing so will not cause overflow.
177: *
178: SCALE = SCALEL*SCALEU
179: IF( SCALE.NE.ONE ) THEN
180: IX = IZAMAX( N, WORK, 1 )
181: IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
182: $ GO TO 20
183: CALL ZDRSCL( N, SCALE, WORK, 1 )
184: END IF
185: GO TO 10
186: END IF
187: *
188: * Compute the estimate of the reciprocal condition number.
189: *
190: IF( AINVNM.NE.ZERO )
191: $ RCOND = ( ONE / AINVNM ) / ANORM
192: *
193: 20 CONTINUE
194: *
195: RETURN
196: *
197: * End of ZPBCON
198: *
199: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zpbcon.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack