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

◆ 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 REAL SROUNDUP_LWORK
335 EXTERNAL lsame, ilaenv, sroundup_lwork
336* ..
337* .. External Subroutines ..
338 EXTERNAL sgemm, slacpy, slaset, sscal, stgsy2, xerbla
339* ..
340* .. Intrinsic Functions ..
341 INTRINSIC 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, 'T' ) ) 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 ) = sroundup_lwork(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( 'STGSYL', -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, 'STGSYL', trans, m, n, -1, -1 )
417 nb = ilaenv( 5, 'STGSYL', 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 slaset( 'F', m, n, zero, zero, c, ldc )
425 CALL slaset( 'F', m, n, zero, zero, 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 DO 30 iround = 1, isolve
435*
436* Use unblocked Level 2 solver
437*
438 dscale = zero
439 dsum = one
440 pq = 0
441 CALL stgsy2( trans, ifunc, m, n, a, lda, b, ldb, c, ldc, d,
442 $ ldd, e, lde, f, ldf, scale, dsum, dscale,
443 $ iwork, pq, info )
444 IF( dscale.NE.zero ) THEN
445 IF( ijob.EQ.1 .OR. ijob.EQ.3 ) THEN
446 dif = sqrt( real( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
447 ELSE
448 dif = sqrt( real( pq ) ) / ( dscale*sqrt( dsum ) )
449 END IF
450 END IF
451*
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 slacpy( 'F', m, n, c, ldc, work, m )
458 CALL slacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
459 CALL slaset( 'F', m, n, zero, zero, c, ldc )
460 CALL slaset( 'F', m, n, zero, zero, f, ldf )
461 ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
462 CALL slacpy( 'F', m, n, work, m, c, ldc )
463 CALL slacpy( 'F', m, n, work( m*n+1 ), m, f, ldf )
464 scale = scale2
465 END IF
466 30 CONTINUE
467*
468 RETURN
469 END IF
470*
471* Determine block structure of A
472*
473 p = 0
474 i = 1
475 40 CONTINUE
476 IF( i.GT.m )
477 $ GO TO 50
478 p = p + 1
479 iwork( p ) = i
480 i = i + mb
481 IF( i.GE.m )
482 $ GO TO 50
483 IF( a( i, i-1 ).NE.zero )
484 $ i = i + 1
485 GO TO 40
486 50 CONTINUE
487*
488 iwork( p+1 ) = m + 1
489 IF( iwork( p ).EQ.iwork( p+1 ) )
490 $ p = p - 1
491*
492* Determine block structure of B
493*
494 q = p + 1
495 j = 1
496 60 CONTINUE
497 IF( j.GT.n )
498 $ GO TO 70
499 q = q + 1
500 iwork( q ) = j
501 j = j + nb
502 IF( j.GE.n )
503 $ GO TO 70
504 IF( b( j, j-1 ).NE.zero )
505 $ j = j + 1
506 GO TO 60
507 70 CONTINUE
508*
509 iwork( q+1 ) = n + 1
510 IF( iwork( q ).EQ.iwork( q+1 ) )
511 $ q = q - 1
512*
513 IF( notran ) THEN
514*
515 DO 150 iround = 1, isolve
516*
517* Solve (I, J)-subsystem
518* A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
519* D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
520* for I = P, P - 1,..., 1; J = 1, 2,..., Q
521*
522 dscale = zero
523 dsum = one
524 pq = 0
525 scale = one
526 DO 130 j = p + 2, q
527 js = iwork( j )
528 je = iwork( j+1 ) - 1
529 nb = je - js + 1
530 DO 120 i = p, 1, -1
531 is = iwork( i )
532 ie = iwork( i+1 ) - 1
533 mb = ie - is + 1
534 ppqq = 0
535 CALL stgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
536 $ b( js, js ), ldb, c( is, js ), ldc,
537 $ d( is, is ), ldd, e( js, js ), lde,
538 $ f( is, js ), ldf, scaloc, dsum, dscale,
539 $ iwork( q+2 ), ppqq, linfo )
540 IF( linfo.GT.0 )
541 $ info = linfo
542*
543 pq = pq + ppqq
544 IF( scaloc.NE.one ) THEN
545 DO 80 k = 1, js - 1
546 CALL sscal( m, scaloc, c( 1, k ), 1 )
547 CALL sscal( m, scaloc, f( 1, k ), 1 )
548 80 CONTINUE
549 DO 90 k = js, je
550 CALL sscal( is-1, scaloc, c( 1, k ), 1 )
551 CALL sscal( is-1, scaloc, f( 1, k ), 1 )
552 90 CONTINUE
553 DO 100 k = js, je
554 CALL sscal( m-ie, scaloc, c( ie+1, k ), 1 )
555 CALL sscal( m-ie, scaloc, f( ie+1, k ), 1 )
556 100 CONTINUE
557 DO 110 k = je + 1, n
558 CALL sscal( m, scaloc, c( 1, k ), 1 )
559 CALL sscal( m, scaloc, f( 1, k ), 1 )
560 110 CONTINUE
561 scale = scale*scaloc
562 END IF
563*
564* Substitute R(I, J) and L(I, J) into remaining
565* equation.
566*
567 IF( i.GT.1 ) THEN
568 CALL sgemm( 'N', 'N', is-1, nb, mb, -one,
569 $ a( 1, is ), lda, c( is, js ), ldc, one,
570 $ c( 1, js ), ldc )
571 CALL sgemm( 'N', 'N', is-1, nb, mb, -one,
572 $ d( 1, is ), ldd, c( is, js ), ldc, one,
573 $ f( 1, js ), ldf )
574 END IF
575 IF( j.LT.q ) THEN
576 CALL sgemm( 'N', 'N', mb, n-je, nb, one,
577 $ f( is, js ), ldf, b( js, je+1 ), ldb,
578 $ one, c( is, je+1 ), ldc )
579 CALL sgemm( 'N', 'N', mb, n-je, nb, one,
580 $ f( is, js ), ldf, e( js, je+1 ), lde,
581 $ one, f( is, je+1 ), ldf )
582 END IF
583 120 CONTINUE
584 130 CONTINUE
585 IF( dscale.NE.zero ) THEN
586 IF( ijob.EQ.1 .OR. ijob.EQ.3 ) THEN
587 dif = sqrt( real( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
588 ELSE
589 dif = sqrt( real( pq ) ) / ( dscale*sqrt( dsum ) )
590 END IF
591 END IF
592 IF( isolve.EQ.2 .AND. iround.EQ.1 ) THEN
593 IF( notran ) THEN
594 ifunc = ijob
595 END IF
596 scale2 = scale
597 CALL slacpy( 'F', m, n, c, ldc, work, m )
598 CALL slacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
599 CALL slaset( 'F', m, n, zero, zero, c, ldc )
600 CALL slaset( 'F', m, n, zero, zero, f, ldf )
601 ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
602 CALL slacpy( 'F', m, n, work, m, c, ldc )
603 CALL slacpy( 'F', m, n, work( m*n+1 ), m, f, ldf )
604 scale = scale2
605 END IF
606 150 CONTINUE
607*
608 ELSE
609*
610* Solve transposed (I, J)-subsystem
611* A(I, I)**T * R(I, J) + D(I, I)**T * L(I, J) = C(I, J)
612* R(I, J) * B(J, J)**T + L(I, J) * E(J, J)**T = -F(I, J)
613* for I = 1,2,..., P; J = Q, Q-1,..., 1
614*
615 scale = one
616 DO 210 i = 1, p
617 is = iwork( i )
618 ie = iwork( i+1 ) - 1
619 mb = ie - is + 1
620 DO 200 j = q, p + 2, -1
621 js = iwork( j )
622 je = iwork( j+1 ) - 1
623 nb = je - js + 1
624 CALL stgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
625 $ b( js, js ), ldb, c( is, js ), ldc,
626 $ d( is, is ), ldd, e( js, js ), lde,
627 $ f( is, js ), ldf, scaloc, dsum, dscale,
628 $ iwork( q+2 ), ppqq, linfo )
629 IF( linfo.GT.0 )
630 $ info = linfo
631 IF( scaloc.NE.one ) THEN
632 DO 160 k = 1, js - 1
633 CALL sscal( m, scaloc, c( 1, k ), 1 )
634 CALL sscal( m, scaloc, f( 1, k ), 1 )
635 160 CONTINUE
636 DO 170 k = js, je
637 CALL sscal( is-1, scaloc, c( 1, k ), 1 )
638 CALL sscal( is-1, scaloc, f( 1, k ), 1 )
639 170 CONTINUE
640 DO 180 k = js, je
641 CALL sscal( m-ie, scaloc, c( ie+1, k ), 1 )
642 CALL sscal( m-ie, scaloc, f( ie+1, k ), 1 )
643 180 CONTINUE
644 DO 190 k = je + 1, n
645 CALL sscal( m, scaloc, c( 1, k ), 1 )
646 CALL sscal( m, scaloc, f( 1, k ), 1 )
647 190 CONTINUE
648 scale = scale*scaloc
649 END IF
650*
651* Substitute R(I, J) and L(I, J) into remaining equation.
652*
653 IF( j.GT.p+2 ) THEN
654 CALL sgemm( 'N', 'T', mb, js-1, nb, one, c( is, js ),
655 $ ldc, b( 1, js ), ldb, one, f( is, 1 ),
656 $ ldf )
657 CALL sgemm( 'N', 'T', mb, js-1, nb, one, f( is, js ),
658 $ ldf, e( 1, js ), lde, one, f( is, 1 ),
659 $ ldf )
660 END IF
661 IF( i.LT.p ) THEN
662 CALL sgemm( 'T', 'N', m-ie, nb, mb, -one,
663 $ a( is, ie+1 ), lda, c( is, js ), ldc, one,
664 $ c( ie+1, js ), ldc )
665 CALL sgemm( 'T', 'N', m-ie, nb, mb, -one,
666 $ d( is, ie+1 ), ldd, f( is, js ), ldf, one,
667 $ c( ie+1, js ), ldc )
668 END IF
669 200 CONTINUE
670 210 CONTINUE
671*
672 END IF
673*
674 work( 1 ) = sroundup_lwork(lwmin)
675*
676 RETURN
677*
678* End of STGSYL
679*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
logical function lde(ri, rj, lr)
Definition dblat2.f:2970
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
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
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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
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
Here is the call graph for this function:
Here is the caller graph for this function: