LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dtgsyl ( character  TRANS,
integer  IJOB,
integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( ldc, * )  C,
integer  LDC,
double precision, dimension( ldd, * )  D,
integer  LDD,
double precision, dimension( lde, * )  E,
integer  LDE,
double precision, dimension( ldf, * )  F,
integer  LDF,
double precision  SCALE,
double precision  DIF,
double precision, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

DTGSYL

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

Purpose:
 DTGSYL 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 real entries. (A, D) and (B, E) must be in
 generalized (real) Schur canonical form, i.e. A, B are upper quasi
 triangular and D, E are upper triangular.

 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**T, Im) ]         (2)
                [ kron(In, D)  -kron(E**T, Im) ].

 Here Ik is the identity matrix of size k and X**T is the transpose of
 X. kron(X, Y) is the Kronecker product between the matrices X and Y.

 If TRANS = 'T', DTGSYL solves the transposed system Z**T*y = scale*b,
 which is equivalent to solve for R and L in

             A**T * R + D**T * L = scale * C           (3)
             R * B**T + L * E**T = scale * -F

 This case (TRANS = 'T') 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 DLACON.

 If IJOB >= 1, DTGSYL 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. See [1-2] for more
 information.

 This is a level 3 BLAS algorithm.
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: 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 IJOB  = 1 is used).
           =4: Only an estimate of Dif[(A,D), (B,E)] is computed.
               ( DGECON on sub-systems is used ).
          Not referenced if TRANS = 'T'.
[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 DOUBLE PRECISION array, dimension (LDA, M)
          The upper quasi triangular matrix A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A. LDA >= max(1, M).
[in]B
          B is DOUBLE PRECISION array, dimension (LDB, N)
          The upper quasi triangular matrix B.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= max(1, N).
[in,out]C
          C is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
          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 = 'T', DIF is not touched.
[out]SCALE
          SCALE is DOUBLE PRECISION
          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, C and F hold the
          solutions R and L, respectively, to the homogeneous system
          with C = F = 0. Normally, SCALE = 1.
[out]WORK
          WORK is DOUBLE PRECISION 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+6)
[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 close eigenvalues.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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 301 of file dtgsyl.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: