00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #ifndef TEMPLATE_LAPACK_LASQ6_HEADER
00038 #define TEMPLATE_LAPACK_LASQ6_HEADER
00039
00040 template<class Treal>
00041 int template_lapack_lasq6(integer *i0, integer *n0, Treal *z__,
00042 integer *pp, Treal *dmin__, Treal *dmin1, Treal *dmin2,
00043 Treal *dn, Treal *dnm1, Treal *dnm2)
00044 {
00045
00046 integer i__1;
00047 Treal d__1, d__2;
00048
00049
00050 Treal d__;
00051 integer j4, j4p2;
00052 Treal emin, temp;
00053 Treal safmin;
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124 --z__;
00125
00126
00127 if (*n0 - *i0 - 1 <= 0) {
00128 return 0;
00129 }
00130
00131 safmin = template_lapack_lamch("Safe minimum", (Treal)0);
00132 j4 = (*i0 << 2) + *pp - 3;
00133 emin = z__[j4 + 4];
00134 d__ = z__[j4];
00135 *dmin__ = d__;
00136
00137 if (*pp == 0) {
00138 i__1 = (*n0 - 3) << 2;
00139 for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
00140 z__[j4 - 2] = d__ + z__[j4 - 1];
00141 if (z__[j4 - 2] == 0.) {
00142 z__[j4] = 0.;
00143 d__ = z__[j4 + 1];
00144 *dmin__ = d__;
00145 emin = 0.;
00146 } else if (safmin * z__[j4 + 1] < z__[j4 - 2] && safmin * z__[j4
00147 - 2] < z__[j4 + 1]) {
00148 temp = z__[j4 + 1] / z__[j4 - 2];
00149 z__[j4] = z__[j4 - 1] * temp;
00150 d__ *= temp;
00151 } else {
00152 z__[j4] = z__[j4 + 1] * (z__[j4 - 1] / z__[j4 - 2]);
00153 d__ = z__[j4 + 1] * (d__ / z__[j4 - 2]);
00154 }
00155 *dmin__ = minMACRO(*dmin__,d__);
00156
00157 d__1 = emin, d__2 = z__[j4];
00158 emin = minMACRO(d__1,d__2);
00159
00160 }
00161 } else {
00162 i__1 = ( *n0 - 3 ) << 2;
00163 for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
00164 z__[j4 - 3] = d__ + z__[j4];
00165 if (z__[j4 - 3] == 0.) {
00166 z__[j4 - 1] = 0.;
00167 d__ = z__[j4 + 2];
00168 *dmin__ = d__;
00169 emin = 0.;
00170 } else if (safmin * z__[j4 + 2] < z__[j4 - 3] && safmin * z__[j4
00171 - 3] < z__[j4 + 2]) {
00172 temp = z__[j4 + 2] / z__[j4 - 3];
00173 z__[j4 - 1] = z__[j4] * temp;
00174 d__ *= temp;
00175 } else {
00176 z__[j4 - 1] = z__[j4 + 2] * (z__[j4] / z__[j4 - 3]);
00177 d__ = z__[j4 + 2] * (d__ / z__[j4 - 3]);
00178 }
00179 *dmin__ = minMACRO(*dmin__,d__);
00180
00181 d__1 = emin, d__2 = z__[j4 - 1];
00182 emin = minMACRO(d__1,d__2);
00183
00184 }
00185 }
00186
00187
00188
00189 *dnm2 = d__;
00190 *dmin2 = *dmin__;
00191 j4 = ( ( *n0 - 2 ) << 2) - *pp;
00192 j4p2 = j4 + (*pp << 1) - 1;
00193 z__[j4 - 2] = *dnm2 + z__[j4p2];
00194 if (z__[j4 - 2] == 0.) {
00195 z__[j4] = 0.;
00196 *dnm1 = z__[j4p2 + 2];
00197 *dmin__ = *dnm1;
00198 emin = 0.;
00199 } else if (safmin * z__[j4p2 + 2] < z__[j4 - 2] && safmin * z__[j4 - 2] <
00200 z__[j4p2 + 2]) {
00201 temp = z__[j4p2 + 2] / z__[j4 - 2];
00202 z__[j4] = z__[j4p2] * temp;
00203 *dnm1 = *dnm2 * temp;
00204 } else {
00205 z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
00206 *dnm1 = z__[j4p2 + 2] * (*dnm2 / z__[j4 - 2]);
00207 }
00208 *dmin__ = minMACRO(*dmin__,*dnm1);
00209
00210 *dmin1 = *dmin__;
00211 j4 += 4;
00212 j4p2 = j4 + (*pp << 1) - 1;
00213 z__[j4 - 2] = *dnm1 + z__[j4p2];
00214 if (z__[j4 - 2] == 0.) {
00215 z__[j4] = 0.;
00216 *dn = z__[j4p2 + 2];
00217 *dmin__ = *dn;
00218 emin = 0.;
00219 } else if (safmin * z__[j4p2 + 2] < z__[j4 - 2] && safmin * z__[j4 - 2] <
00220 z__[j4p2 + 2]) {
00221 temp = z__[j4p2 + 2] / z__[j4 - 2];
00222 z__[j4] = z__[j4p2] * temp;
00223 *dn = *dnm1 * temp;
00224 } else {
00225 z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
00226 *dn = z__[j4p2 + 2] * (*dnm1 / z__[j4 - 2]);
00227 }
00228 *dmin__ = minMACRO(*dmin__,*dn);
00229
00230 z__[j4 + 2] = *dn;
00231 z__[(*n0 << 2) - *pp] = emin;
00232 return 0;
00233
00234
00235
00236 }
00237
00238 #endif