1 DOUBLE PRECISION FUNCTION DLA_PORCOND( UPLO, N, A, LDA, AF, LDAF,
2 $ CMODE, C, INFO, WORK,
3 $ 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 UPLO
17 INTEGER N, LDA, LDAF, INFO, CMODE
18 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * ),
19 $ C( * )
20 * ..
21 * .. Array Arguments ..
22 INTEGER IWORK( * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * DLA_PORCOND 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 * UPLO (input) CHARACTER*1
42 * = 'U': Upper triangle of A is stored;
43 * = 'L': Lower triangle of A is stored.
44 *
45 * N (input) INTEGER
46 * The number of linear equations, i.e., the order of the
47 * matrix A. N >= 0.
48 *
49 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
50 * On entry, the N-by-N matrix A.
51 *
52 * LDA (input) INTEGER
53 * The leading dimension of the array A. LDA >= max(1,N).
54 *
55 * AF (input) DOUBLE PRECISION array, dimension (LDAF,N)
56 * The triangular factor U or L from the Cholesky factorization
57 * A = U**T*U or A = L*L**T, as computed by DPOTRF.
58 *
59 * LDAF (input) INTEGER
60 * The leading dimension of the array AF. LDAF >= max(1,N).
61 *
62 * CMODE (input) INTEGER
63 * Determines op2(C) in the formula op(A) * op2(C) as follows:
64 * CMODE = 1 op2(C) = C
65 * CMODE = 0 op2(C) = I
66 * CMODE = -1 op2(C) = inv(C)
67 *
68 * C (input) DOUBLE PRECISION array, dimension (N)
69 * The vector C in the formula op(A) * op2(C).
70 *
71 * INFO (output) INTEGER
72 * = 0: Successful exit.
73 * i > 0: The ith argument is invalid.
74 *
75 * WORK (input) DOUBLE PRECISION array, dimension (3*N).
76 * Workspace.
77 *
78 * IWORK (input) INTEGER array, dimension (N).
79 * Workspace.
80 *
81 * =====================================================================
82 *
83 * .. Local Scalars ..
84 INTEGER KASE, I, J
85 DOUBLE PRECISION AINVNM, TMP
86 LOGICAL UP
87 * ..
88 * .. Array Arguments ..
89 INTEGER ISAVE( 3 )
90 * ..
91 * .. External Functions ..
92 LOGICAL LSAME
93 INTEGER IDAMAX
94 EXTERNAL LSAME, IDAMAX
95 * ..
96 * .. External Subroutines ..
97 EXTERNAL DLACN2, DPOTRS, XERBLA
98 * ..
99 * .. Intrinsic Functions ..
100 INTRINSIC ABS, MAX
101 * ..
102 * .. Executable Statements ..
103 *
104 DLA_PORCOND = 0.0D+0
105 *
106 INFO = 0
107 IF( N.LT.0 ) THEN
108 INFO = -2
109 END IF
110 IF( INFO.NE.0 ) THEN
111 CALL XERBLA( 'DLA_PORCOND', -INFO )
112 RETURN
113 END IF
114
115 IF( N.EQ.0 ) THEN
116 DLA_PORCOND = 1.0D+0
117 RETURN
118 END IF
119 UP = .FALSE.
120 IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
121 *
122 * Compute the equilibration matrix R such that
123 * inv(R)*A*C has unit 1-norm.
124 *
125 IF ( UP ) THEN
126 DO I = 1, N
127 TMP = 0.0D+0
128 IF ( CMODE .EQ. 1 ) THEN
129 DO J = 1, I
130 TMP = TMP + ABS( A( J, I ) * C( J ) )
131 END DO
132 DO J = I+1, N
133 TMP = TMP + ABS( A( I, J ) * C( J ) )
134 END DO
135 ELSE IF ( CMODE .EQ. 0 ) THEN
136 DO J = 1, I
137 TMP = TMP + ABS( A( J, I ) )
138 END DO
139 DO J = I+1, N
140 TMP = TMP + ABS( A( I, J ) )
141 END DO
142 ELSE
143 DO J = 1, I
144 TMP = TMP + ABS( A( J ,I ) / C( J ) )
145 END DO
146 DO J = I+1, N
147 TMP = TMP + ABS( A( I, J ) / C( J ) )
148 END DO
149 END IF
150 WORK( 2*N+I ) = TMP
151 END DO
152 ELSE
153 DO I = 1, N
154 TMP = 0.0D+0
155 IF ( CMODE .EQ. 1 ) THEN
156 DO J = 1, I
157 TMP = TMP + ABS( A( I, J ) * C( J ) )
158 END DO
159 DO J = I+1, N
160 TMP = TMP + ABS( A( J, I ) * C( J ) )
161 END DO
162 ELSE IF ( CMODE .EQ. 0 ) THEN
163 DO J = 1, I
164 TMP = TMP + ABS( A( I, J ) )
165 END DO
166 DO J = I+1, N
167 TMP = TMP + ABS( A( J, I ) )
168 END DO
169 ELSE
170 DO J = 1, I
171 TMP = TMP + ABS( A( I, J ) / C( J ) )
172 END DO
173 DO J = I+1, N
174 TMP = TMP + ABS( A( J, I ) / C( J ) )
175 END DO
176 END IF
177 WORK( 2*N+I ) = TMP
178 END DO
179 ENDIF
180 *
181 * Estimate the norm of inv(op(A)).
182 *
183 AINVNM = 0.0D+0
184
185 KASE = 0
186 10 CONTINUE
187 CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
188 IF( KASE.NE.0 ) THEN
189 IF( KASE.EQ.2 ) THEN
190 *
191 * Multiply by R.
192 *
193 DO I = 1, N
194 WORK( I ) = WORK( I ) * WORK( 2*N+I )
195 END DO
196
197 IF (UP) THEN
198 CALL DPOTRS( 'Upper', N, 1, AF, LDAF, WORK, N, INFO )
199 ELSE
200 CALL DPOTRS( 'Lower', N, 1, AF, LDAF, WORK, N, INFO )
201 ENDIF
202 *
203 * Multiply by inv(C).
204 *
205 IF ( CMODE .EQ. 1 ) THEN
206 DO I = 1, N
207 WORK( I ) = WORK( I ) / C( I )
208 END DO
209 ELSE IF ( CMODE .EQ. -1 ) THEN
210 DO I = 1, N
211 WORK( I ) = WORK( I ) * C( I )
212 END DO
213 END IF
214 ELSE
215 *
216 * Multiply by inv(C**T).
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
228 IF ( UP ) THEN
229 CALL DPOTRS( 'Upper', N, 1, AF, LDAF, WORK, N, INFO )
230 ELSE
231 CALL DPOTRS( 'Lower', N, 1, AF, LDAF, WORK, N, INFO )
232 ENDIF
233 *
234 * Multiply by R.
235 *
236 DO I = 1, N
237 WORK( I ) = WORK( I ) * WORK( 2*N+I )
238 END DO
239 END IF
240 GO TO 10
241 END IF
242 *
243 * Compute the estimate of the reciprocal condition number.
244 *
245 IF( AINVNM .NE. 0.0D+0 )
246 $ DLA_PORCOND = ( 1.0D+0 / AINVNM )
247 *
248 RETURN
249 *
250 END
2 $ CMODE, C, INFO, WORK,
3 $ 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 UPLO
17 INTEGER N, LDA, LDAF, INFO, CMODE
18 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * ),
19 $ C( * )
20 * ..
21 * .. Array Arguments ..
22 INTEGER IWORK( * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * DLA_PORCOND 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 * UPLO (input) CHARACTER*1
42 * = 'U': Upper triangle of A is stored;
43 * = 'L': Lower triangle of A is stored.
44 *
45 * N (input) INTEGER
46 * The number of linear equations, i.e., the order of the
47 * matrix A. N >= 0.
48 *
49 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
50 * On entry, the N-by-N matrix A.
51 *
52 * LDA (input) INTEGER
53 * The leading dimension of the array A. LDA >= max(1,N).
54 *
55 * AF (input) DOUBLE PRECISION array, dimension (LDAF,N)
56 * The triangular factor U or L from the Cholesky factorization
57 * A = U**T*U or A = L*L**T, as computed by DPOTRF.
58 *
59 * LDAF (input) INTEGER
60 * The leading dimension of the array AF. LDAF >= max(1,N).
61 *
62 * CMODE (input) INTEGER
63 * Determines op2(C) in the formula op(A) * op2(C) as follows:
64 * CMODE = 1 op2(C) = C
65 * CMODE = 0 op2(C) = I
66 * CMODE = -1 op2(C) = inv(C)
67 *
68 * C (input) DOUBLE PRECISION array, dimension (N)
69 * The vector C in the formula op(A) * op2(C).
70 *
71 * INFO (output) INTEGER
72 * = 0: Successful exit.
73 * i > 0: The ith argument is invalid.
74 *
75 * WORK (input) DOUBLE PRECISION array, dimension (3*N).
76 * Workspace.
77 *
78 * IWORK (input) INTEGER array, dimension (N).
79 * Workspace.
80 *
81 * =====================================================================
82 *
83 * .. Local Scalars ..
84 INTEGER KASE, I, J
85 DOUBLE PRECISION AINVNM, TMP
86 LOGICAL UP
87 * ..
88 * .. Array Arguments ..
89 INTEGER ISAVE( 3 )
90 * ..
91 * .. External Functions ..
92 LOGICAL LSAME
93 INTEGER IDAMAX
94 EXTERNAL LSAME, IDAMAX
95 * ..
96 * .. External Subroutines ..
97 EXTERNAL DLACN2, DPOTRS, XERBLA
98 * ..
99 * .. Intrinsic Functions ..
100 INTRINSIC ABS, MAX
101 * ..
102 * .. Executable Statements ..
103 *
104 DLA_PORCOND = 0.0D+0
105 *
106 INFO = 0
107 IF( N.LT.0 ) THEN
108 INFO = -2
109 END IF
110 IF( INFO.NE.0 ) THEN
111 CALL XERBLA( 'DLA_PORCOND', -INFO )
112 RETURN
113 END IF
114
115 IF( N.EQ.0 ) THEN
116 DLA_PORCOND = 1.0D+0
117 RETURN
118 END IF
119 UP = .FALSE.
120 IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
121 *
122 * Compute the equilibration matrix R such that
123 * inv(R)*A*C has unit 1-norm.
124 *
125 IF ( UP ) THEN
126 DO I = 1, N
127 TMP = 0.0D+0
128 IF ( CMODE .EQ. 1 ) THEN
129 DO J = 1, I
130 TMP = TMP + ABS( A( J, I ) * C( J ) )
131 END DO
132 DO J = I+1, N
133 TMP = TMP + ABS( A( I, J ) * C( J ) )
134 END DO
135 ELSE IF ( CMODE .EQ. 0 ) THEN
136 DO J = 1, I
137 TMP = TMP + ABS( A( J, I ) )
138 END DO
139 DO J = I+1, N
140 TMP = TMP + ABS( A( I, J ) )
141 END DO
142 ELSE
143 DO J = 1, I
144 TMP = TMP + ABS( A( J ,I ) / C( J ) )
145 END DO
146 DO J = I+1, N
147 TMP = TMP + ABS( A( I, J ) / C( J ) )
148 END DO
149 END IF
150 WORK( 2*N+I ) = TMP
151 END DO
152 ELSE
153 DO I = 1, N
154 TMP = 0.0D+0
155 IF ( CMODE .EQ. 1 ) THEN
156 DO J = 1, I
157 TMP = TMP + ABS( A( I, J ) * C( J ) )
158 END DO
159 DO J = I+1, N
160 TMP = TMP + ABS( A( J, I ) * C( J ) )
161 END DO
162 ELSE IF ( CMODE .EQ. 0 ) THEN
163 DO J = 1, I
164 TMP = TMP + ABS( A( I, J ) )
165 END DO
166 DO J = I+1, N
167 TMP = TMP + ABS( A( J, I ) )
168 END DO
169 ELSE
170 DO J = 1, I
171 TMP = TMP + ABS( A( I, J ) / C( J ) )
172 END DO
173 DO J = I+1, N
174 TMP = TMP + ABS( A( J, I ) / C( J ) )
175 END DO
176 END IF
177 WORK( 2*N+I ) = TMP
178 END DO
179 ENDIF
180 *
181 * Estimate the norm of inv(op(A)).
182 *
183 AINVNM = 0.0D+0
184
185 KASE = 0
186 10 CONTINUE
187 CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
188 IF( KASE.NE.0 ) THEN
189 IF( KASE.EQ.2 ) THEN
190 *
191 * Multiply by R.
192 *
193 DO I = 1, N
194 WORK( I ) = WORK( I ) * WORK( 2*N+I )
195 END DO
196
197 IF (UP) THEN
198 CALL DPOTRS( 'Upper', N, 1, AF, LDAF, WORK, N, INFO )
199 ELSE
200 CALL DPOTRS( 'Lower', N, 1, AF, LDAF, WORK, N, INFO )
201 ENDIF
202 *
203 * Multiply by inv(C).
204 *
205 IF ( CMODE .EQ. 1 ) THEN
206 DO I = 1, N
207 WORK( I ) = WORK( I ) / C( I )
208 END DO
209 ELSE IF ( CMODE .EQ. -1 ) THEN
210 DO I = 1, N
211 WORK( I ) = WORK( I ) * C( I )
212 END DO
213 END IF
214 ELSE
215 *
216 * Multiply by inv(C**T).
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
228 IF ( UP ) THEN
229 CALL DPOTRS( 'Upper', N, 1, AF, LDAF, WORK, N, INFO )
230 ELSE
231 CALL DPOTRS( 'Lower', N, 1, AF, LDAF, WORK, N, INFO )
232 ENDIF
233 *
234 * Multiply by R.
235 *
236 DO I = 1, N
237 WORK( I ) = WORK( I ) * WORK( 2*N+I )
238 END DO
239 END IF
240 GO TO 10
241 END IF
242 *
243 * Compute the estimate of the reciprocal condition number.
244 *
245 IF( AINVNM .NE. 0.0D+0 )
246 $ DLA_PORCOND = ( 1.0D+0 / AINVNM )
247 *
248 RETURN
249 *
250 END