#include "blaswrap.h" /* -- translated by f2c (version 19990503). You must link the resulting object file with the libraries: -lf2c -lm (in that order) */ #include "f2c.h" /* Common Block Declarations */ struct { doublereal ops, itcnt; } latime_; #define latime_1 latime_ /* Subroutine */ int dlasq5_(integer *i0, integer *n0, doublereal *z__, integer *pp, doublereal *tau, doublereal *dmin__, doublereal *dmin1, doublereal *dmin2, doublereal *dn, doublereal *dnm1, doublereal *dnm2, logical *ieee) { /* System generated locals */ integer i__1; doublereal d__1, d__2; /* Local variables */ static doublereal emin, temp, d__; static integer j4, j4p2; /* -- LAPACK auxiliary routine (instrumented to count ops, version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University May 17, 2000 Purpose ======= DLASQ5 computes one dqds transform in ping-pong form, one version for IEEE machines another for non IEEE machines. Arguments ========= I0 (input) INTEGER First index. N0 (input) INTEGER Last index. Z (input) DOUBLE PRECISION array, dimension ( 4*N ) Z holds the qd array. EMIN is stored in Z(4*N0) to avoid an extra argument. PP (input) INTEGER PP=0 for ping, PP=1 for pong. TAU (input) DOUBLE PRECISION This is the shift. DMIN (output) DOUBLE PRECISION Minimum value of d. DMIN1 (output) DOUBLE PRECISION Minimum value of d, excluding D( N0 ). DMIN2 (output) DOUBLE PRECISION Minimum value of d, excluding D( N0 ) and D( N0-1 ). DN (output) DOUBLE PRECISION d(N0), the last value of d. DNM1 (output) DOUBLE PRECISION d(N0-1). DNM2 (output) DOUBLE PRECISION d(N0-2). IEEE (input) LOGICAL Flag for IEEE or non IEEE arithmetic. ===================================================================== Parameter adjustments */ --z__; /* Function Body */ if (*n0 - *i0 - 1 <= 0) { return 0; } latime_1.ops += 1.; j4 = (*i0 << 2) + *pp - 3; emin = z__[j4 + 4]; d__ = z__[j4] - *tau; *dmin__ = d__; *dmin1 = -z__[j4]; if (*ieee) { /* Code for IEEE arithmetic. */ if (*pp == 0) { i__1 = *n0 - 3 << 2; for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) { latime_1.ops += 5.; z__[j4 - 2] = d__ + z__[j4 - 1]; temp = z__[j4 + 1] / z__[j4 - 2]; d__ = d__ * temp - *tau; *dmin__ = min(*dmin__,d__); z__[j4] = z__[j4 - 1] * temp; /* Computing MIN */ d__1 = z__[j4]; emin = min(d__1,emin); /* L10: */ } } else { i__1 = *n0 - 3 << 2; for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) { latime_1.ops += 5.; z__[j4 - 3] = d__ + z__[j4]; temp = z__[j4 + 2] / z__[j4 - 3]; d__ = d__ * temp - *tau; *dmin__ = min(*dmin__,d__); z__[j4 - 1] = z__[j4] * temp; /* Computing MIN */ d__1 = z__[j4 - 1]; emin = min(d__1,emin); /* L20: */ } } /* Unroll last two steps. */ latime_1.ops += 6.; *dnm2 = d__; *dmin2 = *dmin__; j4 = (*n0 - 2 << 2) - *pp; j4p2 = j4 + (*pp << 1) - 1; z__[j4 - 2] = *dnm2 + z__[j4p2]; z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]); *dnm1 = z__[j4p2 + 2] * (*dnm2 / z__[j4 - 2]) - *tau; *dmin__ = min(*dmin__,*dnm1); latime_1.ops += 6.; *dmin1 = *dmin__; j4 += 4; j4p2 = j4 + (*pp << 1) - 1; z__[j4 - 2] = *dnm1 + z__[j4p2]; z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]); *dn = z__[j4p2 + 2] * (*dnm1 / z__[j4 - 2]) - *tau; *dmin__ = min(*dmin__,*dn); } else { /* Code for non IEEE arithmetic. */ if (*pp == 0) { i__1 = *n0 - 3 << 2; for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) { z__[j4 - 2] = d__ + z__[j4 - 1]; if (d__ < 0.) { return 0; } else { latime_1.ops += 5.; z__[j4] = z__[j4 + 1] * (z__[j4 - 1] / z__[j4 - 2]); d__ = z__[j4 + 1] * (d__ / z__[j4 - 2]) - *tau; } *dmin__ = min(*dmin__,d__); /* Computing MIN */ d__1 = emin, d__2 = z__[j4]; emin = min(d__1,d__2); /* L30: */ } } else { i__1 = *n0 - 3 << 2; for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) { z__[j4 - 3] = d__ + z__[j4]; if (d__ < 0.) { return 0; } else { latime_1.ops += 5.; z__[j4 - 1] = z__[j4 + 2] * (z__[j4] / z__[j4 - 3]); d__ = z__[j4 + 2] * (d__ / z__[j4 - 3]) - *tau; } *dmin__ = min(*dmin__,d__); /* Computing MIN */ d__1 = emin, d__2 = z__[j4 - 1]; emin = min(d__1,d__2); /* L40: */ } } /* Unroll last two steps. */ latime_1.ops += 1.; *dnm2 = d__; *dmin2 = *dmin__; j4 = (*n0 - 2 << 2) - *pp; j4p2 = j4 + (*pp << 1) - 1; z__[j4 - 2] = *dnm2 + z__[j4p2]; if (*dnm2 < 0.) { return 0; } else { latime_1.ops += 5.; z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]); *dnm1 = z__[j4p2 + 2] * (*dnm2 / z__[j4 - 2]) - *tau; } *dmin__ = min(*dmin__,*dnm1); latime_1.ops += 1.; *dmin1 = *dmin__; j4 += 4; j4p2 = j4 + (*pp << 1) - 1; z__[j4 - 2] = *dnm1 + z__[j4p2]; if (*dnm1 < 0.) { return 0; } else { latime_1.ops += 5.; z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]); *dn = z__[j4p2 + 2] * (*dnm1 / z__[j4 - 2]) - *tau; } *dmin__ = min(*dmin__,*dn); } z__[j4 + 2] = *dn; z__[(*n0 << 2) - *pp] = emin; return 0; /* End of DLASQ5 */ } /* dlasq5_ */ .