1 DOUBLE PRECISION FUNCTION ZLA_SYRCOND_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_SYRCOND_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 ZSYTRF.
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 ZSYTRF.
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
80 DOUBLE PRECISION AINVNM, ANORM, TMP
81 INTEGER I, J
82 LOGICAL UP
83 COMPLEX*16 ZDUM
84 * ..
85 * .. Local Arrays ..
86 INTEGER ISAVE( 3 )
87 * ..
88 * .. External Functions ..
89 LOGICAL LSAME
90 EXTERNAL LSAME
91 * ..
92 * .. External Subroutines ..
93 EXTERNAL ZLACN2, ZSYTRS, XERBLA
94 * ..
95 * .. Intrinsic Functions ..
96 INTRINSIC ABS, MAX
97 * ..
98 * .. Statement Functions ..
99 DOUBLE PRECISION CABS1
100 * ..
101 * .. Statement Function Definitions ..
102 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
103 * ..
104 * .. Executable Statements ..
105 *
106 ZLA_SYRCOND_C = 0.0D+0
107 *
108 INFO = 0
109 IF( N.LT.0 ) THEN
110 INFO = -2
111 END IF
112 IF( INFO.NE.0 ) THEN
113 CALL XERBLA( 'ZLA_SYRCOND_C', -INFO )
114 RETURN
115 END IF
116 UP = .FALSE.
117 IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
118 *
119 * Compute norm of op(A)*op2(C).
120 *
121 ANORM = 0.0D+0
122 IF ( UP ) THEN
123 DO I = 1, N
124 TMP = 0.0D+0
125 IF ( CAPPLY ) THEN
126 DO J = 1, I
127 TMP = TMP + CABS1( A( J, I ) ) / C( J )
128 END DO
129 DO J = I+1, N
130 TMP = TMP + CABS1( A( I, J ) ) / C( J )
131 END DO
132 ELSE
133 DO J = 1, I
134 TMP = TMP + CABS1( A( J, I ) )
135 END DO
136 DO J = I+1, N
137 TMP = TMP + CABS1( A( I, J ) )
138 END DO
139 END IF
140 RWORK( I ) = TMP
141 ANORM = MAX( ANORM, TMP )
142 END DO
143 ELSE
144 DO I = 1, N
145 TMP = 0.0D+0
146 IF ( CAPPLY ) THEN
147 DO J = 1, I
148 TMP = TMP + CABS1( A( I, J ) ) / C( J )
149 END DO
150 DO J = I+1, N
151 TMP = TMP + CABS1( A( J, I ) ) / C( J )
152 END DO
153 ELSE
154 DO J = 1, I
155 TMP = TMP + CABS1( A( I, J ) )
156 END DO
157 DO J = I+1, N
158 TMP = TMP + CABS1( A( J, I ) )
159 END DO
160 END IF
161 RWORK( I ) = TMP
162 ANORM = MAX( ANORM, TMP )
163 END DO
164 END IF
165 *
166 * Quick return if possible.
167 *
168 IF( N.EQ.0 ) THEN
169 ZLA_SYRCOND_C = 1.0D+0
170 RETURN
171 ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
172 RETURN
173 END IF
174 *
175 * Estimate the norm of inv(op(A)).
176 *
177 AINVNM = 0.0D+0
178 *
179 KASE = 0
180 10 CONTINUE
181 CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
182 IF( KASE.NE.0 ) THEN
183 IF( KASE.EQ.2 ) THEN
184 *
185 * Multiply by R.
186 *
187 DO I = 1, N
188 WORK( I ) = WORK( I ) * RWORK( I )
189 END DO
190 *
191 IF ( UP ) THEN
192 CALL ZSYTRS( 'U', N, 1, AF, LDAF, IPIV,
193 $ WORK, N, INFO )
194 ELSE
195 CALL ZSYTRS( 'L', N, 1, AF, LDAF, IPIV,
196 $ WORK, N, INFO )
197 ENDIF
198 *
199 * Multiply by inv(C).
200 *
201 IF ( CAPPLY ) THEN
202 DO I = 1, N
203 WORK( I ) = WORK( I ) * C( I )
204 END DO
205 END IF
206 ELSE
207 *
208 * Multiply by inv(C**T).
209 *
210 IF ( CAPPLY ) THEN
211 DO I = 1, N
212 WORK( I ) = WORK( I ) * C( I )
213 END DO
214 END IF
215 *
216 IF ( UP ) THEN
217 CALL ZSYTRS( 'U', N, 1, AF, LDAF, IPIV,
218 $ WORK, N, INFO )
219 ELSE
220 CALL ZSYTRS( 'L', N, 1, AF, LDAF, IPIV,
221 $ WORK, N, INFO )
222 END IF
223 *
224 * Multiply by R.
225 *
226 DO I = 1, N
227 WORK( I ) = WORK( I ) * RWORK( I )
228 END DO
229 END IF
230 GO TO 10
231 END IF
232 *
233 * Compute the estimate of the reciprocal condition number.
234 *
235 IF( AINVNM .NE. 0.0D+0 )
236 $ ZLA_SYRCOND_C = 1.0D+0 / AINVNM
237 *
238 RETURN
239 *
240 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_SYRCOND_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 ZSYTRF.
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 ZSYTRF.
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
80 DOUBLE PRECISION AINVNM, ANORM, TMP
81 INTEGER I, J
82 LOGICAL UP
83 COMPLEX*16 ZDUM
84 * ..
85 * .. Local Arrays ..
86 INTEGER ISAVE( 3 )
87 * ..
88 * .. External Functions ..
89 LOGICAL LSAME
90 EXTERNAL LSAME
91 * ..
92 * .. External Subroutines ..
93 EXTERNAL ZLACN2, ZSYTRS, XERBLA
94 * ..
95 * .. Intrinsic Functions ..
96 INTRINSIC ABS, MAX
97 * ..
98 * .. Statement Functions ..
99 DOUBLE PRECISION CABS1
100 * ..
101 * .. Statement Function Definitions ..
102 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
103 * ..
104 * .. Executable Statements ..
105 *
106 ZLA_SYRCOND_C = 0.0D+0
107 *
108 INFO = 0
109 IF( N.LT.0 ) THEN
110 INFO = -2
111 END IF
112 IF( INFO.NE.0 ) THEN
113 CALL XERBLA( 'ZLA_SYRCOND_C', -INFO )
114 RETURN
115 END IF
116 UP = .FALSE.
117 IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
118 *
119 * Compute norm of op(A)*op2(C).
120 *
121 ANORM = 0.0D+0
122 IF ( UP ) THEN
123 DO I = 1, N
124 TMP = 0.0D+0
125 IF ( CAPPLY ) THEN
126 DO J = 1, I
127 TMP = TMP + CABS1( A( J, I ) ) / C( J )
128 END DO
129 DO J = I+1, N
130 TMP = TMP + CABS1( A( I, J ) ) / C( J )
131 END DO
132 ELSE
133 DO J = 1, I
134 TMP = TMP + CABS1( A( J, I ) )
135 END DO
136 DO J = I+1, N
137 TMP = TMP + CABS1( A( I, J ) )
138 END DO
139 END IF
140 RWORK( I ) = TMP
141 ANORM = MAX( ANORM, TMP )
142 END DO
143 ELSE
144 DO I = 1, N
145 TMP = 0.0D+0
146 IF ( CAPPLY ) THEN
147 DO J = 1, I
148 TMP = TMP + CABS1( A( I, J ) ) / C( J )
149 END DO
150 DO J = I+1, N
151 TMP = TMP + CABS1( A( J, I ) ) / C( J )
152 END DO
153 ELSE
154 DO J = 1, I
155 TMP = TMP + CABS1( A( I, J ) )
156 END DO
157 DO J = I+1, N
158 TMP = TMP + CABS1( A( J, I ) )
159 END DO
160 END IF
161 RWORK( I ) = TMP
162 ANORM = MAX( ANORM, TMP )
163 END DO
164 END IF
165 *
166 * Quick return if possible.
167 *
168 IF( N.EQ.0 ) THEN
169 ZLA_SYRCOND_C = 1.0D+0
170 RETURN
171 ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
172 RETURN
173 END IF
174 *
175 * Estimate the norm of inv(op(A)).
176 *
177 AINVNM = 0.0D+0
178 *
179 KASE = 0
180 10 CONTINUE
181 CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
182 IF( KASE.NE.0 ) THEN
183 IF( KASE.EQ.2 ) THEN
184 *
185 * Multiply by R.
186 *
187 DO I = 1, N
188 WORK( I ) = WORK( I ) * RWORK( I )
189 END DO
190 *
191 IF ( UP ) THEN
192 CALL ZSYTRS( 'U', N, 1, AF, LDAF, IPIV,
193 $ WORK, N, INFO )
194 ELSE
195 CALL ZSYTRS( 'L', N, 1, AF, LDAF, IPIV,
196 $ WORK, N, INFO )
197 ENDIF
198 *
199 * Multiply by inv(C).
200 *
201 IF ( CAPPLY ) THEN
202 DO I = 1, N
203 WORK( I ) = WORK( I ) * C( I )
204 END DO
205 END IF
206 ELSE
207 *
208 * Multiply by inv(C**T).
209 *
210 IF ( CAPPLY ) THEN
211 DO I = 1, N
212 WORK( I ) = WORK( I ) * C( I )
213 END DO
214 END IF
215 *
216 IF ( UP ) THEN
217 CALL ZSYTRS( 'U', N, 1, AF, LDAF, IPIV,
218 $ WORK, N, INFO )
219 ELSE
220 CALL ZSYTRS( 'L', N, 1, AF, LDAF, IPIV,
221 $ WORK, N, INFO )
222 END IF
223 *
224 * Multiply by R.
225 *
226 DO I = 1, N
227 WORK( I ) = WORK( I ) * RWORK( I )
228 END DO
229 END IF
230 GO TO 10
231 END IF
232 *
233 * Compute the estimate of the reciprocal condition number.
234 *
235 IF( AINVNM .NE. 0.0D+0 )
236 $ ZLA_SYRCOND_C = 1.0D+0 / AINVNM
237 *
238 RETURN
239 *
240 END