1: SUBROUTINE DLASY2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR,
2: $ LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO )
3: *
4: * -- LAPACK auxiliary routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: LOGICAL LTRANL, LTRANR
11: INTEGER INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2
12: DOUBLE PRECISION SCALE, XNORM
13: * ..
14: * .. Array Arguments ..
15: DOUBLE PRECISION B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ),
16: $ X( LDX, * )
17: * ..
18: *
19: * Purpose
20: * =======
21: *
22: * DLASY2 solves for the N1 by N2 matrix X, 1 <= N1,N2 <= 2, in
23: *
24: * op(TL)*X + ISGN*X*op(TR) = SCALE*B,
25: *
26: * where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or
27: * -1. op(T) = T or T', where T' denotes the transpose of T.
28: *
29: * Arguments
30: * =========
31: *
32: * LTRANL (input) LOGICAL
33: * On entry, LTRANL specifies the op(TL):
34: * = .FALSE., op(TL) = TL,
35: * = .TRUE., op(TL) = TL'.
36: *
37: * LTRANR (input) LOGICAL
38: * On entry, LTRANR specifies the op(TR):
39: * = .FALSE., op(TR) = TR,
40: * = .TRUE., op(TR) = TR'.
41: *
42: * ISGN (input) INTEGER
43: * On entry, ISGN specifies the sign of the equation
44: * as described before. ISGN may only be 1 or -1.
45: *
46: * N1 (input) INTEGER
47: * On entry, N1 specifies the order of matrix TL.
48: * N1 may only be 0, 1 or 2.
49: *
50: * N2 (input) INTEGER
51: * On entry, N2 specifies the order of matrix TR.
52: * N2 may only be 0, 1 or 2.
53: *
54: * TL (input) DOUBLE PRECISION array, dimension (LDTL,2)
55: * On entry, TL contains an N1 by N1 matrix.
56: *
57: * LDTL (input) INTEGER
58: * The leading dimension of the matrix TL. LDTL >= max(1,N1).
59: *
60: * TR (input) DOUBLE PRECISION array, dimension (LDTR,2)
61: * On entry, TR contains an N2 by N2 matrix.
62: *
63: * LDTR (input) INTEGER
64: * The leading dimension of the matrix TR. LDTR >= max(1,N2).
65: *
66: * B (input) DOUBLE PRECISION array, dimension (LDB,2)
67: * On entry, the N1 by N2 matrix B contains the right-hand
68: * side of the equation.
69: *
70: * LDB (input) INTEGER
71: * The leading dimension of the matrix B. LDB >= max(1,N1).
72: *
73: * SCALE (output) DOUBLE PRECISION
74: * On exit, SCALE contains the scale factor. SCALE is chosen
75: * less than or equal to 1 to prevent the solution overflowing.
76: *
77: * X (output) DOUBLE PRECISION array, dimension (LDX,2)
78: * On exit, X contains the N1 by N2 solution.
79: *
80: * LDX (input) INTEGER
81: * The leading dimension of the matrix X. LDX >= max(1,N1).
82: *
83: * XNORM (output) DOUBLE PRECISION
84: * On exit, XNORM is the infinity-norm of the solution.
85: *
86: * INFO (output) INTEGER
87: * On exit, INFO is set to
88: * 0: successful exit.
89: * 1: TL and TR have too close eigenvalues, so TL or
90: * TR is perturbed to get a nonsingular equation.
91: * NOTE: In the interests of speed, this routine does not
92: * check the inputs for errors.
93: *
94: * =====================================================================
95: *
96: * .. Parameters ..
97: DOUBLE PRECISION ZERO, ONE
98: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
99: DOUBLE PRECISION TWO, HALF, EIGHT
100: PARAMETER ( TWO = 2.0D+0, HALF = 0.5D+0, EIGHT = 8.0D+0 )
101: * ..
102: * .. Local Scalars ..
103: LOGICAL BSWAP, XSWAP
104: INTEGER I, IP, IPIV, IPSV, J, JP, JPSV, K
105: DOUBLE PRECISION BET, EPS, GAM, L21, SGN, SMIN, SMLNUM, TAU1,
106: $ TEMP, U11, U12, U22, XMAX
107: * ..
108: * .. Local Arrays ..
109: LOGICAL BSWPIV( 4 ), XSWPIV( 4 )
110: INTEGER JPIV( 4 ), LOCL21( 4 ), LOCU12( 4 ),
111: $ LOCU22( 4 )
112: DOUBLE PRECISION BTMP( 4 ), T16( 4, 4 ), TMP( 4 ), X2( 2 )
113: * ..
114: * .. External Functions ..
115: INTEGER IDAMAX
116: DOUBLE PRECISION DLAMCH
117: EXTERNAL IDAMAX, DLAMCH
118: * ..
119: * .. External Subroutines ..
120: EXTERNAL DCOPY, DSWAP
121: * ..
122: * .. Intrinsic Functions ..
123: INTRINSIC ABS, MAX
124: * ..
125: * .. Data statements ..
126: DATA LOCU12 / 3, 4, 1, 2 / , LOCL21 / 2, 1, 4, 3 / ,
127: $ LOCU22 / 4, 3, 2, 1 /
128: DATA XSWPIV / .FALSE., .FALSE., .TRUE., .TRUE. /
129: DATA BSWPIV / .FALSE., .TRUE., .FALSE., .TRUE. /
130: * ..
131: * .. Executable Statements ..
132: *
133: * Do not check the input parameters for errors
134: *
135: INFO = 0
136: *
137: * Quick return if possible
138: *
139: IF( N1.EQ.0 .OR. N2.EQ.0 )
140: $ RETURN
141: *
142: * Set constants to control overflow
143: *
144: EPS = DLAMCH( 'P' )
145: SMLNUM = DLAMCH( 'S' ) / EPS
146: SGN = ISGN
147: *
148: K = N1 + N1 + N2 - 2
149: GO TO ( 10, 20, 30, 50 )K
150: *
151: * 1 by 1: TL11*X + SGN*X*TR11 = B11
152: *
153: 10 CONTINUE
154: TAU1 = TL( 1, 1 ) + SGN*TR( 1, 1 )
155: BET = ABS( TAU1 )
156: IF( BET.LE.SMLNUM ) THEN
157: TAU1 = SMLNUM
158: BET = SMLNUM
159: INFO = 1
160: END IF
161: *
162: SCALE = ONE
163: GAM = ABS( B( 1, 1 ) )
164: IF( SMLNUM*GAM.GT.BET )
165: $ SCALE = ONE / GAM
166: *
167: X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / TAU1
168: XNORM = ABS( X( 1, 1 ) )
169: RETURN
170: *
171: * 1 by 2:
172: * TL11*[X11 X12] + ISGN*[X11 X12]*op[TR11 TR12] = [B11 B12]
173: * [TR21 TR22]
174: *
175: 20 CONTINUE
176: *
177: SMIN = MAX( EPS*MAX( ABS( TL( 1, 1 ) ), ABS( TR( 1, 1 ) ),
178: $ ABS( TR( 1, 2 ) ), ABS( TR( 2, 1 ) ), ABS( TR( 2, 2 ) ) ),
179: $ SMLNUM )
180: TMP( 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
181: TMP( 4 ) = TL( 1, 1 ) + SGN*TR( 2, 2 )
182: IF( LTRANR ) THEN
183: TMP( 2 ) = SGN*TR( 2, 1 )
184: TMP( 3 ) = SGN*TR( 1, 2 )
185: ELSE
186: TMP( 2 ) = SGN*TR( 1, 2 )
187: TMP( 3 ) = SGN*TR( 2, 1 )
188: END IF
189: BTMP( 1 ) = B( 1, 1 )
190: BTMP( 2 ) = B( 1, 2 )
191: GO TO 40
192: *
193: * 2 by 1:
194: * op[TL11 TL12]*[X11] + ISGN* [X11]*TR11 = [B11]
195: * [TL21 TL22] [X21] [X21] [B21]
196: *
197: 30 CONTINUE
198: SMIN = MAX( EPS*MAX( ABS( TR( 1, 1 ) ), ABS( TL( 1, 1 ) ),
199: $ ABS( TL( 1, 2 ) ), ABS( TL( 2, 1 ) ), ABS( TL( 2, 2 ) ) ),
200: $ SMLNUM )
201: TMP( 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
202: TMP( 4 ) = TL( 2, 2 ) + SGN*TR( 1, 1 )
203: IF( LTRANL ) THEN
204: TMP( 2 ) = TL( 1, 2 )
205: TMP( 3 ) = TL( 2, 1 )
206: ELSE
207: TMP( 2 ) = TL( 2, 1 )
208: TMP( 3 ) = TL( 1, 2 )
209: END IF
210: BTMP( 1 ) = B( 1, 1 )
211: BTMP( 2 ) = B( 2, 1 )
212: 40 CONTINUE
213: *
214: * Solve 2 by 2 system using complete pivoting.
215: * Set pivots less than SMIN to SMIN.
216: *
217: IPIV = IDAMAX( 4, TMP, 1 )
218: U11 = TMP( IPIV )
219: IF( ABS( U11 ).LE.SMIN ) THEN
220: INFO = 1
221: U11 = SMIN
222: END IF
223: U12 = TMP( LOCU12( IPIV ) )
224: L21 = TMP( LOCL21( IPIV ) ) / U11
225: U22 = TMP( LOCU22( IPIV ) ) - U12*L21
226: XSWAP = XSWPIV( IPIV )
227: BSWAP = BSWPIV( IPIV )
228: IF( ABS( U22 ).LE.SMIN ) THEN
229: INFO = 1
230: U22 = SMIN
231: END IF
232: IF( BSWAP ) THEN
233: TEMP = BTMP( 2 )
234: BTMP( 2 ) = BTMP( 1 ) - L21*TEMP
235: BTMP( 1 ) = TEMP
236: ELSE
237: BTMP( 2 ) = BTMP( 2 ) - L21*BTMP( 1 )
238: END IF
239: SCALE = ONE
240: IF( ( TWO*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( U22 ) .OR.
241: $ ( TWO*SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( U11 ) ) THEN
242: SCALE = HALF / MAX( ABS( BTMP( 1 ) ), ABS( BTMP( 2 ) ) )
243: BTMP( 1 ) = BTMP( 1 )*SCALE
244: BTMP( 2 ) = BTMP( 2 )*SCALE
245: END IF
246: X2( 2 ) = BTMP( 2 ) / U22
247: X2( 1 ) = BTMP( 1 ) / U11 - ( U12 / U11 )*X2( 2 )
248: IF( XSWAP ) THEN
249: TEMP = X2( 2 )
250: X2( 2 ) = X2( 1 )
251: X2( 1 ) = TEMP
252: END IF
253: X( 1, 1 ) = X2( 1 )
254: IF( N1.EQ.1 ) THEN
255: X( 1, 2 ) = X2( 2 )
256: XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) )
257: ELSE
258: X( 2, 1 ) = X2( 2 )
259: XNORM = MAX( ABS( X( 1, 1 ) ), ABS( X( 2, 1 ) ) )
260: END IF
261: RETURN
262: *
263: * 2 by 2:
264: * op[TL11 TL12]*[X11 X12] +ISGN* [X11 X12]*op[TR11 TR12] = [B11 B12]
265: * [TL21 TL22] [X21 X22] [X21 X22] [TR21 TR22] [B21 B22]
266: *
267: * Solve equivalent 4 by 4 system using complete pivoting.
268: * Set pivots less than SMIN to SMIN.
269: *
270: 50 CONTINUE
271: SMIN = MAX( ABS( TR( 1, 1 ) ), ABS( TR( 1, 2 ) ),
272: $ ABS( TR( 2, 1 ) ), ABS( TR( 2, 2 ) ) )
273: SMIN = MAX( SMIN, ABS( TL( 1, 1 ) ), ABS( TL( 1, 2 ) ),
274: $ ABS( TL( 2, 1 ) ), ABS( TL( 2, 2 ) ) )
275: SMIN = MAX( EPS*SMIN, SMLNUM )
276: BTMP( 1 ) = ZERO
277: CALL DCOPY( 16, BTMP, 0, T16, 1 )
278: T16( 1, 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
279: T16( 2, 2 ) = TL( 2, 2 ) + SGN*TR( 1, 1 )
280: T16( 3, 3 ) = TL( 1, 1 ) + SGN*TR( 2, 2 )
281: T16( 4, 4 ) = TL( 2, 2 ) + SGN*TR( 2, 2 )
282: IF( LTRANL ) THEN
283: T16( 1, 2 ) = TL( 2, 1 )
284: T16( 2, 1 ) = TL( 1, 2 )
285: T16( 3, 4 ) = TL( 2, 1 )
286: T16( 4, 3 ) = TL( 1, 2 )
287: ELSE
288: T16( 1, 2 ) = TL( 1, 2 )
289: T16( 2, 1 ) = TL( 2, 1 )
290: T16( 3, 4 ) = TL( 1, 2 )
291: T16( 4, 3 ) = TL( 2, 1 )
292: END IF
293: IF( LTRANR ) THEN
294: T16( 1, 3 ) = SGN*TR( 1, 2 )
295: T16( 2, 4 ) = SGN*TR( 1, 2 )
296: T16( 3, 1 ) = SGN*TR( 2, 1 )
297: T16( 4, 2 ) = SGN*TR( 2, 1 )
298: ELSE
299: T16( 1, 3 ) = SGN*TR( 2, 1 )
300: T16( 2, 4 ) = SGN*TR( 2, 1 )
301: T16( 3, 1 ) = SGN*TR( 1, 2 )
302: T16( 4, 2 ) = SGN*TR( 1, 2 )
303: END IF
304: BTMP( 1 ) = B( 1, 1 )
305: BTMP( 2 ) = B( 2, 1 )
306: BTMP( 3 ) = B( 1, 2 )
307: BTMP( 4 ) = B( 2, 2 )
308: *
309: * Perform elimination
310: *
311: DO 100 I = 1, 3
312: XMAX = ZERO
313: DO 70 IP = I, 4
314: DO 60 JP = I, 4
315: IF( ABS( T16( IP, JP ) ).GE.XMAX ) THEN
316: XMAX = ABS( T16( IP, JP ) )
317: IPSV = IP
318: JPSV = JP
319: END IF
320: 60 CONTINUE
321: 70 CONTINUE
322: IF( IPSV.NE.I ) THEN
323: CALL DSWAP( 4, T16( IPSV, 1 ), 4, T16( I, 1 ), 4 )
324: TEMP = BTMP( I )
325: BTMP( I ) = BTMP( IPSV )
326: BTMP( IPSV ) = TEMP
327: END IF
328: IF( JPSV.NE.I )
329: $ CALL DSWAP( 4, T16( 1, JPSV ), 1, T16( 1, I ), 1 )
330: JPIV( I ) = JPSV
331: IF( ABS( T16( I, I ) ).LT.SMIN ) THEN
332: INFO = 1
333: T16( I, I ) = SMIN
334: END IF
335: DO 90 J = I + 1, 4
336: T16( J, I ) = T16( J, I ) / T16( I, I )
337: BTMP( J ) = BTMP( J ) - T16( J, I )*BTMP( I )
338: DO 80 K = I + 1, 4
339: T16( J, K ) = T16( J, K ) - T16( J, I )*T16( I, K )
340: 80 CONTINUE
341: 90 CONTINUE
342: 100 CONTINUE
343: IF( ABS( T16( 4, 4 ) ).LT.SMIN )
344: $ T16( 4, 4 ) = SMIN
345: SCALE = ONE
346: IF( ( EIGHT*SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( T16( 1, 1 ) ) .OR.
347: $ ( EIGHT*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( T16( 2, 2 ) ) .OR.
348: $ ( EIGHT*SMLNUM )*ABS( BTMP( 3 ) ).GT.ABS( T16( 3, 3 ) ) .OR.
349: $ ( EIGHT*SMLNUM )*ABS( BTMP( 4 ) ).GT.ABS( T16( 4, 4 ) ) ) THEN
350: SCALE = ( ONE / EIGHT ) / MAX( ABS( BTMP( 1 ) ),
351: $ ABS( BTMP( 2 ) ), ABS( BTMP( 3 ) ), ABS( BTMP( 4 ) ) )
352: BTMP( 1 ) = BTMP( 1 )*SCALE
353: BTMP( 2 ) = BTMP( 2 )*SCALE
354: BTMP( 3 ) = BTMP( 3 )*SCALE
355: BTMP( 4 ) = BTMP( 4 )*SCALE
356: END IF
357: DO 120 I = 1, 4
358: K = 5 - I
359: TEMP = ONE / T16( K, K )
360: TMP( K ) = BTMP( K )*TEMP
361: DO 110 J = K + 1, 4
362: TMP( K ) = TMP( K ) - ( TEMP*T16( K, J ) )*TMP( J )
363: 110 CONTINUE
364: 120 CONTINUE
365: DO 130 I = 1, 3
366: IF( JPIV( 4-I ).NE.4-I ) THEN
367: TEMP = TMP( 4-I )
368: TMP( 4-I ) = TMP( JPIV( 4-I ) )
369: TMP( JPIV( 4-I ) ) = TEMP
370: END IF
371: 130 CONTINUE
372: X( 1, 1 ) = TMP( 1 )
373: X( 2, 1 ) = TMP( 2 )
374: X( 1, 2 ) = TMP( 3 )
375: X( 2, 2 ) = TMP( 4 )
376: XNORM = MAX( ABS( TMP( 1 ) )+ABS( TMP( 3 ) ),
377: $ ABS( TMP( 2 ) )+ABS( TMP( 4 ) ) )
378: RETURN
379: *
380: * End of DLASY2
381: *
382: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasy2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack