1 /* dlasq5.f -- translated by f2c (version 20061008).
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
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
10 http://www.netlib.org/f2c/libf2c.zip
16 /* Subroutine */ int dlasq5_(integer *i0, integer *n0, doublereal *z__,
17 integer *pp, doublereal *tau, doublereal *dmin__, doublereal *dmin1,
18 doublereal *dmin2, doublereal *dn, doublereal *dnm1, doublereal *dnm2,
21 /* System generated locals */
23 doublereal d__1, d__2;
28 doublereal emin, temp;
31 /* -- LAPACK routine (version 3.2) -- */
33 /* -- Contributed by Osni Marques of the Lawrence Berkeley National -- */
34 /* -- Laboratory and Beresford Parlett of the Univ. of California at -- */
36 /* -- November 2008 -- */
38 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
39 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
41 /* .. Scalar Arguments .. */
43 /* .. Array Arguments .. */
49 /* DLASQ5 computes one dqds transform in ping-pong form, one */
50 /* version for IEEE machines another for non IEEE machines. */
55 /* I0 (input) INTEGER */
58 /* N0 (input) INTEGER */
61 /* Z (input) DOUBLE PRECISION array, dimension ( 4*N ) */
62 /* Z holds the qd array. EMIN is stored in Z(4*N0) to avoid */
63 /* an extra argument. */
65 /* PP (input) INTEGER */
66 /* PP=0 for ping, PP=1 for pong. */
68 /* TAU (input) DOUBLE PRECISION */
69 /* This is the shift. */
71 /* DMIN (output) DOUBLE PRECISION */
72 /* Minimum value of d. */
74 /* DMIN1 (output) DOUBLE PRECISION */
75 /* Minimum value of d, excluding D( N0 ). */
77 /* DMIN2 (output) DOUBLE PRECISION */
78 /* Minimum value of d, excluding D( N0 ) and D( N0-1 ). */
80 /* DN (output) DOUBLE PRECISION */
81 /* d(N0), the last value of d. */
83 /* DNM1 (output) DOUBLE PRECISION */
86 /* DNM2 (output) DOUBLE PRECISION */
89 /* IEEE (input) LOGICAL */
90 /* Flag for IEEE or non IEEE arithmetic. */
92 /* ===================================================================== */
96 /* .. Local Scalars .. */
98 /* .. Intrinsic Functions .. */
100 /* .. Executable Statements .. */
102 /* Parameter adjustments */
106 if (*n0 - *i0 - 1 <= 0) {
110 j4 = (*i0 << 2) + *pp - 3;
112 d__ = z__[j4] - *tau;
118 /* Code for IEEE arithmetic. */
122 for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
123 z__[j4 - 2] = d__ + z__[j4 - 1];
124 temp = z__[j4 + 1] / z__[j4 - 2];
125 d__ = d__ * temp - *tau;
126 *dmin__ = min(*dmin__,d__);
127 z__[j4] = z__[j4 - 1] * temp;
130 emin = min(d__1,emin);
135 for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
136 z__[j4 - 3] = d__ + z__[j4];
137 temp = z__[j4 + 2] / z__[j4 - 3];
138 d__ = d__ * temp - *tau;
139 *dmin__ = min(*dmin__,d__);
140 z__[j4 - 1] = z__[j4] * temp;
143 emin = min(d__1,emin);
148 /* Unroll last two steps. */
152 j4 = (*n0 - 2 << 2) - *pp;
153 j4p2 = j4 + (*pp << 1) - 1;
154 z__[j4 - 2] = *dnm2 + z__[j4p2];
155 z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
156 *dnm1 = z__[j4p2 + 2] * (*dnm2 / z__[j4 - 2]) - *tau;
157 *dmin__ = min(*dmin__,*dnm1);
161 j4p2 = j4 + (*pp << 1) - 1;
162 z__[j4 - 2] = *dnm1 + z__[j4p2];
163 z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
164 *dn = z__[j4p2 + 2] * (*dnm1 / z__[j4 - 2]) - *tau;
165 *dmin__ = min(*dmin__,*dn);
169 /* Code for non IEEE arithmetic. */
173 for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
174 z__[j4 - 2] = d__ + z__[j4 - 1];
178 z__[j4] = z__[j4 + 1] * (z__[j4 - 1] / z__[j4 - 2]);
179 d__ = z__[j4 + 1] * (d__ / z__[j4 - 2]) - *tau;
181 *dmin__ = min(*dmin__,d__);
183 d__1 = emin, d__2 = z__[j4];
184 emin = min(d__1,d__2);
189 for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
190 z__[j4 - 3] = d__ + z__[j4];
194 z__[j4 - 1] = z__[j4 + 2] * (z__[j4] / z__[j4 - 3]);
195 d__ = z__[j4 + 2] * (d__ / z__[j4 - 3]) - *tau;
197 *dmin__ = min(*dmin__,d__);
199 d__1 = emin, d__2 = z__[j4 - 1];
200 emin = min(d__1,d__2);
205 /* Unroll last two steps. */
209 j4 = (*n0 - 2 << 2) - *pp;
210 j4p2 = j4 + (*pp << 1) - 1;
211 z__[j4 - 2] = *dnm2 + z__[j4p2];
215 z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
216 *dnm1 = z__[j4p2 + 2] * (*dnm2 / z__[j4 - 2]) - *tau;
218 *dmin__ = min(*dmin__,*dnm1);
222 j4p2 = j4 + (*pp << 1) - 1;
223 z__[j4 - 2] = *dnm1 + z__[j4p2];
227 z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
228 *dn = z__[j4p2 + 2] * (*dnm1 / z__[j4 - 2]) - *tau;
230 *dmin__ = min(*dmin__,*dn);
235 z__[(*n0 << 2) - *pp] = emin;