1 /* ./src_f77/cpocon.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
cpocon_(char * uplo,integer * n,complex * a,integer * lda,real * anorm,real * rcond,complex * work,real * rwork,integer * info,ftnlen uplo_len)12 /* Subroutine */ int cpocon_(char *uplo, integer *n, complex *a, integer *lda,
13 real *anorm, real *rcond, complex *work, real *rwork, integer *info,
14 ftnlen uplo_len)
15 {
16 /* System generated locals */
17 integer a_dim1, a_offset, i__1;
18 real r__1, r__2;
19
20 /* Builtin functions */
21 double r_imag(complex *);
22
23 /* Local variables */
24 static integer ix, kase;
25 static real scale;
26 extern logical lsame_(char *, char *, ftnlen, ftnlen);
27 static logical upper;
28 extern /* Subroutine */ int clacon_(integer *, complex *, complex *, real
29 *, integer *);
30 extern integer icamax_(integer *, complex *, integer *);
31 static real scalel;
32 extern doublereal slamch_(char *, ftnlen);
33 static real scaleu;
34 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
35 static real ainvnm;
36 extern /* Subroutine */ int clatrs_(char *, char *, char *, char *,
37 integer *, complex *, integer *, complex *, real *, real *,
38 integer *, ftnlen, ftnlen, ftnlen, ftnlen), csrscl_(integer *,
39 real *, complex *, integer *);
40 static char normin[1];
41 static real smlnum;
42
43
44 /* -- LAPACK routine (version 3.0) -- */
45 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
46 /* Courant Institute, Argonne National Lab, and Rice University */
47 /* March 31, 1993 */
48
49 /* .. Scalar Arguments .. */
50 /* .. */
51 /* .. Array Arguments .. */
52 /* .. */
53
54 /* Purpose */
55 /* ======= */
56
57 /* CPOCON estimates the reciprocal of the condition number (in the */
58 /* 1-norm) of a complex Hermitian positive definite matrix using the */
59 /* Cholesky factorization A = U**H*U or A = L*L**H computed by CPOTRF. */
60
61 /* An estimate is obtained for norm(inv(A)), and the reciprocal of the */
62 /* condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))). */
63
64 /* Arguments */
65 /* ========= */
66
67 /* UPLO (input) CHARACTER*1 */
68 /* = 'U': Upper triangle of A is stored; */
69 /* = 'L': Lower triangle of A is stored. */
70
71 /* N (input) INTEGER */
72 /* The order of the matrix A. N >= 0. */
73
74 /* A (input) COMPLEX array, dimension (LDA,N) */
75 /* The triangular factor U or L from the Cholesky factorization */
76 /* A = U**H*U or A = L*L**H, as computed by CPOTRF. */
77
78 /* LDA (input) INTEGER */
79 /* The leading dimension of the array A. LDA >= max(1,N). */
80
81 /* ANORM (input) REAL */
82 /* The 1-norm (or infinity-norm) of the Hermitian matrix A. */
83
84 /* RCOND (output) REAL */
85 /* The reciprocal of the condition number of the matrix A, */
86 /* computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an */
87 /* estimate of the 1-norm of inv(A) computed in this routine. */
88
89 /* WORK (workspace) COMPLEX array, dimension (2*N) */
90
91 /* RWORK (workspace) REAL array, dimension (N) */
92
93 /* INFO (output) INTEGER */
94 /* = 0: successful exit */
95 /* < 0: if INFO = -i, the i-th argument had an illegal value */
96
97 /* ===================================================================== */
98
99 /* .. Parameters .. */
100 /* .. */
101 /* .. Local Scalars .. */
102 /* .. */
103 /* .. External Functions .. */
104 /* .. */
105 /* .. External Subroutines .. */
106 /* .. */
107 /* .. Intrinsic Functions .. */
108 /* .. */
109 /* .. Statement Functions .. */
110 /* .. */
111 /* .. Statement Function definitions .. */
112 /* .. */
113 /* .. Executable Statements .. */
114
115 /* Test the input parameters. */
116
117 /* Parameter adjustments */
118 a_dim1 = *lda;
119 a_offset = 1 + a_dim1;
120 a -= a_offset;
121 --work;
122 --rwork;
123
124 /* Function Body */
125 *info = 0;
126 upper = lsame_(uplo, "U", (ftnlen)1, (ftnlen)1);
127 if (! upper && ! lsame_(uplo, "L", (ftnlen)1, (ftnlen)1)) {
128 *info = -1;
129 } else if (*n < 0) {
130 *info = -2;
131 } else if (*lda < max(1,*n)) {
132 *info = -4;
133 } else if (*anorm < 0.f) {
134 *info = -5;
135 }
136 if (*info != 0) {
137 i__1 = -(*info);
138 xerbla_("CPOCON", &i__1, (ftnlen)6);
139 return 0;
140 }
141
142 /* Quick return if possible */
143
144 *rcond = 0.f;
145 if (*n == 0) {
146 *rcond = 1.f;
147 return 0;
148 } else if (*anorm == 0.f) {
149 return 0;
150 }
151
152 smlnum = slamch_("Safe minimum", (ftnlen)12);
153
154 /* Estimate the 1-norm of inv(A). */
155
156 kase = 0;
157 *(unsigned char *)normin = 'N';
158 L10:
159 clacon_(n, &work[*n + 1], &work[1], &ainvnm, &kase);
160 if (kase != 0) {
161 if (upper) {
162
163 /* Multiply by inv(U'). */
164
165 clatrs_("Upper", "Conjugate transpose", "Non-unit", normin, n, &a[
166 a_offset], lda, &work[1], &scalel, &rwork[1], info, (
167 ftnlen)5, (ftnlen)19, (ftnlen)8, (ftnlen)1);
168 *(unsigned char *)normin = 'Y';
169
170 /* Multiply by inv(U). */
171
172 clatrs_("Upper", "No transpose", "Non-unit", normin, n, &a[
173 a_offset], lda, &work[1], &scaleu, &rwork[1], info, (
174 ftnlen)5, (ftnlen)12, (ftnlen)8, (ftnlen)1);
175 } else {
176
177 /* Multiply by inv(L). */
178
179 clatrs_("Lower", "No transpose", "Non-unit", normin, n, &a[
180 a_offset], lda, &work[1], &scalel, &rwork[1], info, (
181 ftnlen)5, (ftnlen)12, (ftnlen)8, (ftnlen)1);
182 *(unsigned char *)normin = 'Y';
183
184 /* Multiply by inv(L'). */
185
186 clatrs_("Lower", "Conjugate transpose", "Non-unit", normin, n, &a[
187 a_offset], lda, &work[1], &scaleu, &rwork[1], info, (
188 ftnlen)5, (ftnlen)19, (ftnlen)8, (ftnlen)1);
189 }
190
191 /* Multiply by 1/SCALE if doing so will not cause overflow. */
192
193 scale = scalel * scaleu;
194 if (scale != 1.f) {
195 ix = icamax_(n, &work[1], &c__1);
196 i__1 = ix;
197 if (scale < ((r__1 = work[i__1].r, dabs(r__1)) + (r__2 = r_imag(&
198 work[ix]), dabs(r__2))) * smlnum || scale == 0.f) {
199 goto L20;
200 }
201 csrscl_(n, &scale, &work[1], &c__1);
202 }
203 goto L10;
204 }
205
206 /* Compute the estimate of the reciprocal condition number. */
207
208 if (ainvnm != 0.f) {
209 *rcond = 1.f / ainvnm / *anorm;
210 }
211
212 L20:
213 return 0;
214
215 /* End of CPOCON */
216
217 } /* cpocon_ */
218
219