1: *> \brief \b ZLA_PORCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for Hermitian positive-definite matrices.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZLA_PORCOND_C + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zla_porcond_c.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zla_porcond_c.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zla_porcond_c.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * DOUBLE PRECISION FUNCTION ZLA_PORCOND_C( UPLO, N, A, LDA, AF,
22: * LDAF, C, CAPPLY, INFO,
23: * WORK, RWORK )
24: *
25: * .. Scalar Arguments ..
26: * CHARACTER UPLO
27: * LOGICAL CAPPLY
28: * INTEGER N, LDA, LDAF, INFO
29: * ..
30: * .. Array Arguments ..
31: * COMPLEX*16 A( LDA, * ), AF( LDAF, * ), WORK( * )
32: * DOUBLE PRECISION C( * ), RWORK( * )
33: * ..
34: *
35: *
36: *> \par Purpose:
37: * =============
38: *>
39: *> \verbatim
40: *>
41: *> ZLA_PORCOND_C Computes the infinity norm condition number of
42: *> op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector
43: *> \endverbatim
44: *
45: * Arguments:
46: * ==========
47: *
48: *> \param[in] UPLO
49: *> \verbatim
50: *> UPLO is CHARACTER*1
51: *> = 'U': Upper triangle of A is stored;
52: *> = 'L': Lower triangle of A is stored.
53: *> \endverbatim
54: *>
55: *> \param[in] N
56: *> \verbatim
57: *> N is INTEGER
58: *> The number of linear equations, i.e., the order of the
59: *> matrix A. N >= 0.
60: *> \endverbatim
61: *>
62: *> \param[in] A
63: *> \verbatim
64: *> A is COMPLEX*16 array, dimension (LDA,N)
65: *> On entry, the N-by-N matrix A
66: *> \endverbatim
67: *>
68: *> \param[in] LDA
69: *> \verbatim
70: *> LDA is INTEGER
71: *> The leading dimension of the array A. LDA >= max(1,N).
72: *> \endverbatim
73: *>
74: *> \param[in] AF
75: *> \verbatim
76: *> AF is COMPLEX*16 array, dimension (LDAF,N)
77: *> The triangular factor U or L from the Cholesky factorization
78: *> A = U**H*U or A = L*L**H, as computed by ZPOTRF.
79: *> \endverbatim
80: *>
81: *> \param[in] LDAF
82: *> \verbatim
83: *> LDAF is INTEGER
84: *> The leading dimension of the array AF. LDAF >= max(1,N).
85: *> \endverbatim
86: *>
87: *> \param[in] C
88: *> \verbatim
89: *> C is DOUBLE PRECISION array, dimension (N)
90: *> The vector C in the formula op(A) * inv(diag(C)).
91: *> \endverbatim
92: *>
93: *> \param[in] CAPPLY
94: *> \verbatim
95: *> CAPPLY is LOGICAL
96: *> If .TRUE. then access the vector C in the formula above.
97: *> \endverbatim
98: *>
99: *> \param[out] INFO
100: *> \verbatim
101: *> INFO is INTEGER
102: *> = 0: Successful exit.
103: *> i > 0: The ith argument is invalid.
104: *> \endverbatim
105: *>
106: *> \param[out] WORK
107: *> \verbatim
108: *> WORK is COMPLEX*16 array, dimension (2*N).
109: *> Workspace.
110: *> \endverbatim
111: *>
112: *> \param[out] RWORK
113: *> \verbatim
114: *> RWORK is DOUBLE PRECISION array, dimension (N).
115: *> Workspace.
116: *> \endverbatim
117: *
118: * Authors:
119: * ========
120: *
121: *> \author Univ. of Tennessee
122: *> \author Univ. of California Berkeley
123: *> \author Univ. of Colorado Denver
124: *> \author NAG Ltd.
125: *
126: *> \ingroup complex16POcomputational
127: *
128: * =====================================================================
129: DOUBLE PRECISION FUNCTION ZLA_PORCOND_C( UPLO, N, A, LDA, AF,
130: $ LDAF, C, CAPPLY, INFO,
131: $ WORK, RWORK )
132: *
133: * -- LAPACK computational routine --
134: * -- LAPACK is a software package provided by Univ. of Tennessee, --
135: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
136: *
137: * .. Scalar Arguments ..
138: CHARACTER UPLO
139: LOGICAL CAPPLY
140: INTEGER N, LDA, LDAF, INFO
141: * ..
142: * .. Array Arguments ..
143: COMPLEX*16 A( LDA, * ), AF( LDAF, * ), WORK( * )
144: DOUBLE PRECISION C( * ), RWORK( * )
145: * ..
146: *
147: * =====================================================================
148: *
149: * .. Local Scalars ..
150: INTEGER KASE
151: DOUBLE PRECISION AINVNM, ANORM, TMP
152: INTEGER I, J
153: LOGICAL UP, UPPER
154: COMPLEX*16 ZDUM
155: * ..
156: * .. Local Arrays ..
157: INTEGER ISAVE( 3 )
158: * ..
159: * .. External Functions ..
160: LOGICAL LSAME
161: EXTERNAL LSAME
162: * ..
163: * .. External Subroutines ..
164: EXTERNAL ZLACN2, ZPOTRS, XERBLA
165: * ..
166: * .. Intrinsic Functions ..
167: INTRINSIC ABS, MAX, REAL, DIMAG
168: * ..
169: * .. Statement Functions ..
170: DOUBLE PRECISION CABS1
171: * ..
172: * .. Statement Function Definitions ..
173: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
174: * ..
175: * .. Executable Statements ..
176: *
177: ZLA_PORCOND_C = 0.0D+0
178: *
179: INFO = 0
180: UPPER = LSAME( UPLO, 'U' )
181: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
182: INFO = -1
183: ELSE IF( N.LT.0 ) THEN
184: INFO = -2
185: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
186: INFO = -4
187: ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
188: INFO = -6
189: END IF
190: IF( INFO.NE.0 ) THEN
191: CALL XERBLA( 'ZLA_PORCOND_C', -INFO )
192: RETURN
193: END IF
194: UP = .FALSE.
195: IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
196: *
197: * Compute norm of op(A)*op2(C).
198: *
199: ANORM = 0.0D+0
200: IF ( UP ) THEN
201: DO I = 1, N
202: TMP = 0.0D+0
203: IF ( CAPPLY ) THEN
204: DO J = 1, I
205: TMP = TMP + CABS1( A( J, I ) ) / C( J )
206: END DO
207: DO J = I+1, N
208: TMP = TMP + CABS1( A( I, J ) ) / C( J )
209: END DO
210: ELSE
211: DO J = 1, I
212: TMP = TMP + CABS1( A( J, I ) )
213: END DO
214: DO J = I+1, N
215: TMP = TMP + CABS1( A( I, J ) )
216: END DO
217: END IF
218: RWORK( I ) = TMP
219: ANORM = MAX( ANORM, TMP )
220: END DO
221: ELSE
222: DO I = 1, N
223: TMP = 0.0D+0
224: IF ( CAPPLY ) THEN
225: DO J = 1, I
226: TMP = TMP + CABS1( A( I, J ) ) / C( J )
227: END DO
228: DO J = I+1, N
229: TMP = TMP + CABS1( A( J, I ) ) / C( J )
230: END DO
231: ELSE
232: DO J = 1, I
233: TMP = TMP + CABS1( A( I, J ) )
234: END DO
235: DO J = I+1, N
236: TMP = TMP + CABS1( A( J, I ) )
237: END DO
238: END IF
239: RWORK( I ) = TMP
240: ANORM = MAX( ANORM, TMP )
241: END DO
242: END IF
243: *
244: * Quick return if possible.
245: *
246: IF( N.EQ.0 ) THEN
247: ZLA_PORCOND_C = 1.0D+0
248: RETURN
249: ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
250: RETURN
251: END IF
252: *
253: * Estimate the norm of inv(op(A)).
254: *
255: AINVNM = 0.0D+0
256: *
257: KASE = 0
258: 10 CONTINUE
259: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
260: IF( KASE.NE.0 ) THEN
261: IF( KASE.EQ.2 ) THEN
262: *
263: * Multiply by R.
264: *
265: DO I = 1, N
266: WORK( I ) = WORK( I ) * RWORK( I )
267: END DO
268: *
269: IF ( UP ) THEN
270: CALL ZPOTRS( 'U', N, 1, AF, LDAF,
271: $ WORK, N, INFO )
272: ELSE
273: CALL ZPOTRS( 'L', N, 1, AF, LDAF,
274: $ WORK, N, INFO )
275: ENDIF
276: *
277: * Multiply by inv(C).
278: *
279: IF ( CAPPLY ) THEN
280: DO I = 1, N
281: WORK( I ) = WORK( I ) * C( I )
282: END DO
283: END IF
284: ELSE
285: *
286: * Multiply by inv(C**H).
287: *
288: IF ( CAPPLY ) THEN
289: DO I = 1, N
290: WORK( I ) = WORK( I ) * C( I )
291: END DO
292: END IF
293: *
294: IF ( UP ) THEN
295: CALL ZPOTRS( 'U', N, 1, AF, LDAF,
296: $ WORK, N, INFO )
297: ELSE
298: CALL ZPOTRS( 'L', N, 1, AF, LDAF,
299: $ WORK, N, INFO )
300: END IF
301: *
302: * Multiply by R.
303: *
304: DO I = 1, N
305: WORK( I ) = WORK( I ) * RWORK( I )
306: END DO
307: END IF
308: GO TO 10
309: END IF
310: *
311: * Compute the estimate of the reciprocal condition number.
312: *
313: IF( AINVNM .NE. 0.0D+0 )
314: $ ZLA_PORCOND_C = 1.0D+0 / AINVNM
315: *
316: RETURN
317: *
318: * End of ZLA_PORCOND_C
319: *
320: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zla_porcond_c.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack