LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ slasy2()

 subroutine slasy2 ( logical ltranl, logical ltranr, integer isgn, integer n1, integer n2, real, dimension( ldtl, * ) tl, integer ldtl, real, dimension( ldtr, * ) tr, integer ldtr, real, dimension( ldb, * ) b, integer ldb, real scale, real, dimension( ldx, * ) x, integer ldx, real xnorm, integer info )

SLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2.

Purpose:
``` SLASY2 solves for the N1 by N2 matrix X, 1 <= N1,N2 <= 2, in

op(TL)*X + ISGN*X*op(TR) = SCALE*B,

where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or
-1.  op(T) = T or T**T, where T**T denotes the transpose of T.```
Parameters
 [in] LTRANL ``` LTRANL is LOGICAL On entry, LTRANL specifies the op(TL): = .FALSE., op(TL) = TL, = .TRUE., op(TL) = TL**T.``` [in] LTRANR ``` LTRANR is LOGICAL On entry, LTRANR specifies the op(TR): = .FALSE., op(TR) = TR, = .TRUE., op(TR) = TR**T.``` [in] ISGN ``` ISGN is INTEGER On entry, ISGN specifies the sign of the equation as described before. ISGN may only be 1 or -1.``` [in] N1 ``` N1 is INTEGER On entry, N1 specifies the order of matrix TL. N1 may only be 0, 1 or 2.``` [in] N2 ``` N2 is INTEGER On entry, N2 specifies the order of matrix TR. N2 may only be 0, 1 or 2.``` [in] TL ``` TL is REAL array, dimension (LDTL,2) On entry, TL contains an N1 by N1 matrix.``` [in] LDTL ``` LDTL is INTEGER The leading dimension of the matrix TL. LDTL >= max(1,N1).``` [in] TR ``` TR is REAL array, dimension (LDTR,2) On entry, TR contains an N2 by N2 matrix.``` [in] LDTR ``` LDTR is INTEGER The leading dimension of the matrix TR. LDTR >= max(1,N2).``` [in] B ``` B is REAL array, dimension (LDB,2) On entry, the N1 by N2 matrix B contains the right-hand side of the equation.``` [in] LDB ``` LDB is INTEGER The leading dimension of the matrix B. LDB >= max(1,N1).``` [out] SCALE ``` SCALE is REAL On exit, SCALE contains the scale factor. SCALE is chosen less than or equal to 1 to prevent the solution overflowing.``` [out] X ``` X is REAL array, dimension (LDX,2) On exit, X contains the N1 by N2 solution.``` [in] LDX ``` LDX is INTEGER The leading dimension of the matrix X. LDX >= max(1,N1).``` [out] XNORM ``` XNORM is REAL On exit, XNORM is the infinity-norm of the solution.``` [out] INFO ``` INFO is INTEGER On exit, INFO is set to 0: successful exit. 1: TL and TR have too close eigenvalues, so TL or TR is perturbed to get a nonsingular equation. NOTE: In the interests of speed, this routine does not check the inputs for errors.```

Definition at line 172 of file slasy2.f.

