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