OpenRadioss 2025.1.11
OpenRadioss project
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*> \htmlonly
9*> Download SLAQTR + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaqtr.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaqtr.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaqtr.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SLAQTR( 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* REAL SCALE, W
28* ..
29* .. Array Arguments ..
30* REAL B( * ), T( LDT, * ), WORK( * ), X( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> SLAQTR 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 STRSNA.
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 REAL 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 must
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 REAL 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 REAL
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 REAL
124*> On exit, SCALE is the scale factor.
125*> \endverbatim
126*>
127*> \param[in,out] X
128*> \verbatim
129*> X is REAL 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 REAL 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 SLALN2 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 realOTHERauxiliary
161*
162* =====================================================================
163 SUBROUTINE slaqtr( 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 REAL SCALE, W
174* ..
175* .. Array Arguments ..
176 REAL B( * ), T( LDT, * ), WORK( * ), X( * )
177* ..
178*
179* =====================================================================
180*
181* .. Parameters ..
182 REAL ZERO, ONE
183 parameter( zero = 0.0e+0, one = 1.0e+0 )
184* ..
185* .. Local Scalars ..
186 LOGICAL NOTRAN
187 INTEGER I, IERR, J, J1, J2, JNEXT, K, N1, N2
188 REAL BIGNUM, EPS, REC, SCALOC, SI, SMIN, SMINW,
189 $ smlnum, sr, tjj, tmp, xj, xmax, xnorm, z
190* ..
191* .. Local Arrays ..
192 REAL D( 2, 2 ), V( 2, 2 )
193* ..
194* .. External Functions ..
195 INTEGER ISAMAX
196 REAL SASUM, SDOT, SLAMCH, SLANGE
197 EXTERNAL isamax, sasum, sdot, slamch, slange
198* ..
199* .. External Subroutines ..
200 EXTERNAL saxpy, sladiv, slaln2, sscal
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 = slamch( 'P' )
220 smlnum = slamch( 'S' ) / eps
221 bignum = one / smlnum
222*
223 xnorm = slange( 'M', n, n, t, ldt, d )
224 IF( .NOT.lreal )
225 $ xnorm = max( xnorm, abs( w ), slange( '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 ) = sasum( 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 = isamax( n1, x, 1 )
247 xmax = abs( x( k ) )
248 scale = one
249*
250 IF( xmax.GT.bignum ) THEN
251 scale = bignum / xmax
252 CALL sscal( 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 sscal( 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 sscal( n, rec, x, 1 )
313 scale = scale*rec
314 END IF
315 END IF
316 IF( j1.GT.1 ) THEN
317 CALL saxpy( j1-1, -x( j1 ), t( 1, j1 ), 1, x, 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, 1 )
361 CALL saxpy( j1-1, -x( j2 ), t( 1, j2 ), 1, x, 1 )
362 k = isamax( 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 sscal( n, rec, x, 1 )
400 scale = scale*rec
401 xmax = xmax*rec
402 END IF
403 END IF
404*
405 x( j1 ) = x( j1 ) - sdot( 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 sscal( 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 sscal( 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 ) - sdot( j1-1, t( 1, j1 ), 1, x,
446 $ 1 )
447 d( 2, 1 ) = x( j2 ) - sdot( j1-1, t( 1, j2 ), 1, x,
448 $ 1 )
449*
450 CALL slaln2( .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 sscal( 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 sscal( n2, rec, x, 1 )
514 scale = scale*rec
515 xmax = xmax*rec
516 END IF
517 END IF
518 CALL sladiv( 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 sscal( n2, rec, x, 1 )
530 scale = scale*rec
531 END IF
532 END IF
533*
534 IF( j1.GT.1 ) THEN
535 CALL saxpy( j1-1, -x( j1 ), t( 1, j1 ), 1, x, 1 )
536 CALL saxpy( 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 slaln2( .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 sscal( 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 sscal( 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 saxpy( j1-1, -x( j1 ), t( 1, j1 ), 1, x, 1 )
590 CALL saxpy( j1-1, -x( j2 ), t( 1, j2 ), 1, x, 1 )
591*
592 CALL saxpy( j1-1, -x( n+j1 ), t( 1, j1 ), 1,
593 $ x( n+1 ), 1 )
594 CALL saxpy( 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 sscal( n2, rec, x, 1 )
642 scale = scale*rec
643 xmax = xmax*rec
644 END IF
645 END IF
646*
647 x( j1 ) = x( j1 ) - sdot( j1-1, t( 1, j1 ), 1, x, 1 )
648 x( n+j1 ) = x( n+j1 ) - sdot( 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 sscal( n2, rec, x, 1 )
675 scale = scale*rec
676 xmax = xmax*rec
677 END IF
678 END IF
679 CALL sladiv( 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 sscal( 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 ) - sdot( j1-1, t( 1, j1 ), 1, x,
704 $ 1 )
705 d( 2, 1 ) = x( j2 ) - sdot( j1-1, t( 1, j2 ), 1, x,
706 $ 1 )
707 d( 1, 2 ) = x( n+j1 ) - sdot( j1-1, t( 1, j1 ), 1,
708 $ x( n+1 ), 1 )
709 d( 2, 2 ) = x( n+j2 ) - sdot( 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 slaln2( .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 sscal( 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 SLAQTR
744*
745 END
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:165
sladiv
subroutine sladiv(a, b, c, d, p, q)
SLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition sladiv.f:91
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:218
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
saxpy
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
max
#define max(a, b)
Definition macros.h:21