174*
175* -- LAPACK auxiliary routine --
176* -- LAPACK is a software package provided by Univ. of Tennessee, --
177* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
178*
179* .. Scalar Arguments ..
180 LOGICAL LTRANL, LTRANR
181 INTEGER INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2
182 REAL SCALE, XNORM
183* ..
184* .. Array Arguments ..
185 REAL B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ),
186 \$ X( LDX, * )
187* ..
188*
189* =====================================================================
190*
191* .. Parameters ..
192 REAL ZERO, ONE
193 parameter( zero = 0.0e+0, one = 1.0e+0 )
194 REAL TWO, HALF, EIGHT
195 parameter( two = 2.0e+0, half = 0.5e+0, eight = 8.0e+0 )
196* ..
197* .. Local Scalars ..
198 LOGICAL BSWAP, XSWAP
199 INTEGER I, IP, IPIV, IPSV, J, JP, JPSV, K
200 REAL BET, EPS, GAM, L21, SGN, SMIN, SMLNUM, TAU1,
201 \$ TEMP, U11, U12, U22, XMAX
202* ..
203* .. Local Arrays ..
204 LOGICAL BSWPIV( 4 ), XSWPIV( 4 )
205 INTEGER JPIV( 4 ), LOCL21( 4 ), LOCU12( 4 ),
206 \$ LOCU22( 4 )
207 REAL BTMP( 4 ), T16( 4, 4 ), TMP( 4 ), X2( 2 )
208* ..
209* .. External Functions ..
210 INTEGER ISAMAX
211 REAL SLAMCH
212 EXTERNAL isamax, slamch
213* ..
214* .. External Subroutines ..
215 EXTERNAL scopy, sswap
216* ..
217* .. Intrinsic Functions ..
218 INTRINSIC abs, max
219* ..
220* .. Data statements ..
221 DATA locu12 / 3, 4, 1, 2 / , locl21 / 2, 1, 4, 3 / ,
222 \$ locu22 / 4, 3, 2, 1 /
223 DATA xswpiv / .false., .false., .true., .true. /
224 DATA bswpiv / .false., .true., .false., .true. /
225* ..
226* .. Executable Statements ..
227*
228* Do not check the input parameters for errors
229*
230 info = 0
231*
232* Quick return if possible
233*
234 IF( n1.EQ.0 .OR. n2.EQ.0 )
235 \$ RETURN
236*
237* Set constants to control overflow
238*
239 eps = slamch( 'P' )
240 smlnum = slamch( 'S' ) / eps
241 sgn = isgn
242*
243 k = n1 + n1 + n2 - 2
244 GO TO ( 10, 20, 30, 50 )k
245*
246* 1 by 1: TL11*X + SGN*X*TR11 = B11
247*
248 10 CONTINUE
249 tau1 = tl( 1, 1 ) + sgn*tr( 1, 1 )
250 bet = abs( tau1 )
251 IF( bet.LE.smlnum ) THEN
252 tau1 = smlnum
253 bet = smlnum
254 info = 1
255 END IF
256*
257 scale = one
258 gam = abs( b( 1, 1 ) )
259 IF( smlnum*gam.GT.bet )
260 \$ scale = one / gam
261*
262 x( 1, 1 ) = ( b( 1, 1 )*scale ) / tau1
263 xnorm = abs( x( 1, 1 ) )
264 RETURN
265*
266* 1 by 2:
267* TL11*[X11 X12] + ISGN*[X11 X12]*op[TR11 TR12] = [B11 B12]
268* [TR21 TR22]
269*
270 20 CONTINUE
271*
272 smin = max( eps*max( abs( tl( 1, 1 ) ), abs( tr( 1, 1 ) ),
273 \$ abs( tr( 1, 2 ) ), abs( tr( 2, 1 ) ), abs( tr( 2, 2 ) ) ),
274 \$ smlnum )
275 tmp( 1 ) = tl( 1, 1 ) + sgn*tr( 1, 1 )
276 tmp( 4 ) = tl( 1, 1 ) + sgn*tr( 2, 2 )
277 IF( ltranr ) THEN
278 tmp( 2 ) = sgn*tr( 2, 1 )
279 tmp( 3 ) = sgn*tr( 1, 2 )
280 ELSE
281 tmp( 2 ) = sgn*tr( 1, 2 )
282 tmp( 3 ) = sgn*tr( 2, 1 )
283 END IF
284 btmp( 1 ) = b( 1, 1 )
285 btmp( 2 ) = b( 1, 2 )
286 GO TO 40
287*
288* 2 by 1:
289* op[TL11 TL12]*[X11] + ISGN* [X11]*TR11 = [B11]
290* [TL21 TL22] [X21] [X21] [B21]
291*
292 30 CONTINUE
293 smin = max( eps*max( abs( tr( 1, 1 ) ), abs( tl( 1, 1 ) ),
294 \$ abs( tl( 1, 2 ) ), abs( tl( 2, 1 ) ), abs( tl( 2, 2 ) ) ),
295 \$ smlnum )
296 tmp( 1 ) = tl( 1, 1 ) + sgn*tr( 1, 1 )
297 tmp( 4 ) = tl( 2, 2 ) + sgn*tr( 1, 1 )
298 IF( ltranl ) THEN
299 tmp( 2 ) = tl( 1, 2 )
300 tmp( 3 ) = tl( 2, 1 )
301 ELSE
302 tmp( 2 ) = tl( 2, 1 )
303 tmp( 3 ) = tl( 1, 2 )
304 END IF
305 btmp( 1 ) = b( 1, 1 )
306 btmp( 2 ) = b( 2, 1 )
307 40 CONTINUE
308*
309* Solve 2 by 2 system using complete pivoting.
310* Set pivots less than SMIN to SMIN.
311*
312 ipiv = isamax( 4, tmp, 1 )
313 u11 = tmp( ipiv )
314 IF( abs( u11 ).LE.smin ) THEN
315 info = 1
316 u11 = smin
317 END IF
318 u12 = tmp( locu12( ipiv ) )
319 l21 = tmp( locl21( ipiv ) ) / u11
320 u22 = tmp( locu22( ipiv ) ) - u12*l21
321 xswap = xswpiv( ipiv )
322 bswap = bswpiv( ipiv )
323 IF( abs( u22 ).LE.smin ) THEN
324 info = 1
325 u22 = smin
326 END IF
327 IF( bswap ) THEN
328 temp = btmp( 2 )
329 btmp( 2 ) = btmp( 1 ) - l21*temp
330 btmp( 1 ) = temp
331 ELSE
332 btmp( 2 ) = btmp( 2 ) - l21*btmp( 1 )
333 END IF
334 scale = one
335 IF( ( two*smlnum )*abs( btmp( 2 ) ).GT.abs( u22 ) .OR.
336 \$ ( two*smlnum )*abs( btmp( 1 ) ).GT.abs( u11 ) ) THEN
337 scale = half / max( abs( btmp( 1 ) ), abs( btmp( 2 ) ) )
338 btmp( 1 ) = btmp( 1 )*scale
339 btmp( 2 ) = btmp( 2 )*scale
340 END IF
341 x2( 2 ) = btmp( 2 ) / u22
342 x2( 1 ) = btmp( 1 ) / u11 - ( u12 / u11 )*x2( 2 )
343 IF( xswap ) THEN
344 temp = x2( 2 )
345 x2( 2 ) = x2( 1 )
346 x2( 1 ) = temp
347 END IF
348 x( 1, 1 ) = x2( 1 )
349 IF( n1.EQ.1 ) THEN
350 x( 1, 2 ) = x2( 2 )
351 xnorm = abs( x( 1, 1 ) ) + abs( x( 1, 2 ) )
352 ELSE
353 x( 2, 1 ) = x2( 2 )
354 xnorm = max( abs( x( 1, 1 ) ), abs( x( 2, 1 ) ) )
355 END IF
356 RETURN
357*
358* 2 by 2:
359* op[TL11 TL12]*[X11 X12] +ISGN* [X11 X12]*op[TR11 TR12] = [B11 B12]
360* [TL21 TL22] [X21 X22] [X21 X22] [TR21 TR22] [B21 B22]
361*
362* Solve equivalent 4 by 4 system using complete pivoting.
363* Set pivots less than SMIN to SMIN.
364*
365 50 CONTINUE
366 smin = max( abs( tr( 1, 1 ) ), abs( tr( 1, 2 ) ),
367 \$ abs( tr( 2, 1 ) ), abs( tr( 2, 2 ) ) )
368 smin = max( smin, abs( tl( 1, 1 ) ), abs( tl( 1, 2 ) ),
369 \$ abs( tl( 2, 1 ) ), abs( tl( 2, 2 ) ) )
370 smin = max( eps*smin, smlnum )
371 btmp( 1 ) = zero
372 CALL scopy( 16, btmp, 0, t16, 1 )
373 t16( 1, 1 ) = tl( 1, 1 ) + sgn*tr( 1, 1 )
374 t16( 2, 2 ) = tl( 2, 2 ) + sgn*tr( 1, 1 )
375 t16( 3, 3 ) = tl( 1, 1 ) + sgn*tr( 2, 2 )
376 t16( 4, 4 ) = tl( 2, 2 ) + sgn*tr( 2, 2 )
377 IF( ltranl ) THEN
378 t16( 1, 2 ) = tl( 2, 1 )
379 t16( 2, 1 ) = tl( 1, 2 )
380 t16( 3, 4 ) = tl( 2, 1 )
381 t16( 4, 3 ) = tl( 1, 2 )
382 ELSE
383 t16( 1, 2 ) = tl( 1, 2 )
384 t16( 2, 1 ) = tl( 2, 1 )
385 t16( 3, 4 ) = tl( 1, 2 )
386 t16( 4, 3 ) = tl( 2, 1 )
387 END IF
388 IF( ltranr ) THEN
389 t16( 1, 3 ) = sgn*tr( 1, 2 )
390 t16( 2, 4 ) = sgn*tr( 1, 2 )
391 t16( 3, 1 ) = sgn*tr( 2, 1 )
392 t16( 4, 2 ) = sgn*tr( 2, 1 )
393 ELSE
394 t16( 1, 3 ) = sgn*tr( 2, 1 )
395 t16( 2, 4 ) = sgn*tr( 2, 1 )
396 t16( 3, 1 ) = sgn*tr( 1, 2 )
397 t16( 4, 2 ) = sgn*tr( 1, 2 )
398 END IF
399 btmp( 1 ) = b( 1, 1 )
400 btmp( 2 ) = b( 2, 1 )
401 btmp( 3 ) = b( 1, 2 )
402 btmp( 4 ) = b( 2, 2 )
403*
404* Perform elimination
405*
406 DO 100 i = 1, 3
407 xmax = zero
408 DO 70 ip = i, 4
409 DO 60 jp = i, 4
410 IF( abs( t16( ip, jp ) ).GE.xmax ) THEN
411 xmax = abs( t16( ip, jp ) )
412 ipsv = ip
413 jpsv = jp
414 END IF
415 60 CONTINUE
416 70 CONTINUE
417 IF( ipsv.NE.i ) THEN
418 CALL sswap( 4, t16( ipsv, 1 ), 4, t16( i, 1 ), 4 )
419 temp = btmp( i )
420 btmp( i ) = btmp( ipsv )
421 btmp( ipsv ) = temp
422 END IF
423 IF( jpsv.NE.i )
424 \$ CALL sswap( 4, t16( 1, jpsv ), 1, t16( 1, i ), 1 )
425 jpiv( i ) = jpsv
426 IF( abs( t16( i, i ) ).LT.smin ) THEN
427 info = 1
428 t16( i, i ) = smin
429 END IF
430 DO 90 j = i + 1, 4
431 t16( j, i ) = t16( j, i ) / t16( i, i )
432 btmp( j ) = btmp( j ) - t16( j, i )*btmp( i )
433 DO 80 k = i + 1, 4
434 t16( j, k ) = t16( j, k ) - t16( j, i )*t16( i, k )
435 80 CONTINUE
436 90 CONTINUE
437 100 CONTINUE
438 IF( abs( t16( 4, 4 ) ).LT.smin ) THEN
439 info = 1
440 t16( 4, 4 ) = smin
441 END IF
442 scale = one
443 IF( ( eight*smlnum )*abs( btmp( 1 ) ).GT.abs( t16( 1, 1 ) ) .OR.
444 \$ ( eight*smlnum )*abs( btmp( 2 ) ).GT.abs( t16( 2, 2 ) ) .OR.
445 \$ ( eight*smlnum )*abs( btmp( 3 ) ).GT.abs( t16( 3, 3 ) ) .OR.
446 \$ ( eight*smlnum )*abs( btmp( 4 ) ).GT.abs( t16( 4, 4 ) ) ) THEN
447 scale = ( one / eight ) / max( abs( btmp( 1 ) ),
448 \$ abs( btmp( 2 ) ), abs( btmp( 3 ) ), abs( btmp( 4 ) ) )
449 btmp( 1 ) = btmp( 1 )*scale
450 btmp( 2 ) = btmp( 2 )*scale
451 btmp( 3 ) = btmp( 3 )*scale
452 btmp( 4 ) = btmp( 4 )*scale
453 END IF
454 DO 120 i = 1, 4
455 k = 5 - i
456 temp = one / t16( k, k )
457 tmp( k ) = btmp( k )*temp
458 DO 110 j = k + 1, 4
459 tmp( k ) = tmp( k ) - ( temp*t16( k, j ) )*tmp( j )
460 110 CONTINUE
461 120 CONTINUE
462 DO 130 i = 1, 3
463 IF( jpiv( 4-i ).NE.4-i ) THEN
464 temp = tmp( 4-i )
465 tmp( 4-i ) = tmp( jpiv( 4-i ) )
466 tmp( jpiv( 4-i ) ) = temp
467 END IF
468 130 CONTINUE
469 x( 1, 1 ) = tmp( 1 )
470 x( 2, 1 ) = tmp( 2 )
471 x( 1, 2 ) = tmp( 3 )
472 x( 2, 2 ) = tmp( 4 )
473 xnorm = max( abs( tmp( 1 ) )+abs( tmp( 3 ) ),
474 \$ abs( tmp( 2 ) )+abs( tmp( 4 ) ) )
475 RETURN
476*
477* End of SLASY2
478*
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: