LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 conjuguate 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 256 of file ztgsy2.f.

259 *
260 * -- LAPACK auxiliary routine --
261 * -- LAPACK is a software package provided by Univ. of Tennessee, --
262 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
263 *
264 * .. Scalar Arguments ..
265  CHARACTER TRANS
266  INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N
267  DOUBLE PRECISION RDSCAL, RDSUM, SCALE
268 * ..
269 * .. Array Arguments ..
270  COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ),
271  $ D( LDD, * ), E( LDE, * ), F( LDF, * )
272 * ..
273 *
274 * =====================================================================
275 *
276 * .. Parameters ..
277  DOUBLE PRECISION ZERO, ONE
278  INTEGER LDZ
279  parameter( zero = 0.0d+0, one = 1.0d+0, ldz = 2 )
280 * ..
281 * .. Local Scalars ..
282  LOGICAL NOTRAN
283  INTEGER I, IERR, J, K
284  DOUBLE PRECISION SCALOC
285  COMPLEX*16 ALPHA
286 * ..
287 * .. Local Arrays ..
288  INTEGER IPIV( LDZ ), JPIV( LDZ )
289  COMPLEX*16 RHS( LDZ ), Z( LDZ, LDZ )
290 * ..
291 * .. External Functions ..
292  LOGICAL LSAME
293  EXTERNAL lsame
294 * ..
295 * .. External Subroutines ..
296  EXTERNAL xerbla, zaxpy, zgesc2, zgetc2, zlatdf, 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 ), 1 )
394  CALL zaxpy( i-1, alpha, d( 1, i ), 1, f( 1, j ), 1 )
395  END IF
396  IF( j.LT.n ) THEN
397  CALL zaxpy( n-j, rhs( 2 ), b( j, j+1 ), ldb,
398  $ c( i, j+1 ), ldc )
399  CALL zaxpy( n-j, rhs( 2 ), e( j, j+1 ), lde,
400  $ f( i, j+1 ), ldf )
401  END IF
402 *
403  20 CONTINUE
404  30 CONTINUE
405  ELSE
406 *
407 * Solve transposed (I, J) - system:
408 * A(I, I)**H * R(I, J) + D(I, I)**H * L(J, J) = C(I, J)
409 * R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
410 * for I = 1, 2, ..., M, J = N, N - 1, ..., 1
411 *
412  scale = one
413  scaloc = one
414  DO 80 i = 1, m
415  DO 70 j = n, 1, -1
416 *
417 * Build 2 by 2 system Z**H
418 *
419  z( 1, 1 ) = dconjg( a( i, i ) )
420  z( 2, 1 ) = -dconjg( b( j, j ) )
421  z( 1, 2 ) = dconjg( d( i, i ) )
422  z( 2, 2 ) = -dconjg( e( j, j ) )
423 *
424 *
425 * Set up right hand side(s)
426 *
427  rhs( 1 ) = c( i, j )
428  rhs( 2 ) = f( i, j )
429 *
430 * Solve Z**H * x = RHS
431 *
432  CALL zgetc2( ldz, z, ldz, ipiv, jpiv, ierr )
433  IF( ierr.GT.0 )
434  $ info = ierr
435  CALL zgesc2( ldz, z, ldz, rhs, ipiv, jpiv, scaloc )
436  IF( scaloc.NE.one ) THEN
437  DO 40 k = 1, n
438  CALL zscal( m, dcmplx( scaloc, zero ), c( 1, k ),
439  $ 1 )
440  CALL zscal( m, dcmplx( scaloc, zero ), f( 1, k ),
441  $ 1 )
442  40 CONTINUE
443  scale = scale*scaloc
444  END IF
445 *
446 * Unpack solution vector(s)
447 *
448  c( i, j ) = rhs( 1 )
449  f( i, j ) = rhs( 2 )
450 *
451 * Substitute R(I, J) and L(I, J) into remaining equation.
452 *
453  DO 50 k = 1, j - 1
454  f( i, k ) = f( i, k ) + rhs( 1 )*dconjg( b( k, j ) ) +
455  $ rhs( 2 )*dconjg( e( k, j ) )
456  50 CONTINUE
457  DO 60 k = i + 1, m
458  c( k, j ) = c( k, j ) - dconjg( a( i, k ) )*rhs( 1 ) -
459  $ dconjg( d( i, k ) )*rhs( 2 )
460  60 CONTINUE
461 *
462  70 CONTINUE
463  80 CONTINUE
464  END IF
465  RETURN
466 *
467 * End of ZTGSY2
468 *
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2942
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:88
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78
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:115
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:111
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:169
Here is the call graph for this function:
Here is the caller graph for this function: