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