1 /* ../netlib/dlaln2.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib;
2 on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */
3 #include "FLA_f2c.h" /* > \brief \b DLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form. */
4 /* =========== DOCUMENTATION =========== */
5 /* Online html documentation available at */
6 /* http://www.netlib.org/lapack/explore-html/ */
7 /* > \htmlonly */
8 /* > Download DLALN2 + dependencies */
9 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaln2. f"> */
10 /* > [TGZ]</a> */
11 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaln2. f"> */
12 /* > [ZIP]</a> */
13 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaln2. f"> */
14 /* > [TXT]</a> */
15 /* > \endhtmlonly */
16 /* Definition: */
17 /* =========== */
18 /* SUBROUTINE DLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, */
19 /* LDB, WR, WI, X, LDX, SCALE, XNORM, INFO ) */
20 /* .. Scalar Arguments .. */
21 /* LOGICAL LTRANS */
22 /* INTEGER INFO, LDA, LDB, LDX, NA, NW */
23 /* DOUBLE PRECISION CA, D1, D2, SCALE, SMIN, WI, WR, XNORM */
24 /* .. */
25 /* .. Array Arguments .. */
26 /* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), X( LDX, * ) */
27 /* .. */
28 /* > \par Purpose: */
29 /* ============= */
30 /* > */
31 /* > \verbatim */
32 /* > */
33 /* > DLALN2 solves a system of the form (ca A - w D ) X = s B */
34 /* > or (ca A**T - w D) X = s B with possible scaling ("s") and */
35 /* > perturbation of A. (A**T means A-transpose.) */
36 /* > */
37 /* > A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA */
38 /* > real diagonal matrix, w is a real or complex value, and X and B are */
39 /* > NA x 1 matrices -- real if w is real, complex if w is complex. NA */
40 /* > may be 1 or 2. */
41 /* > */
42 /* > If w is complex, X and B are represented as NA x 2 matrices, */
43 /* > the first column of each being the real part and the second */
44 /* > being the imaginary part. */
45 /* > */
46 /* > "s" is a scaling factor (.LE. 1), computed by DLALN2, which is */
47 /* > so chosen that X can be computed without overflow. X is further */
48 /* > scaled if necessary to assure that norm(ca A - w D)*norm(X) is less */
49 /* > than overflow. */
50 /* > */
51 /* > If both singular values of (ca A - w D) are less than SMIN, */
52 /* > SMIN*identity will be used instead of (ca A - w D). If only one */
53 /* > singular value is less than SMIN, one element of (ca A - w D) will be */
54 /* > perturbed enough to make the smallest singular value roughly SMIN. */
55 /* > If both singular values are at least SMIN, (ca A - w D) will not be */
56 /* > perturbed. In any case, the perturbation will be at most some small */
57 /* > multiple of max( SMIN, ulp*norm(ca A - w D) ). The singular values */
58 /* > are computed by infinity-norm approximations, and thus will only be */
59 /* > correct to a factor of 2 or so. */
60 /* > */
61 /* > Note: all input quantities are assumed to be smaller than overflow */
62 /* > by a reasonable factor. (See BIGNUM.) */
63 /* > \endverbatim */
64 /* Arguments: */
65 /* ========== */
66 /* > \param[in] LTRANS */
67 /* > \verbatim */
68 /* > LTRANS is LOGICAL */
69 /* > =.TRUE.: A-transpose will be used. */
70 /* > =.FALSE.: A will be used (not transposed.) */
71 /* > \endverbatim */
72 /* > */
73 /* > \param[in] NA */
74 /* > \verbatim */
75 /* > NA is INTEGER */
76 /* > The size of the matrix A. It may (only) be 1 or 2. */
77 /* > \endverbatim */
78 /* > */
79 /* > \param[in] NW */
80 /* > \verbatim */
81 /* > NW is INTEGER */
82 /* > 1 if "w" is real, 2 if "w" is complex. It may only be 1 */
83 /* > or 2. */
84 /* > \endverbatim */
85 /* > */
86 /* > \param[in] SMIN */
87 /* > \verbatim */
88 /* > SMIN is DOUBLE PRECISION */
89 /* > The desired lower bound on the singular values of A. This */
90 /* > should be a safe distance away from underflow or overflow, */
91 /* > say, between (underflow/machine precision) and (machine */
92 /* > precision * overflow ). (See BIGNUM and ULP.) */
93 /* > \endverbatim */
94 /* > */
95 /* > \param[in] CA */
96 /* > \verbatim */
97 /* > CA is DOUBLE PRECISION */
98 /* > The coefficient c, which A is multiplied by. */
99 /* > \endverbatim */
100 /* > */
101 /* > \param[in] A */
102 /* > \verbatim */
103 /* > A is DOUBLE PRECISION array, dimension (LDA,NA) */
104 /* > The NA x NA matrix A. */
105 /* > \endverbatim */
106 /* > */
107 /* > \param[in] LDA */
108 /* > \verbatim */
109 /* > LDA is INTEGER */
110 /* > The leading dimension of A. It must be at least NA. */
111 /* > \endverbatim */
112 /* > */
113 /* > \param[in] D1 */
114 /* > \verbatim */
115 /* > D1 is DOUBLE PRECISION */
116 /* > The 1,1 element in the diagonal matrix D. */
117 /* > \endverbatim */
118 /* > */
119 /* > \param[in] D2 */
120 /* > \verbatim */
121 /* > D2 is DOUBLE PRECISION */
122 /* > The 2,2 element in the diagonal matrix D. Not used if NW=1. */
123 /* > \endverbatim */
124 /* > */
125 /* > \param[in] B */
126 /* > \verbatim */
127 /* > B is DOUBLE PRECISION array, dimension (LDB,NW) */
128 /* > The NA x NW matrix B (right-hand side). If NW=2 ("w" is */
129 /* > complex), column 1 contains the real part of B and column 2 */
130 /* > contains the imaginary part. */
131 /* > \endverbatim */
132 /* > */
133 /* > \param[in] LDB */
134 /* > \verbatim */
135 /* > LDB is INTEGER */
136 /* > The leading dimension of B. It must be at least NA. */
137 /* > \endverbatim */
138 /* > */
139 /* > \param[in] WR */
140 /* > \verbatim */
141 /* > WR is DOUBLE PRECISION */
142 /* > The real part of the scalar "w". */
143 /* > \endverbatim */
144 /* > */
145 /* > \param[in] WI */
146 /* > \verbatim */
147 /* > WI is DOUBLE PRECISION */
148 /* > The imaginary part of the scalar "w". Not used if NW=1. */
149 /* > \endverbatim */
150 /* > */
151 /* > \param[out] X */
152 /* > \verbatim */
153 /* > X is DOUBLE PRECISION array, dimension (LDX,NW) */
154 /* > The NA x NW matrix X (unknowns), as computed by DLALN2. */
155 /* > If NW=2 ("w" is complex), on exit, column 1 will contain */
156 /* > the real part of X and column 2 will contain the imaginary */
157 /* > part. */
158 /* > \endverbatim */
159 /* > */
160 /* > \param[in] LDX */
161 /* > \verbatim */
162 /* > LDX is INTEGER */
163 /* > The leading dimension of X. It must be at least NA. */
164 /* > \endverbatim */
165 /* > */
166 /* > \param[out] SCALE */
167 /* > \verbatim */
168 /* > SCALE is DOUBLE PRECISION */
169 /* > The scale factor that B must be multiplied by to insure */
170 /* > that overflow does not occur when computing X. Thus, */
171 /* > (ca A - w D) X will be SCALE*B, not B (ignoring */
172 /* > perturbations of A.) It will be at most 1. */
173 /* > \endverbatim */
174 /* > */
175 /* > \param[out] XNORM */
176 /* > \verbatim */
177 /* > XNORM is DOUBLE PRECISION */
178 /* > The infinity-norm of X, when X is regarded as an NA x NW */
179 /* > real matrix. */
180 /* > \endverbatim */
181 /* > */
182 /* > \param[out] INFO */
183 /* > \verbatim */
184 /* > INFO is INTEGER */
185 /* > An error flag. It will be set to zero if no error occurs, */
186 /* > a negative number if an argument is in error, or a positive */
187 /* > number if ca A - w D had to be perturbed. */
188 /* > The possible values are: */
189 /* > = 0: No error occurred, and (ca A - w D) did not have to be */
190 /* > perturbed. */
191 /* > = 1: (ca A - w D) had to be perturbed to make its smallest */
192 /* > (or only) singular value greater than SMIN. */
193 /* > NOTE: In the interests of speed, this routine does not */
194 /* > check the inputs for errors. */
195 /* > \endverbatim */
196 /* Authors: */
197 /* ======== */
198 /* > \author Univ. of Tennessee */
199 /* > \author Univ. of California Berkeley */
200 /* > \author Univ. of Colorado Denver */
201 /* > \author NAG Ltd. */
202 /* > \date September 2012 */
203 /* > \ingroup doubleOTHERauxiliary */
204 /* ===================================================================== */
205 /* Subroutine */
dlaln2_(logical * ltrans,integer * na,integer * nw,doublereal * smin,doublereal * ca,doublereal * a,integer * lda,doublereal * d1,doublereal * d2,doublereal * b,integer * ldb,doublereal * wr,doublereal * wi,doublereal * x,integer * ldx,doublereal * scale,doublereal * xnorm,integer * info)206 int dlaln2_(logical *ltrans, integer *na, integer *nw, doublereal *smin, doublereal *ca, doublereal *a, integer *lda, doublereal *d1, doublereal *d2, doublereal *b, integer *ldb, doublereal *wr, doublereal *wi, doublereal *x, integer *ldx, doublereal *scale, doublereal *xnorm, integer *info)
207 {
208 /* Initialized data */
209 static logical zswap[4] =
210 {
211 FALSE_,FALSE_,TRUE_,TRUE_
212 }
213 ;
214 static logical rswap[4] =
215 {
216 FALSE_,TRUE_,FALSE_,TRUE_
217 }
218 ;
219 static integer ipivot[16] /* was [4][4] */
220 =
221 {
222 1,2,3,4,2,1,4,3,3,4,1,2, 4,3,2,1
223 }
224 ;
225 /* System generated locals */
226 integer a_dim1, a_offset, b_dim1, b_offset, x_dim1, x_offset;
227 doublereal d__1, d__2, d__3, d__4, d__5, d__6;
228 static doublereal equiv_0[4], equiv_1[4];
229 /* Local variables */
230 integer j;
231 #define ci (equiv_0)
232 #define cr (equiv_1)
233 doublereal bi1, bi2, br1, br2, xi1, xi2, xr1, xr2, ci21, ci22, cr21, cr22, li21, csi, ui11, lr21, ui12, ui22;
234 #define civ (equiv_0)
235 doublereal csr, ur11, ur12, ur22;
236 #define crv (equiv_1)
237 doublereal bbnd, cmax, ui11r, ui12s, temp, ur11r, ur12s, u22abs;
238 integer icmax;
239 doublereal bnorm, cnorm, smini;
240 extern doublereal dlamch_(char *);
241 extern /* Subroutine */
242 int dladiv_(doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *);
243 doublereal bignum, smlnum;
244 /* -- LAPACK auxiliary routine (version 3.4.2) -- */
245 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
246 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
247 /* September 2012 */
248 /* .. Scalar Arguments .. */
249 /* .. */
250 /* .. Array Arguments .. */
251 /* .. */
252 /* ===================================================================== */
253 /* .. Parameters .. */
254 /* .. */
255 /* .. Local Scalars .. */
256 /* .. */
257 /* .. Local Arrays .. */
258 /* .. */
259 /* .. External Functions .. */
260 /* .. */
261 /* .. External Subroutines .. */
262 /* .. */
263 /* .. Intrinsic Functions .. */
264 /* .. */
265 /* .. Equivalences .. */
266 /* .. */
267 /* .. Data statements .. */
268 /* Parameter adjustments */
269 a_dim1 = *lda;
270 a_offset = 1 + a_dim1;
271 a -= a_offset;
272 b_dim1 = *ldb;
273 b_offset = 1 + b_dim1;
274 b -= b_offset;
275 x_dim1 = *ldx;
276 x_offset = 1 + x_dim1;
277 x -= x_offset;
278 /* Function Body */
279 /* .. */
280 /* .. Executable Statements .. */
281 /* Compute BIGNUM */
282 smlnum = 2. * dlamch_("Safe minimum");
283 bignum = 1. / smlnum;
284 smini = max(*smin,smlnum);
285 /* Don't check for input errors */
286 *info = 0;
287 /* Standard Initializations */
288 *scale = 1.;
289 if (*na == 1)
290 {
291 /* 1 x 1 (i.e., scalar) system C X = B */
292 if (*nw == 1)
293 {
294 /* Real 1x1 system. */
295 /* C = ca A - w D */
296 csr = *ca * a[a_dim1 + 1] - *wr * *d1;
297 cnorm = f2c_abs(csr);
298 /* If | C | < SMINI, use C = SMINI */
299 if (cnorm < smini)
300 {
301 csr = smini;
302 cnorm = smini;
303 *info = 1;
304 }
305 /* Check scaling for X = B / C */
306 bnorm = (d__1 = b[b_dim1 + 1], f2c_abs(d__1));
307 if (cnorm < 1. && bnorm > 1.)
308 {
309 if (bnorm > bignum * cnorm)
310 {
311 *scale = 1. / bnorm;
312 }
313 }
314 /* Compute X */
315 x[x_dim1 + 1] = b[b_dim1 + 1] * *scale / csr;
316 *xnorm = (d__1 = x[x_dim1 + 1], f2c_abs(d__1));
317 }
318 else
319 {
320 /* Complex 1x1 system (w is complex) */
321 /* C = ca A - w D */
322 csr = *ca * a[a_dim1 + 1] - *wr * *d1;
323 csi = -(*wi) * *d1;
324 cnorm = f2c_abs(csr) + f2c_abs(csi);
325 /* If | C | < SMINI, use C = SMINI */
326 if (cnorm < smini)
327 {
328 csr = smini;
329 csi = 0.;
330 cnorm = smini;
331 *info = 1;
332 }
333 /* Check scaling for X = B / C */
334 bnorm = (d__1 = b[b_dim1 + 1], f2c_abs(d__1)) + (d__2 = b[(b_dim1 << 1) + 1], f2c_abs(d__2));
335 if (cnorm < 1. && bnorm > 1.)
336 {
337 if (bnorm > bignum * cnorm)
338 {
339 *scale = 1. / bnorm;
340 }
341 }
342 /* Compute X */
343 d__1 = *scale * b[b_dim1 + 1];
344 d__2 = *scale * b[(b_dim1 << 1) + 1];
345 dladiv_(&d__1, &d__2, &csr, &csi, &x[x_dim1 + 1], &x[(x_dim1 << 1) + 1]);
346 *xnorm = (d__1 = x[x_dim1 + 1], f2c_abs(d__1)) + (d__2 = x[(x_dim1 << 1) + 1], f2c_abs(d__2));
347 }
348 }
349 else
350 {
351 /* 2x2 System */
352 /* Compute the real part of C = ca A - w D (or ca A**T - w D ) */
353 cr[0] = *ca * a[a_dim1 + 1] - *wr * *d1;
354 cr[3] = *ca * a[(a_dim1 << 1) + 2] - *wr * *d2;
355 if (*ltrans)
356 {
357 cr[2] = *ca * a[a_dim1 + 2];
358 cr[1] = *ca * a[(a_dim1 << 1) + 1];
359 }
360 else
361 {
362 cr[1] = *ca * a[a_dim1 + 2];
363 cr[2] = *ca * a[(a_dim1 << 1) + 1];
364 }
365 if (*nw == 1)
366 {
367 /* Real 2x2 system (w is real) */
368 /* Find the largest element in C */
369 cmax = 0.;
370 icmax = 0;
371 for (j = 1;
372 j <= 4;
373 ++j)
374 {
375 if ((d__1 = crv[j - 1], f2c_abs(d__1)) > cmax)
376 {
377 cmax = (d__1 = crv[j - 1], f2c_abs(d__1));
378 icmax = j;
379 }
380 /* L10: */
381 }
382 /* If norm(C) < SMINI, use SMINI*identity. */
383 if (cmax < smini)
384 {
385 /* Computing MAX */
386 d__3 = (d__1 = b[b_dim1 + 1], f2c_abs(d__1));
387 d__4 = (d__2 = b[ b_dim1 + 2], f2c_abs(d__2)); // , expr subst
388 bnorm = max(d__3,d__4);
389 if (smini < 1. && bnorm > 1.)
390 {
391 if (bnorm > bignum * smini)
392 {
393 *scale = 1. / bnorm;
394 }
395 }
396 temp = *scale / smini;
397 x[x_dim1 + 1] = temp * b[b_dim1 + 1];
398 x[x_dim1 + 2] = temp * b[b_dim1 + 2];
399 *xnorm = temp * bnorm;
400 *info = 1;
401 return 0;
402 }
403 /* Gaussian elimination with complete pivoting. */
404 ur11 = crv[icmax - 1];
405 cr21 = crv[ipivot[(icmax << 2) - 3] - 1];
406 ur12 = crv[ipivot[(icmax << 2) - 2] - 1];
407 cr22 = crv[ipivot[(icmax << 2) - 1] - 1];
408 ur11r = 1. / ur11;
409 lr21 = ur11r * cr21;
410 ur22 = cr22 - ur12 * lr21;
411 /* If smaller pivot < SMINI, use SMINI */
412 if (f2c_abs(ur22) < smini)
413 {
414 ur22 = smini;
415 *info = 1;
416 }
417 if (rswap[icmax - 1])
418 {
419 br1 = b[b_dim1 + 2];
420 br2 = b[b_dim1 + 1];
421 }
422 else
423 {
424 br1 = b[b_dim1 + 1];
425 br2 = b[b_dim1 + 2];
426 }
427 br2 -= lr21 * br1;
428 /* Computing MAX */
429 d__2 = (d__1 = br1 * (ur22 * ur11r), f2c_abs(d__1));
430 d__3 = f2c_abs(br2); // , expr subst
431 bbnd = max(d__2,d__3);
432 if (bbnd > 1. && f2c_abs(ur22) < 1.)
433 {
434 if (bbnd >= bignum * f2c_abs(ur22))
435 {
436 *scale = 1. / bbnd;
437 }
438 }
439 xr2 = br2 * *scale / ur22;
440 xr1 = *scale * br1 * ur11r - xr2 * (ur11r * ur12);
441 if (zswap[icmax - 1])
442 {
443 x[x_dim1 + 1] = xr2;
444 x[x_dim1 + 2] = xr1;
445 }
446 else
447 {
448 x[x_dim1 + 1] = xr1;
449 x[x_dim1 + 2] = xr2;
450 }
451 /* Computing MAX */
452 d__1 = f2c_abs(xr1);
453 d__2 = f2c_abs(xr2); // , expr subst
454 *xnorm = max(d__1,d__2);
455 /* Further scaling if norm(A) norm(X) > overflow */
456 if (*xnorm > 1. && cmax > 1.)
457 {
458 if (*xnorm > bignum / cmax)
459 {
460 temp = cmax / bignum;
461 x[x_dim1 + 1] = temp * x[x_dim1 + 1];
462 x[x_dim1 + 2] = temp * x[x_dim1 + 2];
463 *xnorm = temp * *xnorm;
464 *scale = temp * *scale;
465 }
466 }
467 }
468 else
469 {
470 /* Complex 2x2 system (w is complex) */
471 /* Find the largest element in C */
472 ci[0] = -(*wi) * *d1;
473 ci[1] = 0.;
474 ci[2] = 0.;
475 ci[3] = -(*wi) * *d2;
476 cmax = 0.;
477 icmax = 0;
478 for (j = 1;
479 j <= 4;
480 ++j)
481 {
482 if ((d__1 = crv[j - 1], f2c_abs(d__1)) + (d__2 = civ[j - 1], f2c_abs( d__2)) > cmax)
483 {
484 cmax = (d__1 = crv[j - 1], f2c_abs(d__1)) + (d__2 = civ[j - 1] , f2c_abs(d__2));
485 icmax = j;
486 }
487 /* L20: */
488 }
489 /* If norm(C) < SMINI, use SMINI*identity. */
490 if (cmax < smini)
491 {
492 /* Computing MAX */
493 d__5 = (d__1 = b[b_dim1 + 1], f2c_abs(d__1)) + (d__2 = b[(b_dim1 << 1) + 1], f2c_abs(d__2));
494 d__6 = (d__3 = b[b_dim1 + 2], f2c_abs(d__3)) + (d__4 = b[(b_dim1 << 1) + 2], f2c_abs(d__4)); // , expr subst
495 bnorm = max(d__5,d__6);
496 if (smini < 1. && bnorm > 1.)
497 {
498 if (bnorm > bignum * smini)
499 {
500 *scale = 1. / bnorm;
501 }
502 }
503 temp = *scale / smini;
504 x[x_dim1 + 1] = temp * b[b_dim1 + 1];
505 x[x_dim1 + 2] = temp * b[b_dim1 + 2];
506 x[(x_dim1 << 1) + 1] = temp * b[(b_dim1 << 1) + 1];
507 x[(x_dim1 << 1) + 2] = temp * b[(b_dim1 << 1) + 2];
508 *xnorm = temp * bnorm;
509 *info = 1;
510 return 0;
511 }
512 /* Gaussian elimination with complete pivoting. */
513 ur11 = crv[icmax - 1];
514 ui11 = civ[icmax - 1];
515 cr21 = crv[ipivot[(icmax << 2) - 3] - 1];
516 ci21 = civ[ipivot[(icmax << 2) - 3] - 1];
517 ur12 = crv[ipivot[(icmax << 2) - 2] - 1];
518 ui12 = civ[ipivot[(icmax << 2) - 2] - 1];
519 cr22 = crv[ipivot[(icmax << 2) - 1] - 1];
520 ci22 = civ[ipivot[(icmax << 2) - 1] - 1];
521 if (icmax == 1 || icmax == 4)
522 {
523 /* Code when off-diagonals of pivoted C are real */
524 if (f2c_abs(ur11) > f2c_abs(ui11))
525 {
526 temp = ui11 / ur11;
527 /* Computing 2nd power */
528 d__1 = temp;
529 ur11r = 1. / (ur11 * (d__1 * d__1 + 1.));
530 ui11r = -temp * ur11r;
531 }
532 else
533 {
534 temp = ur11 / ui11;
535 /* Computing 2nd power */
536 d__1 = temp;
537 ui11r = -1. / (ui11 * (d__1 * d__1 + 1.));
538 ur11r = -temp * ui11r;
539 }
540 lr21 = cr21 * ur11r;
541 li21 = cr21 * ui11r;
542 ur12s = ur12 * ur11r;
543 ui12s = ur12 * ui11r;
544 ur22 = cr22 - ur12 * lr21;
545 ui22 = ci22 - ur12 * li21;
546 }
547 else
548 {
549 /* Code when diagonals of pivoted C are real */
550 ur11r = 1. / ur11;
551 ui11r = 0.;
552 lr21 = cr21 * ur11r;
553 li21 = ci21 * ur11r;
554 ur12s = ur12 * ur11r;
555 ui12s = ui12 * ur11r;
556 ur22 = cr22 - ur12 * lr21 + ui12 * li21;
557 ui22 = -ur12 * li21 - ui12 * lr21;
558 }
559 u22abs = f2c_abs(ur22) + f2c_abs(ui22);
560 /* If smaller pivot < SMINI, use SMINI */
561 if (u22abs < smini)
562 {
563 ur22 = smini;
564 ui22 = 0.;
565 *info = 1;
566 }
567 if (rswap[icmax - 1])
568 {
569 br2 = b[b_dim1 + 1];
570 br1 = b[b_dim1 + 2];
571 bi2 = b[(b_dim1 << 1) + 1];
572 bi1 = b[(b_dim1 << 1) + 2];
573 }
574 else
575 {
576 br1 = b[b_dim1 + 1];
577 br2 = b[b_dim1 + 2];
578 bi1 = b[(b_dim1 << 1) + 1];
579 bi2 = b[(b_dim1 << 1) + 2];
580 }
581 br2 = br2 - lr21 * br1 + li21 * bi1;
582 bi2 = bi2 - li21 * br1 - lr21 * bi1;
583 /* Computing MAX */
584 d__1 = (f2c_abs(br1) + f2c_abs(bi1)) * (u22abs * (f2c_abs(ur11r) + f2c_abs(ui11r)) );
585 d__2 = f2c_abs(br2) + f2c_abs(bi2); // , expr subst
586 bbnd = max(d__1,d__2);
587 if (bbnd > 1. && u22abs < 1.)
588 {
589 if (bbnd >= bignum * u22abs)
590 {
591 *scale = 1. / bbnd;
592 br1 = *scale * br1;
593 bi1 = *scale * bi1;
594 br2 = *scale * br2;
595 bi2 = *scale * bi2;
596 }
597 }
598 dladiv_(&br2, &bi2, &ur22, &ui22, &xr2, &xi2);
599 xr1 = ur11r * br1 - ui11r * bi1 - ur12s * xr2 + ui12s * xi2;
600 xi1 = ui11r * br1 + ur11r * bi1 - ui12s * xr2 - ur12s * xi2;
601 if (zswap[icmax - 1])
602 {
603 x[x_dim1 + 1] = xr2;
604 x[x_dim1 + 2] = xr1;
605 x[(x_dim1 << 1) + 1] = xi2;
606 x[(x_dim1 << 1) + 2] = xi1;
607 }
608 else
609 {
610 x[x_dim1 + 1] = xr1;
611 x[x_dim1 + 2] = xr2;
612 x[(x_dim1 << 1) + 1] = xi1;
613 x[(x_dim1 << 1) + 2] = xi2;
614 }
615 /* Computing MAX */
616 d__1 = f2c_abs(xr1) + f2c_abs(xi1);
617 d__2 = f2c_abs(xr2) + f2c_abs(xi2); // , expr subst
618 *xnorm = max(d__1,d__2);
619 /* Further scaling if norm(A) norm(X) > overflow */
620 if (*xnorm > 1. && cmax > 1.)
621 {
622 if (*xnorm > bignum / cmax)
623 {
624 temp = cmax / bignum;
625 x[x_dim1 + 1] = temp * x[x_dim1 + 1];
626 x[x_dim1 + 2] = temp * x[x_dim1 + 2];
627 x[(x_dim1 << 1) + 1] = temp * x[(x_dim1 << 1) + 1];
628 x[(x_dim1 << 1) + 2] = temp * x[(x_dim1 << 1) + 2];
629 *xnorm = temp * *xnorm;
630 *scale = temp * *scale;
631 }
632 }
633 }
634 }
635 return 0;
636 /* End of DLALN2 */
637 }
638 /* dlaln2_ */
639 #undef crv
640 #undef civ
641 #undef cr
642 #undef ci
643