1 DOUBLE PRECISION FUNCTION ZLA_GBRCOND_C( TRANS, N, KL, KU, AB,
2 $ LDAB, AFB, LDAFB, IPIV,
3 $ C, CAPPLY, INFO, WORK,
4 $ RWORK )
5 *
6 * -- LAPACK routine (version 3.2.1) --
7 * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
8 * -- Jason Riedy of Univ. of California Berkeley. --
9 * -- April 2009 --
10 *
11 * -- LAPACK is a software package provided by Univ. of Tennessee, --
12 * -- Univ. of California Berkeley and NAG Ltd. --
13 *
14 IMPLICIT NONE
15 * ..
16 * .. Scalar Arguments ..
17 CHARACTER TRANS
18 LOGICAL CAPPLY
19 INTEGER N, KL, KU, KD, KE, LDAB, LDAFB, INFO
20 * ..
21 * .. Array Arguments ..
22 INTEGER IPIV( * )
23 COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), WORK( * )
24 DOUBLE PRECISION C( * ), RWORK( * )
25 *
26 *
27 * Purpose
28 * =======
29 *
30 * ZLA_GBRCOND_C Computes the infinity norm condition number of
31 * op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
32 *
33 * Arguments
34 * =========
35 *
36 * TRANS (input) CHARACTER*1
37 * Specifies the form of the system of equations:
38 * = 'N': A * X = B (No transpose)
39 * = 'T': A**T * X = B (Transpose)
40 * = 'C': A**H * X = B (Conjugate Transpose = Transpose)
41 *
42 * N (input) INTEGER
43 * The number of linear equations, i.e., the order of the
44 * matrix A. N >= 0.
45 *
46 * KL (input) INTEGER
47 * The number of subdiagonals within the band of A. KL >= 0.
48 *
49 * KU (input) INTEGER
50 * The number of superdiagonals within the band of A. KU >= 0.
51 *
52 * AB (input) COMPLEX*16 array, dimension (LDAB,N)
53 * On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
54 * The j-th column of A is stored in the j-th column of the
55 * array AB as follows:
56 * AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
57 *
58 * LDAB (input) INTEGER
59 * The leading dimension of the array AB. LDAB >= KL+KU+1.
60 *
61 * AFB (input) COMPLEX*16 array, dimension (LDAFB,N)
62 * Details of the LU factorization of the band matrix A, as
63 * computed by ZGBTRF. U is stored as an upper triangular
64 * band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1,
65 * and the multipliers used during the factorization are stored
66 * in rows KL+KU+2 to 2*KL+KU+1.
67 *
68 * LDAFB (input) INTEGER
69 * The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1.
70 *
71 * IPIV (input) INTEGER array, dimension (N)
72 * The pivot indices from the factorization A = P*L*U
73 * as computed by ZGBTRF; row i of the matrix was interchanged
74 * with row IPIV(i).
75 *
76 * C (input) DOUBLE PRECISION array, dimension (N)
77 * The vector C in the formula op(A) * inv(diag(C)).
78 *
79 * CAPPLY (input) LOGICAL
80 * If .TRUE. then access the vector C in the formula above.
81 *
82 * INFO (output) INTEGER
83 * = 0: Successful exit.
84 * i > 0: The ith argument is invalid.
85 *
86 * WORK (input) COMPLEX*16 array, dimension (2*N).
87 * Workspace.
88 *
89 * RWORK (input) DOUBLE PRECISION array, dimension (N).
90 * Workspace.
91 *
92 * =====================================================================
93 *
94 * .. Local Scalars ..
95 LOGICAL NOTRANS
96 INTEGER KASE, I, J
97 DOUBLE PRECISION AINVNM, ANORM, TMP
98 COMPLEX*16 ZDUM
99 * ..
100 * .. Local Arrays ..
101 INTEGER ISAVE( 3 )
102 * ..
103 * .. External Functions ..
104 LOGICAL LSAME
105 EXTERNAL LSAME
106 * ..
107 * .. External Subroutines ..
108 EXTERNAL ZLACN2, ZGBTRS, XERBLA
109 * ..
110 * .. Intrinsic Functions ..
111 INTRINSIC ABS, MAX
112 * ..
113 * .. Statement Functions ..
114 DOUBLE PRECISION CABS1
115 * ..
116 * .. Statement Function Definitions ..
117 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
118 * ..
119 * .. Executable Statements ..
120 ZLA_GBRCOND_C = 0.0D+0
121 *
122 INFO = 0
123 NOTRANS = LSAME( TRANS, 'N' )
124 IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT.
125 $ LSAME( TRANS, 'C' ) ) THEN
126 INFO = -1
127 ELSE IF( N.LT.0 ) THEN
128 INFO = -2
129 ELSE IF( KL.LT.0 .OR. KL.GT.N-1 ) THEN
130 INFO = -3
131 ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN
132 INFO = -4
133 ELSE IF( LDAB.LT.KL+KU+1 ) THEN
134 INFO = -6
135 ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN
136 INFO = -8
137 END IF
138 IF( INFO.NE.0 ) THEN
139 CALL XERBLA( 'ZLA_GBRCOND_C', -INFO )
140 RETURN
141 END IF
142 *
143 * Compute norm of op(A)*op2(C).
144 *
145 ANORM = 0.0D+0
146 KD = KU + 1
147 KE = KL + 1
148 IF ( NOTRANS ) THEN
149 DO I = 1, N
150 TMP = 0.0D+0
151 IF ( CAPPLY ) THEN
152 DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
153 TMP = TMP + CABS1( AB( KD+I-J, J ) ) / C( J )
154 END DO
155 ELSE
156 DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
157 TMP = TMP + CABS1( AB( KD+I-J, J ) )
158 END DO
159 END IF
160 RWORK( I ) = TMP
161 ANORM = MAX( ANORM, TMP )
162 END DO
163 ELSE
164 DO I = 1, N
165 TMP = 0.0D+0
166 IF ( CAPPLY ) THEN
167 DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
168 TMP = TMP + CABS1( AB( KE-I+J, I ) ) / C( J )
169 END DO
170 ELSE
171 DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
172 TMP = TMP + CABS1( AB( KE-I+J, I ) )
173 END DO
174 END IF
175 RWORK( I ) = TMP
176 ANORM = MAX( ANORM, TMP )
177 END DO
178 END IF
179 *
180 * Quick return if possible.
181 *
182 IF( N.EQ.0 ) THEN
183 ZLA_GBRCOND_C = 1.0D+0
184 RETURN
185 ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
186 RETURN
187 END IF
188 *
189 * Estimate the norm of inv(op(A)).
190 *
191 AINVNM = 0.0D+0
192 *
193 KASE = 0
194 10 CONTINUE
195 CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
196 IF( KASE.NE.0 ) THEN
197 IF( KASE.EQ.2 ) THEN
198 *
199 * Multiply by R.
200 *
201 DO I = 1, N
202 WORK( I ) = WORK( I ) * RWORK( I )
203 END DO
204 *
205 IF ( NOTRANS ) THEN
206 CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB,
207 $ IPIV, WORK, N, INFO )
208 ELSE
209 CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB,
210 $ LDAFB, IPIV, WORK, N, INFO )
211 ENDIF
212 *
213 * Multiply by inv(C).
214 *
215 IF ( CAPPLY ) THEN
216 DO I = 1, N
217 WORK( I ) = WORK( I ) * C( I )
218 END DO
219 END IF
220 ELSE
221 *
222 * Multiply by inv(C**H).
223 *
224 IF ( CAPPLY ) THEN
225 DO I = 1, N
226 WORK( I ) = WORK( I ) * C( I )
227 END DO
228 END IF
229 *
230 IF ( NOTRANS ) THEN
231 CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB,
232 $ LDAFB, IPIV, WORK, N, INFO )
233 ELSE
234 CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB,
235 $ IPIV, WORK, N, INFO )
236 END IF
237 *
238 * Multiply by R.
239 *
240 DO I = 1, N
241 WORK( I ) = WORK( I ) * RWORK( I )
242 END DO
243 END IF
244 GO TO 10
245 END IF
246 *
247 * Compute the estimate of the reciprocal condition number.
248 *
249 IF( AINVNM .NE. 0.0D+0 )
250 $ ZLA_GBRCOND_C = 1.0D+0 / AINVNM
251 *
252 RETURN
253 *
254 END
2 $ LDAB, AFB, LDAFB, IPIV,
3 $ C, CAPPLY, INFO, WORK,
4 $ RWORK )
5 *
6 * -- LAPACK routine (version 3.2.1) --
7 * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
8 * -- Jason Riedy of Univ. of California Berkeley. --
9 * -- April 2009 --
10 *
11 * -- LAPACK is a software package provided by Univ. of Tennessee, --
12 * -- Univ. of California Berkeley and NAG Ltd. --
13 *
14 IMPLICIT NONE
15 * ..
16 * .. Scalar Arguments ..
17 CHARACTER TRANS
18 LOGICAL CAPPLY
19 INTEGER N, KL, KU, KD, KE, LDAB, LDAFB, INFO
20 * ..
21 * .. Array Arguments ..
22 INTEGER IPIV( * )
23 COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), WORK( * )
24 DOUBLE PRECISION C( * ), RWORK( * )
25 *
26 *
27 * Purpose
28 * =======
29 *
30 * ZLA_GBRCOND_C Computes the infinity norm condition number of
31 * op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
32 *
33 * Arguments
34 * =========
35 *
36 * TRANS (input) CHARACTER*1
37 * Specifies the form of the system of equations:
38 * = 'N': A * X = B (No transpose)
39 * = 'T': A**T * X = B (Transpose)
40 * = 'C': A**H * X = B (Conjugate Transpose = Transpose)
41 *
42 * N (input) INTEGER
43 * The number of linear equations, i.e., the order of the
44 * matrix A. N >= 0.
45 *
46 * KL (input) INTEGER
47 * The number of subdiagonals within the band of A. KL >= 0.
48 *
49 * KU (input) INTEGER
50 * The number of superdiagonals within the band of A. KU >= 0.
51 *
52 * AB (input) COMPLEX*16 array, dimension (LDAB,N)
53 * On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
54 * The j-th column of A is stored in the j-th column of the
55 * array AB as follows:
56 * AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
57 *
58 * LDAB (input) INTEGER
59 * The leading dimension of the array AB. LDAB >= KL+KU+1.
60 *
61 * AFB (input) COMPLEX*16 array, dimension (LDAFB,N)
62 * Details of the LU factorization of the band matrix A, as
63 * computed by ZGBTRF. U is stored as an upper triangular
64 * band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1,
65 * and the multipliers used during the factorization are stored
66 * in rows KL+KU+2 to 2*KL+KU+1.
67 *
68 * LDAFB (input) INTEGER
69 * The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1.
70 *
71 * IPIV (input) INTEGER array, dimension (N)
72 * The pivot indices from the factorization A = P*L*U
73 * as computed by ZGBTRF; row i of the matrix was interchanged
74 * with row IPIV(i).
75 *
76 * C (input) DOUBLE PRECISION array, dimension (N)
77 * The vector C in the formula op(A) * inv(diag(C)).
78 *
79 * CAPPLY (input) LOGICAL
80 * If .TRUE. then access the vector C in the formula above.
81 *
82 * INFO (output) INTEGER
83 * = 0: Successful exit.
84 * i > 0: The ith argument is invalid.
85 *
86 * WORK (input) COMPLEX*16 array, dimension (2*N).
87 * Workspace.
88 *
89 * RWORK (input) DOUBLE PRECISION array, dimension (N).
90 * Workspace.
91 *
92 * =====================================================================
93 *
94 * .. Local Scalars ..
95 LOGICAL NOTRANS
96 INTEGER KASE, I, J
97 DOUBLE PRECISION AINVNM, ANORM, TMP
98 COMPLEX*16 ZDUM
99 * ..
100 * .. Local Arrays ..
101 INTEGER ISAVE( 3 )
102 * ..
103 * .. External Functions ..
104 LOGICAL LSAME
105 EXTERNAL LSAME
106 * ..
107 * .. External Subroutines ..
108 EXTERNAL ZLACN2, ZGBTRS, XERBLA
109 * ..
110 * .. Intrinsic Functions ..
111 INTRINSIC ABS, MAX
112 * ..
113 * .. Statement Functions ..
114 DOUBLE PRECISION CABS1
115 * ..
116 * .. Statement Function Definitions ..
117 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
118 * ..
119 * .. Executable Statements ..
120 ZLA_GBRCOND_C = 0.0D+0
121 *
122 INFO = 0
123 NOTRANS = LSAME( TRANS, 'N' )
124 IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT.
125 $ LSAME( TRANS, 'C' ) ) THEN
126 INFO = -1
127 ELSE IF( N.LT.0 ) THEN
128 INFO = -2
129 ELSE IF( KL.LT.0 .OR. KL.GT.N-1 ) THEN
130 INFO = -3
131 ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN
132 INFO = -4
133 ELSE IF( LDAB.LT.KL+KU+1 ) THEN
134 INFO = -6
135 ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN
136 INFO = -8
137 END IF
138 IF( INFO.NE.0 ) THEN
139 CALL XERBLA( 'ZLA_GBRCOND_C', -INFO )
140 RETURN
141 END IF
142 *
143 * Compute norm of op(A)*op2(C).
144 *
145 ANORM = 0.0D+0
146 KD = KU + 1
147 KE = KL + 1
148 IF ( NOTRANS ) THEN
149 DO I = 1, N
150 TMP = 0.0D+0
151 IF ( CAPPLY ) THEN
152 DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
153 TMP = TMP + CABS1( AB( KD+I-J, J ) ) / C( J )
154 END DO
155 ELSE
156 DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
157 TMP = TMP + CABS1( AB( KD+I-J, J ) )
158 END DO
159 END IF
160 RWORK( I ) = TMP
161 ANORM = MAX( ANORM, TMP )
162 END DO
163 ELSE
164 DO I = 1, N
165 TMP = 0.0D+0
166 IF ( CAPPLY ) THEN
167 DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
168 TMP = TMP + CABS1( AB( KE-I+J, I ) ) / C( J )
169 END DO
170 ELSE
171 DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
172 TMP = TMP + CABS1( AB( KE-I+J, I ) )
173 END DO
174 END IF
175 RWORK( I ) = TMP
176 ANORM = MAX( ANORM, TMP )
177 END DO
178 END IF
179 *
180 * Quick return if possible.
181 *
182 IF( N.EQ.0 ) THEN
183 ZLA_GBRCOND_C = 1.0D+0
184 RETURN
185 ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
186 RETURN
187 END IF
188 *
189 * Estimate the norm of inv(op(A)).
190 *
191 AINVNM = 0.0D+0
192 *
193 KASE = 0
194 10 CONTINUE
195 CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
196 IF( KASE.NE.0 ) THEN
197 IF( KASE.EQ.2 ) THEN
198 *
199 * Multiply by R.
200 *
201 DO I = 1, N
202 WORK( I ) = WORK( I ) * RWORK( I )
203 END DO
204 *
205 IF ( NOTRANS ) THEN
206 CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB,
207 $ IPIV, WORK, N, INFO )
208 ELSE
209 CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB,
210 $ LDAFB, IPIV, WORK, N, INFO )
211 ENDIF
212 *
213 * Multiply by inv(C).
214 *
215 IF ( CAPPLY ) THEN
216 DO I = 1, N
217 WORK( I ) = WORK( I ) * C( I )
218 END DO
219 END IF
220 ELSE
221 *
222 * Multiply by inv(C**H).
223 *
224 IF ( CAPPLY ) THEN
225 DO I = 1, N
226 WORK( I ) = WORK( I ) * C( I )
227 END DO
228 END IF
229 *
230 IF ( NOTRANS ) THEN
231 CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB,
232 $ LDAFB, IPIV, WORK, N, INFO )
233 ELSE
234 CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB,
235 $ IPIV, WORK, N, INFO )
236 END IF
237 *
238 * Multiply by R.
239 *
240 DO I = 1, N
241 WORK( I ) = WORK( I ) * RWORK( I )
242 END DO
243 END IF
244 GO TO 10
245 END IF
246 *
247 * Compute the estimate of the reciprocal condition number.
248 *
249 IF( AINVNM .NE. 0.0D+0 )
250 $ ZLA_GBRCOND_C = 1.0D+0 / AINVNM
251 *
252 RETURN
253 *
254 END