LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
slasy2.f
Go to the documentation of this file.
1 *> \brief \b SLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SLASY2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasy2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasy2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasy2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLASY2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR,
22 * LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO )
23 *
24 * .. Scalar Arguments ..
25 * LOGICAL LTRANL, LTRANR
26 * INTEGER INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2
27 * REAL SCALE, XNORM
28 * ..
29 * .. Array Arguments ..
30 * REAL B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ),
31 * $ X( LDX, * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> SLASY2 solves for the N1 by N2 matrix X, 1 <= N1,N2 <= 2, in
41 *>
42 *> op(TL)*X + ISGN*X*op(TR) = SCALE*B,
43 *>
44 *> where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or
45 *> -1. op(T) = T or T**T, where T**T denotes the transpose of T.
46 *> \endverbatim
47 *
48 * Arguments:
49 * ==========
50 *
51 *> \param[in] LTRANL
52 *> \verbatim
53 *> LTRANL is LOGICAL
54 *> On entry, LTRANL specifies the op(TL):
55 *> = .FALSE., op(TL) = TL,
56 *> = .TRUE., op(TL) = TL**T.
57 *> \endverbatim
58 *>
59 *> \param[in] LTRANR
60 *> \verbatim
61 *> LTRANR is LOGICAL
62 *> On entry, LTRANR specifies the op(TR):
63 *> = .FALSE., op(TR) = TR,
64 *> = .TRUE., op(TR) = TR**T.
65 *> \endverbatim
66 *>
67 *> \param[in] ISGN
68 *> \verbatim
69 *> ISGN is INTEGER
70 *> On entry, ISGN specifies the sign of the equation
71 *> as described before. ISGN may only be 1 or -1.
72 *> \endverbatim
73 *>
74 *> \param[in] N1
75 *> \verbatim
76 *> N1 is INTEGER
77 *> On entry, N1 specifies the order of matrix TL.
78 *> N1 may only be 0, 1 or 2.
79 *> \endverbatim
80 *>
81 *> \param[in] N2
82 *> \verbatim
83 *> N2 is INTEGER
84 *> On entry, N2 specifies the order of matrix TR.
85 *> N2 may only be 0, 1 or 2.
86 *> \endverbatim
87 *>
88 *> \param[in] TL
89 *> \verbatim
90 *> TL is REAL array, dimension (LDTL,2)
91 *> On entry, TL contains an N1 by N1 matrix.
92 *> \endverbatim
93 *>
94 *> \param[in] LDTL
95 *> \verbatim
96 *> LDTL is INTEGER
97 *> The leading dimension of the matrix TL. LDTL >= max(1,N1).
98 *> \endverbatim
99 *>
100 *> \param[in] TR
101 *> \verbatim
102 *> TR is REAL array, dimension (LDTR,2)
103 *> On entry, TR contains an N2 by N2 matrix.
104 *> \endverbatim
105 *>
106 *> \param[in] LDTR
107 *> \verbatim
108 *> LDTR is INTEGER
109 *> The leading dimension of the matrix TR. LDTR >= max(1,N2).
110 *> \endverbatim
111 *>
112 *> \param[in] B
113 *> \verbatim
114 *> B is REAL array, dimension (LDB,2)
115 *> On entry, the N1 by N2 matrix B contains the right-hand
116 *> side of the equation.
117 *> \endverbatim
118 *>
119 *> \param[in] LDB
120 *> \verbatim
121 *> LDB is INTEGER
122 *> The leading dimension of the matrix B. LDB >= max(1,N1).
123 *> \endverbatim
124 *>
125 *> \param[out] SCALE
126 *> \verbatim
127 *> SCALE is REAL
128 *> On exit, SCALE contains the scale factor. SCALE is chosen
129 *> less than or equal to 1 to prevent the solution overflowing.
130 *> \endverbatim
131 *>
132 *> \param[out] X
133 *> \verbatim
134 *> X is REAL array, dimension (LDX,2)
135 *> On exit, X contains the N1 by N2 solution.
136 *> \endverbatim
137 *>
138 *> \param[in] LDX
139 *> \verbatim
140 *> LDX is INTEGER
141 *> The leading dimension of the matrix X. LDX >= max(1,N1).
142 *> \endverbatim
143 *>
144 *> \param[out] XNORM
145 *> \verbatim
146 *> XNORM is REAL
147 *> On exit, XNORM is the infinity-norm of the solution.
148 *> \endverbatim
149 *>
150 *> \param[out] INFO
151 *> \verbatim
152 *> INFO is INTEGER
153 *> On exit, INFO is set to
154 *> 0: successful exit.
155 *> 1: TL and TR have too close eigenvalues, so TL or
156 *> TR is perturbed to get a nonsingular equation.
157 *> NOTE: In the interests of speed, this routine does not
158 *> check the inputs for errors.
159 *> \endverbatim
160 *
161 * Authors:
162 * ========
163 *
164 *> \author Univ. of Tennessee
165 *> \author Univ. of California Berkeley
166 *> \author Univ. of Colorado Denver
167 *> \author NAG Ltd.
168 *
169 *> \date June 2016
170 *
171 *> \ingroup realSYauxiliary
172 *
173 * =====================================================================
174  SUBROUTINE slasy2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR,
175  $ ldtr, b, ldb, scale, x, ldx, xnorm, info )
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  REAL SCALE, XNORM
186 * ..
187 * .. Array Arguments ..
188  REAL B( ldb, * ), TL( ldtl, * ), TR( ldtr, * ),
189  $ x( ldx, * )
190 * ..
191 *
192 * =====================================================================
193 *
194 * .. Parameters ..
195  REAL ZERO, ONE
196  parameter ( zero = 0.0e+0, one = 1.0e+0 )
197  REAL TWO, HALF, EIGHT
198  parameter ( two = 2.0e+0, half = 0.5e+0, eight = 8.0e+0 )
199 * ..
200 * .. Local Scalars ..
201  LOGICAL BSWAP, XSWAP
202  INTEGER I, IP, IPIV, IPSV, J, JP, JPSV, K
203  REAL 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  REAL BTMP( 4 ), T16( 4, 4 ), TMP( 4 ), X2( 2 )
211 * ..
212 * .. External Functions ..
213  INTEGER ISAMAX
214  REAL SLAMCH
215  EXTERNAL isamax, slamch
216 * ..
217 * .. External Subroutines ..
218  EXTERNAL scopy, sswap
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 = slamch( 'P' )
243  smlnum = slamch( '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 = isamax( 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 scopy( 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 sswap( 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 sswap( 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 SLASY2
481 *
482  END
slasy2
subroutine slasy2(LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR, LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO)
SLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2.
Definition: slasy2.f:176
sswap
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:53
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53