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