LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dlasy2()

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.

Definition at line 172 of file dlasy2.f.

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*
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: