1: SUBROUTINE ZGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK,
2: $ INFO )
3: *
4: * -- LAPACK routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH.
10: *
11: * .. Scalar Arguments ..
12: CHARACTER NORM
13: INTEGER INFO, LDA, N
14: DOUBLE PRECISION ANORM, RCOND
15: * ..
16: * .. Array Arguments ..
17: DOUBLE PRECISION RWORK( * )
18: COMPLEX*16 A( LDA, * ), WORK( * )
19: * ..
20: *
21: * Purpose
22: * =======
23: *
24: * ZGECON estimates the reciprocal of the condition number of a general
25: * complex matrix A, in either the 1-norm or the infinity-norm, using
26: * the LU factorization computed by ZGETRF.
27: *
28: * An estimate is obtained for norm(inv(A)), and the reciprocal of the
29: * condition number is computed as
30: * RCOND = 1 / ( norm(A) * norm(inv(A)) ).
31: *
32: * Arguments
33: * =========
34: *
35: * NORM (input) CHARACTER*1
36: * Specifies whether the 1-norm condition number or the
37: * infinity-norm condition number is required:
38: * = '1' or 'O': 1-norm;
39: * = 'I': Infinity-norm.
40: *
41: * N (input) INTEGER
42: * The order of the matrix A. N >= 0.
43: *
44: * A (input) COMPLEX*16 array, dimension (LDA,N)
45: * The factors L and U from the factorization A = P*L*U
46: * as computed by ZGETRF.
47: *
48: * LDA (input) INTEGER
49: * The leading dimension of the array A. LDA >= max(1,N).
50: *
51: * ANORM (input) DOUBLE PRECISION
52: * If NORM = '1' or 'O', the 1-norm of the original matrix A.
53: * If NORM = 'I', the infinity-norm of the original matrix A.
54: *
55: * RCOND (output) DOUBLE PRECISION
56: * The reciprocal of the condition number of the matrix A,
57: * computed as RCOND = 1/(norm(A) * norm(inv(A))).
58: *
59: * WORK (workspace) COMPLEX*16 array, dimension (2*N)
60: *
61: * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
62: *
63: * INFO (output) INTEGER
64: * = 0: successful exit
65: * < 0: if INFO = -i, the i-th argument had an illegal value
66: *
67: * =====================================================================
68: *
69: * .. Parameters ..
70: DOUBLE PRECISION ONE, ZERO
71: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
72: * ..
73: * .. Local Scalars ..
74: LOGICAL ONENRM
75: CHARACTER NORMIN
76: INTEGER IX, KASE, KASE1
77: DOUBLE PRECISION AINVNM, SCALE, SL, SMLNUM, SU
78: COMPLEX*16 ZDUM
79: * ..
80: * .. Local Arrays ..
81: INTEGER ISAVE( 3 )
82: * ..
83: * .. External Functions ..
84: LOGICAL LSAME
85: INTEGER IZAMAX
86: DOUBLE PRECISION DLAMCH
87: EXTERNAL LSAME, IZAMAX, DLAMCH
88: * ..
89: * .. External Subroutines ..
90: EXTERNAL XERBLA, ZDRSCL, ZLACN2, ZLATRS
91: * ..
92: * .. Intrinsic Functions ..
93: INTRINSIC ABS, DBLE, DIMAG, MAX
94: * ..
95: * .. Statement Functions ..
96: DOUBLE PRECISION CABS1
97: * ..
98: * .. Statement Function definitions ..
99: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
100: * ..
101: * .. Executable Statements ..
102: *
103: * Test the input parameters.
104: *
105: INFO = 0
106: ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
107: IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
108: INFO = -1
109: ELSE IF( N.LT.0 ) THEN
110: INFO = -2
111: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
112: INFO = -4
113: ELSE IF( ANORM.LT.ZERO ) THEN
114: INFO = -5
115: END IF
116: IF( INFO.NE.0 ) THEN
117: CALL XERBLA( 'ZGECON', -INFO )
118: RETURN
119: END IF
120: *
121: * Quick return if possible
122: *
123: RCOND = ZERO
124: IF( N.EQ.0 ) THEN
125: RCOND = ONE
126: RETURN
127: ELSE IF( ANORM.EQ.ZERO ) THEN
128: RETURN
129: END IF
130: *
131: SMLNUM = DLAMCH( 'Safe minimum' )
132: *
133: * Estimate the norm of inv(A).
134: *
135: AINVNM = ZERO
136: NORMIN = 'N'
137: IF( ONENRM ) THEN
138: KASE1 = 1
139: ELSE
140: KASE1 = 2
141: END IF
142: KASE = 0
143: 10 CONTINUE
144: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
145: IF( KASE.NE.0 ) THEN
146: IF( KASE.EQ.KASE1 ) THEN
147: *
148: * Multiply by inv(L).
149: *
150: CALL ZLATRS( 'Lower', 'No transpose', 'Unit', NORMIN, N, A,
151: $ LDA, WORK, SL, RWORK, INFO )
152: *
153: * Multiply by inv(U).
154: *
155: CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
156: $ A, LDA, WORK, SU, RWORK( N+1 ), INFO )
157: ELSE
158: *
159: * Multiply by inv(U').
160: *
161: CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit',
162: $ NORMIN, N, A, LDA, WORK, SU, RWORK( N+1 ),
163: $ INFO )
164: *
165: * Multiply by inv(L').
166: *
167: CALL ZLATRS( 'Lower', 'Conjugate transpose', 'Unit', NORMIN,
168: $ N, A, LDA, WORK, SL, RWORK, INFO )
169: END IF
170: *
171: * Divide X by 1/(SL*SU) if doing so will not cause overflow.
172: *
173: SCALE = SL*SU
174: NORMIN = 'Y'
175: IF( SCALE.NE.ONE ) THEN
176: IX = IZAMAX( N, WORK, 1 )
177: IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
178: $ GO TO 20
179: CALL ZDRSCL( N, SCALE, WORK, 1 )
180: END IF
181: GO TO 10
182: END IF
183: *
184: * Compute the estimate of the reciprocal condition number.
185: *
186: IF( AINVNM.NE.ZERO )
187: $ RCOND = ( ONE / AINVNM ) / ANORM
188: *
189: 20 CONTINUE
190: RETURN
191: *
192: * End of ZGECON
193: *
194: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgecon.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack