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