1 /* ./src_f77/dlansp.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
dlansp_(char * norm,char * uplo,integer * n,doublereal * ap,doublereal * work,ftnlen norm_len,ftnlen uplo_len)12 doublereal dlansp_(char *norm, char *uplo, integer *n, doublereal *ap,
13 doublereal *work, ftnlen norm_len, ftnlen uplo_len)
14 {
15 /* System generated locals */
16 integer i__1, i__2;
17 doublereal ret_val, d__1, d__2, d__3;
18
19 /* Builtin functions */
20 double sqrt(doublereal);
21
22 /* Local variables */
23 static integer i__, j, k;
24 static doublereal sum, absa, scale;
25 extern logical lsame_(char *, char *, ftnlen, ftnlen);
26 static doublereal value;
27 extern /* Subroutine */ int dlassq_(integer *, doublereal *, integer *,
28 doublereal *, doublereal *);
29
30
31 /* -- LAPACK auxiliary routine (version 3.0) -- */
32 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
33 /* Courant Institute, Argonne National Lab, and Rice University */
34 /* October 31, 1992 */
35
36 /* .. Scalar Arguments .. */
37 /* .. */
38 /* .. Array Arguments .. */
39 /* .. */
40
41 /* Purpose */
42 /* ======= */
43
44 /* DLANSP returns the value of the one norm, or the Frobenius norm, or */
45 /* the infinity norm, or the element of largest absolute value of a */
46 /* real symmetric matrix A, supplied in packed form. */
47
48 /* Description */
49 /* =========== */
50
51 /* DLANSP returns the value */
52
53 /* DLANSP = ( max(abs(A(i,j))), NORM = 'M' or 'm' */
54 /* ( */
55 /* ( norm1(A), NORM = '1', 'O' or 'o' */
56 /* ( */
57 /* ( normI(A), NORM = 'I' or 'i' */
58 /* ( */
59 /* ( normF(A), NORM = 'F', 'f', 'E' or 'e' */
60
61 /* where norm1 denotes the one norm of a matrix (maximum column sum), */
62 /* normI denotes the infinity norm of a matrix (maximum row sum) and */
63 /* normF denotes the Frobenius norm of a matrix (square root of sum of */
64 /* squares). Note that max(abs(A(i,j))) is not a matrix norm. */
65
66 /* Arguments */
67 /* ========= */
68
69 /* NORM (input) CHARACTER*1 */
70 /* Specifies the value to be returned in DLANSP as described */
71 /* above. */
72
73 /* UPLO (input) CHARACTER*1 */
74 /* Specifies whether the upper or lower triangular part of the */
75 /* symmetric matrix A is supplied. */
76 /* = 'U': Upper triangular part of A is supplied */
77 /* = 'L': Lower triangular part of A is supplied */
78
79 /* N (input) INTEGER */
80 /* The order of the matrix A. N >= 0. When N = 0, DLANSP is */
81 /* set to zero. */
82
83 /* AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2) */
84 /* The upper or lower triangle of the symmetric matrix A, packed */
85 /* columnwise in a linear array. The j-th column of A is stored */
86 /* in the array AP as follows: */
87 /* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */
88 /* if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. */
89
90 /* WORK (workspace) DOUBLE PRECISION array, dimension (LWORK), */
91 /* where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, */
92 /* WORK is not referenced. */
93
94 /* ===================================================================== */
95
96 /* .. Parameters .. */
97 /* .. */
98 /* .. Local Scalars .. */
99 /* .. */
100 /* .. External Subroutines .. */
101 /* .. */
102 /* .. External Functions .. */
103 /* .. */
104 /* .. Intrinsic Functions .. */
105 /* .. */
106 /* .. Executable Statements .. */
107
108 /* Parameter adjustments */
109 --work;
110 --ap;
111
112 /* Function Body */
113 if (*n == 0) {
114 value = 0.;
115 } else if (lsame_(norm, "M", (ftnlen)1, (ftnlen)1)) {
116
117 /* Find max(abs(A(i,j))). */
118
119 value = 0.;
120 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
121 k = 1;
122 i__1 = *n;
123 for (j = 1; j <= i__1; ++j) {
124 i__2 = k + j - 1;
125 for (i__ = k; i__ <= i__2; ++i__) {
126 /* Computing MAX */
127 d__2 = value, d__3 = (d__1 = ap[i__], abs(d__1));
128 value = max(d__2,d__3);
129 /* L10: */
130 }
131 k += j;
132 /* L20: */
133 }
134 } else {
135 k = 1;
136 i__1 = *n;
137 for (j = 1; j <= i__1; ++j) {
138 i__2 = k + *n - j;
139 for (i__ = k; i__ <= i__2; ++i__) {
140 /* Computing MAX */
141 d__2 = value, d__3 = (d__1 = ap[i__], abs(d__1));
142 value = max(d__2,d__3);
143 /* L30: */
144 }
145 k = k + *n - j + 1;
146 /* L40: */
147 }
148 }
149 } else if (lsame_(norm, "I", (ftnlen)1, (ftnlen)1) || lsame_(norm, "O", (
150 ftnlen)1, (ftnlen)1) || *(unsigned char *)norm == '1') {
151
152 /* Find normI(A) ( = norm1(A), since A is symmetric). */
153
154 value = 0.;
155 k = 1;
156 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
157 i__1 = *n;
158 for (j = 1; j <= i__1; ++j) {
159 sum = 0.;
160 i__2 = j - 1;
161 for (i__ = 1; i__ <= i__2; ++i__) {
162 absa = (d__1 = ap[k], abs(d__1));
163 sum += absa;
164 work[i__] += absa;
165 ++k;
166 /* L50: */
167 }
168 work[j] = sum + (d__1 = ap[k], abs(d__1));
169 ++k;
170 /* L60: */
171 }
172 i__1 = *n;
173 for (i__ = 1; i__ <= i__1; ++i__) {
174 /* Computing MAX */
175 d__1 = value, d__2 = work[i__];
176 value = max(d__1,d__2);
177 /* L70: */
178 }
179 } else {
180 i__1 = *n;
181 for (i__ = 1; i__ <= i__1; ++i__) {
182 work[i__] = 0.;
183 /* L80: */
184 }
185 i__1 = *n;
186 for (j = 1; j <= i__1; ++j) {
187 sum = work[j] + (d__1 = ap[k], abs(d__1));
188 ++k;
189 i__2 = *n;
190 for (i__ = j + 1; i__ <= i__2; ++i__) {
191 absa = (d__1 = ap[k], abs(d__1));
192 sum += absa;
193 work[i__] += absa;
194 ++k;
195 /* L90: */
196 }
197 value = max(value,sum);
198 /* L100: */
199 }
200 }
201 } else if (lsame_(norm, "F", (ftnlen)1, (ftnlen)1) || lsame_(norm, "E", (
202 ftnlen)1, (ftnlen)1)) {
203
204 /* Find normF(A). */
205
206 scale = 0.;
207 sum = 1.;
208 k = 2;
209 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
210 i__1 = *n;
211 for (j = 2; j <= i__1; ++j) {
212 i__2 = j - 1;
213 dlassq_(&i__2, &ap[k], &c__1, &scale, &sum);
214 k += j;
215 /* L110: */
216 }
217 } else {
218 i__1 = *n - 1;
219 for (j = 1; j <= i__1; ++j) {
220 i__2 = *n - j;
221 dlassq_(&i__2, &ap[k], &c__1, &scale, &sum);
222 k = k + *n - j + 1;
223 /* L120: */
224 }
225 }
226 sum *= 2;
227 k = 1;
228 i__1 = *n;
229 for (i__ = 1; i__ <= i__1; ++i__) {
230 if (ap[k] != 0.) {
231 absa = (d__1 = ap[k], abs(d__1));
232 if (scale < absa) {
233 /* Computing 2nd power */
234 d__1 = scale / absa;
235 sum = sum * (d__1 * d__1) + 1.;
236 scale = absa;
237 } else {
238 /* Computing 2nd power */
239 d__1 = absa / scale;
240 sum += d__1 * d__1;
241 }
242 }
243 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
244 k = k + i__ + 1;
245 } else {
246 k = k + *n - i__ + 1;
247 }
248 /* L130: */
249 }
250 value = scale * sqrt(sum);
251 }
252
253 ret_val = value;
254 return ret_val;
255
256 /* End of DLANSP */
257
258 } /* dlansp_ */
259
260