LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> Download SLASY2 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasy2.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasy2.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasy2.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE SLASY2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR,
20* LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO )
21*
22* .. Scalar Arguments ..
23* LOGICAL LTRANL, LTRANR
24* INTEGER INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2
25* REAL SCALE, XNORM
26* ..
27* .. Array Arguments ..
28* REAL B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ),
29* $ X( LDX, * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> SLASY2 solves for the N1 by N2 matrix X, 1 <= N1,N2 <= 2, in
39*>
40*> op(TL)*X + ISGN*X*op(TR) = SCALE*B,
41*>
42*> where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or
43*> -1. op(T) = T or T**T, where T**T denotes the transpose of T.
44*> \endverbatim
45*
46* Arguments:
47* ==========
48*
49*> \param[in] LTRANL
50*> \verbatim
51*> LTRANL is LOGICAL
52*> On entry, LTRANL specifies the op(TL):
53*> = .FALSE., op(TL) = TL,
54*> = .TRUE., op(TL) = TL**T.
55*> \endverbatim
56*>
57*> \param[in] LTRANR
58*> \verbatim
59*> LTRANR is LOGICAL
60*> On entry, LTRANR specifies the op(TR):
61*> = .FALSE., op(TR) = TR,
62*> = .TRUE., op(TR) = TR**T.
63*> \endverbatim
64*>
65*> \param[in] ISGN
66*> \verbatim
67*> ISGN is INTEGER
68*> On entry, ISGN specifies the sign of the equation
69*> as described before. ISGN may only be 1 or -1.
70*> \endverbatim
71*>
72*> \param[in] N1
73*> \verbatim
74*> N1 is INTEGER
75*> On entry, N1 specifies the order of matrix TL.
76*> N1 may only be 0, 1 or 2.
77*> \endverbatim
78*>
79*> \param[in] N2
80*> \verbatim
81*> N2 is INTEGER
82*> On entry, N2 specifies the order of matrix TR.
83*> N2 may only be 0, 1 or 2.
84*> \endverbatim
85*>
86*> \param[in] TL
87*> \verbatim
88*> TL is REAL array, dimension (LDTL,2)
89*> On entry, TL contains an N1 by N1 matrix.
90*> \endverbatim
91*>
92*> \param[in] LDTL
93*> \verbatim
94*> LDTL is INTEGER
95*> The leading dimension of the matrix TL. LDTL >= max(1,N1).
96*> \endverbatim
97*>
98*> \param[in] TR
99*> \verbatim
100*> TR is REAL array, dimension (LDTR,2)
101*> On entry, TR contains an N2 by N2 matrix.
102*> \endverbatim
103*>
104*> \param[in] LDTR
105*> \verbatim
106*> LDTR is INTEGER
107*> The leading dimension of the matrix TR. LDTR >= max(1,N2).
108*> \endverbatim
109*>
110*> \param[in] B
111*> \verbatim
112*> B is REAL array, dimension (LDB,2)
113*> On entry, the N1 by N2 matrix B contains the right-hand
114*> side of the equation.
115*> \endverbatim
116*>
117*> \param[in] LDB
118*> \verbatim
119*> LDB is INTEGER
120*> The leading dimension of the matrix B. LDB >= max(1,N1).
121*> \endverbatim
122*>
123*> \param[out] SCALE
124*> \verbatim
125*> SCALE is REAL
126*> On exit, SCALE contains the scale factor. SCALE is chosen
127*> less than or equal to 1 to prevent the solution overflowing.
128*> \endverbatim
129*>
130*> \param[out] X
131*> \verbatim
132*> X is REAL array, dimension (LDX,2)
133*> On exit, X contains the N1 by N2 solution.
134*> \endverbatim
135*>
136*> \param[in] LDX
137*> \verbatim
138*> LDX is INTEGER
139*> The leading dimension of the matrix X. LDX >= max(1,N1).
140*> \endverbatim
141*>
142*> \param[out] XNORM
143*> \verbatim
144*> XNORM is REAL
145*> On exit, XNORM is the infinity-norm of the solution.
146*> \endverbatim
147*>
148*> \param[out] INFO
149*> \verbatim
150*> INFO is INTEGER
151*> On exit, INFO is set to
152*> 0: successful exit.
153*> 1: TL and TR have too close eigenvalues, so TL or
154*> TR is perturbed to get a nonsingular equation.
155*> NOTE: In the interests of speed, this routine does not
156*> check the inputs for errors.
157*> \endverbatim
158*
159* Authors:
160* ========
161*
162*> \author Univ. of Tennessee
163*> \author Univ. of California Berkeley
164*> \author Univ. of Colorado Denver
165*> \author NAG Ltd.
166*
167*> \ingroup lasy2
168*
169* =====================================================================
170 SUBROUTINE slasy2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR,
171 $ LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO )
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*
477 END
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
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:172
sswap
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82