LAPACK  3.10.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 *> \ingroup realSYauxiliary
170 *
171 * =====================================================================
172  SUBROUTINE slasy2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR,
173  $ LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO )
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 *
479  END
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:174
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82