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_LALN2_HEADER
00038 #define TEMPLATE_LAPACK_LALN2_HEADER
00039
00040
00041 template<class Treal>
00042 int template_lapack_laln2(const logical *ltrans, const integer *na, const integer *nw,
00043 const Treal *smin, const Treal *ca, const Treal *a, const integer *lda,
00044 const Treal *d1, const Treal *d2, const Treal *b, const integer *ldb,
00045 const Treal *wr, const Treal *wi, Treal *x, const integer *ldx,
00046 Treal *scale, Treal *xnorm, integer *info)
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
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172 logical zswap[4] = { FALSE_,FALSE_,TRUE_,TRUE_ };
00173 logical rswap[4] = { FALSE_,TRUE_,FALSE_,TRUE_ };
00174 integer ipivot[16] = { 1,2,3,4,2,1,4,3,3,4,1,2,
00175 4,3,2,1 };
00176
00177 integer a_dim1, a_offset, b_dim1, b_offset, x_dim1, x_offset;
00178 Treal d__1, d__2, d__3, d__4, d__5, d__6;
00179 Treal equiv_0[4], equiv_1[4];
00180
00181 Treal bbnd, cmax, ui11r, ui12s, temp, ur11r, ur12s;
00182 integer j;
00183 Treal u22abs;
00184 integer icmax;
00185 Treal bnorm, cnorm, smini;
00186 #define ci (equiv_0)
00187 #define cr (equiv_1)
00188 Treal bignum, bi1, bi2, br1, br2, smlnum, xi1, xi2, xr1, xr2,
00189 ci21, ci22, cr21, cr22, li21, csi, ui11, lr21, ui12, ui22;
00190 #define civ (equiv_0)
00191 Treal csr, ur11, ur12, ur22;
00192 #define crv (equiv_1)
00193 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00194 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
00195 #define x_ref(a_1,a_2) x[(a_2)*x_dim1 + a_1]
00196 #define ci_ref(a_1,a_2) ci[(a_2)*2 + a_1 - 3]
00197 #define cr_ref(a_1,a_2) cr[(a_2)*2 + a_1 - 3]
00198 #define ipivot_ref(a_1,a_2) ipivot[(a_2)*4 + a_1 - 5]
00199
00200 a_dim1 = *lda;
00201 a_offset = 1 + a_dim1 * 1;
00202 a -= a_offset;
00203 b_dim1 = *ldb;
00204 b_offset = 1 + b_dim1 * 1;
00205 b -= b_offset;
00206 x_dim1 = *ldx;
00207 x_offset = 1 + x_dim1 * 1;
00208 x -= x_offset;
00209
00210
00211
00212
00213
00214 smlnum = 2. * template_lapack_lamch("Safe minimum", (Treal)0);
00215 bignum = 1. / smlnum;
00216 smini = maxMACRO(*smin,smlnum);
00217
00218
00219
00220 *info = 0;
00221
00222
00223
00224 *scale = 1.;
00225
00226 if (*na == 1) {
00227
00228
00229
00230 if (*nw == 1) {
00231
00232
00233
00234
00235
00236 csr = *ca * a_ref(1, 1) - *wr * *d1;
00237 cnorm = absMACRO(csr);
00238
00239
00240
00241 if (cnorm < smini) {
00242 csr = smini;
00243 cnorm = smini;
00244 *info = 1;
00245 }
00246
00247
00248
00249 bnorm = (d__1 = b_ref(1, 1), absMACRO(d__1));
00250 if (cnorm < 1. && bnorm > 1.) {
00251 if (bnorm > bignum * cnorm) {
00252 *scale = 1. / bnorm;
00253 }
00254 }
00255
00256
00257
00258 x_ref(1, 1) = b_ref(1, 1) * *scale / csr;
00259 *xnorm = (d__1 = x_ref(1, 1), absMACRO(d__1));
00260 } else {
00261
00262
00263
00264
00265
00266 csr = *ca * a_ref(1, 1) - *wr * *d1;
00267 csi = -(*wi) * *d1;
00268 cnorm = absMACRO(csr) + absMACRO(csi);
00269
00270
00271
00272 if (cnorm < smini) {
00273 csr = smini;
00274 csi = 0.;
00275 cnorm = smini;
00276 *info = 1;
00277 }
00278
00279
00280
00281 bnorm = (d__1 = b_ref(1, 1), absMACRO(d__1)) + (d__2 = b_ref(1, 2),
00282 absMACRO(d__2));
00283 if (cnorm < 1. && bnorm > 1.) {
00284 if (bnorm > bignum * cnorm) {
00285 *scale = 1. / bnorm;
00286 }
00287 }
00288
00289
00290
00291 d__1 = *scale * b_ref(1, 1);
00292 d__2 = *scale * b_ref(1, 2);
00293 template_lapack_ladiv(&d__1, &d__2, &csr, &csi, &x_ref(1, 1), &x_ref(1, 2));
00294 *xnorm = (d__1 = x_ref(1, 1), absMACRO(d__1)) + (d__2 = x_ref(1, 2),
00295 absMACRO(d__2));
00296 }
00297
00298 } else {
00299
00300
00301
00302
00303
00304 cr_ref(1, 1) = *ca * a_ref(1, 1) - *wr * *d1;
00305 cr_ref(2, 2) = *ca * a_ref(2, 2) - *wr * *d2;
00306 if (*ltrans) {
00307 cr_ref(1, 2) = *ca * a_ref(2, 1);
00308 cr_ref(2, 1) = *ca * a_ref(1, 2);
00309 } else {
00310 cr_ref(2, 1) = *ca * a_ref(2, 1);
00311 cr_ref(1, 2) = *ca * a_ref(1, 2);
00312 }
00313
00314 if (*nw == 1) {
00315
00316
00317
00318
00319
00320 cmax = 0.;
00321 icmax = 0;
00322
00323 for (j = 1; j <= 4; ++j) {
00324 if ((d__1 = crv[j - 1], absMACRO(d__1)) > cmax) {
00325 cmax = (d__1 = crv[j - 1], absMACRO(d__1));
00326 icmax = j;
00327 }
00328
00329 }
00330
00331
00332
00333 if (cmax < smini) {
00334
00335 d__3 = (d__1 = b_ref(1, 1), absMACRO(d__1)), d__4 = (d__2 = b_ref(
00336 2, 1), absMACRO(d__2));
00337 bnorm = maxMACRO(d__3,d__4);
00338 if (smini < 1. && bnorm > 1.) {
00339 if (bnorm > bignum * smini) {
00340 *scale = 1. / bnorm;
00341 }
00342 }
00343 temp = *scale / smini;
00344 x_ref(1, 1) = temp * b_ref(1, 1);
00345 x_ref(2, 1) = temp * b_ref(2, 1);
00346 *xnorm = temp * bnorm;
00347 *info = 1;
00348 return 0;
00349 }
00350
00351
00352
00353 ur11 = crv[icmax - 1];
00354 cr21 = crv[ipivot_ref(2, icmax) - 1];
00355 ur12 = crv[ipivot_ref(3, icmax) - 1];
00356 cr22 = crv[ipivot_ref(4, icmax) - 1];
00357 ur11r = 1. / ur11;
00358 lr21 = ur11r * cr21;
00359 ur22 = cr22 - ur12 * lr21;
00360
00361
00362
00363 if (absMACRO(ur22) < smini) {
00364 ur22 = smini;
00365 *info = 1;
00366 }
00367 if (rswap[icmax - 1]) {
00368 br1 = b_ref(2, 1);
00369 br2 = b_ref(1, 1);
00370 } else {
00371 br1 = b_ref(1, 1);
00372 br2 = b_ref(2, 1);
00373 }
00374 br2 -= lr21 * br1;
00375
00376 d__2 = (d__1 = br1 * (ur22 * ur11r), absMACRO(d__1)), d__3 = absMACRO(br2);
00377 bbnd = maxMACRO(d__2,d__3);
00378 if (bbnd > 1. && absMACRO(ur22) < 1.) {
00379 if (bbnd >= bignum * absMACRO(ur22)) {
00380 *scale = 1. / bbnd;
00381 }
00382 }
00383
00384 xr2 = br2 * *scale / ur22;
00385 xr1 = *scale * br1 * ur11r - xr2 * (ur11r * ur12);
00386 if (zswap[icmax - 1]) {
00387 x_ref(1, 1) = xr2;
00388 x_ref(2, 1) = xr1;
00389 } else {
00390 x_ref(1, 1) = xr1;
00391 x_ref(2, 1) = xr2;
00392 }
00393
00394 d__1 = absMACRO(xr1), d__2 = absMACRO(xr2);
00395 *xnorm = maxMACRO(d__1,d__2);
00396
00397
00398
00399 if (*xnorm > 1. && cmax > 1.) {
00400 if (*xnorm > bignum / cmax) {
00401 temp = cmax / bignum;
00402 x_ref(1, 1) = temp * x_ref(1, 1);
00403 x_ref(2, 1) = temp * x_ref(2, 1);
00404 *xnorm = temp * *xnorm;
00405 *scale = temp * *scale;
00406 }
00407 }
00408 } else {
00409
00410
00411
00412
00413
00414 ci_ref(1, 1) = -(*wi) * *d1;
00415 ci_ref(2, 1) = 0.;
00416 ci_ref(1, 2) = 0.;
00417 ci_ref(2, 2) = -(*wi) * *d2;
00418 cmax = 0.;
00419 icmax = 0;
00420
00421 for (j = 1; j <= 4; ++j) {
00422 if ((d__1 = crv[j - 1], absMACRO(d__1)) + (d__2 = civ[j - 1], absMACRO(
00423 d__2)) > cmax) {
00424 cmax = (d__1 = crv[j - 1], absMACRO(d__1)) + (d__2 = civ[j - 1]
00425 , absMACRO(d__2));
00426 icmax = j;
00427 }
00428
00429 }
00430
00431
00432
00433 if (cmax < smini) {
00434
00435 d__5 = (d__1 = b_ref(1, 1), absMACRO(d__1)) + (d__2 = b_ref(1, 2),
00436 absMACRO(d__2)), d__6 = (d__3 = b_ref(2, 1), absMACRO(d__3)) + (
00437 d__4 = b_ref(2, 2), absMACRO(d__4));
00438 bnorm = maxMACRO(d__5,d__6);
00439 if (smini < 1. && bnorm > 1.) {
00440 if (bnorm > bignum * smini) {
00441 *scale = 1. / bnorm;
00442 }
00443 }
00444 temp = *scale / smini;
00445 x_ref(1, 1) = temp * b_ref(1, 1);
00446 x_ref(2, 1) = temp * b_ref(2, 1);
00447 x_ref(1, 2) = temp * b_ref(1, 2);
00448 x_ref(2, 2) = temp * b_ref(2, 2);
00449 *xnorm = temp * bnorm;
00450 *info = 1;
00451 return 0;
00452 }
00453
00454
00455
00456 ur11 = crv[icmax - 1];
00457 ui11 = civ[icmax - 1];
00458 cr21 = crv[ipivot_ref(2, icmax) - 1];
00459 ci21 = civ[ipivot_ref(2, icmax) - 1];
00460 ur12 = crv[ipivot_ref(3, icmax) - 1];
00461 ui12 = civ[ipivot_ref(3, icmax) - 1];
00462 cr22 = crv[ipivot_ref(4, icmax) - 1];
00463 ci22 = civ[ipivot_ref(4, icmax) - 1];
00464 if (icmax == 1 || icmax == 4) {
00465
00466
00467
00468 if (absMACRO(ur11) > absMACRO(ui11)) {
00469 temp = ui11 / ur11;
00470
00471 d__1 = temp;
00472 ur11r = 1. / (ur11 * (d__1 * d__1 + 1.));
00473 ui11r = -temp * ur11r;
00474 } else {
00475 temp = ur11 / ui11;
00476
00477 d__1 = temp;
00478 ui11r = -1. / (ui11 * (d__1 * d__1 + 1.));
00479 ur11r = -temp * ui11r;
00480 }
00481 lr21 = cr21 * ur11r;
00482 li21 = cr21 * ui11r;
00483 ur12s = ur12 * ur11r;
00484 ui12s = ur12 * ui11r;
00485 ur22 = cr22 - ur12 * lr21;
00486 ui22 = ci22 - ur12 * li21;
00487 } else {
00488
00489
00490
00491 ur11r = 1. / ur11;
00492 ui11r = 0.;
00493 lr21 = cr21 * ur11r;
00494 li21 = ci21 * ur11r;
00495 ur12s = ur12 * ur11r;
00496 ui12s = ui12 * ur11r;
00497 ur22 = cr22 - ur12 * lr21 + ui12 * li21;
00498 ui22 = -ur12 * li21 - ui12 * lr21;
00499 }
00500 u22abs = absMACRO(ur22) + absMACRO(ui22);
00501
00502
00503
00504 if (u22abs < smini) {
00505 ur22 = smini;
00506 ui22 = 0.;
00507 *info = 1;
00508 }
00509 if (rswap[icmax - 1]) {
00510 br2 = b_ref(1, 1);
00511 br1 = b_ref(2, 1);
00512 bi2 = b_ref(1, 2);
00513 bi1 = b_ref(2, 2);
00514 } else {
00515 br1 = b_ref(1, 1);
00516 br2 = b_ref(2, 1);
00517 bi1 = b_ref(1, 2);
00518 bi2 = b_ref(2, 2);
00519 }
00520 br2 = br2 - lr21 * br1 + li21 * bi1;
00521 bi2 = bi2 - li21 * br1 - lr21 * bi1;
00522
00523 d__1 = (absMACRO(br1) + absMACRO(bi1)) * (u22abs * (absMACRO(ur11r) + absMACRO(ui11r))
00524 ), d__2 = absMACRO(br2) + absMACRO(bi2);
00525 bbnd = maxMACRO(d__1,d__2);
00526 if (bbnd > 1. && u22abs < 1.) {
00527 if (bbnd >= bignum * u22abs) {
00528 *scale = 1. / bbnd;
00529 br1 = *scale * br1;
00530 bi1 = *scale * bi1;
00531 br2 = *scale * br2;
00532 bi2 = *scale * bi2;
00533 }
00534 }
00535
00536 template_lapack_ladiv(&br2, &bi2, &ur22, &ui22, &xr2, &xi2);
00537 xr1 = ur11r * br1 - ui11r * bi1 - ur12s * xr2 + ui12s * xi2;
00538 xi1 = ui11r * br1 + ur11r * bi1 - ui12s * xr2 - ur12s * xi2;
00539 if (zswap[icmax - 1]) {
00540 x_ref(1, 1) = xr2;
00541 x_ref(2, 1) = xr1;
00542 x_ref(1, 2) = xi2;
00543 x_ref(2, 2) = xi1;
00544 } else {
00545 x_ref(1, 1) = xr1;
00546 x_ref(2, 1) = xr2;
00547 x_ref(1, 2) = xi1;
00548 x_ref(2, 2) = xi2;
00549 }
00550
00551 d__1 = absMACRO(xr1) + absMACRO(xi1), d__2 = absMACRO(xr2) + absMACRO(xi2);
00552 *xnorm = maxMACRO(d__1,d__2);
00553
00554
00555
00556 if (*xnorm > 1. && cmax > 1.) {
00557 if (*xnorm > bignum / cmax) {
00558 temp = cmax / bignum;
00559 x_ref(1, 1) = temp * x_ref(1, 1);
00560 x_ref(2, 1) = temp * x_ref(2, 1);
00561 x_ref(1, 2) = temp * x_ref(1, 2);
00562 x_ref(2, 2) = temp * x_ref(2, 2);
00563 *xnorm = temp * *xnorm;
00564 *scale = temp * *scale;
00565 }
00566 }
00567 }
00568 }
00569
00570 return 0;
00571
00572
00573
00574 }
00575
00576 #undef ipivot_ref
00577 #undef cr_ref
00578 #undef ci_ref
00579 #undef x_ref
00580 #undef b_ref
00581 #undef a_ref
00582 #undef crv
00583 #undef civ
00584 #undef cr
00585 #undef ci
00586
00587
00588 #endif