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