1 SUBROUTINE DLAQTR( LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK,
2 $ INFO )
3 *
4 * -- LAPACK auxiliary routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 LOGICAL LREAL, LTRAN
11 INTEGER INFO, LDT, N
12 DOUBLE PRECISION SCALE, W
13 * ..
14 * .. Array Arguments ..
15 DOUBLE PRECISION B( * ), T( LDT, * ), WORK( * ), X( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLAQTR solves the real quasi-triangular system
22 *
23 * op(T)*p = scale*c, if LREAL = .TRUE.
24 *
25 * or the complex quasi-triangular systems
26 *
27 * op(T + iB)*(p+iq) = scale*(c+id), if LREAL = .FALSE.
28 *
29 * in real arithmetic, where T is upper quasi-triangular.
30 * If LREAL = .FALSE., then the first diagonal block of T must be
31 * 1 by 1, B is the specially structured matrix
32 *
33 * B = [ b(1) b(2) ... b(n) ]
34 * [ w ]
35 * [ w ]
36 * [ . ]
37 * [ w ]
38 *
39 * op(A) = A or A**T, A**T denotes the transpose of
40 * matrix A.
41 *
42 * On input, X = [ c ]. On output, X = [ p ].
43 * [ d ] [ q ]
44 *
45 * This subroutine is designed for the condition number estimation
46 * in routine DTRSNA.
47 *
48 * Arguments
49 * =========
50 *
51 * LTRAN (input) LOGICAL
52 * On entry, LTRAN specifies the option of conjugate transpose:
53 * = .FALSE., op(T+i*B) = T+i*B,
54 * = .TRUE., op(T+i*B) = (T+i*B)**T.
55 *
56 * LREAL (input) LOGICAL
57 * On entry, LREAL specifies the input matrix structure:
58 * = .FALSE., the input is complex
59 * = .TRUE., the input is real
60 *
61 * N (input) INTEGER
62 * On entry, N specifies the order of T+i*B. N >= 0.
63 *
64 * T (input) DOUBLE PRECISION array, dimension (LDT,N)
65 * On entry, T contains a matrix in Schur canonical form.
66 * If LREAL = .FALSE., then the first diagonal block of T mu
67 * be 1 by 1.
68 *
69 * LDT (input) INTEGER
70 * The leading dimension of the matrix T. LDT >= max(1,N).
71 *
72 * B (input) DOUBLE PRECISION array, dimension (N)
73 * On entry, B contains the elements to form the matrix
74 * B as described above.
75 * If LREAL = .TRUE., B is not referenced.
76 *
77 * W (input) DOUBLE PRECISION
78 * On entry, W is the diagonal element of the matrix B.
79 * If LREAL = .TRUE., W is not referenced.
80 *
81 * SCALE (output) DOUBLE PRECISION
82 * On exit, SCALE is the scale factor.
83 *
84 * X (input/output) DOUBLE PRECISION array, dimension (2*N)
85 * On entry, X contains the right hand side of the system.
86 * On exit, X is overwritten by the solution.
87 *
88 * WORK (workspace) DOUBLE PRECISION array, dimension (N)
89 *
90 * INFO (output) INTEGER
91 * On exit, INFO is set to
92 * 0: successful exit.
93 * 1: the some diagonal 1 by 1 block has been perturbed by
94 * a small number SMIN to keep nonsingularity.
95 * 2: the some diagonal 2 by 2 block has been perturbed by
96 * a small number in DLALN2 to keep nonsingularity.
97 * NOTE: In the interests of speed, this routine does not
98 * check the inputs for errors.
99 *
100 * =====================================================================
101 *
102 * .. Parameters ..
103 DOUBLE PRECISION ZERO, ONE
104 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
105 * ..
106 * .. Local Scalars ..
107 LOGICAL NOTRAN
108 INTEGER I, IERR, J, J1, J2, JNEXT, K, N1, N2
109 DOUBLE PRECISION BIGNUM, EPS, REC, SCALOC, SI, SMIN, SMINW,
110 $ SMLNUM, SR, TJJ, TMP, XJ, XMAX, XNORM, Z
111 * ..
112 * .. Local Arrays ..
113 DOUBLE PRECISION D( 2, 2 ), V( 2, 2 )
114 * ..
115 * .. External Functions ..
116 INTEGER IDAMAX
117 DOUBLE PRECISION DASUM, DDOT, DLAMCH, DLANGE
118 EXTERNAL IDAMAX, DASUM, DDOT, DLAMCH, DLANGE
119 * ..
120 * .. External Subroutines ..
121 EXTERNAL DAXPY, DLADIV, DLALN2, DSCAL
122 * ..
123 * .. Intrinsic Functions ..
124 INTRINSIC ABS, MAX
125 * ..
126 * .. Executable Statements ..
127 *
128 * Do not test the input parameters for errors
129 *
130 NOTRAN = .NOT.LTRAN
131 INFO = 0
132 *
133 * Quick return if possible
134 *
135 IF( N.EQ.0 )
136 $ RETURN
137 *
138 * Set constants to control overflow
139 *
140 EPS = DLAMCH( 'P' )
141 SMLNUM = DLAMCH( 'S' ) / EPS
142 BIGNUM = ONE / SMLNUM
143 *
144 XNORM = DLANGE( 'M', N, N, T, LDT, D )
145 IF( .NOT.LREAL )
146 $ XNORM = MAX( XNORM, ABS( W ), DLANGE( 'M', N, 1, B, N, D ) )
147 SMIN = MAX( SMLNUM, EPS*XNORM )
148 *
149 * Compute 1-norm of each column of strictly upper triangular
150 * part of T to control overflow in triangular solver.
151 *
152 WORK( 1 ) = ZERO
153 DO 10 J = 2, N
154 WORK( J ) = DASUM( J-1, T( 1, J ), 1 )
155 10 CONTINUE
156 *
157 IF( .NOT.LREAL ) THEN
158 DO 20 I = 2, N
159 WORK( I ) = WORK( I ) + ABS( B( I ) )
160 20 CONTINUE
161 END IF
162 *
163 N2 = 2*N
164 N1 = N
165 IF( .NOT.LREAL )
166 $ N1 = N2
167 K = IDAMAX( N1, X, 1 )
168 XMAX = ABS( X( K ) )
169 SCALE = ONE
170 *
171 IF( XMAX.GT.BIGNUM ) THEN
172 SCALE = BIGNUM / XMAX
173 CALL DSCAL( N1, SCALE, X, 1 )
174 XMAX = BIGNUM
175 END IF
176 *
177 IF( LREAL ) THEN
178 *
179 IF( NOTRAN ) THEN
180 *
181 * Solve T*p = scale*c
182 *
183 JNEXT = N
184 DO 30 J = N, 1, -1
185 IF( J.GT.JNEXT )
186 $ GO TO 30
187 J1 = J
188 J2 = J
189 JNEXT = J - 1
190 IF( J.GT.1 ) THEN
191 IF( T( J, J-1 ).NE.ZERO ) THEN
192 J1 = J - 1
193 JNEXT = J - 2
194 END IF
195 END IF
196 *
197 IF( J1.EQ.J2 ) THEN
198 *
199 * Meet 1 by 1 diagonal block
200 *
201 * Scale to avoid overflow when computing
202 * x(j) = b(j)/T(j,j)
203 *
204 XJ = ABS( X( J1 ) )
205 TJJ = ABS( T( J1, J1 ) )
206 TMP = T( J1, J1 )
207 IF( TJJ.LT.SMIN ) THEN
208 TMP = SMIN
209 TJJ = SMIN
210 INFO = 1
211 END IF
212 *
213 IF( XJ.EQ.ZERO )
214 $ GO TO 30
215 *
216 IF( TJJ.LT.ONE ) THEN
217 IF( XJ.GT.BIGNUM*TJJ ) THEN
218 REC = ONE / XJ
219 CALL DSCAL( N, REC, X, 1 )
220 SCALE = SCALE*REC
221 XMAX = XMAX*REC
222 END IF
223 END IF
224 X( J1 ) = X( J1 ) / TMP
225 XJ = ABS( X( J1 ) )
226 *
227 * Scale x if necessary to avoid overflow when adding a
228 * multiple of column j1 of T.
229 *
230 IF( XJ.GT.ONE ) THEN
231 REC = ONE / XJ
232 IF( WORK( J1 ).GT.( BIGNUM-XMAX )*REC ) THEN
233 CALL DSCAL( N, REC, X, 1 )
234 SCALE = SCALE*REC
235 END IF
236 END IF
237 IF( J1.GT.1 ) THEN
238 CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
239 K = IDAMAX( J1-1, X, 1 )
240 XMAX = ABS( X( K ) )
241 END IF
242 *
243 ELSE
244 *
245 * Meet 2 by 2 diagonal block
246 *
247 * Call 2 by 2 linear system solve, to take
248 * care of possible overflow by scaling factor.
249 *
250 D( 1, 1 ) = X( J1 )
251 D( 2, 1 ) = X( J2 )
252 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, T( J1, J1 ),
253 $ LDT, ONE, ONE, D, 2, ZERO, ZERO, V, 2,
254 $ SCALOC, XNORM, IERR )
255 IF( IERR.NE.0 )
256 $ INFO = 2
257 *
258 IF( SCALOC.NE.ONE ) THEN
259 CALL DSCAL( N, SCALOC, X, 1 )
260 SCALE = SCALE*SCALOC
261 END IF
262 X( J1 ) = V( 1, 1 )
263 X( J2 ) = V( 2, 1 )
264 *
265 * Scale V(1,1) (= X(J1)) and/or V(2,1) (=X(J2))
266 * to avoid overflow in updating right-hand side.
267 *
268 XJ = MAX( ABS( V( 1, 1 ) ), ABS( V( 2, 1 ) ) )
269 IF( XJ.GT.ONE ) THEN
270 REC = ONE / XJ
271 IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
272 $ ( BIGNUM-XMAX )*REC ) THEN
273 CALL DSCAL( N, REC, X, 1 )
274 SCALE = SCALE*REC
275 END IF
276 END IF
277 *
278 * Update right-hand side
279 *
280 IF( J1.GT.1 ) THEN
281 CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
282 CALL DAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 )
283 K = IDAMAX( J1-1, X, 1 )
284 XMAX = ABS( X( K ) )
285 END IF
286 *
287 END IF
288 *
289 30 CONTINUE
290 *
291 ELSE
292 *
293 * Solve T**T*p = scale*c
294 *
295 JNEXT = 1
296 DO 40 J = 1, N
297 IF( J.LT.JNEXT )
298 $ GO TO 40
299 J1 = J
300 J2 = J
301 JNEXT = J + 1
302 IF( J.LT.N ) THEN
303 IF( T( J+1, J ).NE.ZERO ) THEN
304 J2 = J + 1
305 JNEXT = J + 2
306 END IF
307 END IF
308 *
309 IF( J1.EQ.J2 ) THEN
310 *
311 * 1 by 1 diagonal block
312 *
313 * Scale if necessary to avoid overflow in forming the
314 * right-hand side element by inner product.
315 *
316 XJ = ABS( X( J1 ) )
317 IF( XMAX.GT.ONE ) THEN
318 REC = ONE / XMAX
319 IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN
320 CALL DSCAL( N, REC, X, 1 )
321 SCALE = SCALE*REC
322 XMAX = XMAX*REC
323 END IF
324 END IF
325 *
326 X( J1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X, 1 )
327 *
328 XJ = ABS( X( J1 ) )
329 TJJ = ABS( T( J1, J1 ) )
330 TMP = T( J1, J1 )
331 IF( TJJ.LT.SMIN ) THEN
332 TMP = SMIN
333 TJJ = SMIN
334 INFO = 1
335 END IF
336 *
337 IF( TJJ.LT.ONE ) THEN
338 IF( XJ.GT.BIGNUM*TJJ ) THEN
339 REC = ONE / XJ
340 CALL DSCAL( N, REC, X, 1 )
341 SCALE = SCALE*REC
342 XMAX = XMAX*REC
343 END IF
344 END IF
345 X( J1 ) = X( J1 ) / TMP
346 XMAX = MAX( XMAX, ABS( X( J1 ) ) )
347 *
348 ELSE
349 *
350 * 2 by 2 diagonal block
351 *
352 * Scale if necessary to avoid overflow in forming the
353 * right-hand side elements by inner product.
354 *
355 XJ = MAX( ABS( X( J1 ) ), ABS( X( J2 ) ) )
356 IF( XMAX.GT.ONE ) THEN
357 REC = ONE / XMAX
358 IF( MAX( WORK( J2 ), WORK( J1 ) ).GT.( BIGNUM-XJ )*
359 $ REC ) THEN
360 CALL DSCAL( N, REC, X, 1 )
361 SCALE = SCALE*REC
362 XMAX = XMAX*REC
363 END IF
364 END IF
365 *
366 D( 1, 1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X,
367 $ 1 )
368 D( 2, 1 ) = X( J2 ) - DDOT( J1-1, T( 1, J2 ), 1, X,
369 $ 1 )
370 *
371 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J1, J1 ),
372 $ LDT, ONE, ONE, D, 2, ZERO, ZERO, V, 2,
373 $ SCALOC, XNORM, IERR )
374 IF( IERR.NE.0 )
375 $ INFO = 2
376 *
377 IF( SCALOC.NE.ONE ) THEN
378 CALL DSCAL( N, SCALOC, X, 1 )
379 SCALE = SCALE*SCALOC
380 END IF
381 X( J1 ) = V( 1, 1 )
382 X( J2 ) = V( 2, 1 )
383 XMAX = MAX( ABS( X( J1 ) ), ABS( X( J2 ) ), XMAX )
384 *
385 END IF
386 40 CONTINUE
387 END IF
388 *
389 ELSE
390 *
391 SMINW = MAX( EPS*ABS( W ), SMIN )
392 IF( NOTRAN ) THEN
393 *
394 * Solve (T + iB)*(p+iq) = c+id
395 *
396 JNEXT = N
397 DO 70 J = N, 1, -1
398 IF( J.GT.JNEXT )
399 $ GO TO 70
400 J1 = J
401 J2 = J
402 JNEXT = J - 1
403 IF( J.GT.1 ) THEN
404 IF( T( J, J-1 ).NE.ZERO ) THEN
405 J1 = J - 1
406 JNEXT = J - 2
407 END IF
408 END IF
409 *
410 IF( J1.EQ.J2 ) THEN
411 *
412 * 1 by 1 diagonal block
413 *
414 * Scale if necessary to avoid overflow in division
415 *
416 Z = W
417 IF( J1.EQ.1 )
418 $ Z = B( 1 )
419 XJ = ABS( X( J1 ) ) + ABS( X( N+J1 ) )
420 TJJ = ABS( T( J1, J1 ) ) + ABS( Z )
421 TMP = T( J1, J1 )
422 IF( TJJ.LT.SMINW ) THEN
423 TMP = SMINW
424 TJJ = SMINW
425 INFO = 1
426 END IF
427 *
428 IF( XJ.EQ.ZERO )
429 $ GO TO 70
430 *
431 IF( TJJ.LT.ONE ) THEN
432 IF( XJ.GT.BIGNUM*TJJ ) THEN
433 REC = ONE / XJ
434 CALL DSCAL( N2, REC, X, 1 )
435 SCALE = SCALE*REC
436 XMAX = XMAX*REC
437 END IF
438 END IF
439 CALL DLADIV( X( J1 ), X( N+J1 ), TMP, Z, SR, SI )
440 X( J1 ) = SR
441 X( N+J1 ) = SI
442 XJ = ABS( X( J1 ) ) + ABS( X( N+J1 ) )
443 *
444 * Scale x if necessary to avoid overflow when adding a
445 * multiple of column j1 of T.
446 *
447 IF( XJ.GT.ONE ) THEN
448 REC = ONE / XJ
449 IF( WORK( J1 ).GT.( BIGNUM-XMAX )*REC ) THEN
450 CALL DSCAL( N2, REC, X, 1 )
451 SCALE = SCALE*REC
452 END IF
453 END IF
454 *
455 IF( J1.GT.1 ) THEN
456 CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
457 CALL DAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1,
458 $ X( N+1 ), 1 )
459 *
460 X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 )
461 X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 )
462 *
463 XMAX = ZERO
464 DO 50 K = 1, J1 - 1
465 XMAX = MAX( XMAX, ABS( X( K ) )+
466 $ ABS( X( K+N ) ) )
467 50 CONTINUE
468 END IF
469 *
470 ELSE
471 *
472 * Meet 2 by 2 diagonal block
473 *
474 D( 1, 1 ) = X( J1 )
475 D( 2, 1 ) = X( J2 )
476 D( 1, 2 ) = X( N+J1 )
477 D( 2, 2 ) = X( N+J2 )
478 CALL DLALN2( .FALSE., 2, 2, SMINW, ONE, T( J1, J1 ),
479 $ LDT, ONE, ONE, D, 2, ZERO, -W, V, 2,
480 $ SCALOC, XNORM, IERR )
481 IF( IERR.NE.0 )
482 $ INFO = 2
483 *
484 IF( SCALOC.NE.ONE ) THEN
485 CALL DSCAL( 2*N, SCALOC, X, 1 )
486 SCALE = SCALOC*SCALE
487 END IF
488 X( J1 ) = V( 1, 1 )
489 X( J2 ) = V( 2, 1 )
490 X( N+J1 ) = V( 1, 2 )
491 X( N+J2 ) = V( 2, 2 )
492 *
493 * Scale X(J1), .... to avoid overflow in
494 * updating right hand side.
495 *
496 XJ = MAX( ABS( V( 1, 1 ) )+ABS( V( 1, 2 ) ),
497 $ ABS( V( 2, 1 ) )+ABS( V( 2, 2 ) ) )
498 IF( XJ.GT.ONE ) THEN
499 REC = ONE / XJ
500 IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
501 $ ( BIGNUM-XMAX )*REC ) THEN
502 CALL DSCAL( N2, REC, X, 1 )
503 SCALE = SCALE*REC
504 END IF
505 END IF
506 *
507 * Update the right-hand side.
508 *
509 IF( J1.GT.1 ) THEN
510 CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
511 CALL DAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 )
512 *
513 CALL DAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1,
514 $ X( N+1 ), 1 )
515 CALL DAXPY( J1-1, -X( N+J2 ), T( 1, J2 ), 1,
516 $ X( N+1 ), 1 )
517 *
518 X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 ) +
519 $ B( J2 )*X( N+J2 )
520 X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 ) -
521 $ B( J2 )*X( J2 )
522 *
523 XMAX = ZERO
524 DO 60 K = 1, J1 - 1
525 XMAX = MAX( ABS( X( K ) )+ABS( X( K+N ) ),
526 $ XMAX )
527 60 CONTINUE
528 END IF
529 *
530 END IF
531 70 CONTINUE
532 *
533 ELSE
534 *
535 * Solve (T + iB)**T*(p+iq) = c+id
536 *
537 JNEXT = 1
538 DO 80 J = 1, N
539 IF( J.LT.JNEXT )
540 $ GO TO 80
541 J1 = J
542 J2 = J
543 JNEXT = J + 1
544 IF( J.LT.N ) THEN
545 IF( T( J+1, J ).NE.ZERO ) THEN
546 J2 = J + 1
547 JNEXT = J + 2
548 END IF
549 END IF
550 *
551 IF( J1.EQ.J2 ) THEN
552 *
553 * 1 by 1 diagonal block
554 *
555 * Scale if necessary to avoid overflow in forming the
556 * right-hand side element by inner product.
557 *
558 XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) )
559 IF( XMAX.GT.ONE ) THEN
560 REC = ONE / XMAX
561 IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN
562 CALL DSCAL( N2, REC, X, 1 )
563 SCALE = SCALE*REC
564 XMAX = XMAX*REC
565 END IF
566 END IF
567 *
568 X( J1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X, 1 )
569 X( N+J1 ) = X( N+J1 ) - DDOT( J1-1, T( 1, J1 ), 1,
570 $ X( N+1 ), 1 )
571 IF( J1.GT.1 ) THEN
572 X( J1 ) = X( J1 ) - B( J1 )*X( N+1 )
573 X( N+J1 ) = X( N+J1 ) + B( J1 )*X( 1 )
574 END IF
575 XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) )
576 *
577 Z = W
578 IF( J1.EQ.1 )
579 $ Z = B( 1 )
580 *
581 * Scale if necessary to avoid overflow in
582 * complex division
583 *
584 TJJ = ABS( T( J1, J1 ) ) + ABS( Z )
585 TMP = T( J1, J1 )
586 IF( TJJ.LT.SMINW ) THEN
587 TMP = SMINW
588 TJJ = SMINW
589 INFO = 1
590 END IF
591 *
592 IF( TJJ.LT.ONE ) THEN
593 IF( XJ.GT.BIGNUM*TJJ ) THEN
594 REC = ONE / XJ
595 CALL DSCAL( N2, REC, X, 1 )
596 SCALE = SCALE*REC
597 XMAX = XMAX*REC
598 END IF
599 END IF
600 CALL DLADIV( X( J1 ), X( N+J1 ), TMP, -Z, SR, SI )
601 X( J1 ) = SR
602 X( J1+N ) = SI
603 XMAX = MAX( ABS( X( J1 ) )+ABS( X( J1+N ) ), XMAX )
604 *
605 ELSE
606 *
607 * 2 by 2 diagonal block
608 *
609 * Scale if necessary to avoid overflow in forming the
610 * right-hand side element by inner product.
611 *
612 XJ = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ),
613 $ ABS( X( J2 ) )+ABS( X( N+J2 ) ) )
614 IF( XMAX.GT.ONE ) THEN
615 REC = ONE / XMAX
616 IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
617 $ ( BIGNUM-XJ ) / XMAX ) THEN
618 CALL DSCAL( N2, REC, X, 1 )
619 SCALE = SCALE*REC
620 XMAX = XMAX*REC
621 END IF
622 END IF
623 *
624 D( 1, 1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X,
625 $ 1 )
626 D( 2, 1 ) = X( J2 ) - DDOT( J1-1, T( 1, J2 ), 1, X,
627 $ 1 )
628 D( 1, 2 ) = X( N+J1 ) - DDOT( J1-1, T( 1, J1 ), 1,
629 $ X( N+1 ), 1 )
630 D( 2, 2 ) = X( N+J2 ) - DDOT( J1-1, T( 1, J2 ), 1,
631 $ X( N+1 ), 1 )
632 D( 1, 1 ) = D( 1, 1 ) - B( J1 )*X( N+1 )
633 D( 2, 1 ) = D( 2, 1 ) - B( J2 )*X( N+1 )
634 D( 1, 2 ) = D( 1, 2 ) + B( J1 )*X( 1 )
635 D( 2, 2 ) = D( 2, 2 ) + B( J2 )*X( 1 )
636 *
637 CALL DLALN2( .TRUE., 2, 2, SMINW, ONE, T( J1, J1 ),
638 $ LDT, ONE, ONE, D, 2, ZERO, W, V, 2,
639 $ SCALOC, XNORM, IERR )
640 IF( IERR.NE.0 )
641 $ INFO = 2
642 *
643 IF( SCALOC.NE.ONE ) THEN
644 CALL DSCAL( N2, SCALOC, X, 1 )
645 SCALE = SCALOC*SCALE
646 END IF
647 X( J1 ) = V( 1, 1 )
648 X( J2 ) = V( 2, 1 )
649 X( N+J1 ) = V( 1, 2 )
650 X( N+J2 ) = V( 2, 2 )
651 XMAX = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ),
652 $ ABS( X( J2 ) )+ABS( X( N+J2 ) ), XMAX )
653 *
654 END IF
655 *
656 80 CONTINUE
657 *
658 END IF
659 *
660 END IF
661 *
662 RETURN
663 *
664 * End of DLAQTR
665 *
666 END
2 $ INFO )
3 *
4 * -- LAPACK auxiliary routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 LOGICAL LREAL, LTRAN
11 INTEGER INFO, LDT, N
12 DOUBLE PRECISION SCALE, W
13 * ..
14 * .. Array Arguments ..
15 DOUBLE PRECISION B( * ), T( LDT, * ), WORK( * ), X( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLAQTR solves the real quasi-triangular system
22 *
23 * op(T)*p = scale*c, if LREAL = .TRUE.
24 *
25 * or the complex quasi-triangular systems
26 *
27 * op(T + iB)*(p+iq) = scale*(c+id), if LREAL = .FALSE.
28 *
29 * in real arithmetic, where T is upper quasi-triangular.
30 * If LREAL = .FALSE., then the first diagonal block of T must be
31 * 1 by 1, B is the specially structured matrix
32 *
33 * B = [ b(1) b(2) ... b(n) ]
34 * [ w ]
35 * [ w ]
36 * [ . ]
37 * [ w ]
38 *
39 * op(A) = A or A**T, A**T denotes the transpose of
40 * matrix A.
41 *
42 * On input, X = [ c ]. On output, X = [ p ].
43 * [ d ] [ q ]
44 *
45 * This subroutine is designed for the condition number estimation
46 * in routine DTRSNA.
47 *
48 * Arguments
49 * =========
50 *
51 * LTRAN (input) LOGICAL
52 * On entry, LTRAN specifies the option of conjugate transpose:
53 * = .FALSE., op(T+i*B) = T+i*B,
54 * = .TRUE., op(T+i*B) = (T+i*B)**T.
55 *
56 * LREAL (input) LOGICAL
57 * On entry, LREAL specifies the input matrix structure:
58 * = .FALSE., the input is complex
59 * = .TRUE., the input is real
60 *
61 * N (input) INTEGER
62 * On entry, N specifies the order of T+i*B. N >= 0.
63 *
64 * T (input) DOUBLE PRECISION array, dimension (LDT,N)
65 * On entry, T contains a matrix in Schur canonical form.
66 * If LREAL = .FALSE., then the first diagonal block of T mu
67 * be 1 by 1.
68 *
69 * LDT (input) INTEGER
70 * The leading dimension of the matrix T. LDT >= max(1,N).
71 *
72 * B (input) DOUBLE PRECISION array, dimension (N)
73 * On entry, B contains the elements to form the matrix
74 * B as described above.
75 * If LREAL = .TRUE., B is not referenced.
76 *
77 * W (input) DOUBLE PRECISION
78 * On entry, W is the diagonal element of the matrix B.
79 * If LREAL = .TRUE., W is not referenced.
80 *
81 * SCALE (output) DOUBLE PRECISION
82 * On exit, SCALE is the scale factor.
83 *
84 * X (input/output) DOUBLE PRECISION array, dimension (2*N)
85 * On entry, X contains the right hand side of the system.
86 * On exit, X is overwritten by the solution.
87 *
88 * WORK (workspace) DOUBLE PRECISION array, dimension (N)
89 *
90 * INFO (output) INTEGER
91 * On exit, INFO is set to
92 * 0: successful exit.
93 * 1: the some diagonal 1 by 1 block has been perturbed by
94 * a small number SMIN to keep nonsingularity.
95 * 2: the some diagonal 2 by 2 block has been perturbed by
96 * a small number in DLALN2 to keep nonsingularity.
97 * NOTE: In the interests of speed, this routine does not
98 * check the inputs for errors.
99 *
100 * =====================================================================
101 *
102 * .. Parameters ..
103 DOUBLE PRECISION ZERO, ONE
104 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
105 * ..
106 * .. Local Scalars ..
107 LOGICAL NOTRAN
108 INTEGER I, IERR, J, J1, J2, JNEXT, K, N1, N2
109 DOUBLE PRECISION BIGNUM, EPS, REC, SCALOC, SI, SMIN, SMINW,
110 $ SMLNUM, SR, TJJ, TMP, XJ, XMAX, XNORM, Z
111 * ..
112 * .. Local Arrays ..
113 DOUBLE PRECISION D( 2, 2 ), V( 2, 2 )
114 * ..
115 * .. External Functions ..
116 INTEGER IDAMAX
117 DOUBLE PRECISION DASUM, DDOT, DLAMCH, DLANGE
118 EXTERNAL IDAMAX, DASUM, DDOT, DLAMCH, DLANGE
119 * ..
120 * .. External Subroutines ..
121 EXTERNAL DAXPY, DLADIV, DLALN2, DSCAL
122 * ..
123 * .. Intrinsic Functions ..
124 INTRINSIC ABS, MAX
125 * ..
126 * .. Executable Statements ..
127 *
128 * Do not test the input parameters for errors
129 *
130 NOTRAN = .NOT.LTRAN
131 INFO = 0
132 *
133 * Quick return if possible
134 *
135 IF( N.EQ.0 )
136 $ RETURN
137 *
138 * Set constants to control overflow
139 *
140 EPS = DLAMCH( 'P' )
141 SMLNUM = DLAMCH( 'S' ) / EPS
142 BIGNUM = ONE / SMLNUM
143 *
144 XNORM = DLANGE( 'M', N, N, T, LDT, D )
145 IF( .NOT.LREAL )
146 $ XNORM = MAX( XNORM, ABS( W ), DLANGE( 'M', N, 1, B, N, D ) )
147 SMIN = MAX( SMLNUM, EPS*XNORM )
148 *
149 * Compute 1-norm of each column of strictly upper triangular
150 * part of T to control overflow in triangular solver.
151 *
152 WORK( 1 ) = ZERO
153 DO 10 J = 2, N
154 WORK( J ) = DASUM( J-1, T( 1, J ), 1 )
155 10 CONTINUE
156 *
157 IF( .NOT.LREAL ) THEN
158 DO 20 I = 2, N
159 WORK( I ) = WORK( I ) + ABS( B( I ) )
160 20 CONTINUE
161 END IF
162 *
163 N2 = 2*N
164 N1 = N
165 IF( .NOT.LREAL )
166 $ N1 = N2
167 K = IDAMAX( N1, X, 1 )
168 XMAX = ABS( X( K ) )
169 SCALE = ONE
170 *
171 IF( XMAX.GT.BIGNUM ) THEN
172 SCALE = BIGNUM / XMAX
173 CALL DSCAL( N1, SCALE, X, 1 )
174 XMAX = BIGNUM
175 END IF
176 *
177 IF( LREAL ) THEN
178 *
179 IF( NOTRAN ) THEN
180 *
181 * Solve T*p = scale*c
182 *
183 JNEXT = N
184 DO 30 J = N, 1, -1
185 IF( J.GT.JNEXT )
186 $ GO TO 30
187 J1 = J
188 J2 = J
189 JNEXT = J - 1
190 IF( J.GT.1 ) THEN
191 IF( T( J, J-1 ).NE.ZERO ) THEN
192 J1 = J - 1
193 JNEXT = J - 2
194 END IF
195 END IF
196 *
197 IF( J1.EQ.J2 ) THEN
198 *
199 * Meet 1 by 1 diagonal block
200 *
201 * Scale to avoid overflow when computing
202 * x(j) = b(j)/T(j,j)
203 *
204 XJ = ABS( X( J1 ) )
205 TJJ = ABS( T( J1, J1 ) )
206 TMP = T( J1, J1 )
207 IF( TJJ.LT.SMIN ) THEN
208 TMP = SMIN
209 TJJ = SMIN
210 INFO = 1
211 END IF
212 *
213 IF( XJ.EQ.ZERO )
214 $ GO TO 30
215 *
216 IF( TJJ.LT.ONE ) THEN
217 IF( XJ.GT.BIGNUM*TJJ ) THEN
218 REC = ONE / XJ
219 CALL DSCAL( N, REC, X, 1 )
220 SCALE = SCALE*REC
221 XMAX = XMAX*REC
222 END IF
223 END IF
224 X( J1 ) = X( J1 ) / TMP
225 XJ = ABS( X( J1 ) )
226 *
227 * Scale x if necessary to avoid overflow when adding a
228 * multiple of column j1 of T.
229 *
230 IF( XJ.GT.ONE ) THEN
231 REC = ONE / XJ
232 IF( WORK( J1 ).GT.( BIGNUM-XMAX )*REC ) THEN
233 CALL DSCAL( N, REC, X, 1 )
234 SCALE = SCALE*REC
235 END IF
236 END IF
237 IF( J1.GT.1 ) THEN
238 CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
239 K = IDAMAX( J1-1, X, 1 )
240 XMAX = ABS( X( K ) )
241 END IF
242 *
243 ELSE
244 *
245 * Meet 2 by 2 diagonal block
246 *
247 * Call 2 by 2 linear system solve, to take
248 * care of possible overflow by scaling factor.
249 *
250 D( 1, 1 ) = X( J1 )
251 D( 2, 1 ) = X( J2 )
252 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, T( J1, J1 ),
253 $ LDT, ONE, ONE, D, 2, ZERO, ZERO, V, 2,
254 $ SCALOC, XNORM, IERR )
255 IF( IERR.NE.0 )
256 $ INFO = 2
257 *
258 IF( SCALOC.NE.ONE ) THEN
259 CALL DSCAL( N, SCALOC, X, 1 )
260 SCALE = SCALE*SCALOC
261 END IF
262 X( J1 ) = V( 1, 1 )
263 X( J2 ) = V( 2, 1 )
264 *
265 * Scale V(1,1) (= X(J1)) and/or V(2,1) (=X(J2))
266 * to avoid overflow in updating right-hand side.
267 *
268 XJ = MAX( ABS( V( 1, 1 ) ), ABS( V( 2, 1 ) ) )
269 IF( XJ.GT.ONE ) THEN
270 REC = ONE / XJ
271 IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
272 $ ( BIGNUM-XMAX )*REC ) THEN
273 CALL DSCAL( N, REC, X, 1 )
274 SCALE = SCALE*REC
275 END IF
276 END IF
277 *
278 * Update right-hand side
279 *
280 IF( J1.GT.1 ) THEN
281 CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
282 CALL DAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 )
283 K = IDAMAX( J1-1, X, 1 )
284 XMAX = ABS( X( K ) )
285 END IF
286 *
287 END IF
288 *
289 30 CONTINUE
290 *
291 ELSE
292 *
293 * Solve T**T*p = scale*c
294 *
295 JNEXT = 1
296 DO 40 J = 1, N
297 IF( J.LT.JNEXT )
298 $ GO TO 40
299 J1 = J
300 J2 = J
301 JNEXT = J + 1
302 IF( J.LT.N ) THEN
303 IF( T( J+1, J ).NE.ZERO ) THEN
304 J2 = J + 1
305 JNEXT = J + 2
306 END IF
307 END IF
308 *
309 IF( J1.EQ.J2 ) THEN
310 *
311 * 1 by 1 diagonal block
312 *
313 * Scale if necessary to avoid overflow in forming the
314 * right-hand side element by inner product.
315 *
316 XJ = ABS( X( J1 ) )
317 IF( XMAX.GT.ONE ) THEN
318 REC = ONE / XMAX
319 IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN
320 CALL DSCAL( N, REC, X, 1 )
321 SCALE = SCALE*REC
322 XMAX = XMAX*REC
323 END IF
324 END IF
325 *
326 X( J1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X, 1 )
327 *
328 XJ = ABS( X( J1 ) )
329 TJJ = ABS( T( J1, J1 ) )
330 TMP = T( J1, J1 )
331 IF( TJJ.LT.SMIN ) THEN
332 TMP = SMIN
333 TJJ = SMIN
334 INFO = 1
335 END IF
336 *
337 IF( TJJ.LT.ONE ) THEN
338 IF( XJ.GT.BIGNUM*TJJ ) THEN
339 REC = ONE / XJ
340 CALL DSCAL( N, REC, X, 1 )
341 SCALE = SCALE*REC
342 XMAX = XMAX*REC
343 END IF
344 END IF
345 X( J1 ) = X( J1 ) / TMP
346 XMAX = MAX( XMAX, ABS( X( J1 ) ) )
347 *
348 ELSE
349 *
350 * 2 by 2 diagonal block
351 *
352 * Scale if necessary to avoid overflow in forming the
353 * right-hand side elements by inner product.
354 *
355 XJ = MAX( ABS( X( J1 ) ), ABS( X( J2 ) ) )
356 IF( XMAX.GT.ONE ) THEN
357 REC = ONE / XMAX
358 IF( MAX( WORK( J2 ), WORK( J1 ) ).GT.( BIGNUM-XJ )*
359 $ REC ) THEN
360 CALL DSCAL( N, REC, X, 1 )
361 SCALE = SCALE*REC
362 XMAX = XMAX*REC
363 END IF
364 END IF
365 *
366 D( 1, 1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X,
367 $ 1 )
368 D( 2, 1 ) = X( J2 ) - DDOT( J1-1, T( 1, J2 ), 1, X,
369 $ 1 )
370 *
371 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J1, J1 ),
372 $ LDT, ONE, ONE, D, 2, ZERO, ZERO, V, 2,
373 $ SCALOC, XNORM, IERR )
374 IF( IERR.NE.0 )
375 $ INFO = 2
376 *
377 IF( SCALOC.NE.ONE ) THEN
378 CALL DSCAL( N, SCALOC, X, 1 )
379 SCALE = SCALE*SCALOC
380 END IF
381 X( J1 ) = V( 1, 1 )
382 X( J2 ) = V( 2, 1 )
383 XMAX = MAX( ABS( X( J1 ) ), ABS( X( J2 ) ), XMAX )
384 *
385 END IF
386 40 CONTINUE
387 END IF
388 *
389 ELSE
390 *
391 SMINW = MAX( EPS*ABS( W ), SMIN )
392 IF( NOTRAN ) THEN
393 *
394 * Solve (T + iB)*(p+iq) = c+id
395 *
396 JNEXT = N
397 DO 70 J = N, 1, -1
398 IF( J.GT.JNEXT )
399 $ GO TO 70
400 J1 = J
401 J2 = J
402 JNEXT = J - 1
403 IF( J.GT.1 ) THEN
404 IF( T( J, J-1 ).NE.ZERO ) THEN
405 J1 = J - 1
406 JNEXT = J - 2
407 END IF
408 END IF
409 *
410 IF( J1.EQ.J2 ) THEN
411 *
412 * 1 by 1 diagonal block
413 *
414 * Scale if necessary to avoid overflow in division
415 *
416 Z = W
417 IF( J1.EQ.1 )
418 $ Z = B( 1 )
419 XJ = ABS( X( J1 ) ) + ABS( X( N+J1 ) )
420 TJJ = ABS( T( J1, J1 ) ) + ABS( Z )
421 TMP = T( J1, J1 )
422 IF( TJJ.LT.SMINW ) THEN
423 TMP = SMINW
424 TJJ = SMINW
425 INFO = 1
426 END IF
427 *
428 IF( XJ.EQ.ZERO )
429 $ GO TO 70
430 *
431 IF( TJJ.LT.ONE ) THEN
432 IF( XJ.GT.BIGNUM*TJJ ) THEN
433 REC = ONE / XJ
434 CALL DSCAL( N2, REC, X, 1 )
435 SCALE = SCALE*REC
436 XMAX = XMAX*REC
437 END IF
438 END IF
439 CALL DLADIV( X( J1 ), X( N+J1 ), TMP, Z, SR, SI )
440 X( J1 ) = SR
441 X( N+J1 ) = SI
442 XJ = ABS( X( J1 ) ) + ABS( X( N+J1 ) )
443 *
444 * Scale x if necessary to avoid overflow when adding a
445 * multiple of column j1 of T.
446 *
447 IF( XJ.GT.ONE ) THEN
448 REC = ONE / XJ
449 IF( WORK( J1 ).GT.( BIGNUM-XMAX )*REC ) THEN
450 CALL DSCAL( N2, REC, X, 1 )
451 SCALE = SCALE*REC
452 END IF
453 END IF
454 *
455 IF( J1.GT.1 ) THEN
456 CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
457 CALL DAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1,
458 $ X( N+1 ), 1 )
459 *
460 X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 )
461 X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 )
462 *
463 XMAX = ZERO
464 DO 50 K = 1, J1 - 1
465 XMAX = MAX( XMAX, ABS( X( K ) )+
466 $ ABS( X( K+N ) ) )
467 50 CONTINUE
468 END IF
469 *
470 ELSE
471 *
472 * Meet 2 by 2 diagonal block
473 *
474 D( 1, 1 ) = X( J1 )
475 D( 2, 1 ) = X( J2 )
476 D( 1, 2 ) = X( N+J1 )
477 D( 2, 2 ) = X( N+J2 )
478 CALL DLALN2( .FALSE., 2, 2, SMINW, ONE, T( J1, J1 ),
479 $ LDT, ONE, ONE, D, 2, ZERO, -W, V, 2,
480 $ SCALOC, XNORM, IERR )
481 IF( IERR.NE.0 )
482 $ INFO = 2
483 *
484 IF( SCALOC.NE.ONE ) THEN
485 CALL DSCAL( 2*N, SCALOC, X, 1 )
486 SCALE = SCALOC*SCALE
487 END IF
488 X( J1 ) = V( 1, 1 )
489 X( J2 ) = V( 2, 1 )
490 X( N+J1 ) = V( 1, 2 )
491 X( N+J2 ) = V( 2, 2 )
492 *
493 * Scale X(J1), .... to avoid overflow in
494 * updating right hand side.
495 *
496 XJ = MAX( ABS( V( 1, 1 ) )+ABS( V( 1, 2 ) ),
497 $ ABS( V( 2, 1 ) )+ABS( V( 2, 2 ) ) )
498 IF( XJ.GT.ONE ) THEN
499 REC = ONE / XJ
500 IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
501 $ ( BIGNUM-XMAX )*REC ) THEN
502 CALL DSCAL( N2, REC, X, 1 )
503 SCALE = SCALE*REC
504 END IF
505 END IF
506 *
507 * Update the right-hand side.
508 *
509 IF( J1.GT.1 ) THEN
510 CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
511 CALL DAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 )
512 *
513 CALL DAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1,
514 $ X( N+1 ), 1 )
515 CALL DAXPY( J1-1, -X( N+J2 ), T( 1, J2 ), 1,
516 $ X( N+1 ), 1 )
517 *
518 X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 ) +
519 $ B( J2 )*X( N+J2 )
520 X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 ) -
521 $ B( J2 )*X( J2 )
522 *
523 XMAX = ZERO
524 DO 60 K = 1, J1 - 1
525 XMAX = MAX( ABS( X( K ) )+ABS( X( K+N ) ),
526 $ XMAX )
527 60 CONTINUE
528 END IF
529 *
530 END IF
531 70 CONTINUE
532 *
533 ELSE
534 *
535 * Solve (T + iB)**T*(p+iq) = c+id
536 *
537 JNEXT = 1
538 DO 80 J = 1, N
539 IF( J.LT.JNEXT )
540 $ GO TO 80
541 J1 = J
542 J2 = J
543 JNEXT = J + 1
544 IF( J.LT.N ) THEN
545 IF( T( J+1, J ).NE.ZERO ) THEN
546 J2 = J + 1
547 JNEXT = J + 2
548 END IF
549 END IF
550 *
551 IF( J1.EQ.J2 ) THEN
552 *
553 * 1 by 1 diagonal block
554 *
555 * Scale if necessary to avoid overflow in forming the
556 * right-hand side element by inner product.
557 *
558 XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) )
559 IF( XMAX.GT.ONE ) THEN
560 REC = ONE / XMAX
561 IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN
562 CALL DSCAL( N2, REC, X, 1 )
563 SCALE = SCALE*REC
564 XMAX = XMAX*REC
565 END IF
566 END IF
567 *
568 X( J1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X, 1 )
569 X( N+J1 ) = X( N+J1 ) - DDOT( J1-1, T( 1, J1 ), 1,
570 $ X( N+1 ), 1 )
571 IF( J1.GT.1 ) THEN
572 X( J1 ) = X( J1 ) - B( J1 )*X( N+1 )
573 X( N+J1 ) = X( N+J1 ) + B( J1 )*X( 1 )
574 END IF
575 XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) )
576 *
577 Z = W
578 IF( J1.EQ.1 )
579 $ Z = B( 1 )
580 *
581 * Scale if necessary to avoid overflow in
582 * complex division
583 *
584 TJJ = ABS( T( J1, J1 ) ) + ABS( Z )
585 TMP = T( J1, J1 )
586 IF( TJJ.LT.SMINW ) THEN
587 TMP = SMINW
588 TJJ = SMINW
589 INFO = 1
590 END IF
591 *
592 IF( TJJ.LT.ONE ) THEN
593 IF( XJ.GT.BIGNUM*TJJ ) THEN
594 REC = ONE / XJ
595 CALL DSCAL( N2, REC, X, 1 )
596 SCALE = SCALE*REC
597 XMAX = XMAX*REC
598 END IF
599 END IF
600 CALL DLADIV( X( J1 ), X( N+J1 ), TMP, -Z, SR, SI )
601 X( J1 ) = SR
602 X( J1+N ) = SI
603 XMAX = MAX( ABS( X( J1 ) )+ABS( X( J1+N ) ), XMAX )
604 *
605 ELSE
606 *
607 * 2 by 2 diagonal block
608 *
609 * Scale if necessary to avoid overflow in forming the
610 * right-hand side element by inner product.
611 *
612 XJ = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ),
613 $ ABS( X( J2 ) )+ABS( X( N+J2 ) ) )
614 IF( XMAX.GT.ONE ) THEN
615 REC = ONE / XMAX
616 IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
617 $ ( BIGNUM-XJ ) / XMAX ) THEN
618 CALL DSCAL( N2, REC, X, 1 )
619 SCALE = SCALE*REC
620 XMAX = XMAX*REC
621 END IF
622 END IF
623 *
624 D( 1, 1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X,
625 $ 1 )
626 D( 2, 1 ) = X( J2 ) - DDOT( J1-1, T( 1, J2 ), 1, X,
627 $ 1 )
628 D( 1, 2 ) = X( N+J1 ) - DDOT( J1-1, T( 1, J1 ), 1,
629 $ X( N+1 ), 1 )
630 D( 2, 2 ) = X( N+J2 ) - DDOT( J1-1, T( 1, J2 ), 1,
631 $ X( N+1 ), 1 )
632 D( 1, 1 ) = D( 1, 1 ) - B( J1 )*X( N+1 )
633 D( 2, 1 ) = D( 2, 1 ) - B( J2 )*X( N+1 )
634 D( 1, 2 ) = D( 1, 2 ) + B( J1 )*X( 1 )
635 D( 2, 2 ) = D( 2, 2 ) + B( J2 )*X( 1 )
636 *
637 CALL DLALN2( .TRUE., 2, 2, SMINW, ONE, T( J1, J1 ),
638 $ LDT, ONE, ONE, D, 2, ZERO, W, V, 2,
639 $ SCALOC, XNORM, IERR )
640 IF( IERR.NE.0 )
641 $ INFO = 2
642 *
643 IF( SCALOC.NE.ONE ) THEN
644 CALL DSCAL( N2, SCALOC, X, 1 )
645 SCALE = SCALOC*SCALE
646 END IF
647 X( J1 ) = V( 1, 1 )
648 X( J2 ) = V( 2, 1 )
649 X( N+J1 ) = V( 1, 2 )
650 X( N+J2 ) = V( 2, 2 )
651 XMAX = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ),
652 $ ABS( X( J2 ) )+ABS( X( N+J2 ) ), XMAX )
653 *
654 END IF
655 *
656 80 CONTINUE
657 *
658 END IF
659 *
660 END IF
661 *
662 RETURN
663 *
664 * End of DLAQTR
665 *
666 END