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