LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
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.

Download SLASY2 + dependencies [TGZ] [ZIP] [TXT]

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.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 170 of file slasy2.f.

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