1 /* ./src_f77/clantp.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
clantp_(char * norm,char * uplo,char * diag,integer * n,complex * ap,real * work,ftnlen norm_len,ftnlen uplo_len,ftnlen diag_len)12 doublereal clantp_(char *norm, char *uplo, char *diag, integer *n, complex *
13 ap, real *work, ftnlen norm_len, ftnlen uplo_len, ftnlen diag_len)
14 {
15 /* System generated locals */
16 integer i__1, i__2;
17 real ret_val, r__1, r__2;
18
19 /* Builtin functions */
20 double c_abs(complex *), sqrt(doublereal);
21
22 /* Local variables */
23 static integer i__, j, k;
24 static real sum, scale;
25 static logical udiag;
26 extern logical lsame_(char *, char *, ftnlen, ftnlen);
27 static real value;
28 extern /* Subroutine */ int classq_(integer *, complex *, integer *, real
29 *, real *);
30
31
32 /* -- LAPACK auxiliary routine (version 3.0) -- */
33 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
34 /* Courant Institute, Argonne National Lab, and Rice University */
35 /* October 31, 1992 */
36
37 /* .. Scalar Arguments .. */
38 /* .. */
39 /* .. Array Arguments .. */
40 /* .. */
41
42 /* Purpose */
43 /* ======= */
44
45 /* CLANTP returns the value of the one norm, or the Frobenius norm, or */
46 /* the infinity norm, or the element of largest absolute value of a */
47 /* triangular matrix A, supplied in packed form. */
48
49 /* Description */
50 /* =========== */
51
52 /* CLANTP returns the value */
53
54 /* CLANTP = ( max(abs(A(i,j))), NORM = 'M' or 'm' */
55 /* ( */
56 /* ( norm1(A), NORM = '1', 'O' or 'o' */
57 /* ( */
58 /* ( normI(A), NORM = 'I' or 'i' */
59 /* ( */
60 /* ( normF(A), NORM = 'F', 'f', 'E' or 'e' */
61
62 /* where norm1 denotes the one norm of a matrix (maximum column sum), */
63 /* normI denotes the infinity norm of a matrix (maximum row sum) and */
64 /* normF denotes the Frobenius norm of a matrix (square root of sum of */
65 /* squares). Note that max(abs(A(i,j))) is not a matrix norm. */
66
67 /* Arguments */
68 /* ========= */
69
70 /* NORM (input) CHARACTER*1 */
71 /* Specifies the value to be returned in CLANTP as described */
72 /* above. */
73
74 /* UPLO (input) CHARACTER*1 */
75 /* Specifies whether the matrix A is upper or lower triangular. */
76 /* = 'U': Upper triangular */
77 /* = 'L': Lower triangular */
78
79 /* DIAG (input) CHARACTER*1 */
80 /* Specifies whether or not the matrix A is unit triangular. */
81 /* = 'N': Non-unit triangular */
82 /* = 'U': Unit triangular */
83
84 /* N (input) INTEGER */
85 /* The order of the matrix A. N >= 0. When N = 0, CLANTP is */
86 /* set to zero. */
87
88 /* AP (input) COMPLEX array, dimension (N*(N+1)/2) */
89 /* The upper or lower triangular matrix A, packed columnwise in */
90 /* a linear array. The j-th column of A is stored in the array */
91 /* AP as follows: */
92 /* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */
93 /* if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. */
94 /* Note that when DIAG = 'U', the elements of the array AP */
95 /* corresponding to the diagonal elements of the matrix A are */
96 /* not referenced, but are assumed to be one. */
97
98 /* WORK (workspace) REAL array, dimension (LWORK), */
99 /* where LWORK >= N when NORM = 'I'; otherwise, WORK is not */
100 /* referenced. */
101
102 /* ===================================================================== */
103
104 /* .. Parameters .. */
105 /* .. */
106 /* .. Local Scalars .. */
107 /* .. */
108 /* .. External Functions .. */
109 /* .. */
110 /* .. External Subroutines .. */
111 /* .. */
112 /* .. Intrinsic Functions .. */
113 /* .. */
114 /* .. Executable Statements .. */
115
116 /* Parameter adjustments */
117 --work;
118 --ap;
119
120 /* Function Body */
121 if (*n == 0) {
122 value = 0.f;
123 } else if (lsame_(norm, "M", (ftnlen)1, (ftnlen)1)) {
124
125 /* Find max(abs(A(i,j))). */
126
127 k = 1;
128 if (lsame_(diag, "U", (ftnlen)1, (ftnlen)1)) {
129 value = 1.f;
130 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
131 i__1 = *n;
132 for (j = 1; j <= i__1; ++j) {
133 i__2 = k + j - 2;
134 for (i__ = k; i__ <= i__2; ++i__) {
135 /* Computing MAX */
136 r__1 = value, r__2 = c_abs(&ap[i__]);
137 value = dmax(r__1,r__2);
138 /* L10: */
139 }
140 k += j;
141 /* L20: */
142 }
143 } else {
144 i__1 = *n;
145 for (j = 1; j <= i__1; ++j) {
146 i__2 = k + *n - j;
147 for (i__ = k + 1; i__ <= i__2; ++i__) {
148 /* Computing MAX */
149 r__1 = value, r__2 = c_abs(&ap[i__]);
150 value = dmax(r__1,r__2);
151 /* L30: */
152 }
153 k = k + *n - j + 1;
154 /* L40: */
155 }
156 }
157 } else {
158 value = 0.f;
159 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
160 i__1 = *n;
161 for (j = 1; j <= i__1; ++j) {
162 i__2 = k + j - 1;
163 for (i__ = k; i__ <= i__2; ++i__) {
164 /* Computing MAX */
165 r__1 = value, r__2 = c_abs(&ap[i__]);
166 value = dmax(r__1,r__2);
167 /* L50: */
168 }
169 k += j;
170 /* L60: */
171 }
172 } else {
173 i__1 = *n;
174 for (j = 1; j <= i__1; ++j) {
175 i__2 = k + *n - j;
176 for (i__ = k; i__ <= i__2; ++i__) {
177 /* Computing MAX */
178 r__1 = value, r__2 = c_abs(&ap[i__]);
179 value = dmax(r__1,r__2);
180 /* L70: */
181 }
182 k = k + *n - j + 1;
183 /* L80: */
184 }
185 }
186 }
187 } else if (lsame_(norm, "O", (ftnlen)1, (ftnlen)1) || *(unsigned char *)
188 norm == '1') {
189
190 /* Find norm1(A). */
191
192 value = 0.f;
193 k = 1;
194 udiag = lsame_(diag, "U", (ftnlen)1, (ftnlen)1);
195 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
196 i__1 = *n;
197 for (j = 1; j <= i__1; ++j) {
198 if (udiag) {
199 sum = 1.f;
200 i__2 = k + j - 2;
201 for (i__ = k; i__ <= i__2; ++i__) {
202 sum += c_abs(&ap[i__]);
203 /* L90: */
204 }
205 } else {
206 sum = 0.f;
207 i__2 = k + j - 1;
208 for (i__ = k; i__ <= i__2; ++i__) {
209 sum += c_abs(&ap[i__]);
210 /* L100: */
211 }
212 }
213 k += j;
214 value = dmax(value,sum);
215 /* L110: */
216 }
217 } else {
218 i__1 = *n;
219 for (j = 1; j <= i__1; ++j) {
220 if (udiag) {
221 sum = 1.f;
222 i__2 = k + *n - j;
223 for (i__ = k + 1; i__ <= i__2; ++i__) {
224 sum += c_abs(&ap[i__]);
225 /* L120: */
226 }
227 } else {
228 sum = 0.f;
229 i__2 = k + *n - j;
230 for (i__ = k; i__ <= i__2; ++i__) {
231 sum += c_abs(&ap[i__]);
232 /* L130: */
233 }
234 }
235 k = k + *n - j + 1;
236 value = dmax(value,sum);
237 /* L140: */
238 }
239 }
240 } else if (lsame_(norm, "I", (ftnlen)1, (ftnlen)1)) {
241
242 /* Find normI(A). */
243
244 k = 1;
245 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
246 if (lsame_(diag, "U", (ftnlen)1, (ftnlen)1)) {
247 i__1 = *n;
248 for (i__ = 1; i__ <= i__1; ++i__) {
249 work[i__] = 1.f;
250 /* L150: */
251 }
252 i__1 = *n;
253 for (j = 1; j <= i__1; ++j) {
254 i__2 = j - 1;
255 for (i__ = 1; i__ <= i__2; ++i__) {
256 work[i__] += c_abs(&ap[k]);
257 ++k;
258 /* L160: */
259 }
260 ++k;
261 /* L170: */
262 }
263 } else {
264 i__1 = *n;
265 for (i__ = 1; i__ <= i__1; ++i__) {
266 work[i__] = 0.f;
267 /* L180: */
268 }
269 i__1 = *n;
270 for (j = 1; j <= i__1; ++j) {
271 i__2 = j;
272 for (i__ = 1; i__ <= i__2; ++i__) {
273 work[i__] += c_abs(&ap[k]);
274 ++k;
275 /* L190: */
276 }
277 /* L200: */
278 }
279 }
280 } else {
281 if (lsame_(diag, "U", (ftnlen)1, (ftnlen)1)) {
282 i__1 = *n;
283 for (i__ = 1; i__ <= i__1; ++i__) {
284 work[i__] = 1.f;
285 /* L210: */
286 }
287 i__1 = *n;
288 for (j = 1; j <= i__1; ++j) {
289 ++k;
290 i__2 = *n;
291 for (i__ = j + 1; i__ <= i__2; ++i__) {
292 work[i__] += c_abs(&ap[k]);
293 ++k;
294 /* L220: */
295 }
296 /* L230: */
297 }
298 } else {
299 i__1 = *n;
300 for (i__ = 1; i__ <= i__1; ++i__) {
301 work[i__] = 0.f;
302 /* L240: */
303 }
304 i__1 = *n;
305 for (j = 1; j <= i__1; ++j) {
306 i__2 = *n;
307 for (i__ = j; i__ <= i__2; ++i__) {
308 work[i__] += c_abs(&ap[k]);
309 ++k;
310 /* L250: */
311 }
312 /* L260: */
313 }
314 }
315 }
316 value = 0.f;
317 i__1 = *n;
318 for (i__ = 1; i__ <= i__1; ++i__) {
319 /* Computing MAX */
320 r__1 = value, r__2 = work[i__];
321 value = dmax(r__1,r__2);
322 /* L270: */
323 }
324 } else if (lsame_(norm, "F", (ftnlen)1, (ftnlen)1) || lsame_(norm, "E", (
325 ftnlen)1, (ftnlen)1)) {
326
327 /* Find normF(A). */
328
329 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
330 if (lsame_(diag, "U", (ftnlen)1, (ftnlen)1)) {
331 scale = 1.f;
332 sum = (real) (*n);
333 k = 2;
334 i__1 = *n;
335 for (j = 2; j <= i__1; ++j) {
336 i__2 = j - 1;
337 classq_(&i__2, &ap[k], &c__1, &scale, &sum);
338 k += j;
339 /* L280: */
340 }
341 } else {
342 scale = 0.f;
343 sum = 1.f;
344 k = 1;
345 i__1 = *n;
346 for (j = 1; j <= i__1; ++j) {
347 classq_(&j, &ap[k], &c__1, &scale, &sum);
348 k += j;
349 /* L290: */
350 }
351 }
352 } else {
353 if (lsame_(diag, "U", (ftnlen)1, (ftnlen)1)) {
354 scale = 1.f;
355 sum = (real) (*n);
356 k = 2;
357 i__1 = *n - 1;
358 for (j = 1; j <= i__1; ++j) {
359 i__2 = *n - j;
360 classq_(&i__2, &ap[k], &c__1, &scale, &sum);
361 k = k + *n - j + 1;
362 /* L300: */
363 }
364 } else {
365 scale = 0.f;
366 sum = 1.f;
367 k = 1;
368 i__1 = *n;
369 for (j = 1; j <= i__1; ++j) {
370 i__2 = *n - j + 1;
371 classq_(&i__2, &ap[k], &c__1, &scale, &sum);
372 k = k + *n - j + 1;
373 /* L310: */
374 }
375 }
376 }
377 value = scale * sqrt(sum);
378 }
379
380 ret_val = value;
381 return ret_val;
382
383 /* End of CLANTP */
384
385 } /* clantp_ */
386
387