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