1 /* ./src_f77/dlaqtr.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static integer c__1 = 1;
11 static logical c_false = FALSE_;
12 static integer c__2 = 2;
13 static doublereal c_b21 = 1.;
14 static doublereal c_b25 = 0.;
15 static logical c_true = TRUE_;
16
dlaqtr_(logical * ltran,logical * lreal,integer * n,doublereal * t,integer * ldt,doublereal * b,doublereal * w,doublereal * scale,doublereal * x,doublereal * work,integer * info)17 /* Subroutine */ int dlaqtr_(logical *ltran, logical *lreal, integer *n,
18 doublereal *t, integer *ldt, doublereal *b, doublereal *w, doublereal
19 *scale, doublereal *x, doublereal *work, integer *info)
20 {
21 /* System generated locals */
22 integer t_dim1, t_offset, i__1, i__2;
23 doublereal d__1, d__2, d__3, d__4, d__5, d__6;
24
25 /* Local variables */
26 static doublereal d__[4] /* was [2][2] */;
27 static integer i__, j, k;
28 static doublereal v[4] /* was [2][2] */, z__;
29 static integer j1, j2, n1, n2;
30 static doublereal si, xj, sr, rec, eps, tjj, tmp;
31 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
32 integer *);
33 static integer ierr;
34 static doublereal smin, xmax;
35 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
36 integer *);
37 extern doublereal dasum_(integer *, doublereal *, integer *);
38 extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *,
39 integer *, doublereal *, integer *);
40 static integer jnext;
41 static doublereal sminw, xnorm;
42 extern /* Subroutine */ int dlaln2_(logical *, integer *, integer *,
43 doublereal *, doublereal *, doublereal *, integer *, doublereal *,
44 doublereal *, doublereal *, integer *, doublereal *, doublereal *
45 , doublereal *, integer *, doublereal *, doublereal *, integer *);
46 extern doublereal dlamch_(char *, ftnlen), dlange_(char *, integer *,
47 integer *, doublereal *, integer *, doublereal *, ftnlen);
48 extern integer idamax_(integer *, doublereal *, integer *);
49 static doublereal scaloc;
50 extern /* Subroutine */ int dladiv_(doublereal *, doublereal *,
51 doublereal *, doublereal *, doublereal *, doublereal *);
52 static doublereal bignum;
53 static logical notran;
54 static doublereal smlnum;
55
56
57 /* -- LAPACK auxiliary routine (version 3.0) -- */
58 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
59 /* Courant Institute, Argonne National Lab, and Rice University */
60 /* June 30, 1999 */
61
62 /* .. Scalar Arguments .. */
63 /* .. */
64 /* .. Array Arguments .. */
65 /* .. */
66
67 /* Purpose */
68 /* ======= */
69
70 /* DLAQTR solves the real quasi-triangular system */
71
72 /* op(T)*p = scale*c, if LREAL = .TRUE. */
73
74 /* or the complex quasi-triangular systems */
75
76 /* op(T + iB)*(p+iq) = scale*(c+id), if LREAL = .FALSE. */
77
78 /* in real arithmetic, where T is upper quasi-triangular. */
79 /* If LREAL = .FALSE., then the first diagonal block of T must be */
80 /* 1 by 1, B is the specially structured matrix */
81
82 /* B = [ b(1) b(2) ... b(n) ] */
83 /* [ w ] */
84 /* [ w ] */
85 /* [ . ] */
86 /* [ w ] */
87
88 /* op(A) = A or A', A' denotes the conjugate transpose of */
89 /* matrix A. */
90
91 /* On input, X = [ c ]. On output, X = [ p ]. */
92 /* [ d ] [ q ] */
93
94 /* This subroutine is designed for the condition number estimation */
95 /* in routine DTRSNA. */
96
97 /* Arguments */
98 /* ========= */
99
100 /* LTRAN (input) LOGICAL */
101 /* On entry, LTRAN specifies the option of conjugate transpose: */
102 /* = .FALSE., op(T+i*B) = T+i*B, */
103 /* = .TRUE., op(T+i*B) = (T+i*B)'. */
104
105 /* LREAL (input) LOGICAL */
106 /* On entry, LREAL specifies the input matrix structure: */
107 /* = .FALSE., the input is complex */
108 /* = .TRUE., the input is real */
109
110 /* N (input) INTEGER */
111 /* On entry, N specifies the order of T+i*B. N >= 0. */
112
113 /* T (input) DOUBLE PRECISION array, dimension (LDT,N) */
114 /* On entry, T contains a matrix in Schur canonical form. */
115 /* If LREAL = .FALSE., then the first diagonal block of T mu */
116 /* be 1 by 1. */
117
118 /* LDT (input) INTEGER */
119 /* The leading dimension of the matrix T. LDT >= max(1,N). */
120
121 /* B (input) DOUBLE PRECISION array, dimension (N) */
122 /* On entry, B contains the elements to form the matrix */
123 /* B as described above. */
124 /* If LREAL = .TRUE., B is not referenced. */
125
126 /* W (input) DOUBLE PRECISION */
127 /* On entry, W is the diagonal element of the matrix B. */
128 /* If LREAL = .TRUE., W is not referenced. */
129
130 /* SCALE (output) DOUBLE PRECISION */
131 /* On exit, SCALE is the scale factor. */
132
133 /* X (input/output) DOUBLE PRECISION array, dimension (2*N) */
134 /* On entry, X contains the right hand side of the system. */
135 /* On exit, X is overwritten by the solution. */
136
137 /* WORK (workspace) DOUBLE PRECISION array, dimension (N) */
138
139 /* INFO (output) INTEGER */
140 /* On exit, INFO is set to */
141 /* 0: successful exit. */
142 /* 1: the some diagonal 1 by 1 block has been perturbed by */
143 /* a small number SMIN to keep nonsingularity. */
144 /* 2: the some diagonal 2 by 2 block has been perturbed by */
145 /* a small number in DLALN2 to keep nonsingularity. */
146 /* NOTE: In the interests of speed, this routine does not */
147 /* check the inputs for errors. */
148
149 /* ===================================================================== */
150
151 /* .. Parameters .. */
152 /* .. */
153 /* .. Local Scalars .. */
154 /* .. */
155 /* .. Local Arrays .. */
156 /* .. */
157 /* .. External Functions .. */
158 /* .. */
159 /* .. External Subroutines .. */
160 /* .. */
161 /* .. Intrinsic Functions .. */
162 /* .. */
163 /* .. Executable Statements .. */
164
165 /* Do not test the input parameters for errors */
166
167 /* Parameter adjustments */
168 t_dim1 = *ldt;
169 t_offset = 1 + t_dim1;
170 t -= t_offset;
171 --b;
172 --x;
173 --work;
174
175 /* Function Body */
176 notran = ! (*ltran);
177 *info = 0;
178
179 /* Quick return if possible */
180
181 if (*n == 0) {
182 return 0;
183 }
184
185 /* Set constants to control overflow */
186
187 eps = dlamch_("P", (ftnlen)1);
188 smlnum = dlamch_("S", (ftnlen)1) / eps;
189 bignum = 1. / smlnum;
190
191 xnorm = dlange_("M", n, n, &t[t_offset], ldt, d__, (ftnlen)1);
192 if (! (*lreal)) {
193 /* Computing MAX */
194 d__1 = xnorm, d__2 = abs(*w), d__1 = max(d__1,d__2), d__2 = dlange_(
195 "M", n, &c__1, &b[1], n, d__, (ftnlen)1);
196 xnorm = max(d__1,d__2);
197 }
198 /* Computing MAX */
199 d__1 = smlnum, d__2 = eps * xnorm;
200 smin = max(d__1,d__2);
201
202 /* Compute 1-norm of each column of strictly upper triangular */
203 /* part of T to control overflow in triangular solver. */
204
205 work[1] = 0.;
206 i__1 = *n;
207 for (j = 2; j <= i__1; ++j) {
208 i__2 = j - 1;
209 work[j] = dasum_(&i__2, &t[j * t_dim1 + 1], &c__1);
210 /* L10: */
211 }
212
213 if (! (*lreal)) {
214 i__1 = *n;
215 for (i__ = 2; i__ <= i__1; ++i__) {
216 work[i__] += (d__1 = b[i__], abs(d__1));
217 /* L20: */
218 }
219 }
220
221 n2 = *n << 1;
222 n1 = *n;
223 if (! (*lreal)) {
224 n1 = n2;
225 }
226 k = idamax_(&n1, &x[1], &c__1);
227 xmax = (d__1 = x[k], abs(d__1));
228 *scale = 1.;
229
230 if (xmax > bignum) {
231 *scale = bignum / xmax;
232 dscal_(&n1, scale, &x[1], &c__1);
233 xmax = bignum;
234 }
235
236 if (*lreal) {
237
238 if (notran) {
239
240 /* Solve T*p = scale*c */
241
242 jnext = *n;
243 for (j = *n; j >= 1; --j) {
244 if (j > jnext) {
245 goto L30;
246 }
247 j1 = j;
248 j2 = j;
249 jnext = j - 1;
250 if (j > 1) {
251 if (t[j + (j - 1) * t_dim1] != 0.) {
252 j1 = j - 1;
253 jnext = j - 2;
254 }
255 }
256
257 if (j1 == j2) {
258
259 /* Meet 1 by 1 diagonal block */
260
261 /* Scale to avoid overflow when computing */
262 /* x(j) = b(j)/T(j,j) */
263
264 xj = (d__1 = x[j1], abs(d__1));
265 tjj = (d__1 = t[j1 + j1 * t_dim1], abs(d__1));
266 tmp = t[j1 + j1 * t_dim1];
267 if (tjj < smin) {
268 tmp = smin;
269 tjj = smin;
270 *info = 1;
271 }
272
273 if (xj == 0.) {
274 goto L30;
275 }
276
277 if (tjj < 1.) {
278 if (xj > bignum * tjj) {
279 rec = 1. / xj;
280 dscal_(n, &rec, &x[1], &c__1);
281 *scale *= rec;
282 xmax *= rec;
283 }
284 }
285 x[j1] /= tmp;
286 xj = (d__1 = x[j1], abs(d__1));
287
288 /* Scale x if necessary to avoid overflow when adding a */
289 /* multiple of column j1 of T. */
290
291 if (xj > 1.) {
292 rec = 1. / xj;
293 if (work[j1] > (bignum - xmax) * rec) {
294 dscal_(n, &rec, &x[1], &c__1);
295 *scale *= rec;
296 }
297 }
298 if (j1 > 1) {
299 i__1 = j1 - 1;
300 d__1 = -x[j1];
301 daxpy_(&i__1, &d__1, &t[j1 * t_dim1 + 1], &c__1, &x[1]
302 , &c__1);
303 i__1 = j1 - 1;
304 k = idamax_(&i__1, &x[1], &c__1);
305 xmax = (d__1 = x[k], abs(d__1));
306 }
307
308 } else {
309
310 /* Meet 2 by 2 diagonal block */
311
312 /* Call 2 by 2 linear system solve, to take */
313 /* care of possible overflow by scaling factor. */
314
315 d__[0] = x[j1];
316 d__[1] = x[j2];
317 dlaln2_(&c_false, &c__2, &c__1, &smin, &c_b21, &t[j1 + j1
318 * t_dim1], ldt, &c_b21, &c_b21, d__, &c__2, &
319 c_b25, &c_b25, v, &c__2, &scaloc, &xnorm, &ierr);
320 if (ierr != 0) {
321 *info = 2;
322 }
323
324 if (scaloc != 1.) {
325 dscal_(n, &scaloc, &x[1], &c__1);
326 *scale *= scaloc;
327 }
328 x[j1] = v[0];
329 x[j2] = v[1];
330
331 /* Scale V(1,1) (= X(J1)) and/or V(2,1) (=X(J2)) */
332 /* to avoid overflow in updating right-hand side. */
333
334 /* Computing MAX */
335 d__1 = abs(v[0]), d__2 = abs(v[1]);
336 xj = max(d__1,d__2);
337 if (xj > 1.) {
338 rec = 1. / xj;
339 /* Computing MAX */
340 d__1 = work[j1], d__2 = work[j2];
341 if (max(d__1,d__2) > (bignum - xmax) * rec) {
342 dscal_(n, &rec, &x[1], &c__1);
343 *scale *= rec;
344 }
345 }
346
347 /* Update right-hand side */
348
349 if (j1 > 1) {
350 i__1 = j1 - 1;
351 d__1 = -x[j1];
352 daxpy_(&i__1, &d__1, &t[j1 * t_dim1 + 1], &c__1, &x[1]
353 , &c__1);
354 i__1 = j1 - 1;
355 d__1 = -x[j2];
356 daxpy_(&i__1, &d__1, &t[j2 * t_dim1 + 1], &c__1, &x[1]
357 , &c__1);
358 i__1 = j1 - 1;
359 k = idamax_(&i__1, &x[1], &c__1);
360 xmax = (d__1 = x[k], abs(d__1));
361 }
362
363 }
364
365 L30:
366 ;
367 }
368
369 } else {
370
371 /* Solve T'*p = scale*c */
372
373 jnext = 1;
374 i__1 = *n;
375 for (j = 1; j <= i__1; ++j) {
376 if (j < jnext) {
377 goto L40;
378 }
379 j1 = j;
380 j2 = j;
381 jnext = j + 1;
382 if (j < *n) {
383 if (t[j + 1 + j * t_dim1] != 0.) {
384 j2 = j + 1;
385 jnext = j + 2;
386 }
387 }
388
389 if (j1 == j2) {
390
391 /* 1 by 1 diagonal block */
392
393 /* Scale if necessary to avoid overflow in forming the */
394 /* right-hand side element by inner product. */
395
396 xj = (d__1 = x[j1], abs(d__1));
397 if (xmax > 1.) {
398 rec = 1. / xmax;
399 if (work[j1] > (bignum - xj) * rec) {
400 dscal_(n, &rec, &x[1], &c__1);
401 *scale *= rec;
402 xmax *= rec;
403 }
404 }
405
406 i__2 = j1 - 1;
407 x[j1] -= ddot_(&i__2, &t[j1 * t_dim1 + 1], &c__1, &x[1], &
408 c__1);
409
410 xj = (d__1 = x[j1], abs(d__1));
411 tjj = (d__1 = t[j1 + j1 * t_dim1], abs(d__1));
412 tmp = t[j1 + j1 * t_dim1];
413 if (tjj < smin) {
414 tmp = smin;
415 tjj = smin;
416 *info = 1;
417 }
418
419 if (tjj < 1.) {
420 if (xj > bignum * tjj) {
421 rec = 1. / xj;
422 dscal_(n, &rec, &x[1], &c__1);
423 *scale *= rec;
424 xmax *= rec;
425 }
426 }
427 x[j1] /= tmp;
428 /* Computing MAX */
429 d__2 = xmax, d__3 = (d__1 = x[j1], abs(d__1));
430 xmax = max(d__2,d__3);
431
432 } else {
433
434 /* 2 by 2 diagonal block */
435
436 /* Scale if necessary to avoid overflow in forming the */
437 /* right-hand side elements by inner product. */
438
439 /* Computing MAX */
440 d__3 = (d__1 = x[j1], abs(d__1)), d__4 = (d__2 = x[j2],
441 abs(d__2));
442 xj = max(d__3,d__4);
443 if (xmax > 1.) {
444 rec = 1. / xmax;
445 /* Computing MAX */
446 d__1 = work[j2], d__2 = work[j1];
447 if (max(d__1,d__2) > (bignum - xj) * rec) {
448 dscal_(n, &rec, &x[1], &c__1);
449 *scale *= rec;
450 xmax *= rec;
451 }
452 }
453
454 i__2 = j1 - 1;
455 d__[0] = x[j1] - ddot_(&i__2, &t[j1 * t_dim1 + 1], &c__1,
456 &x[1], &c__1);
457 i__2 = j1 - 1;
458 d__[1] = x[j2] - ddot_(&i__2, &t[j2 * t_dim1 + 1], &c__1,
459 &x[1], &c__1);
460
461 dlaln2_(&c_true, &c__2, &c__1, &smin, &c_b21, &t[j1 + j1 *
462 t_dim1], ldt, &c_b21, &c_b21, d__, &c__2, &c_b25,
463 &c_b25, v, &c__2, &scaloc, &xnorm, &ierr);
464 if (ierr != 0) {
465 *info = 2;
466 }
467
468 if (scaloc != 1.) {
469 dscal_(n, &scaloc, &x[1], &c__1);
470 *scale *= scaloc;
471 }
472 x[j1] = v[0];
473 x[j2] = v[1];
474 /* Computing MAX */
475 d__3 = (d__1 = x[j1], abs(d__1)), d__4 = (d__2 = x[j2],
476 abs(d__2)), d__3 = max(d__3,d__4);
477 xmax = max(d__3,xmax);
478
479 }
480 L40:
481 ;
482 }
483 }
484
485 } else {
486
487 /* Computing MAX */
488 d__1 = eps * abs(*w);
489 sminw = max(d__1,smin);
490 if (notran) {
491
492 /* Solve (T + iB)*(p+iq) = c+id */
493
494 jnext = *n;
495 for (j = *n; j >= 1; --j) {
496 if (j > jnext) {
497 goto L70;
498 }
499 j1 = j;
500 j2 = j;
501 jnext = j - 1;
502 if (j > 1) {
503 if (t[j + (j - 1) * t_dim1] != 0.) {
504 j1 = j - 1;
505 jnext = j - 2;
506 }
507 }
508
509 if (j1 == j2) {
510
511 /* 1 by 1 diagonal block */
512
513 /* Scale if necessary to avoid overflow in division */
514
515 z__ = *w;
516 if (j1 == 1) {
517 z__ = b[1];
518 }
519 xj = (d__1 = x[j1], abs(d__1)) + (d__2 = x[*n + j1], abs(
520 d__2));
521 tjj = (d__1 = t[j1 + j1 * t_dim1], abs(d__1)) + abs(z__);
522 tmp = t[j1 + j1 * t_dim1];
523 if (tjj < sminw) {
524 tmp = sminw;
525 tjj = sminw;
526 *info = 1;
527 }
528
529 if (xj == 0.) {
530 goto L70;
531 }
532
533 if (tjj < 1.) {
534 if (xj > bignum * tjj) {
535 rec = 1. / xj;
536 dscal_(&n2, &rec, &x[1], &c__1);
537 *scale *= rec;
538 xmax *= rec;
539 }
540 }
541 dladiv_(&x[j1], &x[*n + j1], &tmp, &z__, &sr, &si);
542 x[j1] = sr;
543 x[*n + j1] = si;
544 xj = (d__1 = x[j1], abs(d__1)) + (d__2 = x[*n + j1], abs(
545 d__2));
546
547 /* Scale x if necessary to avoid overflow when adding a */
548 /* multiple of column j1 of T. */
549
550 if (xj > 1.) {
551 rec = 1. / xj;
552 if (work[j1] > (bignum - xmax) * rec) {
553 dscal_(&n2, &rec, &x[1], &c__1);
554 *scale *= rec;
555 }
556 }
557
558 if (j1 > 1) {
559 i__1 = j1 - 1;
560 d__1 = -x[j1];
561 daxpy_(&i__1, &d__1, &t[j1 * t_dim1 + 1], &c__1, &x[1]
562 , &c__1);
563 i__1 = j1 - 1;
564 d__1 = -x[*n + j1];
565 daxpy_(&i__1, &d__1, &t[j1 * t_dim1 + 1], &c__1, &x[*
566 n + 1], &c__1);
567
568 x[1] += b[j1] * x[*n + j1];
569 x[*n + 1] -= b[j1] * x[j1];
570
571 xmax = 0.;
572 i__1 = j1 - 1;
573 for (k = 1; k <= i__1; ++k) {
574 /* Computing MAX */
575 d__3 = xmax, d__4 = (d__1 = x[k], abs(d__1)) + (
576 d__2 = x[k + *n], abs(d__2));
577 xmax = max(d__3,d__4);
578 /* L50: */
579 }
580 }
581
582 } else {
583
584 /* Meet 2 by 2 diagonal block */
585
586 d__[0] = x[j1];
587 d__[1] = x[j2];
588 d__[2] = x[*n + j1];
589 d__[3] = x[*n + j2];
590 d__1 = -(*w);
591 dlaln2_(&c_false, &c__2, &c__2, &sminw, &c_b21, &t[j1 +
592 j1 * t_dim1], ldt, &c_b21, &c_b21, d__, &c__2, &
593 c_b25, &d__1, v, &c__2, &scaloc, &xnorm, &ierr);
594 if (ierr != 0) {
595 *info = 2;
596 }
597
598 if (scaloc != 1.) {
599 i__1 = *n << 1;
600 dscal_(&i__1, &scaloc, &x[1], &c__1);
601 *scale = scaloc * *scale;
602 }
603 x[j1] = v[0];
604 x[j2] = v[1];
605 x[*n + j1] = v[2];
606 x[*n + j2] = v[3];
607
608 /* Scale X(J1), .... to avoid overflow in */
609 /* updating right hand side. */
610
611 /* Computing MAX */
612 d__1 = abs(v[0]) + abs(v[2]), d__2 = abs(v[1]) + abs(v[3])
613 ;
614 xj = max(d__1,d__2);
615 if (xj > 1.) {
616 rec = 1. / xj;
617 /* Computing MAX */
618 d__1 = work[j1], d__2 = work[j2];
619 if (max(d__1,d__2) > (bignum - xmax) * rec) {
620 dscal_(&n2, &rec, &x[1], &c__1);
621 *scale *= rec;
622 }
623 }
624
625 /* Update the right-hand side. */
626
627 if (j1 > 1) {
628 i__1 = j1 - 1;
629 d__1 = -x[j1];
630 daxpy_(&i__1, &d__1, &t[j1 * t_dim1 + 1], &c__1, &x[1]
631 , &c__1);
632 i__1 = j1 - 1;
633 d__1 = -x[j2];
634 daxpy_(&i__1, &d__1, &t[j2 * t_dim1 + 1], &c__1, &x[1]
635 , &c__1);
636
637 i__1 = j1 - 1;
638 d__1 = -x[*n + j1];
639 daxpy_(&i__1, &d__1, &t[j1 * t_dim1 + 1], &c__1, &x[*
640 n + 1], &c__1);
641 i__1 = j1 - 1;
642 d__1 = -x[*n + j2];
643 daxpy_(&i__1, &d__1, &t[j2 * t_dim1 + 1], &c__1, &x[*
644 n + 1], &c__1);
645
646 x[1] = x[1] + b[j1] * x[*n + j1] + b[j2] * x[*n + j2];
647 x[*n + 1] = x[*n + 1] - b[j1] * x[j1] - b[j2] * x[j2];
648
649 xmax = 0.;
650 i__1 = j1 - 1;
651 for (k = 1; k <= i__1; ++k) {
652 /* Computing MAX */
653 d__3 = (d__1 = x[k], abs(d__1)) + (d__2 = x[k + *
654 n], abs(d__2));
655 xmax = max(d__3,xmax);
656 /* L60: */
657 }
658 }
659
660 }
661 L70:
662 ;
663 }
664
665 } else {
666
667 /* Solve (T + iB)'*(p+iq) = c+id */
668
669 jnext = 1;
670 i__1 = *n;
671 for (j = 1; j <= i__1; ++j) {
672 if (j < jnext) {
673 goto L80;
674 }
675 j1 = j;
676 j2 = j;
677 jnext = j + 1;
678 if (j < *n) {
679 if (t[j + 1 + j * t_dim1] != 0.) {
680 j2 = j + 1;
681 jnext = j + 2;
682 }
683 }
684
685 if (j1 == j2) {
686
687 /* 1 by 1 diagonal block */
688
689 /* Scale if necessary to avoid overflow in forming the */
690 /* right-hand side element by inner product. */
691
692 xj = (d__1 = x[j1], abs(d__1)) + (d__2 = x[j1 + *n], abs(
693 d__2));
694 if (xmax > 1.) {
695 rec = 1. / xmax;
696 if (work[j1] > (bignum - xj) * rec) {
697 dscal_(&n2, &rec, &x[1], &c__1);
698 *scale *= rec;
699 xmax *= rec;
700 }
701 }
702
703 i__2 = j1 - 1;
704 x[j1] -= ddot_(&i__2, &t[j1 * t_dim1 + 1], &c__1, &x[1], &
705 c__1);
706 i__2 = j1 - 1;
707 x[*n + j1] -= ddot_(&i__2, &t[j1 * t_dim1 + 1], &c__1, &x[
708 *n + 1], &c__1);
709 if (j1 > 1) {
710 x[j1] -= b[j1] * x[*n + 1];
711 x[*n + j1] += b[j1] * x[1];
712 }
713 xj = (d__1 = x[j1], abs(d__1)) + (d__2 = x[j1 + *n], abs(
714 d__2));
715
716 z__ = *w;
717 if (j1 == 1) {
718 z__ = b[1];
719 }
720
721 /* Scale if necessary to avoid overflow in */
722 /* complex division */
723
724 tjj = (d__1 = t[j1 + j1 * t_dim1], abs(d__1)) + abs(z__);
725 tmp = t[j1 + j1 * t_dim1];
726 if (tjj < sminw) {
727 tmp = sminw;
728 tjj = sminw;
729 *info = 1;
730 }
731
732 if (tjj < 1.) {
733 if (xj > bignum * tjj) {
734 rec = 1. / xj;
735 dscal_(&n2, &rec, &x[1], &c__1);
736 *scale *= rec;
737 xmax *= rec;
738 }
739 }
740 d__1 = -z__;
741 dladiv_(&x[j1], &x[*n + j1], &tmp, &d__1, &sr, &si);
742 x[j1] = sr;
743 x[j1 + *n] = si;
744 /* Computing MAX */
745 d__3 = (d__1 = x[j1], abs(d__1)) + (d__2 = x[j1 + *n],
746 abs(d__2));
747 xmax = max(d__3,xmax);
748
749 } else {
750
751 /* 2 by 2 diagonal block */
752
753 /* Scale if necessary to avoid overflow in forming the */
754 /* right-hand side element by inner product. */
755
756 /* Computing MAX */
757 d__5 = (d__1 = x[j1], abs(d__1)) + (d__2 = x[*n + j1],
758 abs(d__2)), d__6 = (d__3 = x[j2], abs(d__3)) + (
759 d__4 = x[*n + j2], abs(d__4));
760 xj = max(d__5,d__6);
761 if (xmax > 1.) {
762 rec = 1. / xmax;
763 /* Computing MAX */
764 d__1 = work[j1], d__2 = work[j2];
765 if (max(d__1,d__2) > (bignum - xj) / xmax) {
766 dscal_(&n2, &rec, &x[1], &c__1);
767 *scale *= rec;
768 xmax *= rec;
769 }
770 }
771
772 i__2 = j1 - 1;
773 d__[0] = x[j1] - ddot_(&i__2, &t[j1 * t_dim1 + 1], &c__1,
774 &x[1], &c__1);
775 i__2 = j1 - 1;
776 d__[1] = x[j2] - ddot_(&i__2, &t[j2 * t_dim1 + 1], &c__1,
777 &x[1], &c__1);
778 i__2 = j1 - 1;
779 d__[2] = x[*n + j1] - ddot_(&i__2, &t[j1 * t_dim1 + 1], &
780 c__1, &x[*n + 1], &c__1);
781 i__2 = j1 - 1;
782 d__[3] = x[*n + j2] - ddot_(&i__2, &t[j2 * t_dim1 + 1], &
783 c__1, &x[*n + 1], &c__1);
784 d__[0] -= b[j1] * x[*n + 1];
785 d__[1] -= b[j2] * x[*n + 1];
786 d__[2] += b[j1] * x[1];
787 d__[3] += b[j2] * x[1];
788
789 dlaln2_(&c_true, &c__2, &c__2, &sminw, &c_b21, &t[j1 + j1
790 * t_dim1], ldt, &c_b21, &c_b21, d__, &c__2, &
791 c_b25, w, v, &c__2, &scaloc, &xnorm, &ierr);
792 if (ierr != 0) {
793 *info = 2;
794 }
795
796 if (scaloc != 1.) {
797 dscal_(&n2, &scaloc, &x[1], &c__1);
798 *scale = scaloc * *scale;
799 }
800 x[j1] = v[0];
801 x[j2] = v[1];
802 x[*n + j1] = v[2];
803 x[*n + j2] = v[3];
804 /* Computing MAX */
805 d__5 = (d__1 = x[j1], abs(d__1)) + (d__2 = x[*n + j1],
806 abs(d__2)), d__6 = (d__3 = x[j2], abs(d__3)) + (
807 d__4 = x[*n + j2], abs(d__4)), d__5 = max(d__5,
808 d__6);
809 xmax = max(d__5,xmax);
810
811 }
812
813 L80:
814 ;
815 }
816
817 }
818
819 }
820
821 return 0;
822
823 /* End of DLAQTR */
824
825 } /* dlaqtr_ */
826
827