1 DOUBLE PRECISION FUNCTION ZLA_HERCOND_C( UPLO, N, A, LDA, AF,
2 $ LDAF, IPIV, C, CAPPLY,
3 $ 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 UPLO
17 LOGICAL CAPPLY
18 INTEGER N, LDA, LDAF, INFO
19 * ..
20 * .. Array Arguments ..
21 INTEGER IPIV( * )
22 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), WORK( * )
23 DOUBLE PRECISION C ( * ), RWORK( * )
24 * ..
25 *
26 * Purpose
27 * =======
28 *
29 * ZLA_HERCOND_C computes the infinity norm condition number of
30 * op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
31 *
32 * Arguments
33 * =========
34 *
35 * UPLO (input) CHARACTER*1
36 * = 'U': Upper triangle of A is stored;
37 * = 'L': Lower triangle of A is stored.
38 *
39 * N (input) INTEGER
40 * The number of linear equations, i.e., the order of the
41 * matrix A. N >= 0.
42 *
43 * A (input) COMPLEX*16 array, dimension (LDA,N)
44 * On entry, the N-by-N matrix A
45 *
46 * LDA (input) INTEGER
47 * The leading dimension of the array A. LDA >= max(1,N).
48 *
49 * AF (input) COMPLEX*16 array, dimension (LDAF,N)
50 * The block diagonal matrix D and the multipliers used to
51 * obtain the factor U or L as computed by ZHETRF.
52 *
53 * LDAF (input) INTEGER
54 * The leading dimension of the array AF. LDAF >= max(1,N).
55 *
56 * IPIV (input) INTEGER array, dimension (N)
57 * Details of the interchanges and the block structure of D
58 * as determined by CHETRF.
59 *
60 * C (input) DOUBLE PRECISION array, dimension (N)
61 * The vector C in the formula op(A) * inv(diag(C)).
62 *
63 * CAPPLY (input) LOGICAL
64 * If .TRUE. then access the vector C in the formula above.
65 *
66 * INFO (output) INTEGER
67 * = 0: Successful exit.
68 * i > 0: The ith argument is invalid.
69 *
70 * WORK (input) COMPLEX*16 array, dimension (2*N).
71 * Workspace.
72 *
73 * RWORK (input) DOUBLE PRECISION array, dimension (N).
74 * Workspace.
75 *
76 * =====================================================================
77 *
78 * .. Local Scalars ..
79 INTEGER KASE, I, J
80 DOUBLE PRECISION AINVNM, ANORM, TMP
81 LOGICAL UP
82 COMPLEX*16 ZDUM
83 * ..
84 * .. Local Arrays ..
85 INTEGER ISAVE( 3 )
86 * ..
87 * .. External Functions ..
88 LOGICAL LSAME
89 EXTERNAL LSAME
90 * ..
91 * .. External Subroutines ..
92 EXTERNAL ZLACN2, ZHETRS, XERBLA
93 * ..
94 * .. Intrinsic Functions ..
95 INTRINSIC ABS, MAX
96 * ..
97 * .. Statement Functions ..
98 DOUBLE PRECISION CABS1
99 * ..
100 * .. Statement Function Definitions ..
101 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
102 * ..
103 * .. Executable Statements ..
104 *
105 ZLA_HERCOND_C = 0.0D+0
106 *
107 INFO = 0
108 IF( N.LT.0 ) THEN
109 INFO = -2
110 END IF
111 IF( INFO.NE.0 ) THEN
112 CALL XERBLA( 'ZLA_HERCOND_C', -INFO )
113 RETURN
114 END IF
115 UP = .FALSE.
116 IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
117 *
118 * Compute norm of op(A)*op2(C).
119 *
120 ANORM = 0.0D+0
121 IF ( UP ) THEN
122 DO I = 1, N
123 TMP = 0.0D+0
124 IF ( CAPPLY ) THEN
125 DO J = 1, I
126 TMP = TMP + CABS1( A( J, I ) ) / C( J )
127 END DO
128 DO J = I+1, N
129 TMP = TMP + CABS1( A( I, J ) ) / C( J )
130 END DO
131 ELSE
132 DO J = 1, I
133 TMP = TMP + CABS1( A( J, I ) )
134 END DO
135 DO J = I+1, N
136 TMP = TMP + CABS1( A( I, J ) )
137 END DO
138 END IF
139 RWORK( I ) = TMP
140 ANORM = MAX( ANORM, TMP )
141 END DO
142 ELSE
143 DO I = 1, N
144 TMP = 0.0D+0
145 IF ( CAPPLY ) THEN
146 DO J = 1, I
147 TMP = TMP + CABS1( A( I, J ) ) / C( J )
148 END DO
149 DO J = I+1, N
150 TMP = TMP + CABS1( A( J, I ) ) / C( J )
151 END DO
152 ELSE
153 DO J = 1, I
154 TMP = TMP + CABS1( A( I, J ) )
155 END DO
156 DO J = I+1, N
157 TMP = TMP + CABS1( A( J, I ) )
158 END DO
159 END IF
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_HERCOND_C = 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 ( UP ) THEN
191 CALL ZHETRS( 'U', N, 1, AF, LDAF, IPIV,
192 $ WORK, N, INFO )
193 ELSE
194 CALL ZHETRS( 'L', N, 1, AF, LDAF, IPIV,
195 $ WORK, N, INFO )
196 ENDIF
197 *
198 * Multiply by inv(C).
199 *
200 IF ( CAPPLY ) THEN
201 DO I = 1, N
202 WORK( I ) = WORK( I ) * C( I )
203 END DO
204 END IF
205 ELSE
206 *
207 * Multiply by inv(C**H).
208 *
209 IF ( CAPPLY ) THEN
210 DO I = 1, N
211 WORK( I ) = WORK( I ) * C( I )
212 END DO
213 END IF
214 *
215 IF ( UP ) THEN
216 CALL ZHETRS( 'U', N, 1, AF, LDAF, IPIV,
217 $ WORK, N, INFO )
218 ELSE
219 CALL ZHETRS( 'L', N, 1, AF, LDAF, IPIV,
220 $ WORK, N, INFO )
221 END IF
222 *
223 * Multiply by R.
224 *
225 DO I = 1, N
226 WORK( I ) = WORK( I ) * RWORK( I )
227 END DO
228 END IF
229 GO TO 10
230 END IF
231 *
232 * Compute the estimate of the reciprocal condition number.
233 *
234 IF( AINVNM .NE. 0.0D+0 )
235 $ ZLA_HERCOND_C = 1.0D+0 / AINVNM
236 *
237 RETURN
238 *
239 END
2 $ LDAF, IPIV, C, CAPPLY,
3 $ 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 UPLO
17 LOGICAL CAPPLY
18 INTEGER N, LDA, LDAF, INFO
19 * ..
20 * .. Array Arguments ..
21 INTEGER IPIV( * )
22 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), WORK( * )
23 DOUBLE PRECISION C ( * ), RWORK( * )
24 * ..
25 *
26 * Purpose
27 * =======
28 *
29 * ZLA_HERCOND_C computes the infinity norm condition number of
30 * op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
31 *
32 * Arguments
33 * =========
34 *
35 * UPLO (input) CHARACTER*1
36 * = 'U': Upper triangle of A is stored;
37 * = 'L': Lower triangle of A is stored.
38 *
39 * N (input) INTEGER
40 * The number of linear equations, i.e., the order of the
41 * matrix A. N >= 0.
42 *
43 * A (input) COMPLEX*16 array, dimension (LDA,N)
44 * On entry, the N-by-N matrix A
45 *
46 * LDA (input) INTEGER
47 * The leading dimension of the array A. LDA >= max(1,N).
48 *
49 * AF (input) COMPLEX*16 array, dimension (LDAF,N)
50 * The block diagonal matrix D and the multipliers used to
51 * obtain the factor U or L as computed by ZHETRF.
52 *
53 * LDAF (input) INTEGER
54 * The leading dimension of the array AF. LDAF >= max(1,N).
55 *
56 * IPIV (input) INTEGER array, dimension (N)
57 * Details of the interchanges and the block structure of D
58 * as determined by CHETRF.
59 *
60 * C (input) DOUBLE PRECISION array, dimension (N)
61 * The vector C in the formula op(A) * inv(diag(C)).
62 *
63 * CAPPLY (input) LOGICAL
64 * If .TRUE. then access the vector C in the formula above.
65 *
66 * INFO (output) INTEGER
67 * = 0: Successful exit.
68 * i > 0: The ith argument is invalid.
69 *
70 * WORK (input) COMPLEX*16 array, dimension (2*N).
71 * Workspace.
72 *
73 * RWORK (input) DOUBLE PRECISION array, dimension (N).
74 * Workspace.
75 *
76 * =====================================================================
77 *
78 * .. Local Scalars ..
79 INTEGER KASE, I, J
80 DOUBLE PRECISION AINVNM, ANORM, TMP
81 LOGICAL UP
82 COMPLEX*16 ZDUM
83 * ..
84 * .. Local Arrays ..
85 INTEGER ISAVE( 3 )
86 * ..
87 * .. External Functions ..
88 LOGICAL LSAME
89 EXTERNAL LSAME
90 * ..
91 * .. External Subroutines ..
92 EXTERNAL ZLACN2, ZHETRS, XERBLA
93 * ..
94 * .. Intrinsic Functions ..
95 INTRINSIC ABS, MAX
96 * ..
97 * .. Statement Functions ..
98 DOUBLE PRECISION CABS1
99 * ..
100 * .. Statement Function Definitions ..
101 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
102 * ..
103 * .. Executable Statements ..
104 *
105 ZLA_HERCOND_C = 0.0D+0
106 *
107 INFO = 0
108 IF( N.LT.0 ) THEN
109 INFO = -2
110 END IF
111 IF( INFO.NE.0 ) THEN
112 CALL XERBLA( 'ZLA_HERCOND_C', -INFO )
113 RETURN
114 END IF
115 UP = .FALSE.
116 IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
117 *
118 * Compute norm of op(A)*op2(C).
119 *
120 ANORM = 0.0D+0
121 IF ( UP ) THEN
122 DO I = 1, N
123 TMP = 0.0D+0
124 IF ( CAPPLY ) THEN
125 DO J = 1, I
126 TMP = TMP + CABS1( A( J, I ) ) / C( J )
127 END DO
128 DO J = I+1, N
129 TMP = TMP + CABS1( A( I, J ) ) / C( J )
130 END DO
131 ELSE
132 DO J = 1, I
133 TMP = TMP + CABS1( A( J, I ) )
134 END DO
135 DO J = I+1, N
136 TMP = TMP + CABS1( A( I, J ) )
137 END DO
138 END IF
139 RWORK( I ) = TMP
140 ANORM = MAX( ANORM, TMP )
141 END DO
142 ELSE
143 DO I = 1, N
144 TMP = 0.0D+0
145 IF ( CAPPLY ) THEN
146 DO J = 1, I
147 TMP = TMP + CABS1( A( I, J ) ) / C( J )
148 END DO
149 DO J = I+1, N
150 TMP = TMP + CABS1( A( J, I ) ) / C( J )
151 END DO
152 ELSE
153 DO J = 1, I
154 TMP = TMP + CABS1( A( I, J ) )
155 END DO
156 DO J = I+1, N
157 TMP = TMP + CABS1( A( J, I ) )
158 END DO
159 END IF
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_HERCOND_C = 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 ( UP ) THEN
191 CALL ZHETRS( 'U', N, 1, AF, LDAF, IPIV,
192 $ WORK, N, INFO )
193 ELSE
194 CALL ZHETRS( 'L', N, 1, AF, LDAF, IPIV,
195 $ WORK, N, INFO )
196 ENDIF
197 *
198 * Multiply by inv(C).
199 *
200 IF ( CAPPLY ) THEN
201 DO I = 1, N
202 WORK( I ) = WORK( I ) * C( I )
203 END DO
204 END IF
205 ELSE
206 *
207 * Multiply by inv(C**H).
208 *
209 IF ( CAPPLY ) THEN
210 DO I = 1, N
211 WORK( I ) = WORK( I ) * C( I )
212 END DO
213 END IF
214 *
215 IF ( UP ) THEN
216 CALL ZHETRS( 'U', N, 1, AF, LDAF, IPIV,
217 $ WORK, N, INFO )
218 ELSE
219 CALL ZHETRS( 'L', N, 1, AF, LDAF, IPIV,
220 $ WORK, N, INFO )
221 END IF
222 *
223 * Multiply by R.
224 *
225 DO I = 1, N
226 WORK( I ) = WORK( I ) * RWORK( I )
227 END DO
228 END IF
229 GO TO 10
230 END IF
231 *
232 * Compute the estimate of the reciprocal condition number.
233 *
234 IF( AINVNM .NE. 0.0D+0 )
235 $ ZLA_HERCOND_C = 1.0D+0 / AINVNM
236 *
237 RETURN
238 *
239 END