LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ 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 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 *
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:71
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: