1 DOUBLE PRECISION FUNCTION DLA_GERCOND ( TRANS, N, A, LDA, AF,
2 $ LDAF, IPIV, CMODE, C,
3 $ INFO, WORK, IWORK )
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, LDA, LDAF, INFO, CMODE
18 * ..
19 * .. Array Arguments ..
20 INTEGER IPIV( * ), IWORK( * )
21 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * ),
22 $ C( * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * DLA_GERCOND 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 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
52 * On entry, the N-by-N matrix A.
53 *
54 * LDA (input) INTEGER
55 * The leading dimension of the array A. LDA >= max(1,N).
56 *
57 * AF (input) DOUBLE PRECISION array, dimension (LDAF,N)
58 * The factors L and U from the factorization
59 * A = P*L*U as computed by DGETRF.
60 *
61 * LDAF (input) INTEGER
62 * The leading dimension of the array AF. LDAF >= max(1,N).
63 *
64 * IPIV (input) INTEGER array, dimension (N)
65 * The pivot indices from the factorization A = P*L*U
66 * as computed by DGETRF; row i of the matrix was interchanged
67 * with row IPIV(i).
68 *
69 * CMODE (input) INTEGER
70 * Determines op2(C) in the formula op(A) * op2(C) as follows:
71 * CMODE = 1 op2(C) = C
72 * CMODE = 0 op2(C) = I
73 * CMODE = -1 op2(C) = inv(C)
74 *
75 * C (input) DOUBLE PRECISION array, dimension (N)
76 * The vector C in the formula op(A) * op2(C).
77 *
78 * INFO (output) INTEGER
79 * = 0: Successful exit.
80 * i > 0: The ith argument is invalid.
81 *
82 * WORK (input) DOUBLE PRECISION array, dimension (3*N).
83 * Workspace.
84 *
85 * IWORK (input) INTEGER array, dimension (N).
86 * Workspace.
87 *
88 * =====================================================================
89 *
90 * .. Local Scalars ..
91 LOGICAL NOTRANS
92 INTEGER KASE, I, J
93 DOUBLE PRECISION AINVNM, TMP
94 * ..
95 * .. Local Arrays ..
96 INTEGER ISAVE( 3 )
97 * ..
98 * .. External Functions ..
99 LOGICAL LSAME
100 EXTERNAL LSAME
101 * ..
102 * .. External Subroutines ..
103 EXTERNAL DLACN2, DGETRS, XERBLA
104 * ..
105 * .. Intrinsic Functions ..
106 INTRINSIC ABS, MAX
107 * ..
108 * .. Executable Statements ..
109 *
110 DLA_GERCOND = 0.0D+0
111 *
112 INFO = 0
113 NOTRANS = LSAME( TRANS, 'N' )
114 IF ( .NOT. NOTRANS .AND. .NOT. LSAME(TRANS, 'T')
115 $ .AND. .NOT. LSAME(TRANS, 'C') ) THEN
116 INFO = -1
117 ELSE IF( N.LT.0 ) THEN
118 INFO = -2
119 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
120 INFO = -4
121 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
122 INFO = -6
123 END IF
124 IF( INFO.NE.0 ) THEN
125 CALL XERBLA( 'DLA_GERCOND', -INFO )
126 RETURN
127 END IF
128 IF( N.EQ.0 ) THEN
129 DLA_GERCOND = 1.0D+0
130 RETURN
131 END IF
132 *
133 * Compute the equilibration matrix R such that
134 * inv(R)*A*C has unit 1-norm.
135 *
136 IF (NOTRANS) THEN
137 DO I = 1, N
138 TMP = 0.0D+0
139 IF ( CMODE .EQ. 1 ) THEN
140 DO J = 1, N
141 TMP = TMP + ABS( A( I, J ) * C( J ) )
142 END DO
143 ELSE IF ( CMODE .EQ. 0 ) THEN
144 DO J = 1, N
145 TMP = TMP + ABS( A( I, J ) )
146 END DO
147 ELSE
148 DO J = 1, N
149 TMP = TMP + ABS( A( I, J ) / C( J ) )
150 END DO
151 END IF
152 WORK( 2*N+I ) = TMP
153 END DO
154 ELSE
155 DO I = 1, N
156 TMP = 0.0D+0
157 IF ( CMODE .EQ. 1 ) THEN
158 DO J = 1, N
159 TMP = TMP + ABS( A( J, I ) * C( J ) )
160 END DO
161 ELSE IF ( CMODE .EQ. 0 ) THEN
162 DO J = 1, N
163 TMP = TMP + ABS( A( J, I ) )
164 END DO
165 ELSE
166 DO J = 1, N
167 TMP = TMP + ABS( A( J, I ) / C( J ) )
168 END DO
169 END IF
170 WORK( 2*N+I ) = TMP
171 END DO
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 DLACN2( N, WORK( N+1 ), WORK, IWORK, 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) * WORK(2*N+I)
188 END DO
189
190 IF (NOTRANS) THEN
191 CALL DGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
192 $ WORK, N, INFO )
193 ELSE
194 CALL DGETRS( 'Transpose', N, 1, AF, LDAF, IPIV,
195 $ WORK, N, INFO )
196 END IF
197 *
198 * Multiply by inv(C).
199 *
200 IF ( CMODE .EQ. 1 ) THEN
201 DO I = 1, N
202 WORK( I ) = WORK( I ) / C( I )
203 END DO
204 ELSE IF ( CMODE .EQ. -1 ) THEN
205 DO I = 1, N
206 WORK( I ) = WORK( I ) * C( I )
207 END DO
208 END IF
209 ELSE
210 *
211 * Multiply by inv(C**T).
212 *
213 IF ( CMODE .EQ. 1 ) THEN
214 DO I = 1, N
215 WORK( I ) = WORK( I ) / C( I )
216 END DO
217 ELSE IF ( CMODE .EQ. -1 ) THEN
218 DO I = 1, N
219 WORK( I ) = WORK( I ) * C( I )
220 END DO
221 END IF
222
223 IF (NOTRANS) THEN
224 CALL DGETRS( 'Transpose', N, 1, AF, LDAF, IPIV,
225 $ WORK, N, INFO )
226 ELSE
227 CALL DGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
228 $ WORK, N, INFO )
229 END IF
230 *
231 * Multiply by R.
232 *
233 DO I = 1, N
234 WORK( I ) = WORK( I ) * WORK( 2*N+I )
235 END DO
236 END IF
237 GO TO 10
238 END IF
239 *
240 * Compute the estimate of the reciprocal condition number.
241 *
242 IF( AINVNM .NE. 0.0D+0 )
243 $ DLA_GERCOND = ( 1.0D+0 / AINVNM )
244 *
245 RETURN
246 *
247 END
2 $ LDAF, IPIV, CMODE, C,
3 $ INFO, WORK, IWORK )
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, LDA, LDAF, INFO, CMODE
18 * ..
19 * .. Array Arguments ..
20 INTEGER IPIV( * ), IWORK( * )
21 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * ),
22 $ C( * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * DLA_GERCOND 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 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
52 * On entry, the N-by-N matrix A.
53 *
54 * LDA (input) INTEGER
55 * The leading dimension of the array A. LDA >= max(1,N).
56 *
57 * AF (input) DOUBLE PRECISION array, dimension (LDAF,N)
58 * The factors L and U from the factorization
59 * A = P*L*U as computed by DGETRF.
60 *
61 * LDAF (input) INTEGER
62 * The leading dimension of the array AF. LDAF >= max(1,N).
63 *
64 * IPIV (input) INTEGER array, dimension (N)
65 * The pivot indices from the factorization A = P*L*U
66 * as computed by DGETRF; row i of the matrix was interchanged
67 * with row IPIV(i).
68 *
69 * CMODE (input) INTEGER
70 * Determines op2(C) in the formula op(A) * op2(C) as follows:
71 * CMODE = 1 op2(C) = C
72 * CMODE = 0 op2(C) = I
73 * CMODE = -1 op2(C) = inv(C)
74 *
75 * C (input) DOUBLE PRECISION array, dimension (N)
76 * The vector C in the formula op(A) * op2(C).
77 *
78 * INFO (output) INTEGER
79 * = 0: Successful exit.
80 * i > 0: The ith argument is invalid.
81 *
82 * WORK (input) DOUBLE PRECISION array, dimension (3*N).
83 * Workspace.
84 *
85 * IWORK (input) INTEGER array, dimension (N).
86 * Workspace.
87 *
88 * =====================================================================
89 *
90 * .. Local Scalars ..
91 LOGICAL NOTRANS
92 INTEGER KASE, I, J
93 DOUBLE PRECISION AINVNM, TMP
94 * ..
95 * .. Local Arrays ..
96 INTEGER ISAVE( 3 )
97 * ..
98 * .. External Functions ..
99 LOGICAL LSAME
100 EXTERNAL LSAME
101 * ..
102 * .. External Subroutines ..
103 EXTERNAL DLACN2, DGETRS, XERBLA
104 * ..
105 * .. Intrinsic Functions ..
106 INTRINSIC ABS, MAX
107 * ..
108 * .. Executable Statements ..
109 *
110 DLA_GERCOND = 0.0D+0
111 *
112 INFO = 0
113 NOTRANS = LSAME( TRANS, 'N' )
114 IF ( .NOT. NOTRANS .AND. .NOT. LSAME(TRANS, 'T')
115 $ .AND. .NOT. LSAME(TRANS, 'C') ) THEN
116 INFO = -1
117 ELSE IF( N.LT.0 ) THEN
118 INFO = -2
119 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
120 INFO = -4
121 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
122 INFO = -6
123 END IF
124 IF( INFO.NE.0 ) THEN
125 CALL XERBLA( 'DLA_GERCOND', -INFO )
126 RETURN
127 END IF
128 IF( N.EQ.0 ) THEN
129 DLA_GERCOND = 1.0D+0
130 RETURN
131 END IF
132 *
133 * Compute the equilibration matrix R such that
134 * inv(R)*A*C has unit 1-norm.
135 *
136 IF (NOTRANS) THEN
137 DO I = 1, N
138 TMP = 0.0D+0
139 IF ( CMODE .EQ. 1 ) THEN
140 DO J = 1, N
141 TMP = TMP + ABS( A( I, J ) * C( J ) )
142 END DO
143 ELSE IF ( CMODE .EQ. 0 ) THEN
144 DO J = 1, N
145 TMP = TMP + ABS( A( I, J ) )
146 END DO
147 ELSE
148 DO J = 1, N
149 TMP = TMP + ABS( A( I, J ) / C( J ) )
150 END DO
151 END IF
152 WORK( 2*N+I ) = TMP
153 END DO
154 ELSE
155 DO I = 1, N
156 TMP = 0.0D+0
157 IF ( CMODE .EQ. 1 ) THEN
158 DO J = 1, N
159 TMP = TMP + ABS( A( J, I ) * C( J ) )
160 END DO
161 ELSE IF ( CMODE .EQ. 0 ) THEN
162 DO J = 1, N
163 TMP = TMP + ABS( A( J, I ) )
164 END DO
165 ELSE
166 DO J = 1, N
167 TMP = TMP + ABS( A( J, I ) / C( J ) )
168 END DO
169 END IF
170 WORK( 2*N+I ) = TMP
171 END DO
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 DLACN2( N, WORK( N+1 ), WORK, IWORK, 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) * WORK(2*N+I)
188 END DO
189
190 IF (NOTRANS) THEN
191 CALL DGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
192 $ WORK, N, INFO )
193 ELSE
194 CALL DGETRS( 'Transpose', N, 1, AF, LDAF, IPIV,
195 $ WORK, N, INFO )
196 END IF
197 *
198 * Multiply by inv(C).
199 *
200 IF ( CMODE .EQ. 1 ) THEN
201 DO I = 1, N
202 WORK( I ) = WORK( I ) / C( I )
203 END DO
204 ELSE IF ( CMODE .EQ. -1 ) THEN
205 DO I = 1, N
206 WORK( I ) = WORK( I ) * C( I )
207 END DO
208 END IF
209 ELSE
210 *
211 * Multiply by inv(C**T).
212 *
213 IF ( CMODE .EQ. 1 ) THEN
214 DO I = 1, N
215 WORK( I ) = WORK( I ) / C( I )
216 END DO
217 ELSE IF ( CMODE .EQ. -1 ) THEN
218 DO I = 1, N
219 WORK( I ) = WORK( I ) * C( I )
220 END DO
221 END IF
222
223 IF (NOTRANS) THEN
224 CALL DGETRS( 'Transpose', N, 1, AF, LDAF, IPIV,
225 $ WORK, N, INFO )
226 ELSE
227 CALL DGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
228 $ WORK, N, INFO )
229 END IF
230 *
231 * Multiply by R.
232 *
233 DO I = 1, N
234 WORK( I ) = WORK( I ) * WORK( 2*N+I )
235 END DO
236 END IF
237 GO TO 10
238 END IF
239 *
240 * Compute the estimate of the reciprocal condition number.
241 *
242 IF( AINVNM .NE. 0.0D+0 )
243 $ DLA_GERCOND = ( 1.0D+0 / AINVNM )
244 *
245 RETURN
246 *
247 END