1 /* ../netlib/dlacn2.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib;
2 on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */
3 #include "FLA_f2c.h" /* Table of constant values */
4 static integer c__1 = 1;
5 static doublereal c_b11 = 1.;
6 /* > \brief \b DLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matr ix-vector products. */
7 /* =========== DOCUMENTATION =========== */
8 /* Online html documentation available at */
9 /* http://www.netlib.org/lapack/explore-html/ */
10 /* > \htmlonly */
11 /* > Download DLACN2 + dependencies */
12 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlacn2. f"> */
13 /* > [TGZ]</a> */
14 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlacn2. f"> */
15 /* > [ZIP]</a> */
16 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlacn2. f"> */
17 /* > [TXT]</a> */
18 /* > \endhtmlonly */
19 /* Definition: */
20 /* =========== */
21 /* SUBROUTINE DLACN2( N, V, X, ISGN, EST, KASE, ISAVE ) */
22 /* .. Scalar Arguments .. */
23 /* INTEGER KASE, N */
24 /* DOUBLE PRECISION EST */
25 /* .. */
26 /* .. Array Arguments .. */
27 /* INTEGER ISGN( * ), ISAVE( 3 ) */
28 /* DOUBLE PRECISION V( * ), X( * ) */
29 /* .. */
30 /* > \par Purpose: */
31 /* ============= */
32 /* > */
33 /* > \verbatim */
34 /* > */
35 /* > DLACN2 estimates the 1-norm of a square, real matrix A. */
36 /* > Reverse communication is used for evaluating matrix-vector products. */
37 /* > \endverbatim */
38 /* Arguments: */
39 /* ========== */
40 /* > \param[in] N */
41 /* > \verbatim */
42 /* > N is INTEGER */
43 /* > The order of the matrix. N >= 1. */
44 /* > \endverbatim */
45 /* > */
46 /* > \param[out] V */
47 /* > \verbatim */
48 /* > V is DOUBLE PRECISION array, dimension (N) */
49 /* > On the final return, V = A*W, where EST = norm(V)/norm(W) */
50 /* > (W is not returned). */
51 /* > \endverbatim */
52 /* > */
53 /* > \param[in,out] X */
54 /* > \verbatim */
55 /* > X is DOUBLE PRECISION array, dimension (N) */
56 /* > On an intermediate return, X should be overwritten by */
57 /* > A * X, if KASE=1, */
58 /* > A**T * X, if KASE=2, */
59 /* > and DLACN2 must be re-called with all the other parameters */
60 /* > unchanged. */
61 /* > \endverbatim */
62 /* > */
63 /* > \param[out] ISGN */
64 /* > \verbatim */
65 /* > ISGN is INTEGER array, dimension (N) */
66 /* > \endverbatim */
67 /* > */
68 /* > \param[in,out] EST */
69 /* > \verbatim */
70 /* > EST is DOUBLE PRECISION */
71 /* > On entry with KASE = 1 or 2 and ISAVE(1) = 3, EST should be */
72 /* > unchanged from the previous call to DLACN2. */
73 /* > On exit, EST is an estimate (a lower bound) for norm(A). */
74 /* > \endverbatim */
75 /* > */
76 /* > \param[in,out] KASE */
77 /* > \verbatim */
78 /* > KASE is INTEGER */
79 /* > On the initial call to DLACN2, KASE should be 0. */
80 /* > On an intermediate return, KASE will be 1 or 2, indicating */
81 /* > whether X should be overwritten by A * X or A**T * X. */
82 /* > On the final return from DLACN2, KASE will again be 0. */
83 /* > \endverbatim */
84 /* > */
85 /* > \param[in,out] ISAVE */
86 /* > \verbatim */
87 /* > ISAVE is INTEGER array, dimension (3) */
88 /* > ISAVE is used to save variables between calls to DLACN2 */
89 /* > \endverbatim */
90 /* Authors: */
91 /* ======== */
92 /* > \author Univ. of Tennessee */
93 /* > \author Univ. of California Berkeley */
94 /* > \author Univ. of Colorado Denver */
95 /* > \author NAG Ltd. */
96 /* > \date September 2012 */
97 /* > \ingroup doubleOTHERauxiliary */
98 /* > \par Further Details: */
99 /* ===================== */
100 /* > */
101 /* > \verbatim */
102 /* > */
103 /* > Originally named SONEST, dated March 16, 1988. */
104 /* > */
105 /* > This is a thread safe version of DLACON, which uses the array ISAVE */
106 /* > in place of a SAVE statement, as follows: */
107 /* > */
108 /* > DLACON DLACN2 */
109 /* > JUMP ISAVE(1) */
110 /* > J ISAVE(2) */
111 /* > ITER ISAVE(3) */
112 /* > \endverbatim */
113 /* > \par Contributors: */
114 /* ================== */
115 /* > */
116 /* > Nick Higham, University of Manchester */
117 /* > \par References: */
118 /* ================ */
119 /* > */
120 /* > N.J. Higham, "FORTRAN codes for estimating the one-norm of */
121 /* > a real or complex matrix, with applications to condition estimation", */
122 /* > ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988. */
123 /* > */
124 /* ===================================================================== */
125 /* Subroutine */
dlacn2_(integer * n,doublereal * v,doublereal * x,integer * isgn,doublereal * est,integer * kase,integer * isave)126 int dlacn2_(integer *n, doublereal *v, doublereal *x, integer *isgn, doublereal *est, integer *kase, integer *isave)
127 {
128 /* System generated locals */
129 integer i__1;
130 doublereal d__1;
131 /* Builtin functions */
132 double d_sign(doublereal *, doublereal *);
133 integer i_dnnt(doublereal *);
134 /* Local variables */
135 integer i__;
136 doublereal temp;
137 extern doublereal dasum_(integer *, doublereal *, integer *);
138 integer jlast;
139 extern /* Subroutine */
140 int dcopy_(integer *, doublereal *, integer *, doublereal *, integer *);
141 extern integer idamax_(integer *, doublereal *, integer *);
142 doublereal altsgn, estold;
143 /* -- LAPACK auxiliary routine (version 3.4.2) -- */
144 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
145 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
146 /* September 2012 */
147 /* .. Scalar Arguments .. */
148 /* .. */
149 /* .. Array Arguments .. */
150 /* .. */
151 /* ===================================================================== */
152 /* .. Parameters .. */
153 /* .. */
154 /* .. Local Scalars .. */
155 /* .. */
156 /* .. External Functions .. */
157 /* .. */
158 /* .. External Subroutines .. */
159 /* .. */
160 /* .. Intrinsic Functions .. */
161 /* .. */
162 /* .. Executable Statements .. */
163 /* Parameter adjustments */
164 --isave;
165 --isgn;
166 --x;
167 --v;
168 /* Function Body */
169 if (*kase == 0)
170 {
171 i__1 = *n;
172 for (i__ = 1;
173 i__ <= i__1;
174 ++i__)
175 {
176 x[i__] = 1. / (doublereal) (*n);
177 /* L10: */
178 }
179 *kase = 1;
180 isave[1] = 1;
181 return 0;
182 }
183 switch (isave[1])
184 {
185 case 1:
186 goto L20;
187 case 2:
188 goto L40;
189 case 3:
190 goto L70;
191 case 4:
192 goto L110;
193 case 5:
194 goto L140;
195 }
196 /* ................ ENTRY (ISAVE( 1 ) = 1) */
197 /* FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X. */
198 L20:
199 if (*n == 1)
200 {
201 v[1] = x[1];
202 *est = f2c_abs(v[1]);
203 /* ... QUIT */
204 goto L150;
205 }
206 *est = dasum_(n, &x[1], &c__1);
207 i__1 = *n;
208 for (i__ = 1;
209 i__ <= i__1;
210 ++i__)
211 {
212 x[i__] = d_sign(&c_b11, &x[i__]);
213 isgn[i__] = i_dnnt(&x[i__]);
214 /* L30: */
215 }
216 *kase = 2;
217 isave[1] = 2;
218 return 0;
219 /* ................ ENTRY (ISAVE( 1 ) = 2) */
220 /* FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. */
221 L40:
222 isave[2] = idamax_(n, &x[1], &c__1);
223 isave[3] = 2;
224 /* MAIN LOOP - ITERATIONS 2,3,...,ITMAX. */
225 L50:
226 i__1 = *n;
227 for (i__ = 1;
228 i__ <= i__1;
229 ++i__)
230 {
231 x[i__] = 0.;
232 /* L60: */
233 }
234 x[isave[2]] = 1.;
235 *kase = 1;
236 isave[1] = 3;
237 return 0;
238 /* ................ ENTRY (ISAVE( 1 ) = 3) */
239 /* X HAS BEEN OVERWRITTEN BY A*X. */
240 L70:
241 dcopy_(n, &x[1], &c__1, &v[1], &c__1);
242 estold = *est;
243 *est = dasum_(n, &v[1], &c__1);
244 i__1 = *n;
245 for (i__ = 1;
246 i__ <= i__1;
247 ++i__)
248 {
249 d__1 = d_sign(&c_b11, &x[i__]);
250 if (i_dnnt(&d__1) != isgn[i__])
251 {
252 goto L90;
253 }
254 /* L80: */
255 }
256 /* REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED. */
257 goto L120;
258 L90: /* TEST FOR CYCLING. */
259 if (*est <= estold)
260 {
261 goto L120;
262 }
263 i__1 = *n;
264 for (i__ = 1;
265 i__ <= i__1;
266 ++i__)
267 {
268 x[i__] = d_sign(&c_b11, &x[i__]);
269 isgn[i__] = i_dnnt(&x[i__]);
270 /* L100: */
271 }
272 *kase = 2;
273 isave[1] = 4;
274 return 0;
275 /* ................ ENTRY (ISAVE( 1 ) = 4) */
276 /* X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. */
277 L110:
278 jlast = isave[2];
279 isave[2] = idamax_(n, &x[1], &c__1);
280 if (x[jlast] != (d__1 = x[isave[2]], f2c_abs(d__1)) && isave[3] < 5)
281 {
282 ++isave[3];
283 goto L50;
284 }
285 /* ITERATION COMPLETE. FINAL STAGE. */
286 L120:
287 altsgn = 1.;
288 i__1 = *n;
289 for (i__ = 1;
290 i__ <= i__1;
291 ++i__)
292 {
293 x[i__] = altsgn * ((doublereal) (i__ - 1) / (doublereal) (*n - 1) + 1.);
294 altsgn = -altsgn;
295 /* L130: */
296 }
297 *kase = 1;
298 isave[1] = 5;
299 return 0;
300 /* ................ ENTRY (ISAVE( 1 ) = 5) */
301 /* X HAS BEEN OVERWRITTEN BY A*X. */
302 L140:
303 temp = dasum_(n, &x[1], &c__1) / (doublereal) (*n * 3) * 2.;
304 if (temp > *est)
305 {
306 dcopy_(n, &x[1], &c__1, &v[1], &c__1);
307 *est = temp;
308 }
309 L150:
310 *kase = 0;
311 return 0;
312 /* End of DLACN2 */
313 }
314 /* dlacn2_ */
315