LAPACK  3.10.1
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.
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 296 of file stgsyl.f.

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