dlasq6.c
Go to the documentation of this file.
00001 /* dlasq6.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Subroutine */ int dlasq6_(integer *i0, integer *n0, doublereal *z__, 
00017         integer *pp, doublereal *dmin__, doublereal *dmin1, doublereal *dmin2, 
00018          doublereal *dn, doublereal *dnm1, doublereal *dnm2)
00019 {
00020     /* System generated locals */
00021     integer i__1;
00022     doublereal d__1, d__2;
00023 
00024     /* Local variables */
00025     doublereal d__;
00026     integer j4, j4p2;
00027     doublereal emin, temp;
00028     extern doublereal dlamch_(char *);
00029     doublereal safmin;
00030 
00031 
00032 /*  -- LAPACK routine (version 3.2)                                    -- */
00033 
00034 /*  -- Contributed by Osni Marques of the Lawrence Berkeley National   -- */
00035 /*  -- Laboratory and Beresford Parlett of the Univ. of California at  -- */
00036 /*  -- Berkeley                                                        -- */
00037 /*  -- November 2008                                                   -- */
00038 
00039 /*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
00040 /*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
00041 
00042 /*     .. Scalar Arguments .. */
00043 /*     .. */
00044 /*     .. Array Arguments .. */
00045 /*     .. */
00046 
00047 /*  Purpose */
00048 /*  ======= */
00049 
00050 /*  DLASQ6 computes one dqd (shift equal to zero) transform in */
00051 /*  ping-pong form, with protection against underflow and overflow. */
00052 
00053 /*  Arguments */
00054 /*  ========= */
00055 
00056 /*  I0    (input) INTEGER */
00057 /*        First index. */
00058 
00059 /*  N0    (input) INTEGER */
00060 /*        Last index. */
00061 
00062 /*  Z     (input) DOUBLE PRECISION array, dimension ( 4*N ) */
00063 /*        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid */
00064 /*        an extra argument. */
00065 
00066 /*  PP    (input) INTEGER */
00067 /*        PP=0 for ping, PP=1 for pong. */
00068 
00069 /*  DMIN  (output) DOUBLE PRECISION */
00070 /*        Minimum value of d. */
00071 
00072 /*  DMIN1 (output) DOUBLE PRECISION */
00073 /*        Minimum value of d, excluding D( N0 ). */
00074 
00075 /*  DMIN2 (output) DOUBLE PRECISION */
00076 /*        Minimum value of d, excluding D( N0 ) and D( N0-1 ). */
00077 
00078 /*  DN    (output) DOUBLE PRECISION */
00079 /*        d(N0), the last value of d. */
00080 
00081 /*  DNM1  (output) DOUBLE PRECISION */
00082 /*        d(N0-1). */
00083 
00084 /*  DNM2  (output) DOUBLE PRECISION */
00085 /*        d(N0-2). */
00086 
00087 /*  ===================================================================== */
00088 
00089 /*     .. Parameter .. */
00090 /*     .. */
00091 /*     .. Local Scalars .. */
00092 /*     .. */
00093 /*     .. External Function .. */
00094 /*     .. */
00095 /*     .. Intrinsic Functions .. */
00096 /*     .. */
00097 /*     .. Executable Statements .. */
00098 
00099     /* Parameter adjustments */
00100     --z__;
00101 
00102     /* Function Body */
00103     if (*n0 - *i0 - 1 <= 0) {
00104         return 0;
00105     }
00106 
00107     safmin = dlamch_("Safe minimum");
00108     j4 = (*i0 << 2) + *pp - 3;
00109     emin = z__[j4 + 4];
00110     d__ = z__[j4];
00111     *dmin__ = d__;
00112 
00113     if (*pp == 0) {
00114         i__1 = *n0 - 3 << 2;
00115         for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
00116             z__[j4 - 2] = d__ + z__[j4 - 1];
00117             if (z__[j4 - 2] == 0.) {
00118                 z__[j4] = 0.;
00119                 d__ = z__[j4 + 1];
00120                 *dmin__ = d__;
00121                 emin = 0.;
00122             } else if (safmin * z__[j4 + 1] < z__[j4 - 2] && safmin * z__[j4 
00123                     - 2] < z__[j4 + 1]) {
00124                 temp = z__[j4 + 1] / z__[j4 - 2];
00125                 z__[j4] = z__[j4 - 1] * temp;
00126                 d__ *= temp;
00127             } else {
00128                 z__[j4] = z__[j4 + 1] * (z__[j4 - 1] / z__[j4 - 2]);
00129                 d__ = z__[j4 + 1] * (d__ / z__[j4 - 2]);
00130             }
00131             *dmin__ = min(*dmin__,d__);
00132 /* Computing MIN */
00133             d__1 = emin, d__2 = z__[j4];
00134             emin = min(d__1,d__2);
00135 /* L10: */
00136         }
00137     } else {
00138         i__1 = *n0 - 3 << 2;
00139         for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
00140             z__[j4 - 3] = d__ + z__[j4];
00141             if (z__[j4 - 3] == 0.) {
00142                 z__[j4 - 1] = 0.;
00143                 d__ = z__[j4 + 2];
00144                 *dmin__ = d__;
00145                 emin = 0.;
00146             } else if (safmin * z__[j4 + 2] < z__[j4 - 3] && safmin * z__[j4 
00147                     - 3] < z__[j4 + 2]) {
00148                 temp = z__[j4 + 2] / z__[j4 - 3];
00149                 z__[j4 - 1] = z__[j4] * temp;
00150                 d__ *= temp;
00151             } else {
00152                 z__[j4 - 1] = z__[j4 + 2] * (z__[j4] / z__[j4 - 3]);
00153                 d__ = z__[j4 + 2] * (d__ / z__[j4 - 3]);
00154             }
00155             *dmin__ = min(*dmin__,d__);
00156 /* Computing MIN */
00157             d__1 = emin, d__2 = z__[j4 - 1];
00158             emin = min(d__1,d__2);
00159 /* L20: */
00160         }
00161     }
00162 
00163 /*     Unroll last two steps. */
00164 
00165     *dnm2 = d__;
00166     *dmin2 = *dmin__;
00167     j4 = (*n0 - 2 << 2) - *pp;
00168     j4p2 = j4 + (*pp << 1) - 1;
00169     z__[j4 - 2] = *dnm2 + z__[j4p2];
00170     if (z__[j4 - 2] == 0.) {
00171         z__[j4] = 0.;
00172         *dnm1 = z__[j4p2 + 2];
00173         *dmin__ = *dnm1;
00174         emin = 0.;
00175     } else if (safmin * z__[j4p2 + 2] < z__[j4 - 2] && safmin * z__[j4 - 2] < 
00176             z__[j4p2 + 2]) {
00177         temp = z__[j4p2 + 2] / z__[j4 - 2];
00178         z__[j4] = z__[j4p2] * temp;
00179         *dnm1 = *dnm2 * temp;
00180     } else {
00181         z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
00182         *dnm1 = z__[j4p2 + 2] * (*dnm2 / z__[j4 - 2]);
00183     }
00184     *dmin__ = min(*dmin__,*dnm1);
00185 
00186     *dmin1 = *dmin__;
00187     j4 += 4;
00188     j4p2 = j4 + (*pp << 1) - 1;
00189     z__[j4 - 2] = *dnm1 + z__[j4p2];
00190     if (z__[j4 - 2] == 0.) {
00191         z__[j4] = 0.;
00192         *dn = z__[j4p2 + 2];
00193         *dmin__ = *dn;
00194         emin = 0.;
00195     } else if (safmin * z__[j4p2 + 2] < z__[j4 - 2] && safmin * z__[j4 - 2] < 
00196             z__[j4p2 + 2]) {
00197         temp = z__[j4p2 + 2] / z__[j4 - 2];
00198         z__[j4] = z__[j4p2] * temp;
00199         *dn = *dnm1 * temp;
00200     } else {
00201         z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
00202         *dn = z__[j4p2 + 2] * (*dnm1 / z__[j4 - 2]);
00203     }
00204     *dmin__ = min(*dmin__,*dn);
00205 
00206     z__[j4 + 2] = *dn;
00207     z__[(*n0 << 2) - *pp] = emin;
00208     return 0;
00209 
00210 /*     End of DLASQ6 */
00211 
00212 } /* dlasq6_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:46