1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
|
/* iparmq.f -- translated by f2c (version 20061008).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
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
*/
#include "f2c.h"
#include "blaswrap.h"
integer iparmq_(integer *ispec, char *name__, char *opts, integer *n, integer
*ilo, integer *ihi, integer *lwork)
{
/* System generated locals */
integer ret_val, i__1, i__2;
real r__1;
/* Builtin functions */
double log(doublereal);
integer i_nint(real *);
/* Local variables */
integer nh, ns;
/* -- LAPACK auxiliary routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* Purpose */
/* ======= */
/* This program sets problem and machine dependent parameters */
/* useful for xHSEQR and its subroutines. It is called whenever */
/* ILAENV is called with 12 <= ISPEC <= 16 */
/* Arguments */
/* ========= */
/* ISPEC (input) integer scalar */
/* ISPEC specifies which tunable parameter IPARMQ should */
/* return. */
/* ISPEC=12: (INMIN) Matrices of order nmin or less */
/* are sent directly to xLAHQR, the implicit */
/* double shift QR algorithm. NMIN must be */
/* at least 11. */
/* ISPEC=13: (INWIN) Size of the deflation window. */
/* This is best set greater than or equal to */
/* the number of simultaneous shifts NS. */
/* Larger matrices benefit from larger deflation */
/* windows. */
/* ISPEC=14: (INIBL) Determines when to stop nibbling and */
/* invest in an (expensive) multi-shift QR sweep. */
/* If the aggressive early deflation subroutine */
/* finds LD converged eigenvalues from an order */
/* NW deflation window and LD.GT.(NW*NIBBLE)/100, */
/* then the next QR sweep is skipped and early */
/* deflation is applied immediately to the */
/* remaining active diagonal block. Setting */
/* IPARMQ(ISPEC=14) = 0 causes TTQRE to skip a */
/* multi-shift QR sweep whenever early deflation */
/* finds a converged eigenvalue. Setting */
/* IPARMQ(ISPEC=14) greater than or equal to 100 */
/* prevents TTQRE from skipping a multi-shift */
/* QR sweep. */
/* ISPEC=15: (NSHFTS) The number of simultaneous shifts in */
/* a multi-shift QR iteration. */
/* ISPEC=16: (IACC22) IPARMQ is set to 0, 1 or 2 with the */
/* following meanings. */
/* 0: During the multi-shift QR sweep, */
/* xLAQR5 does not accumulate reflections and */
/* does not use matrix-matrix multiply to */
/* update the far-from-diagonal matrix */
/* entries. */
/* 1: During the multi-shift QR sweep, */
/* xLAQR5 and/or xLAQRaccumulates reflections and uses */
/* matrix-matrix multiply to update the */
/* far-from-diagonal matrix entries. */
/* 2: During the multi-shift QR sweep. */
/* xLAQR5 accumulates reflections and takes */
/* advantage of 2-by-2 block structure during */
/* matrix-matrix multiplies. */
/* (If xTRMM is slower than xGEMM, then */
/* IPARMQ(ISPEC=16)=1 may be more efficient than */
/* IPARMQ(ISPEC=16)=2 despite the greater level of */
/* arithmetic work implied by the latter choice.) */
/* NAME (input) character string */
/* Name of the calling subroutine */
/* OPTS (input) character string */
/* This is a concatenation of the string arguments to */
/* TTQRE. */
/* N (input) integer scalar */
/* N is the order of the Hessenberg matrix H. */
/* ILO (input) INTEGER */
/* IHI (input) INTEGER */
/* It is assumed that H is already upper triangular */
/* in rows and columns 1:ILO-1 and IHI+1:N. */
/* LWORK (input) integer scalar */
/* The amount of workspace available. */
/* Further Details */
/* =============== */
/* Little is known about how best to choose these parameters. */
/* It is possible to use different values of the parameters */
/* for each of CHSEQR, DHSEQR, SHSEQR and ZHSEQR. */
/* It is probably best to choose different parameters for */
/* different matrices and different parameters at different */
/* times during the iteration, but this has not been */
/* implemented --- yet. */
/* The best choices of most of the parameters depend */
/* in an ill-understood way on the relative execution */
/* rate of xLAQR3 and xLAQR5 and on the nature of each */
/* particular eigenvalue problem. Experiment may be the */
/* only practical way to determine which choices are most */
/* effective. */
/* Following is a list of default values supplied by IPARMQ. */
/* These defaults may be adjusted in order to attain better */
/* performance in any particular computational environment. */
/* IPARMQ(ISPEC=12) The xLAHQR vs xLAQR0 crossover point. */
/* Default: 75. (Must be at least 11.) */
/* IPARMQ(ISPEC=13) Recommended deflation window size. */
/* This depends on ILO, IHI and NS, the */
/* number of simultaneous shifts returned */
/* by IPARMQ(ISPEC=15). The default for */
/* (IHI-ILO+1).LE.500 is NS. The default */
/* for (IHI-ILO+1).GT.500 is 3*NS/2. */
/* IPARMQ(ISPEC=14) Nibble crossover point. Default: 14. */
/* IPARMQ(ISPEC=15) Number of simultaneous shifts, NS. */
/* a multi-shift QR iteration. */
/* If IHI-ILO+1 is ... */
/* greater than ...but less ... the */
/* or equal to ... than default is */
/* 0 30 NS = 2+ */
/* 30 60 NS = 4+ */
/* 60 150 NS = 10 */
/* 150 590 NS = ** */
/* 590 3000 NS = 64 */
/* 3000 6000 NS = 128 */
/* 6000 infinity NS = 256 */
/* (+) By default matrices of this order are */
/* passed to the implicit double shift routine */
/* xLAHQR. See IPARMQ(ISPEC=12) above. These */
/* values of NS are used only in case of a rare */
/* xLAHQR failure. */
/* (**) The asterisks (**) indicate an ad-hoc */
/* function increasing from 10 to 64. */
/* IPARMQ(ISPEC=16) Select structured matrix multiply. */
/* (See ISPEC=16 above for details.) */
/* Default: 3. */
/* ================================================================ */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
if (*ispec == 15 || *ispec == 13 || *ispec == 16) {
/* ==== Set the number simultaneous shifts ==== */
nh = *ihi - *ilo + 1;
ns = 2;
if (nh >= 30) {
ns = 4;
}
if (nh >= 60) {
ns = 10;
}
if (nh >= 150) {
/* Computing MAX */
r__1 = log((real) nh) / log(2.f);
i__1 = 10, i__2 = nh / i_nint(&r__1);
ns = max(i__1,i__2);
}
if (nh >= 590) {
ns = 64;
}
if (nh >= 3000) {
ns = 128;
}
if (nh >= 6000) {
ns = 256;
}
/* Computing MAX */
i__1 = 2, i__2 = ns - ns % 2;
ns = max(i__1,i__2);
}
if (*ispec == 12) {
/* ===== Matrices of order smaller than NMIN get sent */
/* . to xLAHQR, the classic double shift algorithm. */
/* . This must be at least 11. ==== */
ret_val = 75;
} else if (*ispec == 14) {
/* ==== INIBL: skip a multi-shift qr iteration and */
/* . whenever aggressive early deflation finds */
/* . at least (NIBBLE*(window size)/100) deflations. ==== */
ret_val = 14;
} else if (*ispec == 15) {
/* ==== NSHFTS: The number of simultaneous shifts ===== */
ret_val = ns;
} else if (*ispec == 13) {
/* ==== NW: deflation window size. ==== */
if (nh <= 500) {
ret_val = ns;
} else {
ret_val = ns * 3 / 2;
}
} else if (*ispec == 16) {
/* ==== IACC22: Whether to accumulate reflections */
/* . before updating the far-from-diagonal elements */
/* . and whether to use 2-by-2 block structure while */
/* . doing it. A small amount of work could be saved */
/* . by making this choice dependent also upon the */
/* . NH=IHI-ILO+1. */
ret_val = 0;
if (ns >= 14) {
ret_val = 1;
}
if (ns >= 14) {
ret_val = 2;
}
} else {
/* ===== invalid value of ispec ===== */
ret_val = -1;
}
/* ==== End of IPARMQ ==== */
return ret_val;
} /* iparmq_ */
|
