1 SUBROUTINE ZGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK,
2 $ INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
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**H).
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**H).
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
2 $ INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
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**H).
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**H).
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