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