1 DOUBLE PRECISION FUNCTION ZLA_GERCOND_C( TRANS, 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 Aguments ..
16 CHARACTER TRANS
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_GERCOND_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 * TRANS (input) CHARACTER*1
36 * Specifies the form of the system of equations:
37 * = 'N': A * X = B (No transpose)
38 * = 'T': A**T * X = B (Transpose)
39 * = 'C': A**H * X = B (Conjugate Transpose = Transpose)
40 *
41 * N (input) INTEGER
42 * The number of linear equations, i.e., the order of the
43 * matrix A. N >= 0.
44 *
45 * A (input) COMPLEX*16 array, dimension (LDA,N)
46 * On entry, the N-by-N matrix A
47 *
48 * LDA (input) INTEGER
49 * The leading dimension of the array A. LDA >= max(1,N).
50 *
51 * AF (input) COMPLEX*16 array, dimension (LDAF,N)
52 * The factors L and U from the factorization
53 * A = P*L*U as computed by ZGETRF.
54 *
55 * LDAF (input) INTEGER
56 * The leading dimension of the array AF. LDAF >= max(1,N).
57 *
58 * IPIV (input) INTEGER array, dimension (N)
59 * The pivot indices from the factorization A = P*L*U
60 * as computed by ZGETRF; row i of the matrix was interchanged
61 * with row IPIV(i).
62 *
63 * C (input) DOUBLE PRECISION array, dimension (N)
64 * The vector C in the formula op(A) * inv(diag(C)).
65 *
66 * CAPPLY (input) LOGICAL
67 * If .TRUE. then access the vector C in the formula above.
68 *
69 * INFO (output) INTEGER
70 * = 0: Successful exit.
71 * i > 0: The ith argument is invalid.
72 *
73 * WORK (input) COMPLEX*16 array, dimension (2*N).
74 * Workspace.
75 *
76 * RWORK (input) DOUBLE PRECISION array, dimension (N).
77 * Workspace.
78 *
79 * =====================================================================
80 *
81 * .. Local Scalars ..
82 LOGICAL NOTRANS
83 INTEGER KASE, I, J
84 DOUBLE PRECISION AINVNM, ANORM, TMP
85 COMPLEX*16 ZDUM
86 * ..
87 * .. Local Arrays ..
88 INTEGER ISAVE( 3 )
89 * ..
90 * .. External Functions ..
91 LOGICAL LSAME
92 EXTERNAL LSAME
93 * ..
94 * .. External Subroutines ..
95 EXTERNAL ZLACN2, ZGETRS, XERBLA
96 * ..
97 * .. Intrinsic Functions ..
98 INTRINSIC ABS, MAX, REAL, DIMAG
99 * ..
100 * .. Statement Functions ..
101 DOUBLE PRECISION CABS1
102 * ..
103 * .. Statement Function Definitions ..
104 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
105 * ..
106 * .. Executable Statements ..
107 ZLA_GERCOND_C = 0.0D+0
108 *
109 INFO = 0
110 NOTRANS = LSAME( TRANS, 'N' )
111 IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT.
112 $ LSAME( TRANS, 'C' ) ) THEN
113 ELSE IF( N.LT.0 ) THEN
114 INFO = -2
115 END IF
116 IF( INFO.NE.0 ) THEN
117 CALL XERBLA( 'ZLA_GERCOND_C', -INFO )
118 RETURN
119 END IF
120 *
121 * Compute norm of op(A)*op2(C).
122 *
123 ANORM = 0.0D+0
124 IF ( NOTRANS ) THEN
125 DO I = 1, N
126 TMP = 0.0D+0
127 IF ( CAPPLY ) THEN
128 DO J = 1, N
129 TMP = TMP + CABS1( A( I, J ) ) / C( J )
130 END DO
131 ELSE
132 DO J = 1, N
133 TMP = TMP + CABS1( A( I, J ) )
134 END DO
135 END IF
136 RWORK( I ) = TMP
137 ANORM = MAX( ANORM, TMP )
138 END DO
139 ELSE
140 DO I = 1, N
141 TMP = 0.0D+0
142 IF ( CAPPLY ) THEN
143 DO J = 1, N
144 TMP = TMP + CABS1( A( J, I ) ) / C( J )
145 END DO
146 ELSE
147 DO J = 1, N
148 TMP = TMP + CABS1( A( J, I ) )
149 END DO
150 END IF
151 RWORK( I ) = TMP
152 ANORM = MAX( ANORM, TMP )
153 END DO
154 END IF
155 *
156 * Quick return if possible.
157 *
158 IF( N.EQ.0 ) THEN
159 ZLA_GERCOND_C = 1.0D+0
160 RETURN
161 ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
162 RETURN
163 END IF
164 *
165 * Estimate the norm of inv(op(A)).
166 *
167 AINVNM = 0.0D+0
168 *
169 KASE = 0
170 10 CONTINUE
171 CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
172 IF( KASE.NE.0 ) THEN
173 IF( KASE.EQ.2 ) THEN
174 *
175 * Multiply by R.
176 *
177 DO I = 1, N
178 WORK( I ) = WORK( I ) * RWORK( I )
179 END DO
180 *
181 IF (NOTRANS) THEN
182 CALL ZGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
183 $ WORK, N, INFO )
184 ELSE
185 CALL ZGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV,
186 $ WORK, N, INFO )
187 ENDIF
188 *
189 * Multiply by inv(C).
190 *
191 IF ( CAPPLY ) THEN
192 DO I = 1, N
193 WORK( I ) = WORK( I ) * C( I )
194 END DO
195 END IF
196 ELSE
197 *
198 * Multiply by inv(C**H).
199 *
200 IF ( CAPPLY ) THEN
201 DO I = 1, N
202 WORK( I ) = WORK( I ) * C( I )
203 END DO
204 END IF
205 *
206 IF ( NOTRANS ) THEN
207 CALL ZGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV,
208 $ WORK, N, INFO )
209 ELSE
210 CALL ZGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
211 $ WORK, N, INFO )
212 END IF
213 *
214 * Multiply by R.
215 *
216 DO I = 1, N
217 WORK( I ) = WORK( I ) * RWORK( I )
218 END DO
219 END IF
220 GO TO 10
221 END IF
222 *
223 * Compute the estimate of the reciprocal condition number.
224 *
225 IF( AINVNM .NE. 0.0D+0 )
226 $ ZLA_GERCOND_C = 1.0D+0 / AINVNM
227 *
228 RETURN
229 *
230 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 Aguments ..
16 CHARACTER TRANS
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_GERCOND_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 * TRANS (input) CHARACTER*1
36 * Specifies the form of the system of equations:
37 * = 'N': A * X = B (No transpose)
38 * = 'T': A**T * X = B (Transpose)
39 * = 'C': A**H * X = B (Conjugate Transpose = Transpose)
40 *
41 * N (input) INTEGER
42 * The number of linear equations, i.e., the order of the
43 * matrix A. N >= 0.
44 *
45 * A (input) COMPLEX*16 array, dimension (LDA,N)
46 * On entry, the N-by-N matrix A
47 *
48 * LDA (input) INTEGER
49 * The leading dimension of the array A. LDA >= max(1,N).
50 *
51 * AF (input) COMPLEX*16 array, dimension (LDAF,N)
52 * The factors L and U from the factorization
53 * A = P*L*U as computed by ZGETRF.
54 *
55 * LDAF (input) INTEGER
56 * The leading dimension of the array AF. LDAF >= max(1,N).
57 *
58 * IPIV (input) INTEGER array, dimension (N)
59 * The pivot indices from the factorization A = P*L*U
60 * as computed by ZGETRF; row i of the matrix was interchanged
61 * with row IPIV(i).
62 *
63 * C (input) DOUBLE PRECISION array, dimension (N)
64 * The vector C in the formula op(A) * inv(diag(C)).
65 *
66 * CAPPLY (input) LOGICAL
67 * If .TRUE. then access the vector C in the formula above.
68 *
69 * INFO (output) INTEGER
70 * = 0: Successful exit.
71 * i > 0: The ith argument is invalid.
72 *
73 * WORK (input) COMPLEX*16 array, dimension (2*N).
74 * Workspace.
75 *
76 * RWORK (input) DOUBLE PRECISION array, dimension (N).
77 * Workspace.
78 *
79 * =====================================================================
80 *
81 * .. Local Scalars ..
82 LOGICAL NOTRANS
83 INTEGER KASE, I, J
84 DOUBLE PRECISION AINVNM, ANORM, TMP
85 COMPLEX*16 ZDUM
86 * ..
87 * .. Local Arrays ..
88 INTEGER ISAVE( 3 )
89 * ..
90 * .. External Functions ..
91 LOGICAL LSAME
92 EXTERNAL LSAME
93 * ..
94 * .. External Subroutines ..
95 EXTERNAL ZLACN2, ZGETRS, XERBLA
96 * ..
97 * .. Intrinsic Functions ..
98 INTRINSIC ABS, MAX, REAL, DIMAG
99 * ..
100 * .. Statement Functions ..
101 DOUBLE PRECISION CABS1
102 * ..
103 * .. Statement Function Definitions ..
104 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
105 * ..
106 * .. Executable Statements ..
107 ZLA_GERCOND_C = 0.0D+0
108 *
109 INFO = 0
110 NOTRANS = LSAME( TRANS, 'N' )
111 IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT.
112 $ LSAME( TRANS, 'C' ) ) THEN
113 ELSE IF( N.LT.0 ) THEN
114 INFO = -2
115 END IF
116 IF( INFO.NE.0 ) THEN
117 CALL XERBLA( 'ZLA_GERCOND_C', -INFO )
118 RETURN
119 END IF
120 *
121 * Compute norm of op(A)*op2(C).
122 *
123 ANORM = 0.0D+0
124 IF ( NOTRANS ) THEN
125 DO I = 1, N
126 TMP = 0.0D+0
127 IF ( CAPPLY ) THEN
128 DO J = 1, N
129 TMP = TMP + CABS1( A( I, J ) ) / C( J )
130 END DO
131 ELSE
132 DO J = 1, N
133 TMP = TMP + CABS1( A( I, J ) )
134 END DO
135 END IF
136 RWORK( I ) = TMP
137 ANORM = MAX( ANORM, TMP )
138 END DO
139 ELSE
140 DO I = 1, N
141 TMP = 0.0D+0
142 IF ( CAPPLY ) THEN
143 DO J = 1, N
144 TMP = TMP + CABS1( A( J, I ) ) / C( J )
145 END DO
146 ELSE
147 DO J = 1, N
148 TMP = TMP + CABS1( A( J, I ) )
149 END DO
150 END IF
151 RWORK( I ) = TMP
152 ANORM = MAX( ANORM, TMP )
153 END DO
154 END IF
155 *
156 * Quick return if possible.
157 *
158 IF( N.EQ.0 ) THEN
159 ZLA_GERCOND_C = 1.0D+0
160 RETURN
161 ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
162 RETURN
163 END IF
164 *
165 * Estimate the norm of inv(op(A)).
166 *
167 AINVNM = 0.0D+0
168 *
169 KASE = 0
170 10 CONTINUE
171 CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
172 IF( KASE.NE.0 ) THEN
173 IF( KASE.EQ.2 ) THEN
174 *
175 * Multiply by R.
176 *
177 DO I = 1, N
178 WORK( I ) = WORK( I ) * RWORK( I )
179 END DO
180 *
181 IF (NOTRANS) THEN
182 CALL ZGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
183 $ WORK, N, INFO )
184 ELSE
185 CALL ZGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV,
186 $ WORK, N, INFO )
187 ENDIF
188 *
189 * Multiply by inv(C).
190 *
191 IF ( CAPPLY ) THEN
192 DO I = 1, N
193 WORK( I ) = WORK( I ) * C( I )
194 END DO
195 END IF
196 ELSE
197 *
198 * Multiply by inv(C**H).
199 *
200 IF ( CAPPLY ) THEN
201 DO I = 1, N
202 WORK( I ) = WORK( I ) * C( I )
203 END DO
204 END IF
205 *
206 IF ( NOTRANS ) THEN
207 CALL ZGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV,
208 $ WORK, N, INFO )
209 ELSE
210 CALL ZGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
211 $ WORK, N, INFO )
212 END IF
213 *
214 * Multiply by R.
215 *
216 DO I = 1, N
217 WORK( I ) = WORK( I ) * RWORK( I )
218 END DO
219 END IF
220 GO TO 10
221 END IF
222 *
223 * Compute the estimate of the reciprocal condition number.
224 *
225 IF( AINVNM .NE. 0.0D+0 )
226 $ ZLA_GERCOND_C = 1.0D+0 / AINVNM
227 *
228 RETURN
229 *
230 END