LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine ctgsy2 ( character  TRANS,
integer  IJOB,
integer  M,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( ldb, * )  B,
integer  LDB,
complex, dimension( ldc, * )  C,
integer  LDC,
complex, dimension( ldd, * )  D,
integer  LDD,
complex, dimension( lde, * )  E,
integer  LDE,
complex, dimension( ldf, * )  F,
integer  LDF,
real  SCALE,
real  RDSUM,
real  RDSCAL,
integer  INFO 
)

CTGSY2 solves the generalized Sylvester equation (unblocked algorithm).

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

Purpose:
 CTGSY2 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 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 communicaton with CLACON.

 CTGSY2 also (IJOB >= 1) contributes to the computation in CTGSYL
 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
 CTGSYL.
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. (SGECON 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 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 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 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 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 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 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 REAL
          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 REAL
          On entry, the sum of squares of computed contributions to
          the Dif-estimate under computation by CTGSYL, 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 CTGSY2 is called by
          CTGSYL.
[in,out]RDSCAL
          RDSCAL is REAL
          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 CTGSY2 is called by
          CTGSYL.
[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.
Date
November 2015
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.

Definition at line 261 of file ctgsy2.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: