LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages

◆ ztgsy2()

subroutine ztgsy2 ( character trans,
integer ijob,
integer m,
integer n,
complex*16, dimension( lda, * ) a,
integer lda,
complex*16, dimension( ldb, * ) b,
integer ldb,
complex*16, dimension( ldc, * ) c,
integer ldc,
complex*16, dimension( ldd, * ) d,
integer ldd,
complex*16, dimension( lde, * ) e,
integer lde,
complex*16, dimension( ldf, * ) f,
integer ldf,
double precision scale,
double precision rdsum,
double precision rdscal,
integer info )

ZTGSY2 solves the generalized Sylvester equation (unblocked algorithm).

Download ZTGSY2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZTGSY2 solves the generalized Sylvester equation
!>
!>             A * R - L * B = scale * C               (1)
!>             D * R - L * E = scale * F
!>
!> using Level 1 and 2 BLAS, where R and L are unknown M-by-N matrices,
!> (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
!> N-by-N and M-by-N, respectively. A, B, D and E are upper triangular
!> (i.e., (A,D) and (B,E) in generalized Schur form).
!>
!> The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output
!> scaling factor chosen to avoid overflow.
!>
!> In matrix notation solving equation (1) corresponds to solve
!> Zx = scale * b, where Z is defined as
!>
!>        Z = [ kron(In, A)  -kron(B**H, Im) ]             (2)
!>            [ kron(In, D)  -kron(E**H, Im) ],
!>
!> Ik is the identity matrix of size k and X**H is the conjugate transpose of X.
!> kron(X, Y) is the Kronecker product between the matrices X and Y.
!>
!> If TRANS = 'C', y in the conjugate transposed system Z**H*y = scale*b
!> is solved for, which is equivalent to solve for R and L in
!>
!>             A**H * R  + D**H * L   = scale * C           (3)
!>             R  * B**H + L  * E**H  = scale * -F
!>
!> This case is used to compute an estimate of Dif[(A, D), (B, E)] =
!> = sigma_min(Z) using reverse communication with ZLACON.
!>
!> ZTGSY2 also (IJOB >= 1) contributes to the computation in ZTGSYL
!> of an upper bound on the separation between to matrix pairs. Then
!> the input (A, D), (B, E) are sub-pencils of two matrix pairs in
!> ZTGSYL.
!>
Parameters
[in]TRANS
!>          TRANS is CHARACTER*1
!>          = 'N': solve the generalized Sylvester equation (1).
!>          = 'T': solve the 'transposed' system (3).
!>
[in]IJOB
!>          IJOB is INTEGER
!>          Specifies what kind of functionality to be performed.
!>          =0: solve (1) only.
!>          =1: A contribution from this subsystem to a Frobenius
!>              norm-based estimate of the separation between two matrix
!>              pairs is computed. (look ahead strategy is used).
!>          =2: A contribution from this subsystem to a Frobenius
!>              norm-based estimate of the separation between two matrix
!>              pairs is computed. (DGECON on sub-systems is used.)
!>          Not referenced if TRANS = 'T'.
!>
[in]M
!>          M is INTEGER
!>          On entry, M specifies the order of A and D, and the row
!>          dimension of C, F, R and L.
!>
[in]N
!>          N is INTEGER
!>          On entry, N specifies the order of B and E, and the column
!>          dimension of C, F, R and L.
!>
[in]A
!>          A is COMPLEX*16 array, dimension (LDA, M)
!>          On entry, A contains an upper triangular matrix.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the matrix A. LDA >= max(1, M).
!>
[in]B
!>          B is COMPLEX*16 array, dimension (LDB, N)
!>          On entry, B contains an upper triangular matrix.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the matrix B. LDB >= max(1, N).
!>
[in,out]C
!>          C is COMPLEX*16 array, dimension (LDC, N)
!>          On entry, C contains the right-hand-side of the first matrix
!>          equation in (1).
!>          On exit, if IJOB = 0, C has been overwritten by the solution
!>          R.
!>
[in]LDC
!>          LDC is INTEGER
!>          The leading dimension of the matrix C. LDC >= max(1, M).
!>
[in]D
!>          D is COMPLEX*16 array, dimension (LDD, M)
!>          On entry, D contains an upper triangular matrix.
!>
[in]LDD
!>          LDD is INTEGER
!>          The leading dimension of the matrix D. LDD >= max(1, M).
!>
[in]E
!>          E is COMPLEX*16 array, dimension (LDE, N)
!>          On entry, E contains an upper triangular matrix.
!>
[in]LDE
!>          LDE is INTEGER
!>          The leading dimension of the matrix E. LDE >= max(1, N).
!>
[in,out]F
!>          F is COMPLEX*16 array, dimension (LDF, N)
!>          On entry, F contains the right-hand-side of the second matrix
!>          equation in (1).
!>          On exit, if IJOB = 0, F has been overwritten by the solution
!>          L.
!>
[in]LDF
!>          LDF is INTEGER
!>          The leading dimension of the matrix F. LDF >= max(1, M).
!>
[out]SCALE
!>          SCALE is DOUBLE PRECISION
!>          On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions
!>          R and L (C and F on entry) will hold the solutions to a
!>          slightly perturbed system but the input matrices A, B, D and
!>          E have not been changed. If SCALE = 0, R and L will hold the
!>          solutions to the homogeneous system with C = F = 0.
!>          Normally, SCALE = 1.
!>
[in,out]RDSUM
!>          RDSUM is DOUBLE PRECISION
!>          On entry, the sum of squares of computed contributions to
!>          the Dif-estimate under computation by ZTGSYL, where the
!>          scaling factor RDSCAL (see below) has been factored out.
!>          On exit, the corresponding sum of squares updated with the
!>          contributions from the current sub-system.
!>          If TRANS = 'T' RDSUM is not touched.
!>          NOTE: RDSUM only makes sense when ZTGSY2 is called by
!>          ZTGSYL.
!>
[in,out]RDSCAL
!>          RDSCAL is DOUBLE PRECISION
!>          On entry, scaling factor used to prevent overflow in RDSUM.
!>          On exit, RDSCAL is updated w.r.t. the current contributions
!>          in RDSUM.
!>          If TRANS = 'T', RDSCAL is not touched.
!>          NOTE: RDSCAL only makes sense when ZTGSY2 is called by
!>          ZTGSYL.
!>
[out]INFO
!>          INFO is INTEGER
!>          On exit, if INFO is set to
!>            =0: Successful exit
!>            <0: If INFO = -i, input argument number i is illegal.
!>            >0: The matrix pairs (A, D) and (B, E) have common or very
!>                close eigenvalues.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.

Definition at line 254 of file ztgsy2.f.

258*
259* -- LAPACK auxiliary routine --
260* -- LAPACK is a software package provided by Univ. of Tennessee, --
261* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
262*
263* .. Scalar Arguments ..
264 CHARACTER TRANS
265 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N
266 DOUBLE PRECISION RDSCAL, RDSUM, SCALE
267* ..
268* .. Array Arguments ..
269 COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ),
270 $ D( LDD, * ), E( LDE, * ), F( LDF, * )
271* ..
272*
273* =====================================================================
274*
275* .. Parameters ..
276 DOUBLE PRECISION ZERO, ONE
277 INTEGER LDZ
278 parameter( zero = 0.0d+0, one = 1.0d+0, ldz = 2 )
279* ..
280* .. Local Scalars ..
281 LOGICAL NOTRAN
282 INTEGER I, IERR, J, K
283 DOUBLE PRECISION SCALOC
284 COMPLEX*16 ALPHA
285* ..
286* .. Local Arrays ..
287 INTEGER IPIV( LDZ ), JPIV( LDZ )
288 COMPLEX*16 RHS( LDZ ), Z( LDZ, LDZ )
289* ..
290* .. External Functions ..
291 LOGICAL LSAME
292 EXTERNAL lsame
293* ..
294* .. External Subroutines ..
295 EXTERNAL xerbla, zaxpy, zgesc2, zgetc2, zlatdf,
296 $ zscal
297* ..
298* .. Intrinsic Functions ..
299 INTRINSIC dcmplx, dconjg, max
300* ..
301* .. Executable Statements ..
302*
303* Decode and test input parameters
304*
305 info = 0
306 ierr = 0
307 notran = lsame( trans, 'N' )
308 IF( .NOT.notran .AND. .NOT.lsame( trans, 'C' ) ) THEN
309 info = -1
310 ELSE IF( notran ) THEN
311 IF( ( ijob.LT.0 ) .OR. ( ijob.GT.2 ) ) THEN
312 info = -2
313 END IF
314 END IF
315 IF( info.EQ.0 ) THEN
316 IF( m.LE.0 ) THEN
317 info = -3
318 ELSE IF( n.LE.0 ) THEN
319 info = -4
320 ELSE IF( lda.LT.max( 1, m ) ) THEN
321 info = -6
322 ELSE IF( ldb.LT.max( 1, n ) ) THEN
323 info = -8
324 ELSE IF( ldc.LT.max( 1, m ) ) THEN
325 info = -10
326 ELSE IF( ldd.LT.max( 1, m ) ) THEN
327 info = -12
328 ELSE IF( lde.LT.max( 1, n ) ) THEN
329 info = -14
330 ELSE IF( ldf.LT.max( 1, m ) ) THEN
331 info = -16
332 END IF
333 END IF
334 IF( info.NE.0 ) THEN
335 CALL xerbla( 'ZTGSY2', -info )
336 RETURN
337 END IF
338*
339 IF( notran ) THEN
340*
341* Solve (I, J) - system
342* A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
343* D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
344* for I = M, M - 1, ..., 1; J = 1, 2, ..., N
345*
346 scale = one
347 scaloc = one
348 DO 30 j = 1, n
349 DO 20 i = m, 1, -1
350*
351* Build 2 by 2 system
352*
353 z( 1, 1 ) = a( i, i )
354 z( 2, 1 ) = d( i, i )
355 z( 1, 2 ) = -b( j, j )
356 z( 2, 2 ) = -e( j, j )
357*
358* Set up right hand side(s)
359*
360 rhs( 1 ) = c( i, j )
361 rhs( 2 ) = f( i, j )
362*
363* Solve Z * x = RHS
364*
365 CALL zgetc2( ldz, z, ldz, ipiv, jpiv, ierr )
366 IF( ierr.GT.0 )
367 $ info = ierr
368 IF( ijob.EQ.0 ) THEN
369 CALL zgesc2( ldz, z, ldz, rhs, ipiv, jpiv, scaloc )
370 IF( scaloc.NE.one ) THEN
371 DO 10 k = 1, n
372 CALL zscal( m, dcmplx( scaloc, zero ),
373 $ c( 1, k ), 1 )
374 CALL zscal( m, dcmplx( scaloc, zero ),
375 $ f( 1, k ), 1 )
376 10 CONTINUE
377 scale = scale*scaloc
378 END IF
379 ELSE
380 CALL zlatdf( ijob, ldz, z, ldz, rhs, rdsum, rdscal,
381 $ ipiv, jpiv )
382 END IF
383*
384* Unpack solution vector(s)
385*
386 c( i, j ) = rhs( 1 )
387 f( i, j ) = rhs( 2 )
388*
389* Substitute R(I, J) and L(I, J) into remaining equation.
390*
391 IF( i.GT.1 ) THEN
392 alpha = -rhs( 1 )
393 CALL zaxpy( i-1, alpha, a( 1, i ), 1, c( 1, j ),
394 $ 1 )
395 CALL zaxpy( i-1, alpha, d( 1, i ), 1, f( 1, j ),
396 $ 1 )
397 END IF
398 IF( j.LT.n ) THEN
399 CALL zaxpy( n-j, rhs( 2 ), b( j, j+1 ), ldb,
400 $ c( i, j+1 ), ldc )
401 CALL zaxpy( n-j, rhs( 2 ), e( j, j+1 ), lde,
402 $ f( i, j+1 ), ldf )
403 END IF
404*
405 20 CONTINUE
406 30 CONTINUE
407 ELSE
408*
409* Solve transposed (I, J) - system:
410* A(I, I)**H * R(I, J) + D(I, I)**H * L(J, J) = C(I, J)
411* R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
412* for I = 1, 2, ..., M, J = N, N - 1, ..., 1
413*
414 scale = one
415 scaloc = one
416 DO 80 i = 1, m
417 DO 70 j = n, 1, -1
418*
419* Build 2 by 2 system Z**H
420*
421 z( 1, 1 ) = dconjg( a( i, i ) )
422 z( 2, 1 ) = -dconjg( b( j, j ) )
423 z( 1, 2 ) = dconjg( d( i, i ) )
424 z( 2, 2 ) = -dconjg( e( j, j ) )
425*
426*
427* Set up right hand side(s)
428*
429 rhs( 1 ) = c( i, j )
430 rhs( 2 ) = f( i, j )
431*
432* Solve Z**H * x = RHS
433*
434 CALL zgetc2( ldz, z, ldz, ipiv, jpiv, ierr )
435 IF( ierr.GT.0 )
436 $ info = ierr
437 CALL zgesc2( ldz, z, ldz, rhs, ipiv, jpiv, scaloc )
438 IF( scaloc.NE.one ) THEN
439 DO 40 k = 1, n
440 CALL zscal( m, dcmplx( scaloc, zero ), c( 1,
441 $ k ),
442 $ 1 )
443 CALL zscal( m, dcmplx( scaloc, zero ), f( 1,
444 $ k ),
445 $ 1 )
446 40 CONTINUE
447 scale = scale*scaloc
448 END IF
449*
450* Unpack solution vector(s)
451*
452 c( i, j ) = rhs( 1 )
453 f( i, j ) = rhs( 2 )
454*
455* Substitute R(I, J) and L(I, J) into remaining equation.
456*
457 DO 50 k = 1, j - 1
458 f( i, k ) = f( i, k ) + rhs( 1 )*dconjg( b( k, j ) ) +
459 $ rhs( 2 )*dconjg( e( k, j ) )
460 50 CONTINUE
461 DO 60 k = i + 1, m
462 c( k, j ) = c( k, j ) - dconjg( a( i, k ) )*rhs( 1 ) -
463 $ dconjg( d( i, k ) )*rhs( 2 )
464 60 CONTINUE
465*
466 70 CONTINUE
467 80 CONTINUE
468 END IF
469 RETURN
470*
471* End of ZTGSY2
472*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
lde
logical function lde(ri, rj, lr)
Definition dblat2.f:2970
zaxpy
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
Definition zaxpy.f:88
zgesc2
subroutine zgesc2(n, a, lda, rhs, ipiv, jpiv, scale)
ZGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed...
Definition zgesc2.f:113
zgetc2
subroutine zgetc2(n, a, lda, ipiv, jpiv, info)
ZGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix.
Definition zgetc2.f:109
zlatdf
subroutine zlatdf(ijob, n, z, ldz, rhs, rdsum, rdscal, ipiv, jpiv)
ZLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution ...
Definition zlatdf.f:167
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
zscal
subroutine zscal(n, za, zx, incx)
ZSCAL
Definition zscal.f:78
Here is the call graph for this function:
Here is the caller graph for this function: