1 /* ./src_f77/zgecon.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static integer c__1 = 1;
11
zgecon_(char * norm,integer * n,doublecomplex * a,integer * lda,doublereal * anorm,doublereal * rcond,doublecomplex * work,doublereal * rwork,integer * info,ftnlen norm_len)12 /* Subroutine */ int zgecon_(char *norm, integer *n, doublecomplex *a,
13 integer *lda, doublereal *anorm, doublereal *rcond, doublecomplex *
14 work, doublereal *rwork, integer *info, ftnlen norm_len)
15 {
16 /* System generated locals */
17 integer a_dim1, a_offset, i__1;
18 doublereal d__1, d__2;
19
20 /* Builtin functions */
21 double d_imag(doublecomplex *);
22
23 /* Local variables */
24 static doublereal sl;
25 static integer ix;
26 static doublereal su;
27 static integer kase, kase1;
28 static doublereal scale;
29 extern logical lsame_(char *, char *, ftnlen, ftnlen);
30 extern doublereal dlamch_(char *, ftnlen);
31 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen), zlacon_(
32 integer *, doublecomplex *, doublecomplex *, doublereal *,
33 integer *);
34 static doublereal ainvnm;
35 extern integer izamax_(integer *, doublecomplex *, integer *);
36 static logical onenrm;
37 extern /* Subroutine */ int zdrscl_(integer *, doublereal *,
38 doublecomplex *, integer *);
39 static char normin[1];
40 static doublereal smlnum;
41 extern /* Subroutine */ int zlatrs_(char *, char *, char *, char *,
42 integer *, doublecomplex *, integer *, doublecomplex *,
43 doublereal *, doublereal *, integer *, ftnlen, ftnlen, ftnlen,
44 ftnlen);
45
46
47 /* -- LAPACK routine (version 3.0) -- */
48 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
49 /* Courant Institute, Argonne National Lab, and Rice University */
50 /* March 31, 1993 */
51
52 /* .. Scalar Arguments .. */
53 /* .. */
54 /* .. Array Arguments .. */
55 /* .. */
56
57 /* Purpose */
58 /* ======= */
59
60 /* ZGECON estimates the reciprocal of the condition number of a general */
61 /* complex matrix A, in either the 1-norm or the infinity-norm, using */
62 /* the LU factorization computed by ZGETRF. */
63
64 /* An estimate is obtained for norm(inv(A)), and the reciprocal of the */
65 /* condition number is computed as */
66 /* RCOND = 1 / ( norm(A) * norm(inv(A)) ). */
67
68 /* Arguments */
69 /* ========= */
70
71 /* NORM (input) CHARACTER*1 */
72 /* Specifies whether the 1-norm condition number or the */
73 /* infinity-norm condition number is required: */
74 /* = '1' or 'O': 1-norm; */
75 /* = 'I': Infinity-norm. */
76
77 /* N (input) INTEGER */
78 /* The order of the matrix A. N >= 0. */
79
80 /* A (input) COMPLEX*16 array, dimension (LDA,N) */
81 /* The factors L and U from the factorization A = P*L*U */
82 /* as computed by ZGETRF. */
83
84 /* LDA (input) INTEGER */
85 /* The leading dimension of the array A. LDA >= max(1,N). */
86
87 /* ANORM (input) DOUBLE PRECISION */
88 /* If NORM = '1' or 'O', the 1-norm of the original matrix A. */
89 /* If NORM = 'I', the infinity-norm of the original matrix A. */
90
91 /* RCOND (output) DOUBLE PRECISION */
92 /* The reciprocal of the condition number of the matrix A, */
93 /* computed as RCOND = 1/(norm(A) * norm(inv(A))). */
94
95 /* WORK (workspace) COMPLEX*16 array, dimension (2*N) */
96
97 /* RWORK (workspace) DOUBLE PRECISION array, dimension (2*N) */
98
99 /* INFO (output) INTEGER */
100 /* = 0: successful exit */
101 /* < 0: if INFO = -i, the i-th argument had an illegal value */
102
103 /* ===================================================================== */
104
105 /* .. Parameters .. */
106 /* .. */
107 /* .. Local Scalars .. */
108 /* .. */
109 /* .. External Functions .. */
110 /* .. */
111 /* .. External Subroutines .. */
112 /* .. */
113 /* .. Intrinsic Functions .. */
114 /* .. */
115 /* .. Statement Functions .. */
116 /* .. */
117 /* .. Statement Function definitions .. */
118 /* .. */
119 /* .. Executable Statements .. */
120
121 /* Test the input parameters. */
122
123 /* Parameter adjustments */
124 a_dim1 = *lda;
125 a_offset = 1 + a_dim1;
126 a -= a_offset;
127 --work;
128 --rwork;
129
130 /* Function Body */
131 *info = 0;
132 onenrm = *(unsigned char *)norm == '1' || lsame_(norm, "O", (ftnlen)1, (
133 ftnlen)1);
134 if (! onenrm && ! lsame_(norm, "I", (ftnlen)1, (ftnlen)1)) {
135 *info = -1;
136 } else if (*n < 0) {
137 *info = -2;
138 } else if (*lda < max(1,*n)) {
139 *info = -4;
140 } else if (*anorm < 0.) {
141 *info = -5;
142 }
143 if (*info != 0) {
144 i__1 = -(*info);
145 xerbla_("ZGECON", &i__1, (ftnlen)6);
146 return 0;
147 }
148
149 /* Quick return if possible */
150
151 *rcond = 0.;
152 if (*n == 0) {
153 *rcond = 1.;
154 return 0;
155 } else if (*anorm == 0.) {
156 return 0;
157 }
158
159 smlnum = dlamch_("Safe minimum", (ftnlen)12);
160
161 /* Estimate the norm of inv(A). */
162
163 ainvnm = 0.;
164 *(unsigned char *)normin = 'N';
165 if (onenrm) {
166 kase1 = 1;
167 } else {
168 kase1 = 2;
169 }
170 kase = 0;
171 L10:
172 zlacon_(n, &work[*n + 1], &work[1], &ainvnm, &kase);
173 if (kase != 0) {
174 if (kase == kase1) {
175
176 /* Multiply by inv(L). */
177
178 zlatrs_("Lower", "No transpose", "Unit", normin, n, &a[a_offset],
179 lda, &work[1], &sl, &rwork[1], info, (ftnlen)5, (ftnlen)
180 12, (ftnlen)4, (ftnlen)1);
181
182 /* Multiply by inv(U). */
183
184 zlatrs_("Upper", "No transpose", "Non-unit", normin, n, &a[
185 a_offset], lda, &work[1], &su, &rwork[*n + 1], info, (
186 ftnlen)5, (ftnlen)12, (ftnlen)8, (ftnlen)1);
187 } else {
188
189 /* Multiply by inv(U'). */
190
191 zlatrs_("Upper", "Conjugate transpose", "Non-unit", normin, n, &a[
192 a_offset], lda, &work[1], &su, &rwork[*n + 1], info, (
193 ftnlen)5, (ftnlen)19, (ftnlen)8, (ftnlen)1);
194
195 /* Multiply by inv(L'). */
196
197 zlatrs_("Lower", "Conjugate transpose", "Unit", normin, n, &a[
198 a_offset], lda, &work[1], &sl, &rwork[1], info, (ftnlen)5,
199 (ftnlen)19, (ftnlen)4, (ftnlen)1);
200 }
201
202 /* Divide X by 1/(SL*SU) if doing so will not cause overflow. */
203
204 scale = sl * su;
205 *(unsigned char *)normin = 'Y';
206 if (scale != 1.) {
207 ix = izamax_(n, &work[1], &c__1);
208 i__1 = ix;
209 if (scale < ((d__1 = work[i__1].r, abs(d__1)) + (d__2 = d_imag(&
210 work[ix]), abs(d__2))) * smlnum || scale == 0.) {
211 goto L20;
212 }
213 zdrscl_(n, &scale, &work[1], &c__1);
214 }
215 goto L10;
216 }
217
218 /* Compute the estimate of the reciprocal condition number. */
219
220 if (ainvnm != 0.) {
221 *rcond = 1. / ainvnm / *anorm;
222 }
223
224 L20:
225 return 0;
226
227 /* End of ZGECON */
228
229 } /* zgecon_ */
230
231