001:       SUBROUTINE DLASY2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR,
002:      $                   LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO )
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       LOGICAL            LTRANL, LTRANR
011:       INTEGER            INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2
012:       DOUBLE PRECISION   SCALE, XNORM
013: *     ..
014: *     .. Array Arguments ..
015:       DOUBLE PRECISION   B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ),
016:      $                   X( LDX, * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  DLASY2 solves for the N1 by N2 matrix X, 1 <= N1,N2 <= 2, in
023: *
024: *         op(TL)*X + ISGN*X*op(TR) = SCALE*B,
025: *
026: *  where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or
027: *  -1.  op(T) = T or T', where T' denotes the transpose of T.
028: *
029: *  Arguments
030: *  =========
031: *
032: *  LTRANL  (input) LOGICAL
033: *          On entry, LTRANL specifies the op(TL):
034: *             = .FALSE., op(TL) = TL,
035: *             = .TRUE., op(TL) = TL'.
036: *
037: *  LTRANR  (input) LOGICAL
038: *          On entry, LTRANR specifies the op(TR):
039: *            = .FALSE., op(TR) = TR,
040: *            = .TRUE., op(TR) = TR'.
041: *
042: *  ISGN    (input) INTEGER
043: *          On entry, ISGN specifies the sign of the equation
044: *          as described before. ISGN may only be 1 or -1.
045: *
046: *  N1      (input) INTEGER
047: *          On entry, N1 specifies the order of matrix TL.
048: *          N1 may only be 0, 1 or 2.
049: *
050: *  N2      (input) INTEGER
051: *          On entry, N2 specifies the order of matrix TR.
052: *          N2 may only be 0, 1 or 2.
053: *
054: *  TL      (input) DOUBLE PRECISION array, dimension (LDTL,2)
055: *          On entry, TL contains an N1 by N1 matrix.
056: *
057: *  LDTL    (input) INTEGER
058: *          The leading dimension of the matrix TL. LDTL >= max(1,N1).
059: *
060: *  TR      (input) DOUBLE PRECISION array, dimension (LDTR,2)
061: *          On entry, TR contains an N2 by N2 matrix.
062: *
063: *  LDTR    (input) INTEGER
064: *          The leading dimension of the matrix TR. LDTR >= max(1,N2).
065: *
066: *  B       (input) DOUBLE PRECISION array, dimension (LDB,2)
067: *          On entry, the N1 by N2 matrix B contains the right-hand
068: *          side of the equation.
069: *
070: *  LDB     (input) INTEGER
071: *          The leading dimension of the matrix B. LDB >= max(1,N1).
072: *
073: *  SCALE   (output) DOUBLE PRECISION
074: *          On exit, SCALE contains the scale factor. SCALE is chosen
075: *          less than or equal to 1 to prevent the solution overflowing.
076: *
077: *  X       (output) DOUBLE PRECISION array, dimension (LDX,2)
078: *          On exit, X contains the N1 by N2 solution.
079: *
080: *  LDX     (input) INTEGER
081: *          The leading dimension of the matrix X. LDX >= max(1,N1).
082: *
083: *  XNORM   (output) DOUBLE PRECISION
084: *          On exit, XNORM is the infinity-norm of the solution.
085: *
086: *  INFO    (output) INTEGER
087: *          On exit, INFO is set to
088: *             0: successful exit.
089: *             1: TL and TR have too close eigenvalues, so TL or
090: *                TR is perturbed to get a nonsingular equation.
091: *          NOTE: In the interests of speed, this routine does not
092: *                check the inputs for errors.
093: *
094: * =====================================================================
095: *
096: *     .. Parameters ..
097:       DOUBLE PRECISION   ZERO, ONE
098:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
099:       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
383: