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

◆ ctgsy2()

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 communication 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.
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.

Definition at line 254 of file ctgsy2.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 REAL RDSCAL, RDSUM, SCALE
267* ..
268* .. Array Arguments ..
269 COMPLEX A( LDA, * ), B( LDB, * ), C( LDC, * ),
270 $ D( LDD, * ), E( LDE, * ), F( LDF, * )
271* ..
272*
273* =====================================================================
274*
275* .. Parameters ..
276 REAL ZERO, ONE
277 INTEGER LDZ
278 parameter( zero = 0.0e+0, one = 1.0e+0, ldz = 2 )
279* ..
280* .. Local Scalars ..
281 LOGICAL NOTRAN
282 INTEGER I, IERR, J, K
283 REAL SCALOC
284 COMPLEX ALPHA
285* ..
286* .. Local Arrays ..
287 INTEGER IPIV( LDZ ), JPIV( LDZ )
288 COMPLEX RHS( LDZ ), Z( LDZ, LDZ )
289* ..
290* .. External Functions ..
291 LOGICAL LSAME
292 EXTERNAL lsame
293* ..
294* .. External Subroutines ..
295 EXTERNAL caxpy, cgesc2, cgetc2, cscal, clatdf,
296 $ xerbla
297* ..
298* .. Intrinsic Functions ..
299 INTRINSIC cmplx, conjg, 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( 'CTGSY2', -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 cgetc2( ldz, z, ldz, ipiv, jpiv, ierr )
366 IF( ierr.GT.0 )
367 $ info = ierr
368 IF( ijob.EQ.0 ) THEN
369 CALL cgesc2( ldz, z, ldz, rhs, ipiv, jpiv, scaloc )
370 IF( scaloc.NE.one ) THEN
371 DO 10 k = 1, n
372 CALL cscal( m, cmplx( scaloc, zero ), c( 1,
373 $ k ),
374 $ 1 )
375 CALL cscal( m, cmplx( scaloc, zero ), f( 1,
376 $ k ),
377 $ 1 )
378 10 CONTINUE
379 scale = scale*scaloc
380 END IF
381 ELSE
382 CALL clatdf( ijob, ldz, z, ldz, rhs, rdsum, rdscal,
383 $ ipiv, jpiv )
384 END IF
385*
386* Unpack solution vector(s)
387*
388 c( i, j ) = rhs( 1 )
389 f( i, j ) = rhs( 2 )
390*
391* Substitute R(I, J) and L(I, J) into remaining equation.
392*
393 IF( i.GT.1 ) THEN
394 alpha = -rhs( 1 )
395 CALL caxpy( i-1, alpha, a( 1, i ), 1, c( 1, j ),
396 $ 1 )
397 CALL caxpy( i-1, alpha, d( 1, i ), 1, f( 1, j ),
398 $ 1 )
399 END IF
400 IF( j.LT.n ) THEN
401 CALL caxpy( n-j, rhs( 2 ), b( j, j+1 ), ldb,
402 $ c( i, j+1 ), ldc )
403 CALL caxpy( n-j, rhs( 2 ), e( j, j+1 ), lde,
404 $ f( i, j+1 ), ldf )
405 END IF
406*
407 20 CONTINUE
408 30 CONTINUE
409 ELSE
410*
411* Solve transposed (I, J) - system:
412* A(I, I)**H * R(I, J) + D(I, I)**H * L(J, J) = C(I, J)
413* R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
414* for I = 1, 2, ..., M, J = N, N - 1, ..., 1
415*
416 scale = one
417 scaloc = one
418 DO 80 i = 1, m
419 DO 70 j = n, 1, -1
420*
421* Build 2 by 2 system Z**H
422*
423 z( 1, 1 ) = conjg( a( i, i ) )
424 z( 2, 1 ) = -conjg( b( j, j ) )
425 z( 1, 2 ) = conjg( d( i, i ) )
426 z( 2, 2 ) = -conjg( e( j, j ) )
427*
428*
429* Set up right hand side(s)
430*
431 rhs( 1 ) = c( i, j )
432 rhs( 2 ) = f( i, j )
433*
434* Solve Z**H * x = RHS
435*
436 CALL cgetc2( ldz, z, ldz, ipiv, jpiv, ierr )
437 IF( ierr.GT.0 )
438 $ info = ierr
439 CALL cgesc2( ldz, z, ldz, rhs, ipiv, jpiv, scaloc )
440 IF( scaloc.NE.one ) THEN
441 DO 40 k = 1, n
442 CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
443 $ 1 )
444 CALL cscal( m, cmplx( scaloc, zero ), f( 1, 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 )*conjg( b( k, j ) ) +
459 $ rhs( 2 )*conjg( e( k, j ) )
460 50 CONTINUE
461 DO 60 k = i + 1, m
462 c( k, j ) = c( k, j ) - conjg( a( i, k ) )*rhs( 1 ) -
463 $ conjg( d( i, k ) )*rhs( 2 )
464 60 CONTINUE
465*
466 70 CONTINUE
467 80 CONTINUE
468 END IF
469 RETURN
470*
471* End of CTGSY2
472*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
lde
logical function lde(ri, rj, lr)
Definition dblat2.f:2970
caxpy
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
cgesc2
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:113
cgetc2
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:109
clatdf
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:167
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
cscal
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
Here is the call graph for this function:
Here is the caller graph for this function: