LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ stgsyl()

subroutine stgsyl ( character  TRANS,
integer  IJOB,
integer  M,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( ldb, * )  B,
integer  LDB,
real, dimension( ldc, * )  C,
integer  LDC,
real, dimension( ldd, * )  D,
integer  LDD,
real, dimension( lde, * )  E,
integer  LDE,
real, dimension( ldf, * )  F,
integer  LDF,
real  SCALE,
real  DIF,
real, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

STGSYL

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

Purpose:
 STGSYL 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', STGSYL 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 SLACON.

 If IJOB >= 1, STGSYL 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.
               ( SGECON 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 = 'T', DIF is not touched.
[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, 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 REAL 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
December 2016
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 stgsyl.f.

301 *
302 * -- LAPACK computational routine (version 3.7.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 * December 2016
306 *
307 * .. Scalar Arguments ..
308  CHARACTER trans
309  INTEGER ijob, info, lda, ldb, ldc, ldd, lde, ldf,
310  $ lwork, m, n
311  REAL dif, scale
312 * ..
313 * .. Array Arguments ..
314  INTEGER iwork( * )
315  REAL a( lda, * ), b( ldb, * ), c( ldc, * ),
316  $ d( ldd, * ), e( lde, * ), f( ldf, * ),
317  $ work( * )
318 * ..
319 *
320 * =====================================================================
321 * Replaced various illegal calls to SCOPY by calls to SLASET.
322 * Sven Hammarling, 1/5/02.
323 *
324 * .. Parameters ..
325  REAL zero, one
326  parameter( zero = 0.0e+0, one = 1.0e+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  REAL dscale, dsum, scale2, scaloc
333 * ..
334 * .. External Functions ..
335  LOGICAL lsame
336  INTEGER ilaenv
337  EXTERNAL lsame, ilaenv
338 * ..
339 * .. External Subroutines ..
340  EXTERNAL sgemm, slacpy, slaset, sscal, stgsy2, xerbla
341 * ..
342 * .. Intrinsic Functions ..
343  INTRINSIC max, REAL, 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( 'STGSYL', -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, 'STGSYL', trans, m, n, -1, -1 )
419  nb = ilaenv( 5, 'STGSYL', 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 slaset( 'F', m, n, zero, zero, c, ldc )
427  CALL slaset( 'F', m, n, zero, zero, f, ldf )
428  ELSE IF( ijob.GE.1 .AND. notran ) 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 stgsy2( 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( REAL( 2*M*N ) ) / ( dscale*sqrt( dsum ) )
449  ELSE
450  dif = sqrt( REAL( 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 slacpy( 'F', m, n, c, ldc, work, m )
460  CALL slacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
461  CALL slaset( 'F', m, n, zero, zero, c, ldc )
462  CALL slaset( 'F', m, n, zero, zero, f, ldf )
463  ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
464  CALL slacpy( 'F', m, n, work, m, c, ldc )
465  CALL slacpy( '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 stgsy2( 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 sscal( m, scaloc, c( 1, k ), 1 )
549  CALL sscal( m, scaloc, f( 1, k ), 1 )
550  80 CONTINUE
551  DO 90 k = js, je
552  CALL sscal( is-1, scaloc, c( 1, k ), 1 )
553  CALL sscal( is-1, scaloc, f( 1, k ), 1 )
554  90 CONTINUE
555  DO 100 k = js, je
556  CALL sscal( m-ie, scaloc, c( ie+1, k ), 1 )
557  CALL sscal( m-ie, scaloc, f( ie+1, k ), 1 )
558  100 CONTINUE
559  DO 110 k = je + 1, n
560  CALL sscal( m, scaloc, c( 1, k ), 1 )
561  CALL sscal( 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 sgemm( 'N', 'N', is-1, nb, mb, -one,
571  $ a( 1, is ), lda, c( is, js ), ldc, one,
572  $ c( 1, js ), ldc )
573  CALL sgemm( '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 sgemm( '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 sgemm( '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( REAL( 2*M*N ) ) / ( dscale*sqrt( dsum ) )
590  ELSE
591  dif = sqrt( REAL( 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 slacpy( 'F', m, n, c, ldc, work, m )
600  CALL slacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
601  CALL slaset( 'F', m, n, zero, zero, c, ldc )
602  CALL slaset( 'F', m, n, zero, zero, f, ldf )
603  ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
604  CALL slacpy( 'F', m, n, work, m, c, ldc )
605  CALL slacpy( '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 stgsy2( 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 sscal( m, scaloc, c( 1, k ), 1 )
636  CALL sscal( m, scaloc, f( 1, k ), 1 )
637  160 CONTINUE
638  DO 170 k = js, je
639  CALL sscal( is-1, scaloc, c( 1, k ), 1 )
640  CALL sscal( is-1, scaloc, f( 1, k ), 1 )
641  170 CONTINUE
642  DO 180 k = js, je
643  CALL sscal( m-ie, scaloc, c( ie+1, k ), 1 )
644  CALL sscal( m-ie, scaloc, f( ie+1, k ), 1 )
645  180 CONTINUE
646  DO 190 k = je + 1, n
647  CALL sscal( m, scaloc, c( 1, k ), 1 )
648  CALL sscal( 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 sgemm( 'N', 'T', mb, js-1, nb, one, c( is, js ),
657  $ ldc, b( 1, js ), ldb, one, f( is, 1 ),
658  $ ldf )
659  CALL sgemm( '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 sgemm( '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 sgemm( '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 STGSYL
681 *
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2945
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:112
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine stgsy2(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, IWORK, PQ, INFO)
STGSY2 solves the generalized Sylvester equation (unblocked algorithm).
Definition: stgsy2.f:276
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:81
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
Here is the call graph for this function:
Here is the caller graph for this function: