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