1 DOUBLE PRECISION FUNCTION ZLA_PORCOND_C( UPLO, N, A, LDA, AF,
2 $ LDAF, C, CAPPLY, 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 LOGICAL CAPPLY
18 INTEGER N, LDA, LDAF, INFO
19 * ..
20 * .. Array Arguments ..
21 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), WORK( * )
22 DOUBLE PRECISION C( * ), RWORK( * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * ZLA_PORCOND_C Computes the infinity norm condition number of
29 * op(A) * inv(diag(C)) where C is a DOUBLE PRECISION 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 triangular factor U or L from the Cholesky factorization
50 * A = U**H*U or A = L*L**H, as computed by ZPOTRF.
51 *
52 * LDAF (input) INTEGER
53 * The leading dimension of the array AF. LDAF >= max(1,N).
54 *
55 * C (input) DOUBLE PRECISION array, dimension (N)
56 * The vector C in the formula op(A) * inv(diag(C)).
57 *
58 * CAPPLY (input) LOGICAL
59 * If .TRUE. then access the vector C in the formula above.
60 *
61 * INFO (output) INTEGER
62 * = 0: Successful exit.
63 * i > 0: The ith argument is invalid.
64 *
65 * WORK (input) COMPLEX*16 array, dimension (2*N).
66 * Workspace.
67 *
68 * RWORK (input) DOUBLE PRECISION array, dimension (N).
69 * Workspace.
70 *
71 * =====================================================================
72 *
73 * .. Local Scalars ..
74 INTEGER KASE
75 DOUBLE PRECISION AINVNM, ANORM, TMP
76 INTEGER I, J
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, ZPOTRS, XERBLA
89 * ..
90 * .. Intrinsic Functions ..
91 INTRINSIC ABS, MAX, REAL, DIMAG
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_PORCOND_C = 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_PORCOND_C', -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 IF ( CAPPLY ) THEN
121 DO J = 1, I
122 TMP = TMP + CABS1( A( J, I ) ) / C( J )
123 END DO
124 DO J = I+1, N
125 TMP = TMP + CABS1( A( I, J ) ) / C( J )
126 END DO
127 ELSE
128 DO J = 1, I
129 TMP = TMP + CABS1( A( J, I ) )
130 END DO
131 DO J = I+1, N
132 TMP = TMP + CABS1( A( I, J ) )
133 END DO
134 END IF
135 RWORK( I ) = TMP
136 ANORM = MAX( ANORM, TMP )
137 END DO
138 ELSE
139 DO I = 1, N
140 TMP = 0.0D+0
141 IF ( CAPPLY ) THEN
142 DO J = 1, I
143 TMP = TMP + CABS1( A( I, J ) ) / C( J )
144 END DO
145 DO J = I+1, N
146 TMP = TMP + CABS1( A( J, I ) ) / C( J )
147 END DO
148 ELSE
149 DO J = 1, I
150 TMP = TMP + CABS1( A( I, J ) )
151 END DO
152 DO J = I+1, N
153 TMP = TMP + CABS1( A( J, I ) )
154 END DO
155 END IF
156 RWORK( I ) = TMP
157 ANORM = MAX( ANORM, TMP )
158 END DO
159 END IF
160 *
161 * Quick return if possible.
162 *
163 IF( N.EQ.0 ) THEN
164 ZLA_PORCOND_C = 1.0D+0
165 RETURN
166 ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
167 RETURN
168 END IF
169 *
170 * Estimate the norm of inv(op(A)).
171 *
172 AINVNM = 0.0D+0
173 *
174 KASE = 0
175 10 CONTINUE
176 CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
177 IF( KASE.NE.0 ) THEN
178 IF( KASE.EQ.2 ) THEN
179 *
180 * Multiply by R.
181 *
182 DO I = 1, N
183 WORK( I ) = WORK( I ) * RWORK( I )
184 END DO
185 *
186 IF ( UP ) THEN
187 CALL ZPOTRS( 'U', N, 1, AF, LDAF,
188 $ WORK, N, INFO )
189 ELSE
190 CALL ZPOTRS( 'L', N, 1, AF, LDAF,
191 $ WORK, N, INFO )
192 ENDIF
193 *
194 * Multiply by inv(C).
195 *
196 IF ( CAPPLY ) THEN
197 DO I = 1, N
198 WORK( I ) = WORK( I ) * C( I )
199 END DO
200 END IF
201 ELSE
202 *
203 * Multiply by inv(C**H).
204 *
205 IF ( CAPPLY ) THEN
206 DO I = 1, N
207 WORK( I ) = WORK( I ) * C( I )
208 END DO
209 END IF
210 *
211 IF ( UP ) THEN
212 CALL ZPOTRS( 'U', N, 1, AF, LDAF,
213 $ WORK, N, INFO )
214 ELSE
215 CALL ZPOTRS( 'L', N, 1, AF, LDAF,
216 $ WORK, N, INFO )
217 END IF
218 *
219 * Multiply by R.
220 *
221 DO I = 1, N
222 WORK( I ) = WORK( I ) * RWORK( I )
223 END DO
224 END IF
225 GO TO 10
226 END IF
227 *
228 * Compute the estimate of the reciprocal condition number.
229 *
230 IF( AINVNM .NE. 0.0D+0 )
231 $ ZLA_PORCOND_C = 1.0D+0 / AINVNM
232 *
233 RETURN
234 *
235 END
2 $ LDAF, C, CAPPLY, 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 LOGICAL CAPPLY
18 INTEGER N, LDA, LDAF, INFO
19 * ..
20 * .. Array Arguments ..
21 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), WORK( * )
22 DOUBLE PRECISION C( * ), RWORK( * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * ZLA_PORCOND_C Computes the infinity norm condition number of
29 * op(A) * inv(diag(C)) where C is a DOUBLE PRECISION 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 triangular factor U or L from the Cholesky factorization
50 * A = U**H*U or A = L*L**H, as computed by ZPOTRF.
51 *
52 * LDAF (input) INTEGER
53 * The leading dimension of the array AF. LDAF >= max(1,N).
54 *
55 * C (input) DOUBLE PRECISION array, dimension (N)
56 * The vector C in the formula op(A) * inv(diag(C)).
57 *
58 * CAPPLY (input) LOGICAL
59 * If .TRUE. then access the vector C in the formula above.
60 *
61 * INFO (output) INTEGER
62 * = 0: Successful exit.
63 * i > 0: The ith argument is invalid.
64 *
65 * WORK (input) COMPLEX*16 array, dimension (2*N).
66 * Workspace.
67 *
68 * RWORK (input) DOUBLE PRECISION array, dimension (N).
69 * Workspace.
70 *
71 * =====================================================================
72 *
73 * .. Local Scalars ..
74 INTEGER KASE
75 DOUBLE PRECISION AINVNM, ANORM, TMP
76 INTEGER I, J
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, ZPOTRS, XERBLA
89 * ..
90 * .. Intrinsic Functions ..
91 INTRINSIC ABS, MAX, REAL, DIMAG
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_PORCOND_C = 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_PORCOND_C', -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 IF ( CAPPLY ) THEN
121 DO J = 1, I
122 TMP = TMP + CABS1( A( J, I ) ) / C( J )
123 END DO
124 DO J = I+1, N
125 TMP = TMP + CABS1( A( I, J ) ) / C( J )
126 END DO
127 ELSE
128 DO J = 1, I
129 TMP = TMP + CABS1( A( J, I ) )
130 END DO
131 DO J = I+1, N
132 TMP = TMP + CABS1( A( I, J ) )
133 END DO
134 END IF
135 RWORK( I ) = TMP
136 ANORM = MAX( ANORM, TMP )
137 END DO
138 ELSE
139 DO I = 1, N
140 TMP = 0.0D+0
141 IF ( CAPPLY ) THEN
142 DO J = 1, I
143 TMP = TMP + CABS1( A( I, J ) ) / C( J )
144 END DO
145 DO J = I+1, N
146 TMP = TMP + CABS1( A( J, I ) ) / C( J )
147 END DO
148 ELSE
149 DO J = 1, I
150 TMP = TMP + CABS1( A( I, J ) )
151 END DO
152 DO J = I+1, N
153 TMP = TMP + CABS1( A( J, I ) )
154 END DO
155 END IF
156 RWORK( I ) = TMP
157 ANORM = MAX( ANORM, TMP )
158 END DO
159 END IF
160 *
161 * Quick return if possible.
162 *
163 IF( N.EQ.0 ) THEN
164 ZLA_PORCOND_C = 1.0D+0
165 RETURN
166 ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
167 RETURN
168 END IF
169 *
170 * Estimate the norm of inv(op(A)).
171 *
172 AINVNM = 0.0D+0
173 *
174 KASE = 0
175 10 CONTINUE
176 CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
177 IF( KASE.NE.0 ) THEN
178 IF( KASE.EQ.2 ) THEN
179 *
180 * Multiply by R.
181 *
182 DO I = 1, N
183 WORK( I ) = WORK( I ) * RWORK( I )
184 END DO
185 *
186 IF ( UP ) THEN
187 CALL ZPOTRS( 'U', N, 1, AF, LDAF,
188 $ WORK, N, INFO )
189 ELSE
190 CALL ZPOTRS( 'L', N, 1, AF, LDAF,
191 $ WORK, N, INFO )
192 ENDIF
193 *
194 * Multiply by inv(C).
195 *
196 IF ( CAPPLY ) THEN
197 DO I = 1, N
198 WORK( I ) = WORK( I ) * C( I )
199 END DO
200 END IF
201 ELSE
202 *
203 * Multiply by inv(C**H).
204 *
205 IF ( CAPPLY ) THEN
206 DO I = 1, N
207 WORK( I ) = WORK( I ) * C( I )
208 END DO
209 END IF
210 *
211 IF ( UP ) THEN
212 CALL ZPOTRS( 'U', N, 1, AF, LDAF,
213 $ WORK, N, INFO )
214 ELSE
215 CALL ZPOTRS( 'L', N, 1, AF, LDAF,
216 $ WORK, N, INFO )
217 END IF
218 *
219 * Multiply by R.
220 *
221 DO I = 1, N
222 WORK( I ) = WORK( I ) * RWORK( I )
223 END DO
224 END IF
225 GO TO 10
226 END IF
227 *
228 * Compute the estimate of the reciprocal condition number.
229 *
230 IF( AINVNM .NE. 0.0D+0 )
231 $ ZLA_PORCOND_C = 1.0D+0 / AINVNM
232 *
233 RETURN
234 *
235 END