LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
 subroutine ctgsyl ( 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 DIF, complex, dimension( * ) WORK, integer LWORK, integer, dimension( * ) IWORK, integer INFO )

CTGSYL

Purpose:
``` CTGSYL solves the generalized Sylvester equation:

A * R - L * B = scale * C            (1)
D * R - L * E = scale * F

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, with complex entries. 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 (1) is equivalent 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) ],

Here Ix is the identity matrix of size x 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 (TRANS = 'C') is used to compute an one-norm-based estimate
of Dif[(A,D), (B,E)], the separation between the matrix pairs (A,D)
and (B,E), using CLACON.

If IJOB >= 1, CTGSYL computes a Frobenius norm-based estimate of
Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the
reciprocal of the smallest singular value of Z.

This is a level-3 BLAS algorithm.```
Parameters
 [in] TRANS ``` TRANS is CHARACTER*1 = 'N': solve the generalized sylvester equation (1). = 'C': solve the "conjugate transposed" system (3).``` [in] IJOB ``` IJOB is INTEGER Specifies what kind of functionality to be performed. =0: solve (1) only. =1: The functionality of 0 and 3. =2: The functionality of 0 and 4. =3: Only an estimate of Dif[(A,D), (B,E)] is computed. (look ahead strategy is used). =4: Only an estimate of Dif[(A,D), (B,E)] is computed. (CGECON on sub-systems is used). Not referenced if TRANS = 'C'.``` [in] M ``` M is INTEGER The order of the matrices A and D, and the row dimension of the matrices C, F, R and L.``` [in] N ``` N is INTEGER The order of the matrices B and E, and the column dimension of the matrices C, F, R and L.``` [in] A ``` A is COMPLEX array, dimension (LDA, M) The upper triangular matrix A.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1, M).``` [in] B ``` B is COMPLEX array, dimension (LDB, N) The upper triangular matrix B.``` [in] LDB ``` LDB is INTEGER The leading dimension of the array 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) or (3). On exit, if IJOB = 0, 1 or 2, C has been overwritten by the solution R. If IJOB = 3 or 4 and TRANS = 'N', C holds R, the solution achieved during the computation of the Dif-estimate.``` [in] LDC ``` LDC is INTEGER The leading dimension of the array C. LDC >= max(1, M).``` [in] D ``` D is COMPLEX array, dimension (LDD, M) The upper triangular matrix D.``` [in] LDD ``` LDD is INTEGER The leading dimension of the array D. LDD >= max(1, M).``` [in] E ``` E is COMPLEX array, dimension (LDE, N) The upper triangular matrix E.``` [in] LDE ``` LDE is INTEGER The leading dimension of the array 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) or (3). On exit, if IJOB = 0, 1 or 2, F has been overwritten by the solution L. If IJOB = 3 or 4 and TRANS = 'N', F holds L, the solution achieved during the computation of the Dif-estimate.``` [in] LDF ``` LDF is INTEGER The leading dimension of the array F. LDF >= max(1, M).``` [out] DIF ``` DIF is REAL On exit DIF is the reciprocal of a lower bound of the reciprocal of the Dif-function, i.e. DIF is an upper bound of Dif[(A,D), (B,E)] = sigma-min(Z), where Z as in (2). IF IJOB = 0 or TRANS = 'C', DIF is not referenced.``` [out] SCALE ``` SCALE is REAL On exit SCALE is the scaling factor in (1) or (3). If 0 < SCALE < 1, C and F hold the solutions R and L, resp., 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 homogenious system with C = F = 0.``` [out] WORK ``` WORK is COMPLEX array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK.``` [in] LWORK ``` LWORK is INTEGER The dimension of the array WORK. LWORK > = 1. If IJOB = 1 or 2 and TRANS = 'N', LWORK >= max(1,2*M*N). If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA.``` [out] IWORK ` IWORK is INTEGER array, dimension (M+N+2)` [out] INFO ``` INFO is INTEGER =0: successful exit <0: If INFO = -i, the i-th argument had an illegal value. >0: (A, D) and (B, E) have common or very close eigenvalues.```
