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 doubleOTHERauxiliary
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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaqtr.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack