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