Date
November 2011
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.
References:
[1] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software for Solving the Generalized Sylvester Equation and Estimating the Separation between Regular Matrix Pairs, Report UMINF - 93.23, Department of Computing Science, Umea University, S-901 87 Umea, Sweden, December 1993, Revised April 1994, Also as LAPACK Working Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996.
[2] B. Kagstrom, A Perturbation Analysis of the Generalized Sylvester Equation (AR - LB, DR - LE ) = (C, F), SIAM J. Matrix Anal. Appl., 15(4):1045-1060, 1994.
[3] B. Kagstrom and L. Westin, Generalized Schur Methods with Condition Estimators for Solving the Generalized Sylvester Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751.

Definition at line 297 of file ctgsyl.f.

297 *
298 * -- LAPACK computational routine (version 3.4.0) --
299 * -- LAPACK is a software package provided by Univ. of Tennessee, --
300 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
301 * November 2011
302 *
303 * .. Scalar Arguments ..
304  CHARACTER trans
305  INTEGER ijob, info, lda, ldb, ldc, ldd, lde, ldf,
306  \$ lwork, m, n
307  REAL dif, scale
308 * ..
309 * .. Array Arguments ..
310  INTEGER iwork( * )
311  COMPLEX a( lda, * ), b( ldb, * ), c( ldc, * ),
312  \$ d( ldd, * ), e( lde, * ), f( ldf, * ),
313  \$ work( * )
314 * ..
315 *
316 * =====================================================================
317 * Replaced various illegal calls to CCOPY by calls to CLASET.
318 * Sven Hammarling, 1/5/02.
319 *
320 * .. Parameters ..
321  REAL zero, one
322  parameter ( zero = 0.0e+0, one = 1.0e+0 )
323  COMPLEX czero
324  parameter ( czero = (0.0e+0, 0.0e+0) )
325 * ..
326 * .. Local Scalars ..
327  LOGICAL lquery, notran
328  INTEGER i, ie, ifunc, iround, is, isolve, j, je, js, k,
329  \$ linfo, lwmin, mb, nb, p, pq, q
330  REAL dscale, dsum, scale2, scaloc
331 * ..
332 * .. External Functions ..
333  LOGICAL lsame
334  INTEGER ilaenv
335  EXTERNAL lsame, ilaenv
336 * ..
337 * .. External Subroutines ..
338  EXTERNAL cgemm, clacpy, claset, cscal, ctgsy2, xerbla
339 * ..
340 * .. Intrinsic Functions ..
341  INTRINSIC cmplx, max, REAL, sqrt
342 * ..
343 * .. Executable Statements ..
344 *
345 * Decode and test input parameters
346 *
347  info = 0
348  notran = lsame( trans, 'N' )
349  lquery = ( lwork.EQ.-1 )
350 *
351  IF( .NOT.notran .AND. .NOT.lsame( trans, 'C' ) ) THEN
352  info = -1
353  ELSE IF( notran ) THEN
354  IF( ( ijob.LT.0 ) .OR. ( ijob.GT.4 ) ) THEN
355  info = -2
356  END IF
357  END IF
358  IF( info.EQ.0 ) THEN
359  IF( m.LE.0 ) THEN
360  info = -3
361  ELSE IF( n.LE.0 ) THEN
362  info = -4
363  ELSE IF( lda.LT.max( 1, m ) ) THEN
364  info = -6
365  ELSE IF( ldb.LT.max( 1, n ) ) THEN
366  info = -8
367  ELSE IF( ldc.LT.max( 1, m ) ) THEN
368  info = -10
369  ELSE IF( ldd.LT.max( 1, m ) ) THEN
370  info = -12
371  ELSE IF( lde.LT.max( 1, n ) ) THEN
372  info = -14
373  ELSE IF( ldf.LT.max( 1, m ) ) THEN
374  info = -16
375  END IF
376  END IF
377 *
378  IF( info.EQ.0 ) THEN
379  IF( notran ) THEN
380  IF( ijob.EQ.1 .OR. ijob.EQ.2 ) THEN
381  lwmin = max( 1, 2*m*n )
382  ELSE
383  lwmin = 1
384  END IF
385  ELSE
386  lwmin = 1
387  END IF
388  work( 1 ) = lwmin
389 *
390  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
391  info = -20
392  END IF
393  END IF
394 *
395  IF( info.NE.0 ) THEN
396  CALL xerbla( 'CTGSYL', -info )
397  RETURN
398  ELSE IF( lquery ) THEN
399  RETURN
400  END IF
401 *
402 * Quick return if possible
403 *
404  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
405  scale = 1
406  IF( notran ) THEN
407  IF( ijob.NE.0 ) THEN
408  dif = 0
409  END IF
410  END IF
411  RETURN
412  END IF
413 *
414 * Determine optimal block sizes MB and NB
415 *
416  mb = ilaenv( 2, 'CTGSYL', trans, m, n, -1, -1 )
417  nb = ilaenv( 5, 'CTGSYL', trans, m, n, -1, -1 )
418 *
419  isolve = 1
420  ifunc = 0
421  IF( notran ) THEN
422  IF( ijob.GE.3 ) THEN
423  ifunc = ijob - 2
424  CALL claset( 'F', m, n, czero, czero, c, ldc )
425  CALL claset( 'F', m, n, czero, czero, f, ldf )
426  ELSE IF( ijob.GE.1 .AND. notran ) THEN
427  isolve = 2
428  END IF
429  END IF
430 *
431  IF( ( mb.LE.1 .AND. nb.LE.1 ) .OR. ( mb.GE.m .AND. nb.GE.n ) )
432  \$ THEN
433 *
434 * Use unblocked Level 2 solver
435 *
436  DO 30 iround = 1, isolve
437 *
438  scale = one
439  dscale = zero
440  dsum = one
441  pq = m*n
442  CALL ctgsy2( trans, ifunc, m, n, a, lda, b, ldb, c, ldc, d,
443  \$ ldd, e, lde, f, ldf, scale, dsum, dscale,
444  \$ info )
445  IF( dscale.NE.zero ) THEN
446  IF( ijob.EQ.1 .OR. ijob.EQ.3 ) THEN
447  dif = sqrt( REAL( 2*M*N ) ) / ( dscale*sqrt( dsum ) )
448  ELSE
449  dif = sqrt( REAL( PQ ) ) / ( dscale*sqrt( dsum ) )
450  END IF
451  END IF
452  IF( isolve.EQ.2 .AND. iround.EQ.1 ) THEN
453  IF( notran ) THEN
454  ifunc = ijob
455  END IF
456  scale2 = scale
457  CALL clacpy( 'F', m, n, c, ldc, work, m )
458  CALL clacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
459  CALL claset( 'F', m, n, czero, czero, c, ldc )
460  CALL claset( 'F', m, n, czero, czero, f, ldf )
461  ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
462  CALL clacpy( 'F', m, n, work, m, c, ldc )
463  CALL clacpy( 'F', m, n, work( m*n+1 ), m, f, ldf )
464  scale = scale2
465  END IF
466  30 CONTINUE
467 *
468  RETURN
469 *
470  END IF
471 *
472 * Determine block structure of A
473 *
474  p = 0
475  i = 1
476  40 CONTINUE
477  IF( i.GT.m )
478  \$ GO TO 50
479  p = p + 1
480  iwork( p ) = i
481  i = i + mb
482  IF( i.GE.m )
483  \$ GO TO 50
484  GO TO 40
485  50 CONTINUE
486  iwork( p+1 ) = m + 1
487  IF( iwork( p ).EQ.iwork( p+1 ) )
488  \$ p = p - 1
489 *
490 * Determine block structure of B
491 *
492  q = p + 1
493  j = 1
494  60 CONTINUE
495  IF( j.GT.n )
496  \$ GO TO 70
497 *
498  q = q + 1
499  iwork( q ) = j
500  j = j + nb
501  IF( j.GE.n )
502  \$ GO TO 70
503  GO TO 60
504 *
505  70 CONTINUE
506  iwork( q+1 ) = n + 1
507  IF( iwork( q ).EQ.iwork( q+1 ) )
508  \$ q = q - 1
509 *
510  IF( notran ) THEN
511  DO 150 iround = 1, isolve
512 *
513 * Solve (I, J) - subsystem
514 * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
515 * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
516 * for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
517 *
518  pq = 0
519  scale = one
520  dscale = zero
521  dsum = one
522  DO 130 j = p + 2, q
523  js = iwork( j )
524  je = iwork( j+1 ) - 1
525  nb = je - js + 1
526  DO 120 i = p, 1, -1
527  is = iwork( i )
528  ie = iwork( i+1 ) - 1
529  mb = ie - is + 1
530  CALL ctgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
531  \$ b( js, js ), ldb, c( is, js ), ldc,
532  \$ d( is, is ), ldd, e( js, js ), lde,
533  \$ f( is, js ), ldf, scaloc, dsum, dscale,
534  \$ linfo )
535  IF( linfo.GT.0 )
536  \$ info = linfo
537  pq = pq + mb*nb
538  IF( scaloc.NE.one ) THEN
539  DO 80 k = 1, js - 1
540  CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
541  \$ 1 )
542  CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
543  \$ 1 )
544  80 CONTINUE
545  DO 90 k = js, je
546  CALL cscal( is-1, cmplx( scaloc, zero ),
547  \$ c( 1, k ), 1 )
548  CALL cscal( is-1, cmplx( scaloc, zero ),
549  \$ f( 1, k ), 1 )
550  90 CONTINUE
551  DO 100 k = js, je
552  CALL cscal( m-ie, cmplx( scaloc, zero ),
553  \$ c( ie+1, k ), 1 )
554  CALL cscal( m-ie, cmplx( scaloc, zero ),
555  \$ f( ie+1, k ), 1 )
556  100 CONTINUE
557  DO 110 k = je + 1, n
558  CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
559  \$ 1 )
560  CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
561  \$ 1 )
562  110 CONTINUE
563  scale = scale*scaloc
564  END IF
565 *
566 * Substitute R(I,J) and L(I,J) into remaining equation.
567 *
568  IF( i.GT.1 ) THEN
569  CALL cgemm( 'N', 'N', is-1, nb, mb,
570  \$ cmplx( -one, zero ), a( 1, is ), lda,
571  \$ c( is, js ), ldc, cmplx( one, zero ),
572  \$ c( 1, js ), ldc )
573  CALL cgemm( 'N', 'N', is-1, nb, mb,
574  \$ cmplx( -one, zero ), d( 1, is ), ldd,
575  \$ c( is, js ), ldc, cmplx( one, zero ),
576  \$ f( 1, js ), ldf )
577  END IF
578  IF( j.LT.q ) THEN
579  CALL cgemm( 'N', 'N', mb, n-je, nb,
580  \$ cmplx( one, zero ), f( is, js ), ldf,
581  \$ b( js, je+1 ), ldb, cmplx( one, zero ),
582  \$ c( is, je+1 ), ldc )
583  CALL cgemm( 'N', 'N', mb, n-je, nb,
584  \$ cmplx( one, zero ), f( is, js ), ldf,
585  \$ e( js, je+1 ), lde, cmplx( one, zero ),
586  \$ f( is, je+1 ), ldf )
587  END IF
588  120 CONTINUE
589  130 CONTINUE
590  IF( dscale.NE.zero ) THEN
591  IF( ijob.EQ.1 .OR. ijob.EQ.3 ) THEN
592  dif = sqrt( REAL( 2*M*N ) ) / ( dscale*sqrt( dsum ) )
593  ELSE
594  dif = sqrt( REAL( PQ ) ) / ( dscale*sqrt( dsum ) )
595  END IF
596  END IF
597  IF( isolve.EQ.2 .AND. iround.EQ.1 ) THEN
598  IF( notran ) THEN
599  ifunc = ijob
600  END IF
601  scale2 = scale
602  CALL clacpy( 'F', m, n, c, ldc, work, m )
603  CALL clacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
604  CALL claset( 'F', m, n, czero, czero, c, ldc )
605  CALL claset( 'F', m, n, czero, czero, f, ldf )
606  ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
607  CALL clacpy( 'F', m, n, work, m, c, ldc )
608  CALL clacpy( 'F', m, n, work( m*n+1 ), m, f, ldf )
609  scale = scale2
610  END IF
611  150 CONTINUE
612  ELSE
613 *
614 * Solve transposed (I, J)-subsystem
615 * A(I, I)**H * R(I, J) + D(I, I)**H * L(I, J) = C(I, J)
616 * R(I, J) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
617 * for I = 1,2,..., P; J = Q, Q-1,..., 1
618 *
619  scale = one
620  DO 210 i = 1, p
621  is = iwork( i )
622  ie = iwork( i+1 ) - 1
623  mb = ie - is + 1
624  DO 200 j = q, p + 2, -1
625  js = iwork( j )
626  je = iwork( j+1 ) - 1
627  nb = je - js + 1
628  CALL ctgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
629  \$ b( js, js ), ldb, c( is, js ), ldc,
630  \$ d( is, is ), ldd, e( js, js ), lde,
631  \$ f( is, js ), ldf, scaloc, dsum, dscale,
632  \$ linfo )
633  IF( linfo.GT.0 )
634  \$ info = linfo
635  IF( scaloc.NE.one ) THEN
636  DO 160 k = 1, js - 1
637  CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
638  \$ 1 )
639  CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
640  \$ 1 )
641  160 CONTINUE
642  DO 170 k = js, je
643  CALL cscal( is-1, cmplx( scaloc, zero ), c( 1, k ),
644  \$ 1 )
645  CALL cscal( is-1, cmplx( scaloc, zero ), f( 1, k ),
646  \$ 1 )
647  170 CONTINUE
648  DO 180 k = js, je
649  CALL cscal( m-ie, cmplx( scaloc, zero ),
650  \$ c( ie+1, k ), 1 )
651  CALL cscal( m-ie, cmplx( scaloc, zero ),
652  \$ f( ie+1, k ), 1 )
653  180 CONTINUE
654  DO 190 k = je + 1, n
655  CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
656  \$ 1 )
657  CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
658  \$ 1 )
659  190 CONTINUE
660  scale = scale*scaloc
661  END IF
662 *
663 * Substitute R(I,J) and L(I,J) into remaining equation.
664 *
665  IF( j.GT.p+2 ) THEN
666  CALL cgemm( 'N', 'C', mb, js-1, nb,
667  \$ cmplx( one, zero ), c( is, js ), ldc,
668  \$ b( 1, js ), ldb, cmplx( one, zero ),
669  \$ f( is, 1 ), ldf )
670  CALL cgemm( 'N', 'C', mb, js-1, nb,
671  \$ cmplx( one, zero ), f( is, js ), ldf,
672  \$ e( 1, js ), lde, cmplx( one, zero ),
673  \$ f( is, 1 ), ldf )
674  END IF
675  IF( i.LT.p ) THEN
676  CALL cgemm( 'C', 'N', m-ie, nb, mb,
677  \$ cmplx( -one, zero ), a( is, ie+1 ), lda,
678  \$ c( is, js ), ldc, cmplx( one, zero ),
679  \$ c( ie+1, js ), ldc )
680  CALL cgemm( 'C', 'N', m-ie, nb, mb,
681  \$ cmplx( -one, zero ), d( is, ie+1 ), ldd,
682  \$ f( is, js ), ldf, cmplx( one, zero ),
683  \$ c( ie+1, js ), ldc )
684  END IF
685  200 CONTINUE
686  210 CONTINUE
687  END IF
688 *
689  work( 1 ) = lwmin
690 *
691  RETURN
692 *
693 * End of CTGSYL
694 *
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 claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: claset.f:108
subroutine ctgsy2(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, INFO)
CTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
Definition: ctgsy2.f:261
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
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: