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: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       LOGICAL            LTRANL, LTRANR
010:       INTEGER            INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2
011:       DOUBLE PRECISION   SCALE, XNORM
012: *     ..
013: *     .. Array Arguments ..
014:       DOUBLE PRECISION   B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ),
015:      $                   X( LDX, * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  DLASY2 solves for the N1 by N2 matrix X, 1 <= N1,N2 <= 2, in
022: *
023: *         op(TL)*X + ISGN*X*op(TR) = SCALE*B,
024: *
025: *  where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or
026: *  -1.  op(T) = T or T', where T' denotes the transpose of T.
027: *
028: *  Arguments
029: *  =========
030: *
031: *  LTRANL  (input) LOGICAL
032: *          On entry, LTRANL specifies the op(TL):
033: *             = .FALSE., op(TL) = TL,
034: *             = .TRUE., op(TL) = TL'.
035: *
036: *  LTRANR  (input) LOGICAL
037: *          On entry, LTRANR specifies the op(TR):
038: *            = .FALSE., op(TR) = TR,
039: *            = .TRUE., op(TR) = TR'.
040: *
041: *  ISGN    (input) INTEGER
042: *          On entry, ISGN specifies the sign of the equation
043: *          as described before. ISGN may only be 1 or -1.
044: *
045: *  N1      (input) INTEGER
046: *          On entry, N1 specifies the order of matrix TL.
047: *          N1 may only be 0, 1 or 2.
048: *
049: *  N2      (input) INTEGER
050: *          On entry, N2 specifies the order of matrix TR.
051: *          N2 may only be 0, 1 or 2.
052: *
053: *  TL      (input) DOUBLE PRECISION array, dimension (LDTL,2)
054: *          On entry, TL contains an N1 by N1 matrix.
055: *
056: *  LDTL    (input) INTEGER
057: *          The leading dimension of the matrix TL. LDTL >= max(1,N1).
058: *
059: *  TR      (input) DOUBLE PRECISION array, dimension (LDTR,2)
060: *          On entry, TR contains an N2 by N2 matrix.
061: *
062: *  LDTR    (input) INTEGER
063: *          The leading dimension of the matrix TR. LDTR >= max(1,N2).
064: *
065: *  B       (input) DOUBLE PRECISION array, dimension (LDB,2)
066: *          On entry, the N1 by N2 matrix B contains the right-hand
067: *          side of the equation.
068: *
069: *  LDB     (input) INTEGER
070: *          The leading dimension of the matrix B. LDB >= max(1,N1).
071: *
072: *  SCALE   (output) DOUBLE PRECISION
073: *          On exit, SCALE contains the scale factor. SCALE is chosen
074: *          less than or equal to 1 to prevent the solution overflowing.
075: *
076: *  X       (output) DOUBLE PRECISION array, dimension (LDX,2)
077: *          On exit, X contains the N1 by N2 solution.
078: *
079: *  LDX     (input) INTEGER
080: *          The leading dimension of the matrix X. LDX >= max(1,N1).
081: *
082: *  XNORM   (output) DOUBLE PRECISION
083: *          On exit, XNORM is the infinity-norm of the solution.
084: *
085: *  INFO    (output) INTEGER
086: *          On exit, INFO is set to
087: *             0: successful exit.
088: *             1: TL and TR have too close eigenvalues, so TL or
089: *                TR is perturbed to get a nonsingular equation.
090: *          NOTE: In the interests of speed, this routine does not
091: *                check the inputs for errors.
092: *
093: * =====================================================================
094: *
095: *     .. Parameters ..
096:       DOUBLE PRECISION   ZERO, ONE
097:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
098:       DOUBLE PRECISION   TWO, HALF, EIGHT
099:       PARAMETER          ( TWO = 2.0D+0, HALF = 0.5D+0, EIGHT = 8.0D+0 )
100: *     ..
101: *     .. Local Scalars ..
102:       LOGICAL            BSWAP, XSWAP
103:       INTEGER            I, IP, IPIV, IPSV, J, JP, JPSV, K
104:       DOUBLE PRECISION   BET, EPS, GAM, L21, SGN, SMIN, SMLNUM, TAU1,
105:      $                   TEMP, U11, U12, U22, XMAX
106: *     ..
107: *     .. Local Arrays ..
108:       LOGICAL            BSWPIV( 4 ), XSWPIV( 4 )
109:       INTEGER            JPIV( 4 ), LOCL21( 4 ), LOCU12( 4 ),
110:      $                   LOCU22( 4 )
111:       DOUBLE PRECISION   BTMP( 4 ), T16( 4, 4 ), TMP( 4 ), X2( 2 )
112: *     ..
113: *     .. External Functions ..
114:       INTEGER            IDAMAX
115:       DOUBLE PRECISION   DLAMCH
116:       EXTERNAL           IDAMAX, DLAMCH
117: *     ..
118: *     .. External Subroutines ..
119:       EXTERNAL           DCOPY, DSWAP
120: *     ..
121: *     .. Intrinsic Functions ..
122:       INTRINSIC          ABS, MAX
123: *     ..
124: *     .. Data statements ..
125:       DATA               LOCU12 / 3, 4, 1, 2 / , LOCL21 / 2, 1, 4, 3 / ,
126:      $                   LOCU22 / 4, 3, 2, 1 /
127:       DATA               XSWPIV / .FALSE., .FALSE., .TRUE., .TRUE. /
128:       DATA               BSWPIV / .FALSE., .TRUE., .FALSE., .TRUE. /
129: *     ..
130: *     .. Executable Statements ..
131: *
132: *     Do not check the input parameters for errors
133: *
134:       INFO = 0
135: *
136: *     Quick return if possible
137: *
138:       IF( N1.EQ.0 .OR. N2.EQ.0 )
139:      $   RETURN
140: *
141: *     Set constants to control overflow
142: *
143:       EPS = DLAMCH( 'P' )
144:       SMLNUM = DLAMCH( 'S' ) / EPS
145:       SGN = ISGN
146: *
147:       K = N1 + N1 + N2 - 2
148:       GO TO ( 10, 20, 30, 50 )K
149: *
150: *     1 by 1: TL11*X + SGN*X*TR11 = B11
151: *
152:    10 CONTINUE
153:       TAU1 = TL( 1, 1 ) + SGN*TR( 1, 1 )
154:       BET = ABS( TAU1 )
155:       IF( BET.LE.SMLNUM ) THEN
156:          TAU1 = SMLNUM
157:          BET = SMLNUM
158:          INFO = 1
159:       END IF
160: *
161:       SCALE = ONE
162:       GAM = ABS( B( 1, 1 ) )
163:       IF( SMLNUM*GAM.GT.BET )
164:      $   SCALE = ONE / GAM
165: *
166:       X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / TAU1
167:       XNORM = ABS( X( 1, 1 ) )
168:       RETURN
169: *
170: *     1 by 2:
171: *     TL11*[X11 X12] + ISGN*[X11 X12]*op[TR11 TR12]  = [B11 B12]
172: *                                       [TR21 TR22]
173: *
174:    20 CONTINUE
175: *
176:       SMIN = MAX( EPS*MAX( ABS( TL( 1, 1 ) ), ABS( TR( 1, 1 ) ),
177:      $       ABS( TR( 1, 2 ) ), ABS( TR( 2, 1 ) ), ABS( TR( 2, 2 ) ) ),
178:      $       SMLNUM )
179:       TMP( 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
180:       TMP( 4 ) = TL( 1, 1 ) + SGN*TR( 2, 2 )
181:       IF( LTRANR ) THEN
182:          TMP( 2 ) = SGN*TR( 2, 1 )
183:          TMP( 3 ) = SGN*TR( 1, 2 )
184:       ELSE
185:          TMP( 2 ) = SGN*TR( 1, 2 )
186:          TMP( 3 ) = SGN*TR( 2, 1 )
187:       END IF
188:       BTMP( 1 ) = B( 1, 1 )
189:       BTMP( 2 ) = B( 1, 2 )
190:       GO TO 40
191: *
192: *     2 by 1:
193: *          op[TL11 TL12]*[X11] + ISGN* [X11]*TR11  = [B11]
194: *            [TL21 TL22] [X21]         [X21]         [B21]
195: *
196:    30 CONTINUE
197:       SMIN = MAX( EPS*MAX( ABS( TR( 1, 1 ) ), ABS( TL( 1, 1 ) ),
198:      $       ABS( TL( 1, 2 ) ), ABS( TL( 2, 1 ) ), ABS( TL( 2, 2 ) ) ),
199:      $       SMLNUM )
200:       TMP( 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
201:       TMP( 4 ) = TL( 2, 2 ) + SGN*TR( 1, 1 )
202:       IF( LTRANL ) THEN
203:          TMP( 2 ) = TL( 1, 2 )
204:          TMP( 3 ) = TL( 2, 1 )
205:       ELSE
206:          TMP( 2 ) = TL( 2, 1 )
207:          TMP( 3 ) = TL( 1, 2 )
208:       END IF
209:       BTMP( 1 ) = B( 1, 1 )
210:       BTMP( 2 ) = B( 2, 1 )
211:    40 CONTINUE
212: *
213: *     Solve 2 by 2 system using complete pivoting.
214: *     Set pivots less than SMIN to SMIN.
215: *
216:       IPIV = IDAMAX( 4, TMP, 1 )
217:       U11 = TMP( IPIV )
218:       IF( ABS( U11 ).LE.SMIN ) THEN
219:          INFO = 1
220:          U11 = SMIN
221:       END IF
222:       U12 = TMP( LOCU12( IPIV ) )
223:       L21 = TMP( LOCL21( IPIV ) ) / U11
224:       U22 = TMP( LOCU22( IPIV ) ) - U12*L21
225:       XSWAP = XSWPIV( IPIV )
226:       BSWAP = BSWPIV( IPIV )
227:       IF( ABS( U22 ).LE.SMIN ) THEN
228:          INFO = 1
229:          U22 = SMIN
230:       END IF
231:       IF( BSWAP ) THEN
232:          TEMP = BTMP( 2 )
233:          BTMP( 2 ) = BTMP( 1 ) - L21*TEMP
234:          BTMP( 1 ) = TEMP
235:       ELSE
236:          BTMP( 2 ) = BTMP( 2 ) - L21*BTMP( 1 )
237:       END IF
238:       SCALE = ONE
239:       IF( ( TWO*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( U22 ) .OR.
240:      $    ( TWO*SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( U11 ) ) THEN
241:          SCALE = HALF / MAX( ABS( BTMP( 1 ) ), ABS( BTMP( 2 ) ) )
242:          BTMP( 1 ) = BTMP( 1 )*SCALE
243:          BTMP( 2 ) = BTMP( 2 )*SCALE
244:       END IF
245:       X2( 2 ) = BTMP( 2 ) / U22
246:       X2( 1 ) = BTMP( 1 ) / U11 - ( U12 / U11 )*X2( 2 )
247:       IF( XSWAP ) THEN
248:          TEMP = X2( 2 )
249:          X2( 2 ) = X2( 1 )
250:          X2( 1 ) = TEMP
251:       END IF
252:       X( 1, 1 ) = X2( 1 )
253:       IF( N1.EQ.1 ) THEN
254:          X( 1, 2 ) = X2( 2 )
255:          XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) )
256:       ELSE
257:          X( 2, 1 ) = X2( 2 )
258:          XNORM = MAX( ABS( X( 1, 1 ) ), ABS( X( 2, 1 ) ) )
259:       END IF
260:       RETURN
261: *
262: *     2 by 2:
263: *     op[TL11 TL12]*[X11 X12] +ISGN* [X11 X12]*op[TR11 TR12] = [B11 B12]
264: *       [TL21 TL22] [X21 X22]        [X21 X22]   [TR21 TR22]   [B21 B22]
265: *
266: *     Solve equivalent 4 by 4 system using complete pivoting.
267: *     Set pivots less than SMIN to SMIN.
268: *
269:    50 CONTINUE
270:       SMIN = MAX( ABS( TR( 1, 1 ) ), ABS( TR( 1, 2 ) ),
271:      $       ABS( TR( 2, 1 ) ), ABS( TR( 2, 2 ) ) )
272:       SMIN = MAX( SMIN, ABS( TL( 1, 1 ) ), ABS( TL( 1, 2 ) ),
273:      $       ABS( TL( 2, 1 ) ), ABS( TL( 2, 2 ) ) )
274:       SMIN = MAX( EPS*SMIN, SMLNUM )
275:       BTMP( 1 ) = ZERO
276:       CALL DCOPY( 16, BTMP, 0, T16, 1 )
277:       T16( 1, 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
278:       T16( 2, 2 ) = TL( 2, 2 ) + SGN*TR( 1, 1 )
279:       T16( 3, 3 ) = TL( 1, 1 ) + SGN*TR( 2, 2 )
280:       T16( 4, 4 ) = TL( 2, 2 ) + SGN*TR( 2, 2 )
281:       IF( LTRANL ) THEN
282:          T16( 1, 2 ) = TL( 2, 1 )
283:          T16( 2, 1 ) = TL( 1, 2 )
284:          T16( 3, 4 ) = TL( 2, 1 )
285:          T16( 4, 3 ) = TL( 1, 2 )
286:       ELSE
287:          T16( 1, 2 ) = TL( 1, 2 )
288:          T16( 2, 1 ) = TL( 2, 1 )
289:          T16( 3, 4 ) = TL( 1, 2 )
290:          T16( 4, 3 ) = TL( 2, 1 )
291:       END IF
292:       IF( LTRANR ) THEN
293:          T16( 1, 3 ) = SGN*TR( 1, 2 )
294:          T16( 2, 4 ) = SGN*TR( 1, 2 )
295:          T16( 3, 1 ) = SGN*TR( 2, 1 )
296:          T16( 4, 2 ) = SGN*TR( 2, 1 )
297:       ELSE
298:          T16( 1, 3 ) = SGN*TR( 2, 1 )
299:          T16( 2, 4 ) = SGN*TR( 2, 1 )
300:          T16( 3, 1 ) = SGN*TR( 1, 2 )
301:          T16( 4, 2 ) = SGN*TR( 1, 2 )
302:       END IF
303:       BTMP( 1 ) = B( 1, 1 )
304:       BTMP( 2 ) = B( 2, 1 )
305:       BTMP( 3 ) = B( 1, 2 )
306:       BTMP( 4 ) = B( 2, 2 )
307: *
308: *     Perform elimination
309: *
310:       DO 100 I = 1, 3
311:          XMAX = ZERO
312:          DO 70 IP = I, 4
313:             DO 60 JP = I, 4
314:                IF( ABS( T16( IP, JP ) ).GE.XMAX ) THEN
315:                   XMAX = ABS( T16( IP, JP ) )
316:                   IPSV = IP
317:                   JPSV = JP
318:                END IF
319:    60       CONTINUE
320:    70    CONTINUE
321:          IF( IPSV.NE.I ) THEN
322:             CALL DSWAP( 4, T16( IPSV, 1 ), 4, T16( I, 1 ), 4 )
323:             TEMP = BTMP( I )
324:             BTMP( I ) = BTMP( IPSV )
325:             BTMP( IPSV ) = TEMP
326:          END IF
327:          IF( JPSV.NE.I )
328:      $      CALL DSWAP( 4, T16( 1, JPSV ), 1, T16( 1, I ), 1 )
329:          JPIV( I ) = JPSV
330:          IF( ABS( T16( I, I ) ).LT.SMIN ) THEN
331:             INFO = 1
332:             T16( I, I ) = SMIN
333:          END IF
334:          DO 90 J = I + 1, 4
335:             T16( J, I ) = T16( J, I ) / T16( I, I )
336:             BTMP( J ) = BTMP( J ) - T16( J, I )*BTMP( I )
337:             DO 80 K = I + 1, 4
338:                T16( J, K ) = T16( J, K ) - T16( J, I )*T16( I, K )
339:    80       CONTINUE
340:    90    CONTINUE
341:   100 CONTINUE
342:       IF( ABS( T16( 4, 4 ) ).LT.SMIN )
343:      $   T16( 4, 4 ) = SMIN
344:       SCALE = ONE
345:       IF( ( EIGHT*SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( T16( 1, 1 ) ) .OR.
346:      $    ( EIGHT*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( T16( 2, 2 ) ) .OR.
347:      $    ( EIGHT*SMLNUM )*ABS( BTMP( 3 ) ).GT.ABS( T16( 3, 3 ) ) .OR.
348:      $    ( EIGHT*SMLNUM )*ABS( BTMP( 4 ) ).GT.ABS( T16( 4, 4 ) ) ) THEN
349:          SCALE = ( ONE / EIGHT ) / MAX( ABS( BTMP( 1 ) ),
350:      $           ABS( BTMP( 2 ) ), ABS( BTMP( 3 ) ), ABS( BTMP( 4 ) ) )
351:          BTMP( 1 ) = BTMP( 1 )*SCALE
352:          BTMP( 2 ) = BTMP( 2 )*SCALE
353:          BTMP( 3 ) = BTMP( 3 )*SCALE
354:          BTMP( 4 ) = BTMP( 4 )*SCALE
355:       END IF
356:       DO 120 I = 1, 4
357:          K = 5 - I
358:          TEMP = ONE / T16( K, K )
359:          TMP( K ) = BTMP( K )*TEMP
360:          DO 110 J = K + 1, 4
361:             TMP( K ) = TMP( K ) - ( TEMP*T16( K, J ) )*TMP( J )
362:   110    CONTINUE
363:   120 CONTINUE
364:       DO 130 I = 1, 3
365:          IF( JPIV( 4-I ).NE.4-I ) THEN
366:             TEMP = TMP( 4-I )
367:             TMP( 4-I ) = TMP( JPIV( 4-I ) )
368:             TMP( JPIV( 4-I ) ) = TEMP
369:          END IF
370:   130 CONTINUE
371:       X( 1, 1 ) = TMP( 1 )
372:       X( 2, 1 ) = TMP( 2 )
373:       X( 1, 2 ) = TMP( 3 )
374:       X( 2, 2 ) = TMP( 4 )
375:       XNORM = MAX( ABS( TMP( 1 ) )+ABS( TMP( 3 ) ),
376:      $        ABS( TMP( 2 ) )+ABS( TMP( 4 ) ) )
377:       RETURN
378: *
379: *     End of DLASY2
380: *
381:       END
382: