LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlasy2.f
Go to the documentation of this file.
1*> \brief \b DLASY2 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 DLASY2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasy2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasy2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasy2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DLASY2( 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* DOUBLE PRECISION SCALE, XNORM
28* ..
29* .. Array Arguments ..
30* DOUBLE PRECISION B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ),
31* $ X( LDX, * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> DLASY2 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 lasy2
170*
171* =====================================================================
172 SUBROUTINE dlasy2( 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 DOUBLE PRECISION SCALE, XNORM
183* ..
184* .. Array Arguments ..
185 DOUBLE PRECISION B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ),
186 $ x( ldx, * )
187* ..
188*
189* =====================================================================
190*
191* .. Parameters ..
192 DOUBLE PRECISION ZERO, ONE
193 parameter( zero = 0.0d+0, one = 1.0d+0 )
194 DOUBLE PRECISION TWO, HALF, EIGHT
195 parameter( two = 2.0d+0, half = 0.5d+0, eight = 8.0d+0 )
196* ..
197* .. Local Scalars ..
198 LOGICAL BSWAP, XSWAP
199 INTEGER I, IP, IPIV, IPSV, J, JP, JPSV, K
200 DOUBLE PRECISION 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 DOUBLE PRECISION BTMP( 4 ), T16( 4, 4 ), TMP( 4 ), X2( 2 )
208* ..
209* .. External Functions ..
210 INTEGER IDAMAX
211 DOUBLE PRECISION DLAMCH
212 EXTERNAL idamax, dlamch
213* ..
214* .. External Subroutines ..
215 EXTERNAL dcopy, dswap
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 = dlamch( 'P' )
240 smlnum = dlamch( '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 = idamax( 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 dcopy( 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 dswap( 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 dswap( 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 DLASY2
478*
479 END
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dlasy2
subroutine dlasy2(ltranl, ltranr, isgn, n1, n2, tl, ldtl, tr, ldtr, b, ldb, scale, x, ldx, xnorm, info)
DLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2.
Definition dlasy2.f:174
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82