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