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

◆ dtgsyl()

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.
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 dtgsyl.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 DOUBLE PRECISION DIF, SCALE
309* ..
310* .. Array Arguments ..
311 INTEGER IWORK( * )
312 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
313 $ D( LDD, * ), E( LDE, * ), F( LDF, * ),
314 $ WORK( * )
315* ..
316*
317* =====================================================================
318* Replaced various illegal calls to DCOPY by calls to DLASET.
319* Sven Hammarling, 1/5/02.
320*
321* .. Parameters ..
322 DOUBLE PRECISION ZERO, ONE
323 parameter( zero = 0.0d+0, one = 1.0d+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 DOUBLE PRECISION DSCALE, DSUM, SCALE2, SCALOC
330* ..
331* .. External Functions ..
332 LOGICAL LSAME
333 INTEGER ILAENV
334 EXTERNAL lsame, ilaenv
335* ..
336* .. External Subroutines ..
337 EXTERNAL dgemm, dlacpy, dlaset, dscal, dtgsy2, xerbla
338* ..
339* .. Intrinsic Functions ..
340 INTRINSIC dble, max, 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( 'DTGSYL', -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, 'DTGSYL', trans, m, n, -1, -1 )
416 nb = ilaenv( 5, 'DTGSYL', 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 dlaset( 'F', m, n, zero, zero, c, ldc )
424 CALL dlaset( 'F', m, n, zero, zero, f, ldf )
425 ELSE IF( ijob.GE.1 ) 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 dtgsy2( 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( dble( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
446 ELSE
447 dif = sqrt( dble( 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 dlacpy( 'F', m, n, c, ldc, work, m )
457 CALL dlacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
458 CALL dlaset( 'F', m, n, zero, zero, c, ldc )
459 CALL dlaset( 'F', m, n, zero, zero, f, ldf )
460 ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
461 CALL dlacpy( 'F', m, n, work, m, c, ldc )
462 CALL dlacpy( '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 dtgsy2( 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 dscal( m, scaloc, c( 1, k ), 1 )
546 CALL dscal( m, scaloc, f( 1, k ), 1 )
547 80 CONTINUE
548 DO 90 k = js, je
549 CALL dscal( is-1, scaloc, c( 1, k ), 1 )
550 CALL dscal( is-1, scaloc, f( 1, k ), 1 )
551 90 CONTINUE
552 DO 100 k = js, je
553 CALL dscal( m-ie, scaloc, c( ie+1, k ), 1 )
554 CALL dscal( m-ie, scaloc, f( ie+1, k ), 1 )
555 100 CONTINUE
556 DO 110 k = je + 1, n
557 CALL dscal( m, scaloc, c( 1, k ), 1 )
558 CALL dscal( 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 dgemm( 'N', 'N', is-1, nb, mb, -one,
568 $ a( 1, is ), lda, c( is, js ), ldc, one,
569 $ c( 1, js ), ldc )
570 CALL dgemm( '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 dgemm( '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 dgemm( '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( dble( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
587 ELSE
588 dif = sqrt( dble( 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 dlacpy( 'F', m, n, c, ldc, work, m )
597 CALL dlacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
598 CALL dlaset( 'F', m, n, zero, zero, c, ldc )
599 CALL dlaset( 'F', m, n, zero, zero, f, ldf )
600 ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
601 CALL dlacpy( 'F', m, n, work, m, c, ldc )
602 CALL dlacpy( '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 dtgsy2( 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 dscal( m, scaloc, c( 1, k ), 1 )
633 CALL dscal( m, scaloc, f( 1, k ), 1 )
634 160 CONTINUE
635 DO 170 k = js, je
636 CALL dscal( is-1, scaloc, c( 1, k ), 1 )
637 CALL dscal( is-1, scaloc, f( 1, k ), 1 )
638 170 CONTINUE
639 DO 180 k = js, je
640 CALL dscal( m-ie, scaloc, c( ie+1, k ), 1 )
641 CALL dscal( m-ie, scaloc, f( ie+1, k ), 1 )
642 180 CONTINUE
643 DO 190 k = je + 1, n
644 CALL dscal( m, scaloc, c( 1, k ), 1 )
645 CALL dscal( 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 dgemm( 'N', 'T', mb, js-1, nb, one, c( is, js ),
654 $ ldc, b( 1, js ), ldb, one, f( is, 1 ),
655 $ ldf )
656 CALL dgemm( '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 dgemm( '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 dgemm( '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 DTGSYL
678*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
logical function lde(ri, rj, lr)
Definition dblat2.f:2970
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
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:110
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
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:274
Here is the call graph for this function:
Here is the caller graph for this function: