1 /* lapack/double/dlanhs.f -- translated by f2c (version 20050501).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11 */
12
13 #ifdef __cplusplus
14 extern "C" {
15 #endif
16 #include "v3p_netlib.h"
17
18 /* Table of constant values */
19
20 static integer c__1 = 1;
21
22 /*< DOUBLE PRECISION FUNCTION DLANHS( NORM, N, A, LDA, WORK ) >*/
dlanhs_(char * norm,integer * n,doublereal * a,integer * lda,doublereal * work,ftnlen norm_len)23 doublereal dlanhs_(char *norm, integer *n, doublereal *a, integer *lda,
24 doublereal *work, ftnlen norm_len)
25 {
26 /* System generated locals */
27 integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
28 doublereal ret_val, d__1, d__2, d__3;
29
30 /* Builtin functions */
31 double sqrt(doublereal);
32
33 /* Local variables */
34 integer i__, j;
35 doublereal sum, scale;
36 extern logical lsame_(const char *, const char *, ftnlen, ftnlen);
37 doublereal value=0;
38 extern /* Subroutine */ int dlassq_(integer *, doublereal *, integer *,
39 doublereal *, doublereal *);
40 (void)norm_len;
41
42 /* -- LAPACK auxiliary routine (version 3.0) -- */
43 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
44 /* Courant Institute, Argonne National Lab, and Rice University */
45 /* October 31, 1992 */
46
47 /* .. Scalar Arguments .. */
48 /*< CHARACTER NORM >*/
49 /*< INTEGER LDA, N >*/
50 /* .. */
51 /* .. Array Arguments .. */
52 /*< DOUBLE PRECISION A( LDA, * ), WORK( * ) >*/
53 /* .. */
54
55 /* Purpose */
56 /* ======= */
57
58 /* DLANHS returns the value of the one norm, or the Frobenius norm, or */
59 /* the infinity norm, or the element of largest absolute value of a */
60 /* Hessenberg matrix A. */
61
62 /* Description */
63 /* =========== */
64
65 /* DLANHS returns the value */
66
67 /* DLANHS = ( max(abs(A(i,j))), NORM = 'M' or 'm' */
68 /* ( */
69 /* ( norm1(A), NORM = '1', 'O' or 'o' */
70 /* ( */
71 /* ( normI(A), NORM = 'I' or 'i' */
72 /* ( */
73 /* ( normF(A), NORM = 'F', 'f', 'E' or 'e' */
74
75 /* where norm1 denotes the one norm of a matrix (maximum column sum), */
76 /* normI denotes the infinity norm of a matrix (maximum row sum) and */
77 /* normF denotes the Frobenius norm of a matrix (square root of sum of */
78 /* squares). Note that max(abs(A(i,j))) is not a matrix norm. */
79
80 /* Arguments */
81 /* ========= */
82
83 /* NORM (input) CHARACTER*1 */
84 /* Specifies the value to be returned in DLANHS as described */
85 /* above. */
86
87 /* N (input) INTEGER */
88 /* The order of the matrix A. N >= 0. When N = 0, DLANHS is */
89 /* set to zero. */
90
91 /* A (input) DOUBLE PRECISION array, dimension (LDA,N) */
92 /* The n by n upper Hessenberg matrix A; the part of A below the */
93 /* first sub-diagonal is not referenced. */
94
95 /* LDA (input) INTEGER */
96 /* The leading dimension of the array A. LDA >= max(N,1). */
97
98 /* WORK (workspace) DOUBLE PRECISION array, dimension (LWORK), */
99 /* where LWORK >= N when NORM = 'I'; otherwise, WORK is not */
100 /* referenced. */
101
102 /* ===================================================================== */
103
104 /* .. Parameters .. */
105 /*< DOUBLE PRECISION ONE, ZERO >*/
106 /*< PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) >*/
107 /* .. */
108 /* .. Local Scalars .. */
109 /*< INTEGER I, J >*/
110 /*< DOUBLE PRECISION SCALE, SUM, VALUE >*/
111 /* .. */
112 /* .. External Subroutines .. */
113 /*< EXTERNAL DLASSQ >*/
114 /* .. */
115 /* .. External Functions .. */
116 /*< LOGICAL LSAME >*/
117 /*< EXTERNAL LSAME >*/
118 /* .. */
119 /* .. Intrinsic Functions .. */
120 /*< INTRINSIC ABS, MAX, MIN, SQRT >*/
121 /* .. */
122 /* .. Executable Statements .. */
123
124 /*< IF( N.EQ.0 ) THEN >*/
125 /* Parameter adjustments */
126 a_dim1 = *lda;
127 a_offset = 1 + a_dim1;
128 a -= a_offset;
129 --work;
130
131 /* Function Body */
132 if (*n == 0) {
133 /*< VALUE = ZERO >*/
134 value = 0.;
135 /*< ELSE IF( LSAME( NORM, 'M' ) ) THEN >*/
136 } else if (lsame_(norm, "M", (ftnlen)1, (ftnlen)1)) {
137
138 /* Find max(abs(A(i,j))). */
139
140 /*< VALUE = ZERO >*/
141 value = 0.;
142 /*< DO 20 J = 1, N >*/
143 i__1 = *n;
144 for (j = 1; j <= i__1; ++j) {
145 /*< DO 10 I = 1, MIN( N, J+1 ) >*/
146 /* Computing MIN */
147 i__3 = *n, i__4 = j + 1;
148 i__2 = min(i__3,i__4);
149 for (i__ = 1; i__ <= i__2; ++i__) {
150 /*< VALUE = MAX( VALUE, ABS( A( I, J ) ) ) >*/
151 /* Computing MAX */
152 d__2 = value, d__3 = (d__1 = a[i__ + j * a_dim1], abs(d__1));
153 value = max(d__2,d__3);
154 /*< 10 CONTINUE >*/
155 /* L10: */
156 }
157 /*< 20 CONTINUE >*/
158 /* L20: */
159 }
160 /*< ELSE IF( ( LSAME( NORM, 'O' ) ) .OR. ( NORM.EQ.'1' ) ) THEN >*/
161 } else if (lsame_(norm, "O", (ftnlen)1, (ftnlen)1) || *(unsigned char *)
162 norm == '1') {
163
164 /* Find norm1(A). */
165
166 /*< VALUE = ZERO >*/
167 value = 0.;
168 /*< DO 40 J = 1, N >*/
169 i__1 = *n;
170 for (j = 1; j <= i__1; ++j) {
171 /*< SUM = ZERO >*/
172 sum = 0.;
173 /*< DO 30 I = 1, MIN( N, J+1 ) >*/
174 /* Computing MIN */
175 i__3 = *n, i__4 = j + 1;
176 i__2 = min(i__3,i__4);
177 for (i__ = 1; i__ <= i__2; ++i__) {
178 /*< SUM = SUM + ABS( A( I, J ) ) >*/
179 sum += (d__1 = a[i__ + j * a_dim1], abs(d__1));
180 /*< 30 CONTINUE >*/
181 /* L30: */
182 }
183 /*< VALUE = MAX( VALUE, SUM ) >*/
184 value = max(value,sum);
185 /*< 40 CONTINUE >*/
186 /* L40: */
187 }
188 /*< ELSE IF( LSAME( NORM, 'I' ) ) THEN >*/
189 } else if (lsame_(norm, "I", (ftnlen)1, (ftnlen)1)) {
190
191 /* Find normI(A). */
192
193 /*< DO 50 I = 1, N >*/
194 i__1 = *n;
195 for (i__ = 1; i__ <= i__1; ++i__) {
196 /*< WORK( I ) = ZERO >*/
197 work[i__] = 0.;
198 /*< 50 CONTINUE >*/
199 /* L50: */
200 }
201 /*< DO 70 J = 1, N >*/
202 i__1 = *n;
203 for (j = 1; j <= i__1; ++j) {
204 /*< DO 60 I = 1, MIN( N, J+1 ) >*/
205 /* Computing MIN */
206 i__3 = *n, i__4 = j + 1;
207 i__2 = min(i__3,i__4);
208 for (i__ = 1; i__ <= i__2; ++i__) {
209 /*< WORK( I ) = WORK( I ) + ABS( A( I, J ) ) >*/
210 work[i__] += (d__1 = a[i__ + j * a_dim1], abs(d__1));
211 /*< 60 CONTINUE >*/
212 /* L60: */
213 }
214 /*< 70 CONTINUE >*/
215 /* L70: */
216 }
217 /*< VALUE = ZERO >*/
218 value = 0.;
219 /*< DO 80 I = 1, N >*/
220 i__1 = *n;
221 for (i__ = 1; i__ <= i__1; ++i__) {
222 /*< VALUE = MAX( VALUE, WORK( I ) ) >*/
223 /* Computing MAX */
224 d__1 = value, d__2 = work[i__];
225 value = max(d__1,d__2);
226 /*< 80 CONTINUE >*/
227 /* L80: */
228 }
229 /*< ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN >*/
230 } else if (lsame_(norm, "F", (ftnlen)1, (ftnlen)1) || lsame_(norm, "E", (
231 ftnlen)1, (ftnlen)1)) {
232
233 /* Find normF(A). */
234
235 /*< SCALE = ZERO >*/
236 scale = 0.;
237 /*< SUM = ONE >*/
238 sum = 1.;
239 /*< DO 90 J = 1, N >*/
240 i__1 = *n;
241 for (j = 1; j <= i__1; ++j) {
242 /*< CALL DLASSQ( MIN( N, J+1 ), A( 1, J ), 1, SCALE, SUM ) >*/
243 /* Computing MIN */
244 i__3 = *n, i__4 = j + 1;
245 i__2 = min(i__3,i__4);
246 dlassq_(&i__2, &a[j * a_dim1 + 1], &c__1, &scale, &sum);
247 /*< 90 CONTINUE >*/
248 /* L90: */
249 }
250 /*< VALUE = SCALE*SQRT( SUM ) >*/
251 value = scale * sqrt(sum);
252 /*< END IF >*/
253 }
254
255 /*< DLANHS = VALUE >*/
256 ret_val = value;
257 /*< RETURN >*/
258 return ret_val;
259
260 /* End of DLANHS */
261
262 /*< END >*/
263 } /* dlanhs_ */
264
265 #ifdef __cplusplus
266 }
267 #endif
268