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