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: *> \date September 2012
161: *
162: *> \ingroup doubleOTHERauxiliary
163: *
164: * =====================================================================
165: SUBROUTINE DLAQTR( LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK,
166: $ INFO )
167: *
168: * -- LAPACK auxiliary routine (version 3.4.2) --
169: * -- LAPACK is a software package provided by Univ. of Tennessee, --
170: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
171: * September 2012
172: *
173: * .. Scalar Arguments ..
174: LOGICAL LREAL, LTRAN
175: INTEGER INFO, LDT, N
176: DOUBLE PRECISION SCALE, W
177: * ..
178: * .. Array Arguments ..
179: DOUBLE PRECISION B( * ), T( LDT, * ), WORK( * ), X( * )
180: * ..
181: *
182: * =====================================================================
183: *
184: * .. Parameters ..
185: DOUBLE PRECISION ZERO, ONE
186: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
187: * ..
188: * .. Local Scalars ..
189: LOGICAL NOTRAN
190: INTEGER I, IERR, J, J1, J2, JNEXT, K, N1, N2
191: DOUBLE PRECISION BIGNUM, EPS, REC, SCALOC, SI, SMIN, SMINW,
192: $ SMLNUM, SR, TJJ, TMP, XJ, XMAX, XNORM, Z
193: * ..
194: * .. Local Arrays ..
195: DOUBLE PRECISION D( 2, 2 ), V( 2, 2 )
196: * ..
197: * .. External Functions ..
198: INTEGER IDAMAX
199: DOUBLE PRECISION DASUM, DDOT, DLAMCH, DLANGE
200: EXTERNAL IDAMAX, DASUM, DDOT, DLAMCH, DLANGE
201: * ..
202: * .. External Subroutines ..
203: EXTERNAL DAXPY, DLADIV, DLALN2, DSCAL
204: * ..
205: * .. Intrinsic Functions ..
206: INTRINSIC ABS, MAX
207: * ..
208: * .. Executable Statements ..
209: *
210: * Do not test the input parameters for errors
211: *
212: NOTRAN = .NOT.LTRAN
213: INFO = 0
214: *
215: * Quick return if possible
216: *
217: IF( N.EQ.0 )
218: $ RETURN
219: *
220: * Set constants to control overflow
221: *
222: EPS = DLAMCH( 'P' )
223: SMLNUM = DLAMCH( 'S' ) / EPS
224: BIGNUM = ONE / SMLNUM
225: *
226: XNORM = DLANGE( 'M', N, N, T, LDT, D )
227: IF( .NOT.LREAL )
228: $ XNORM = MAX( XNORM, ABS( W ), DLANGE( 'M', N, 1, B, N, D ) )
229: SMIN = MAX( SMLNUM, EPS*XNORM )
230: *
231: * Compute 1-norm of each column of strictly upper triangular
232: * part of T to control overflow in triangular solver.
233: *
234: WORK( 1 ) = ZERO
235: DO 10 J = 2, N
236: WORK( J ) = DASUM( J-1, T( 1, J ), 1 )
237: 10 CONTINUE
238: *
239: IF( .NOT.LREAL ) THEN
240: DO 20 I = 2, N
241: WORK( I ) = WORK( I ) + ABS( B( I ) )
242: 20 CONTINUE
243: END IF
244: *
245: N2 = 2*N
246: N1 = N
247: IF( .NOT.LREAL )
248: $ N1 = N2
249: K = IDAMAX( N1, X, 1 )
250: XMAX = ABS( X( K ) )
251: SCALE = ONE
252: *
253: IF( XMAX.GT.BIGNUM ) THEN
254: SCALE = BIGNUM / XMAX
255: CALL DSCAL( N1, SCALE, X, 1 )
256: XMAX = BIGNUM
257: END IF
258: *
259: IF( LREAL ) THEN
260: *
261: IF( NOTRAN ) THEN
262: *
263: * Solve T*p = scale*c
264: *
265: JNEXT = N
266: DO 30 J = N, 1, -1
267: IF( J.GT.JNEXT )
268: $ GO TO 30
269: J1 = J
270: J2 = J
271: JNEXT = J - 1
272: IF( J.GT.1 ) THEN
273: IF( T( J, J-1 ).NE.ZERO ) THEN
274: J1 = J - 1
275: JNEXT = J - 2
276: END IF
277: END IF
278: *
279: IF( J1.EQ.J2 ) THEN
280: *
281: * Meet 1 by 1 diagonal block
282: *
283: * Scale to avoid overflow when computing
284: * x(j) = b(j)/T(j,j)
285: *
286: XJ = ABS( X( J1 ) )
287: TJJ = ABS( T( J1, J1 ) )
288: TMP = T( J1, J1 )
289: IF( TJJ.LT.SMIN ) THEN
290: TMP = SMIN
291: TJJ = SMIN
292: INFO = 1
293: END IF
294: *
295: IF( XJ.EQ.ZERO )
296: $ GO TO 30
297: *
298: IF( TJJ.LT.ONE ) THEN
299: IF( XJ.GT.BIGNUM*TJJ ) THEN
300: REC = ONE / XJ
301: CALL DSCAL( N, REC, X, 1 )
302: SCALE = SCALE*REC
303: XMAX = XMAX*REC
304: END IF
305: END IF
306: X( J1 ) = X( J1 ) / TMP
307: XJ = ABS( X( J1 ) )
308: *
309: * Scale x if necessary to avoid overflow when adding a
310: * multiple of column j1 of T.
311: *
312: IF( XJ.GT.ONE ) THEN
313: REC = ONE / XJ
314: IF( WORK( J1 ).GT.( BIGNUM-XMAX )*REC ) THEN
315: CALL DSCAL( N, REC, X, 1 )
316: SCALE = SCALE*REC
317: END IF
318: END IF
319: IF( J1.GT.1 ) THEN
320: CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
321: K = IDAMAX( J1-1, X, 1 )
322: XMAX = ABS( X( K ) )
323: END IF
324: *
325: ELSE
326: *
327: * Meet 2 by 2 diagonal block
328: *
329: * Call 2 by 2 linear system solve, to take
330: * care of possible overflow by scaling factor.
331: *
332: D( 1, 1 ) = X( J1 )
333: D( 2, 1 ) = X( J2 )
334: CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, T( J1, J1 ),
335: $ LDT, ONE, ONE, D, 2, ZERO, ZERO, V, 2,
336: $ SCALOC, XNORM, IERR )
337: IF( IERR.NE.0 )
338: $ INFO = 2
339: *
340: IF( SCALOC.NE.ONE ) THEN
341: CALL DSCAL( N, SCALOC, X, 1 )
342: SCALE = SCALE*SCALOC
343: END IF
344: X( J1 ) = V( 1, 1 )
345: X( J2 ) = V( 2, 1 )
346: *
347: * Scale V(1,1) (= X(J1)) and/or V(2,1) (=X(J2))
348: * to avoid overflow in updating right-hand side.
349: *
350: XJ = MAX( ABS( V( 1, 1 ) ), ABS( V( 2, 1 ) ) )
351: IF( XJ.GT.ONE ) THEN
352: REC = ONE / XJ
353: IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
354: $ ( BIGNUM-XMAX )*REC ) THEN
355: CALL DSCAL( N, REC, X, 1 )
356: SCALE = SCALE*REC
357: END IF
358: END IF
359: *
360: * Update right-hand side
361: *
362: IF( J1.GT.1 ) THEN
363: CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
364: CALL DAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 )
365: K = IDAMAX( J1-1, X, 1 )
366: XMAX = ABS( X( K ) )
367: END IF
368: *
369: END IF
370: *
371: 30 CONTINUE
372: *
373: ELSE
374: *
375: * Solve T**T*p = scale*c
376: *
377: JNEXT = 1
378: DO 40 J = 1, N
379: IF( J.LT.JNEXT )
380: $ GO TO 40
381: J1 = J
382: J2 = J
383: JNEXT = J + 1
384: IF( J.LT.N ) THEN
385: IF( T( J+1, J ).NE.ZERO ) THEN
386: J2 = J + 1
387: JNEXT = J + 2
388: END IF
389: END IF
390: *
391: IF( J1.EQ.J2 ) THEN
392: *
393: * 1 by 1 diagonal block
394: *
395: * Scale if necessary to avoid overflow in forming the
396: * right-hand side element by inner product.
397: *
398: XJ = ABS( X( J1 ) )
399: IF( XMAX.GT.ONE ) THEN
400: REC = ONE / XMAX
401: IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN
402: CALL DSCAL( N, REC, X, 1 )
403: SCALE = SCALE*REC
404: XMAX = XMAX*REC
405: END IF
406: END IF
407: *
408: X( J1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X, 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 DSCAL( 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 DSCAL( 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 ) - DDOT( J1-1, T( 1, J1 ), 1, X,
449: $ 1 )
450: D( 2, 1 ) = X( J2 ) - DDOT( J1-1, T( 1, J2 ), 1, X,
451: $ 1 )
452: *
453: CALL DLALN2( .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 DSCAL( 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 DSCAL( N2, REC, X, 1 )
517: SCALE = SCALE*REC
518: XMAX = XMAX*REC
519: END IF
520: END IF
521: CALL DLADIV( 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 DSCAL( N2, REC, X, 1 )
533: SCALE = SCALE*REC
534: END IF
535: END IF
536: *
537: IF( J1.GT.1 ) THEN
538: CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
539: CALL DAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1,
540: $ X( N+1 ), 1 )
541: *
542: X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 )
543: X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 )
544: *
545: XMAX = ZERO
546: DO 50 K = 1, J1 - 1
547: XMAX = MAX( XMAX, ABS( X( K ) )+
548: $ ABS( X( K+N ) ) )
549: 50 CONTINUE
550: END IF
551: *
552: ELSE
553: *
554: * Meet 2 by 2 diagonal block
555: *
556: D( 1, 1 ) = X( J1 )
557: D( 2, 1 ) = X( J2 )
558: D( 1, 2 ) = X( N+J1 )
559: D( 2, 2 ) = X( N+J2 )
560: CALL DLALN2( .FALSE., 2, 2, SMINW, ONE, T( J1, J1 ),
561: $ LDT, ONE, ONE, D, 2, ZERO, -W, V, 2,
562: $ SCALOC, XNORM, IERR )
563: IF( IERR.NE.0 )
564: $ INFO = 2
565: *
566: IF( SCALOC.NE.ONE ) THEN
567: CALL DSCAL( 2*N, SCALOC, X, 1 )
568: SCALE = SCALOC*SCALE
569: END IF
570: X( J1 ) = V( 1, 1 )
571: X( J2 ) = V( 2, 1 )
572: X( N+J1 ) = V( 1, 2 )
573: X( N+J2 ) = V( 2, 2 )
574: *
575: * Scale X(J1), .... to avoid overflow in
576: * updating right hand side.
577: *
578: XJ = MAX( ABS( V( 1, 1 ) )+ABS( V( 1, 2 ) ),
579: $ ABS( V( 2, 1 ) )+ABS( V( 2, 2 ) ) )
580: IF( XJ.GT.ONE ) THEN
581: REC = ONE / XJ
582: IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
583: $ ( BIGNUM-XMAX )*REC ) THEN
584: CALL DSCAL( N2, REC, X, 1 )
585: SCALE = SCALE*REC
586: END IF
587: END IF
588: *
589: * Update the right-hand side.
590: *
591: IF( J1.GT.1 ) THEN
592: CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
593: CALL DAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 )
594: *
595: CALL DAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1,
596: $ X( N+1 ), 1 )
597: CALL DAXPY( J1-1, -X( N+J2 ), T( 1, J2 ), 1,
598: $ X( N+1 ), 1 )
599: *
600: X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 ) +
601: $ B( J2 )*X( N+J2 )
602: X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 ) -
603: $ B( J2 )*X( J2 )
604: *
605: XMAX = ZERO
606: DO 60 K = 1, J1 - 1
607: XMAX = MAX( ABS( X( K ) )+ABS( X( K+N ) ),
608: $ XMAX )
609: 60 CONTINUE
610: END IF
611: *
612: END IF
613: 70 CONTINUE
614: *
615: ELSE
616: *
617: * Solve (T + iB)**T*(p+iq) = c+id
618: *
619: JNEXT = 1
620: DO 80 J = 1, N
621: IF( J.LT.JNEXT )
622: $ GO TO 80
623: J1 = J
624: J2 = J
625: JNEXT = J + 1
626: IF( J.LT.N ) THEN
627: IF( T( J+1, J ).NE.ZERO ) THEN
628: J2 = J + 1
629: JNEXT = J + 2
630: END IF
631: END IF
632: *
633: IF( J1.EQ.J2 ) THEN
634: *
635: * 1 by 1 diagonal block
636: *
637: * Scale if necessary to avoid overflow in forming the
638: * right-hand side element by inner product.
639: *
640: XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) )
641: IF( XMAX.GT.ONE ) THEN
642: REC = ONE / XMAX
643: IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN
644: CALL DSCAL( N2, REC, X, 1 )
645: SCALE = SCALE*REC
646: XMAX = XMAX*REC
647: END IF
648: END IF
649: *
650: X( J1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X, 1 )
651: X( N+J1 ) = X( N+J1 ) - DDOT( J1-1, T( 1, J1 ), 1,
652: $ X( N+1 ), 1 )
653: IF( J1.GT.1 ) THEN
654: X( J1 ) = X( J1 ) - B( J1 )*X( N+1 )
655: X( N+J1 ) = X( N+J1 ) + B( J1 )*X( 1 )
656: END IF
657: XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) )
658: *
659: Z = W
660: IF( J1.EQ.1 )
661: $ Z = B( 1 )
662: *
663: * Scale if necessary to avoid overflow in
664: * complex division
665: *
666: TJJ = ABS( T( J1, J1 ) ) + ABS( Z )
667: TMP = T( J1, J1 )
668: IF( TJJ.LT.SMINW ) THEN
669: TMP = SMINW
670: TJJ = SMINW
671: INFO = 1
672: END IF
673: *
674: IF( TJJ.LT.ONE ) THEN
675: IF( XJ.GT.BIGNUM*TJJ ) THEN
676: REC = ONE / XJ
677: CALL DSCAL( N2, REC, X, 1 )
678: SCALE = SCALE*REC
679: XMAX = XMAX*REC
680: END IF
681: END IF
682: CALL DLADIV( X( J1 ), X( N+J1 ), TMP, -Z, SR, SI )
683: X( J1 ) = SR
684: X( J1+N ) = SI
685: XMAX = MAX( ABS( X( J1 ) )+ABS( X( J1+N ) ), XMAX )
686: *
687: ELSE
688: *
689: * 2 by 2 diagonal block
690: *
691: * Scale if necessary to avoid overflow in forming the
692: * right-hand side element by inner product.
693: *
694: XJ = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ),
695: $ ABS( X( J2 ) )+ABS( X( N+J2 ) ) )
696: IF( XMAX.GT.ONE ) THEN
697: REC = ONE / XMAX
698: IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
699: $ ( BIGNUM-XJ ) / XMAX ) THEN
700: CALL DSCAL( N2, REC, X, 1 )
701: SCALE = SCALE*REC
702: XMAX = XMAX*REC
703: END IF
704: END IF
705: *
706: D( 1, 1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X,
707: $ 1 )
708: D( 2, 1 ) = X( J2 ) - DDOT( J1-1, T( 1, J2 ), 1, X,
709: $ 1 )
710: D( 1, 2 ) = X( N+J1 ) - DDOT( J1-1, T( 1, J1 ), 1,
711: $ X( N+1 ), 1 )
712: D( 2, 2 ) = X( N+J2 ) - DDOT( J1-1, T( 1, J2 ), 1,
713: $ X( N+1 ), 1 )
714: D( 1, 1 ) = D( 1, 1 ) - B( J1 )*X( N+1 )
715: D( 2, 1 ) = D( 2, 1 ) - B( J2 )*X( N+1 )
716: D( 1, 2 ) = D( 1, 2 ) + B( J1 )*X( 1 )
717: D( 2, 2 ) = D( 2, 2 ) + B( J2 )*X( 1 )
718: *
719: CALL DLALN2( .TRUE., 2, 2, SMINW, ONE, T( J1, J1 ),
720: $ LDT, ONE, ONE, D, 2, ZERO, W, V, 2,
721: $ SCALOC, XNORM, IERR )
722: IF( IERR.NE.0 )
723: $ INFO = 2
724: *
725: IF( SCALOC.NE.ONE ) THEN
726: CALL DSCAL( N2, SCALOC, X, 1 )
727: SCALE = SCALOC*SCALE
728: END IF
729: X( J1 ) = V( 1, 1 )
730: X( J2 ) = V( 2, 1 )
731: X( N+J1 ) = V( 1, 2 )
732: X( N+J2 ) = V( 2, 2 )
733: XMAX = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ),
734: $ ABS( X( J2 ) )+ABS( X( N+J2 ) ), XMAX )
735: *
736: END IF
737: *
738: 80 CONTINUE
739: *
740: END IF
741: *
742: END IF
743: *
744: RETURN
745: *
746: * End of DLAQTR
747: *
748: 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