LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dlasy2 ( logical  LTRANL,
logical  LTRANR,
integer  ISGN,
integer  N1,
integer  N2,
double precision, dimension( ldtl, * )  TL,
integer  LDTL,
double precision, dimension( ldtr, * )  TR,
integer  LDTR,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision  SCALE,
double precision, dimension( ldx, * )  X,
integer  LDX,
double precision  XNORM,
integer  INFO 
)

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

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

Purpose: 
 DLASY2 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
          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 DOUBLE PRECISION 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 DOUBLE PRECISION
          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.
Date
June 2016

Definition at line 176 of file dlasy2.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: