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 #ifndef TEMPLATE_LAPACK_LAG2_HEADER
00036 #define TEMPLATE_LAPACK_LAG2_HEADER
00037
00038
00039 template<class Treal>
00040 int template_lapack_lag2(const Treal *a, const integer *lda, const Treal *b,
00041 const integer *ldb, const Treal *safmin, Treal *scale1, Treal *
00042 scale2, Treal *wr1, Treal *wr2, Treal *wi)
00043 {
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
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
00125
00126
00127
00128
00129
00130
00131
00132 integer a_dim1, a_offset, b_dim1, b_offset;
00133 Treal d__1, d__2, d__3, d__4, d__5, d__6;
00134
00135 Treal diff, bmin, wbig, wabs, wdet, r__, binv11, binv22,
00136 discr, anorm, bnorm, bsize, shift, c1, c2, c3, c4, c5, rtmin,
00137 rtmax, wsize, s1, s2, a11, a12, a21, a22, b11, b12, b22, ascale,
00138 bscale, pp, qq, ss, wscale, safmax, wsmall, as11, as12, as22, sum,
00139 abi22;
00140 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00141 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
00142
00143 a_dim1 = *lda;
00144 a_offset = 1 + a_dim1 * 1;
00145 a -= a_offset;
00146 b_dim1 = *ldb;
00147 b_offset = 1 + b_dim1 * 1;
00148 b -= b_offset;
00149
00150
00151 rtmin = template_blas_sqrt(*safmin);
00152 rtmax = 1. / rtmin;
00153 safmax = 1. / *safmin;
00154
00155
00156
00157
00158 d__5 = (d__1 = a_ref(1, 1), absMACRO(d__1)) + (d__2 = a_ref(2, 1), absMACRO(d__2)),
00159 d__6 = (d__3 = a_ref(1, 2), absMACRO(d__3)) + (d__4 = a_ref(2, 2), absMACRO(
00160 d__4)), d__5 = maxMACRO(d__5,d__6);
00161 anorm = maxMACRO(d__5,*safmin);
00162 ascale = 1. / anorm;
00163 a11 = ascale * a_ref(1, 1);
00164 a21 = ascale * a_ref(2, 1);
00165 a12 = ascale * a_ref(1, 2);
00166 a22 = ascale * a_ref(2, 2);
00167
00168
00169
00170 b11 = b_ref(1, 1);
00171 b12 = b_ref(1, 2);
00172 b22 = b_ref(2, 2);
00173
00174 d__1 = absMACRO(b11), d__2 = absMACRO(b12), d__1 = maxMACRO(d__1,d__2), d__2 = absMACRO(b22),
00175 d__1 = maxMACRO(d__1,d__2);
00176 bmin = rtmin * maxMACRO(d__1,rtmin);
00177 if (absMACRO(b11) < bmin) {
00178 b11 = template_lapack_d_sign(&bmin, &b11);
00179 }
00180 if (absMACRO(b22) < bmin) {
00181 b22 = template_lapack_d_sign(&bmin, &b22);
00182 }
00183
00184
00185
00186
00187 d__1 = absMACRO(b11), d__2 = absMACRO(b12) + absMACRO(b22), d__1 = maxMACRO(d__1,d__2);
00188 bnorm = maxMACRO(d__1,*safmin);
00189
00190 d__1 = absMACRO(b11), d__2 = absMACRO(b22);
00191 bsize = maxMACRO(d__1,d__2);
00192 bscale = 1. / bsize;
00193 b11 *= bscale;
00194 b12 *= bscale;
00195 b22 *= bscale;
00196
00197
00198
00199
00200
00201 binv11 = 1. / b11;
00202 binv22 = 1. / b22;
00203 s1 = a11 * binv11;
00204 s2 = a22 * binv22;
00205 if (absMACRO(s1) <= absMACRO(s2)) {
00206 as12 = a12 - s1 * b12;
00207 as22 = a22 - s1 * b22;
00208 ss = a21 * (binv11 * binv22);
00209 abi22 = as22 * binv22 - ss * b12;
00210 pp = abi22 * .5;
00211 shift = s1;
00212 } else {
00213 as12 = a12 - s2 * b12;
00214 as11 = a11 - s2 * b11;
00215 ss = a21 * (binv11 * binv22);
00216 abi22 = -ss * b12;
00217 pp = (as11 * binv11 + abi22) * .5;
00218 shift = s2;
00219 }
00220 qq = ss * as12;
00221 if ((d__1 = pp * rtmin, absMACRO(d__1)) >= 1.) {
00222
00223 d__1 = rtmin * pp;
00224 discr = d__1 * d__1 + qq * *safmin;
00225 r__ = template_blas_sqrt((absMACRO(discr))) * rtmax;
00226 } else {
00227
00228 d__1 = pp;
00229 if (d__1 * d__1 + absMACRO(qq) <= *safmin) {
00230
00231 d__1 = rtmax * pp;
00232 discr = d__1 * d__1 + qq * safmax;
00233 r__ = template_blas_sqrt((absMACRO(discr))) * rtmin;
00234 } else {
00235
00236 d__1 = pp;
00237 discr = d__1 * d__1 + qq;
00238 r__ = template_blas_sqrt((absMACRO(discr)));
00239 }
00240 }
00241
00242
00243
00244
00245
00246
00247
00248 if (discr >= 0. || r__ == 0.) {
00249 sum = pp + template_lapack_d_sign(&r__, &pp);
00250 diff = pp - template_lapack_d_sign(&r__, &pp);
00251 wbig = shift + sum;
00252
00253
00254
00255 wsmall = shift + diff;
00256
00257 d__1 = absMACRO(wsmall);
00258 if (absMACRO(wbig) * .5 > maxMACRO(d__1,*safmin)) {
00259 wdet = (a11 * a22 - a12 * a21) * (binv11 * binv22);
00260 wsmall = wdet / wbig;
00261 }
00262
00263
00264
00265
00266 if (pp > abi22) {
00267 *wr1 = minMACRO(wbig,wsmall);
00268 *wr2 = maxMACRO(wbig,wsmall);
00269 } else {
00270 *wr1 = maxMACRO(wbig,wsmall);
00271 *wr2 = minMACRO(wbig,wsmall);
00272 }
00273 *wi = 0.;
00274 } else {
00275
00276
00277
00278 *wr1 = shift + pp;
00279 *wr2 = *wr1;
00280 *wi = r__;
00281 }
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295 c1 = bsize * (*safmin * maxMACRO(1.,ascale));
00296 c2 = *safmin * maxMACRO(1.,bnorm);
00297 c3 = bsize * *safmin;
00298 if (ascale <= 1. && bsize <= 1.) {
00299
00300 d__1 = 1., d__2 = ascale / *safmin * bsize;
00301 c4 = minMACRO(d__1,d__2);
00302 } else {
00303 c4 = 1.;
00304 }
00305 if (ascale <= 1. || bsize <= 1.) {
00306
00307 d__1 = 1., d__2 = ascale * bsize;
00308 c5 = minMACRO(d__1,d__2);
00309 } else {
00310 c5 = 1.;
00311 }
00312
00313
00314
00315 wabs = absMACRO(*wr1) + absMACRO(*wi);
00316
00317
00318 d__3 = c4, d__4 = maxMACRO(wabs,c5) * .5;
00319 d__1 = maxMACRO(*safmin,c1), d__2 = (wabs * c2 + c3) * 1.0000100000000001,
00320 d__1 = maxMACRO(d__1,d__2), d__2 = minMACRO(d__3,d__4);
00321 wsize = maxMACRO(d__1,d__2);
00322 if (wsize != 1.) {
00323 wscale = 1. / wsize;
00324 if (wsize > 1.) {
00325 *scale1 = maxMACRO(ascale,bsize) * wscale * minMACRO(ascale,bsize);
00326 } else {
00327 *scale1 = minMACRO(ascale,bsize) * wscale * maxMACRO(ascale,bsize);
00328 }
00329 *wr1 *= wscale;
00330 if (*wi != 0.) {
00331 *wi *= wscale;
00332 *wr2 = *wr1;
00333 *scale2 = *scale1;
00334 }
00335 } else {
00336 *scale1 = ascale * bsize;
00337 *scale2 = *scale1;
00338 }
00339
00340
00341
00342 if (*wi == 0.) {
00343
00344
00345
00346 d__5 = absMACRO(*wr2);
00347 d__3 = c4, d__4 = maxMACRO(d__5,c5) * .5;
00348 d__1 = maxMACRO(*safmin,c1), d__2 = (absMACRO(*wr2) * c2 + c3) *
00349 1.0000100000000001, d__1 = maxMACRO(d__1,d__2), d__2 = minMACRO(d__3,
00350 d__4);
00351 wsize = maxMACRO(d__1,d__2);
00352 if (wsize != 1.) {
00353 wscale = 1. / wsize;
00354 if (wsize > 1.) {
00355 *scale2 = maxMACRO(ascale,bsize) * wscale * minMACRO(ascale,bsize);
00356 } else {
00357 *scale2 = minMACRO(ascale,bsize) * wscale * maxMACRO(ascale,bsize);
00358 }
00359 *wr2 *= wscale;
00360 } else {
00361 *scale2 = ascale * bsize;
00362 }
00363 }
00364
00365
00366
00367 return 0;
00368 }
00369
00370 #undef b_ref
00371 #undef a_ref
00372
00373
00374 #endif