complex16

Functions
subroutine	zlaesy (a, b, c, rt1, rt2, evscal, cs1, sn1)
	ZLAESY computes the eigenvalues and eigenvectors of a 2-by-2 complex symmetric matrix.
double precision function	zlansy (norm, uplo, n, a, lda, work)
	ZLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex symmetric matrix.
subroutine	zlaqsy (uplo, n, a, lda, s, scond, amax, equed)
	ZLAQSY scales a symmetric/Hermitian matrix, using scaling factors computed by spoequ.
subroutine	zsymv (uplo, n, alpha, a, lda, x, incx, beta, y, incy)
	ZSYMV computes a matrix-vector product for a complex symmetric matrix.
subroutine	zsyr (uplo, n, alpha, x, incx, a, lda)
	ZSYR performs the symmetric rank-1 update of a complex symmetric matrix.
subroutine	zsyswapr (uplo, n, a, lda, i1, i2)
	ZSYSWAPR
subroutine	ztgsy2 (trans, ijob, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, rdsum, rdscal, info)
	ZTGSY2 solves the generalized Sylvester equation (unblocked algorithm).

Detailed Description

This is the group of complex16 auxiliary functions for SY matrices

Function Documentation

zlaesy()

subroutine zlaesy	(	complex*16	a,
		complex*16	b,
		complex*16	c,
		complex*16	rt1,
		complex*16	rt2,
		complex*16	evscal,
		complex*16	cs1,
		complex*16	sn1 )

ZLAESY computes the eigenvalues and eigenvectors of a 2-by-2 complex symmetric matrix.

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

Purpose:

!>
!> ZLAESY computes the eigendecomposition of a 2-by-2 symmetric matrix
!>    ( ( A, B );( B, C ) )
!> provided the norm of the matrix of eigenvectors is larger than
!> some threshold value.
!>
!> RT1 is the eigenvalue of larger absolute value, and RT2 of
!> smaller absolute value.  If the eigenvectors are computed, then
!> on return ( CS1, SN1 ) is the unit eigenvector for RT1, hence
!>
!> [  CS1     SN1   ] . [ A  B ] . [ CS1    -SN1   ] = [ RT1  0  ]
!> [ -SN1     CS1   ]   [ B  C ]   [ SN1     CS1   ]   [  0  RT2 ]
!> 

Parameters

[in]	A	!> A is COMPLEX*16 !> The ( 1, 1 ) element of input matrix. !>
[in]	B	!> B is COMPLEX*16 !> The ( 1, 2 ) element of input matrix. The ( 2, 1 ) element !> is also given by B, since the 2-by-2 matrix is symmetric. !>
[in]	C	!> C is COMPLEX*16 !> The ( 2, 2 ) element of input matrix. !>
[out]	RT1	!> RT1 is COMPLEX*16 !> The eigenvalue of larger modulus. !>
[out]	RT2	!> RT2 is COMPLEX*16 !> The eigenvalue of smaller modulus. !>
[out]	EVSCAL	!> EVSCAL is COMPLEX*16 !> The complex value by which the eigenvector matrix was scaled !> to make it orthonormal. If EVSCAL is zero, the eigenvectors !> were not computed. This means one of two things: the 2-by-2 !> matrix could not be diagonalized, or the norm of the matrix !> of eigenvectors before scaling was larger than the threshold !> value THRESH (set below). !>
[out]	CS1	!> CS1 is COMPLEX*16 !>
[out]	SN1	!> SN1 is COMPLEX*16 !> If EVSCAL .NE. 0, ( CS1, SN1 ) is the unit right eigenvector !> for RT1. !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 114 of file zlaesy.f.

*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      COMPLEX*16         A, B, C, CS1, EVSCAL, RT1, RT2, SN1
*     ..
*
* =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO
      parameter( zero = 0.0d0 )
      DOUBLE PRECISION   ONE
      parameter( one = 1.0d0 )
      COMPLEX*16         CONE
      parameter( cone = ( 1.0d0, 0.0d0 ) )
      DOUBLE PRECISION   HALF
      parameter( half = 0.5d0 )
      DOUBLE PRECISION   THRESH
      parameter( thresh = 0.1d0 )
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION   BABS, EVNORM, TABS, Z
      COMPLEX*16         S, T, TMP
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, sqrt
*     ..
*     .. Executable Statements ..
*
*
*     Special case:  The matrix is actually diagonal.
*     To avoid divide by zero later, we treat this case separately.
*
      IF( abs( b ).EQ.zero ) THEN
         rt1 = a
         rt2 = c
         IF( abs( rt1 ).LT.abs( rt2 ) ) THEN
            tmp = rt1
            rt1 = rt2
            rt2 = tmp
            cs1 = zero
            sn1 = one
         ELSE
            cs1 = one
            sn1 = zero
         END IF
      ELSE
*
*        Compute the eigenvalues and eigenvectors.
*        The characteristic equation is
*           lambda **2 - (A+C) lambda + (A*C - B*B)
*        and we solve it using the quadratic formula.
*
         s = ( a+c )*half
         t = ( a-c )*half
*
*        Take the square root carefully to avoid over/under flow.
*
         babs = abs( b )
         tabs = abs( t )
         z = max( babs, tabs )
         IF( z.GT.zero )
     $      t = z*sqrt( ( t / z )**2+( b / z )**2 )
*
*        Compute the two eigenvalues.  RT1 and RT2 are exchanged
*        if necessary so that RT1 will have the greater magnitude.
*
         rt1 = s + t
         rt2 = s - t
         IF( abs( rt1 ).LT.abs( rt2 ) ) THEN
            tmp = rt1
            rt1 = rt2
            rt2 = tmp
         END IF
*
*        Choose CS1 = 1 and SN1 to satisfy the first equation, then
*        scale the components of this eigenvector so that the matrix
*        of eigenvectors X satisfies  X * X**T = I .  (No scaling is
*        done if the norm of the eigenvalue matrix is less than THRESH.)
*
         sn1 = ( rt1-a ) / b
         tabs = abs( sn1 )
         IF( tabs.GT.one ) THEN
            t = tabs*sqrt( ( one / tabs )**2+( sn1 / tabs )**2 )
         ELSE
            t = sqrt( cone+sn1*sn1 )
         END IF
         evnorm = abs( t )
         IF( evnorm.GE.thresh ) THEN
            evscal = cone / t
            cs1 = evscal
            sn1 = sn1*evscal
         ELSE
            evscal = zero
         END IF
      END IF
      RETURN
*
*     End of ZLAESY
*

zlansy()

double precision function zlansy	(	character	norm,
		character	uplo,
		integer	n,
		complex*16, dimension( lda, * )	a,
		integer	lda,
		double precision, dimension( * )	work )

ZLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex symmetric matrix.

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

Purpose:

!>
!> ZLANSY  returns the value of the one norm,  or the Frobenius norm, or
!> the  infinity norm,  or the  element of  largest absolute value  of a
!> complex symmetric matrix A.
!> 

Returns

ZLANSY

!>
!>    ZLANSY = ( max(abs(A(i,j))), NORM = 'M' or 'm'
!>             (
!>             ( norm1(A),         NORM = '1', 'O' or 'o'
!>             (
!>             ( normI(A),         NORM = 'I' or 'i'
!>             (
!>             ( normF(A),         NORM = 'F', 'f', 'E' or 'e'
!>
!> where  norm1  denotes the  one norm of a matrix (maximum column sum),
!> normI  denotes the  infinity norm  of a matrix  (maximum row sum) and
!> normF  denotes the  Frobenius norm of a matrix (square root of sum of
!> squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.
!> 

Parameters

[in]	NORM	!> NORM is CHARACTER*1 !> Specifies the value to be returned in ZLANSY as described !> above. !>
[in]	UPLO	!> UPLO is CHARACTER*1 !> Specifies whether the upper or lower triangular part of the !> symmetric matrix A is to be referenced. !> = 'U': Upper triangular part of A is referenced !> = 'L': Lower triangular part of A is referenced !>
[in]	N	!> N is INTEGER !> The order of the matrix A. N >= 0. When N = 0, ZLANSY is !> set to zero. !>
[in]	A	!> A is COMPLEX*16 array, dimension (LDA,N) !> The symmetric matrix A. If UPLO = 'U', the leading n by n !> upper triangular part of A contains the upper triangular part !> of the matrix A, and the strictly lower triangular part of A !> is not referenced. If UPLO = 'L', the leading n by n lower !> triangular part of A contains the lower triangular part of !> the matrix A, and the strictly upper triangular part of A is !> not referenced. !>
[in]	LDA	!> LDA is INTEGER !> The leading dimension of the array A. LDA >= max(N,1). !>
[out]	WORK	!> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)), !> where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, !> WORK is not referenced. !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 122 of file zlansy.f.

*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          NORM, UPLO
      INTEGER            LDA, N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   WORK( * )
      COMPLEX*16         A( LDA, * )
*     ..
*
* =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ONE, ZERO
      parameter( one = 1.0d+0, zero = 0.0d+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, J
      DOUBLE PRECISION   ABSA, SCALE, SUM, VALUE
*     ..
*     .. External Functions ..
      LOGICAL            LSAME, DISNAN
      EXTERNAL           lsame, disnan
*     ..
*     .. External Subroutines ..
      EXTERNAL           zlassq
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, sqrt
*     ..
*     .. Executable Statements ..
*
      IF( n.EQ.0 ) THEN
         VALUE = zero
      ELSE IF( lsame( norm, 'M' ) ) THEN
*
*        Find max(abs(A(i,j))).
*
         VALUE = zero
         IF( lsame( uplo, 'U' ) ) THEN
            DO 20 j = 1, n
               DO 10 i = 1, j
                  sum = abs( a( i, j ) )
                  IF( VALUE .LT. sum .OR. disnan( sum ) ) VALUE = sum
   10          CONTINUE
   20       CONTINUE
         ELSE
            DO 40 j = 1, n
               DO 30 i = j, n
                  sum = abs( a( i, j ) )
                  IF( VALUE .LT. sum .OR. disnan( sum ) ) VALUE = sum
   30          CONTINUE
   40       CONTINUE
         END IF
      ELSE IF( ( lsame( norm, 'I' ) ) .OR. ( lsame( norm, 'O' ) ) .OR.
     $         ( norm.EQ.'1' ) ) THEN
*
*        Find normI(A) ( = norm1(A), since A is symmetric).
*
         VALUE = zero
         IF( lsame( uplo, 'U' ) ) THEN
            DO 60 j = 1, n
               sum = zero
               DO 50 i = 1, j - 1
                  absa = abs( a( i, j ) )
                  sum = sum + absa
                  work( i ) = work( i ) + absa
   50          CONTINUE
               work( j ) = sum + abs( a( j, j ) )
   60       CONTINUE
            DO 70 i = 1, n
               sum = work( i )
               IF( VALUE .LT. sum .OR. disnan( sum ) ) VALUE = sum
   70       CONTINUE
         ELSE
            DO 80 i = 1, n
               work( i ) = zero
   80       CONTINUE
            DO 100 j = 1, n
               sum = work( j ) + abs( a( j, j ) )
               DO 90 i = j + 1, n
                  absa = abs( a( i, j ) )
                  sum = sum + absa
                  work( i ) = work( i ) + absa
   90          CONTINUE
               IF( VALUE .LT. sum .OR. disnan( sum ) ) VALUE = sum
  100       CONTINUE
         END IF
      ELSE IF( ( lsame( norm, 'F' ) ) .OR. ( lsame( norm, 'E' ) ) ) THEN
*
*        Find normF(A).
*
         scale = zero
         sum = one
         IF( lsame( uplo, 'U' ) ) THEN
            DO 110 j = 2, n
               CALL zlassq( j-1, a( 1, j ), 1, scale, sum )
  110       CONTINUE
         ELSE
            DO 120 j = 1, n - 1
               CALL zlassq( n-j, a( j+1, j ), 1, scale, sum )
  120       CONTINUE
         END IF
         sum = 2*sum
         CALL zlassq( n, a, lda+1, scale, sum )
         VALUE = scale*sqrt( sum )
      END IF
*
      zlansy = VALUE
      RETURN
*
*     End of ZLANSY
*

zlaqsy()

subroutine zlaqsy	(	character	uplo,
		integer	n,
		complex*16, dimension( lda, * )	a,
		integer	lda,
		double precision, dimension( * )	s,
		double precision	scond,
		double precision	amax,
		character	equed )

ZLAQSY scales a symmetric/Hermitian matrix, using scaling factors computed by spoequ.

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

Purpose:

!>
!> ZLAQSY equilibrates a symmetric matrix A using the scaling factors
!> in the vector S.
!> 

Parameters

[in]	UPLO	!> UPLO is CHARACTER*1 !> Specifies whether the upper or lower triangular part of the !> symmetric matrix A is stored. !> = 'U': Upper triangular !> = 'L': Lower triangular !>
[in]	N	!> N is INTEGER !> The order of the matrix A. N >= 0. !>
[in,out]	A	!> A is COMPLEX*16 array, dimension (LDA,N) !> On entry, the symmetric matrix A. If UPLO = 'U', the leading !> n by n upper triangular part of A contains the upper !> triangular part of the matrix A, and the strictly lower !> triangular part of A is not referenced. If UPLO = 'L', the !> leading n by n lower triangular part of A contains the lower !> triangular part of the matrix A, and the strictly upper !> triangular part of A is not referenced. !> !> On exit, if EQUED = 'Y', the equilibrated matrix: !> diag(S) * A * diag(S). !>
[in]	LDA	!> LDA is INTEGER !> The leading dimension of the array A. LDA >= max(N,1). !>
[in]	S	!> S is DOUBLE PRECISION array, dimension (N) !> The scale factors for A. !>
[in]	SCOND	!> SCOND is DOUBLE PRECISION !> Ratio of the smallest S(i) to the largest S(i). !>
[in]	AMAX	!> AMAX is DOUBLE PRECISION !> Absolute value of largest matrix entry. !>
[out]	EQUED	!> EQUED is CHARACTER*1 !> Specifies whether or not equilibration was done. !> = 'N': No equilibration. !> = 'Y': Equilibration was done, i.e., A has been replaced by !> diag(S) * A * diag(S). !>

Internal Parameters:

!>  THRESH is a threshold value used to decide if scaling should be done
!>  based on the ratio of the scaling factors.  If SCOND < THRESH,
!>  scaling is done.
!>
!>  LARGE and SMALL are threshold values used to decide if scaling should
!>  be done based on the absolute size of the largest matrix element.
!>  If AMAX > LARGE or AMAX < SMALL, scaling is done.
!> 

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 133 of file zlaqsy.f.

*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          EQUED, UPLO
      INTEGER            LDA, N
      DOUBLE PRECISION   AMAX, SCOND
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   S( * )
      COMPLEX*16         A( LDA, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ONE, THRESH
      parameter( one = 1.0d+0, thresh = 0.1d+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, J
      DOUBLE PRECISION   CJ, LARGE, SMALL
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           lsame, dlamch
*     ..
*     .. Executable Statements ..
*
*     Quick return if possible
*
      IF( n.LE.0 ) THEN
         equed = 'N'
         RETURN
      END IF
*
*     Initialize LARGE and SMALL.
*
      small = dlamch( 'Safe minimum' ) / dlamch( 'Precision' )
      large = one / small
*
      IF( scond.GE.thresh .AND. amax.GE.small .AND. amax.LE.large ) THEN
*
*        No equilibration
*
         equed = 'N'
      ELSE
*
*        Replace A by diag(S) * A * diag(S).
*
         IF( lsame( uplo, 'U' ) ) THEN
*
*           Upper triangle of A is stored.
*
            DO 20 j = 1, n
               cj = s( j )
               DO 10 i = 1, j
                  a( i, j ) = cj*s( i )*a( i, j )
   10          CONTINUE
   20       CONTINUE
         ELSE
*
*           Lower triangle of A is stored.
*
            DO 40 j = 1, n
               cj = s( j )
               DO 30 i = j, n
                  a( i, j ) = cj*s( i )*a( i, j )
   30          CONTINUE
   40       CONTINUE
         END IF
         equed = 'Y'
      END IF
*
      RETURN
*
*     End of ZLAQSY
*

zsymv()

subroutine zsymv	(	character	uplo,
		integer	n,
		complex*16	alpha,
		complex*16, dimension( lda, * )	a,
		integer	lda,
		complex*16, dimension( * )	x,
		integer	incx,
		complex*16	beta,
		complex*16, dimension( * )	y,
		integer	incy )

ZSYMV computes a matrix-vector product for a complex symmetric matrix.

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

Purpose:

!>
!> ZSYMV  performs the matrix-vector  operation
!>
!>    y := alpha*A*x + beta*y,
!>
!> where alpha and beta are scalars, x and y are n element vectors and
!> A is an n by n symmetric matrix.
!> 

Parameters

[in]	UPLO	!> UPLO is CHARACTER*1 !> On entry, UPLO specifies whether the upper or lower !> triangular part of the array A is to be referenced as !> follows: !> !> UPLO = 'U' or 'u' Only the upper triangular part of A !> is to be referenced. !> !> UPLO = 'L' or 'l' Only the lower triangular part of A !> is to be referenced. !> !> Unchanged on exit. !>
[in]	N	!> N is INTEGER !> On entry, N specifies the order of the matrix A. !> N must be at least zero. !> Unchanged on exit. !>
[in]	ALPHA	!> ALPHA is COMPLEX*16 !> On entry, ALPHA specifies the scalar alpha. !> Unchanged on exit. !>
[in]	A	!> A is COMPLEX*16 array, dimension ( LDA, N ) !> Before entry, with UPLO = 'U' or 'u', the leading n by n !> upper triangular part of the array A must contain the upper !> triangular part of the symmetric matrix and the strictly !> lower triangular part of A is not referenced. !> Before entry, with UPLO = 'L' or 'l', the leading n by n !> lower triangular part of the array A must contain the lower !> triangular part of the symmetric matrix and the strictly !> upper triangular part of A is not referenced. !> Unchanged on exit. !>
[in]	LDA	!> LDA is INTEGER !> On entry, LDA specifies the first dimension of A as declared !> in the calling (sub) program. LDA must be at least !> max( 1, N ). !> Unchanged on exit. !>
[in]	X	!> X is COMPLEX*16 array, dimension at least !> ( 1 + ( N - 1 )*abs( INCX ) ). !> Before entry, the incremented array X must contain the N- !> element vector x. !> Unchanged on exit. !>
[in]	INCX	!> INCX is INTEGER !> On entry, INCX specifies the increment for the elements of !> X. INCX must not be zero. !> Unchanged on exit. !>
[in]	BETA	!> BETA is COMPLEX*16 !> On entry, BETA specifies the scalar beta. When BETA is !> supplied as zero then Y need not be set on input. !> Unchanged on exit. !>
[in,out]	Y	!> Y is COMPLEX*16 array, dimension at least !> ( 1 + ( N - 1 )*abs( INCY ) ). !> Before entry, the incremented array Y must contain the n !> element vector y. On exit, Y is overwritten by the updated !> vector y. !>
[in]	INCY	!> INCY is INTEGER !> On entry, INCY specifies the increment for the elements of !> Y. INCY must not be zero. !> Unchanged on exit. !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 156 of file zsymv.f.

*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            INCX, INCY, LDA, N
      COMPLEX*16         ALPHA, BETA
*     ..
*     .. Array Arguments ..
      COMPLEX*16         A( LDA, * ), X( * ), Y( * )
*     ..
*
* =====================================================================
*
*     .. Parameters ..
      COMPLEX*16         ONE
      parameter( one = ( 1.0d+0, 0.0d+0 ) )
      COMPLEX*16         ZERO
      parameter( zero = ( 0.0d+0, 0.0d+0 ) )
*     ..
*     .. Local Scalars ..
      INTEGER            I, INFO, IX, IY, J, JX, JY, KX, KY
      COMPLEX*16         TEMP1, TEMP2
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL           xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          max
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      info = 0
      IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
         info = 1
      ELSE IF( n.LT.0 ) THEN
         info = 2
      ELSE IF( lda.LT.max( 1, n ) ) THEN
         info = 5
      ELSE IF( incx.EQ.0 ) THEN
         info = 7
      ELSE IF( incy.EQ.0 ) THEN
         info = 10
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'ZSYMV ', info )
         RETURN
      END IF
*
*     Quick return if possible.
*
      IF( ( n.EQ.0 ) .OR. ( ( alpha.EQ.zero ) .AND. ( beta.EQ.one ) ) )
     $   RETURN
*
*     Set up the start points in  X  and  Y.
*
      IF( incx.GT.0 ) THEN
         kx = 1
      ELSE
         kx = 1 - ( n-1 )*incx
      END IF
      IF( incy.GT.0 ) THEN
         ky = 1
      ELSE
         ky = 1 - ( n-1 )*incy
      END IF
*
*     Start the operations. In this version the elements of A are
*     accessed sequentially with one pass through the triangular part
*     of A.
*
*     First form  y := beta*y.
*
      IF( beta.NE.one ) THEN
         IF( incy.EQ.1 ) THEN
            IF( beta.EQ.zero ) THEN
               DO 10 i = 1, n
                  y( i ) = zero
   10          CONTINUE
            ELSE
               DO 20 i = 1, n
                  y( i ) = beta*y( i )
   20          CONTINUE
            END IF
         ELSE
            iy = ky
            IF( beta.EQ.zero ) THEN
               DO 30 i = 1, n
                  y( iy ) = zero
                  iy = iy + incy
   30          CONTINUE
            ELSE
               DO 40 i = 1, n
                  y( iy ) = beta*y( iy )
                  iy = iy + incy
   40          CONTINUE
            END IF
         END IF
      END IF
      IF( alpha.EQ.zero )
     $   RETURN
      IF( lsame( uplo, 'U' ) ) THEN
*
*        Form  y  when A is stored in upper triangle.
*
         IF( ( incx.EQ.1 ) .AND. ( incy.EQ.1 ) ) THEN
            DO 60 j = 1, n
               temp1 = alpha*x( j )
               temp2 = zero
               DO 50 i = 1, j - 1
                  y( i ) = y( i ) + temp1*a( i, j )
                  temp2 = temp2 + a( i, j )*x( i )
   50          CONTINUE
               y( j ) = y( j ) + temp1*a( j, j ) + alpha*temp2
   60       CONTINUE
         ELSE
            jx = kx
            jy = ky
            DO 80 j = 1, n
               temp1 = alpha*x( jx )
               temp2 = zero
               ix = kx
               iy = ky
               DO 70 i = 1, j - 1
                  y( iy ) = y( iy ) + temp1*a( i, j )
                  temp2 = temp2 + a( i, j )*x( ix )
                  ix = ix + incx
                  iy = iy + incy
   70          CONTINUE
               y( jy ) = y( jy ) + temp1*a( j, j ) + alpha*temp2
               jx = jx + incx
               jy = jy + incy
   80       CONTINUE
         END IF
      ELSE
*
*        Form  y  when A is stored in lower triangle.
*
         IF( ( incx.EQ.1 ) .AND. ( incy.EQ.1 ) ) THEN
            DO 100 j = 1, n
               temp1 = alpha*x( j )
               temp2 = zero
               y( j ) = y( j ) + temp1*a( j, j )
               DO 90 i = j + 1, n
                  y( i ) = y( i ) + temp1*a( i, j )
                  temp2 = temp2 + a( i, j )*x( i )
   90          CONTINUE
               y( j ) = y( j ) + alpha*temp2
  100       CONTINUE
         ELSE
            jx = kx
            jy = ky
            DO 120 j = 1, n
               temp1 = alpha*x( jx )
               temp2 = zero
               y( jy ) = y( jy ) + temp1*a( j, j )
               ix = jx
               iy = jy
               DO 110 i = j + 1, n
                  ix = ix + incx
                  iy = iy + incy
                  y( iy ) = y( iy ) + temp1*a( i, j )
                  temp2 = temp2 + a( i, j )*x( ix )
  110          CONTINUE
               y( jy ) = y( jy ) + alpha*temp2
               jx = jx + incx
               jy = jy + incy
  120       CONTINUE
         END IF
      END IF
*
      RETURN
*
*     End of ZSYMV
*

zsyr()

subroutine zsyr	(	character	uplo,
		integer	n,
		complex*16	alpha,
		complex*16, dimension( * )	x,
		integer	incx,
		complex*16, dimension( lda, * )	a,
		integer	lda )

ZSYR performs the symmetric rank-1 update of a complex symmetric matrix.

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

Purpose:

!>
!> ZSYR   performs the symmetric rank 1 operation
!>
!>    A := alpha*x*x**H + A,
!>
!> where alpha is a complex scalar, x is an n element vector and A is an
!> n by n symmetric matrix.
!> 

Parameters

[in]	UPLO	!> UPLO is CHARACTER*1 !> On entry, UPLO specifies whether the upper or lower !> triangular part of the array A is to be referenced as !> follows: !> !> UPLO = 'U' or 'u' Only the upper triangular part of A !> is to be referenced. !> !> UPLO = 'L' or 'l' Only the lower triangular part of A !> is to be referenced. !> !> Unchanged on exit. !>
[in]	N	!> N is INTEGER !> On entry, N specifies the order of the matrix A. !> N must be at least zero. !> Unchanged on exit. !>
[in]	ALPHA	!> ALPHA is COMPLEX*16 !> On entry, ALPHA specifies the scalar alpha. !> Unchanged on exit. !>
[in]	X	!> X is COMPLEX*16 array, dimension at least !> ( 1 + ( N - 1 )*abs( INCX ) ). !> Before entry, the incremented array X must contain the N- !> element vector x. !> Unchanged on exit. !>
[in]	INCX	!> INCX is INTEGER !> On entry, INCX specifies the increment for the elements of !> X. INCX must not be zero. !> Unchanged on exit. !>
[in,out]	A	!> A is COMPLEX*16 array, dimension ( LDA, N ) !> Before entry, with UPLO = 'U' or 'u', the leading n by n !> upper triangular part of the array A must contain the upper !> triangular part of the symmetric matrix and the strictly !> lower triangular part of A is not referenced. On exit, the !> upper triangular part of the array A is overwritten by the !> upper triangular part of the updated matrix. !> Before entry, with UPLO = 'L' or 'l', the leading n by n !> lower triangular part of the array A must contain the lower !> triangular part of the symmetric matrix and the strictly !> upper triangular part of A is not referenced. On exit, the !> lower triangular part of the array A is overwritten by the !> lower triangular part of the updated matrix. !>
[in]	LDA	!> LDA is INTEGER !> On entry, LDA specifies the first dimension of A as declared !> in the calling (sub) program. LDA must be at least !> max( 1, N ). !> Unchanged on exit. !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 134 of file zsyr.f.

*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            INCX, LDA, N
      COMPLEX*16         ALPHA
*     ..
*     .. Array Arguments ..
      COMPLEX*16         A( LDA, * ), X( * )
*     ..
*
* =====================================================================
*
*     .. Parameters ..
      COMPLEX*16         ZERO
      parameter( zero = ( 0.0d+0, 0.0d+0 ) )
*     ..
*     .. Local Scalars ..
      INTEGER            I, INFO, IX, J, JX, KX
      COMPLEX*16         TEMP
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL           xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          max
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      info = 0
      IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
         info = 1
      ELSE IF( n.LT.0 ) THEN
         info = 2
      ELSE IF( incx.EQ.0 ) THEN
         info = 5
      ELSE IF( lda.LT.max( 1, n ) ) THEN
         info = 7
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'ZSYR  ', info )
         RETURN
      END IF
*
*     Quick return if possible.
*
      IF( ( n.EQ.0 ) .OR. ( alpha.EQ.zero ) )
     $   RETURN
*
*     Set the start point in X if the increment is not unity.
*
      IF( incx.LE.0 ) THEN
         kx = 1 - ( n-1 )*incx
      ELSE IF( incx.NE.1 ) THEN
         kx = 1
      END IF
*
*     Start the operations. In this version the elements of A are
*     accessed sequentially with one pass through the triangular part
*     of A.
*
      IF( lsame( uplo, 'U' ) ) THEN
*
*        Form  A  when A is stored in upper triangle.
*
         IF( incx.EQ.1 ) THEN
            DO 20 j = 1, n
               IF( x( j ).NE.zero ) THEN
                  temp = alpha*x( j )
                  DO 10 i = 1, j
                     a( i, j ) = a( i, j ) + x( i )*temp
   10             CONTINUE
               END IF
   20       CONTINUE
         ELSE
            jx = kx
            DO 40 j = 1, n
               IF( x( jx ).NE.zero ) THEN
                  temp = alpha*x( jx )
                  ix = kx
                  DO 30 i = 1, j
                     a( i, j ) = a( i, j ) + x( ix )*temp
                     ix = ix + incx
   30             CONTINUE
               END IF
               jx = jx + incx
   40       CONTINUE
         END IF
      ELSE
*
*        Form  A  when A is stored in lower triangle.
*
         IF( incx.EQ.1 ) THEN
            DO 60 j = 1, n
               IF( x( j ).NE.zero ) THEN
                  temp = alpha*x( j )
                  DO 50 i = j, n
                     a( i, j ) = a( i, j ) + x( i )*temp
   50             CONTINUE
               END IF
   60       CONTINUE
         ELSE
            jx = kx
            DO 80 j = 1, n
               IF( x( jx ).NE.zero ) THEN
                  temp = alpha*x( jx )
                  ix = jx
                  DO 70 i = j, n
                     a( i, j ) = a( i, j ) + x( ix )*temp
                     ix = ix + incx
   70             CONTINUE
               END IF
               jx = jx + incx
   80       CONTINUE
         END IF
      END IF
*
      RETURN
*
*     End of ZSYR
*

zsyswapr()

subroutine zsyswapr	(	character	uplo,
		integer	n,
		complex*16, dimension( lda, n )	a,
		integer	lda,
		integer	i1,
		integer	i2 )

ZSYSWAPR

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

Purpose:

!>
!> ZSYSWAPR applies an elementary permutation on the rows and the columns of
!> a symmetric matrix.
!> 

Parameters

[in]	UPLO	!> UPLO is CHARACTER*1 !> Specifies whether the details of the factorization are stored !> as an upper or lower triangular matrix. !> = 'U': Upper triangular, form is A = U*D*U**T; !> = 'L': Lower triangular, form is A = L*D*L**T. !>
[in]	N	!> N is INTEGER !> The order of the matrix A. N >= 0. !>
[in,out]	A	!> A is COMPLEX*16 array, dimension (LDA,N) !> On entry, the NB diagonal matrix D and the multipliers !> used to obtain the factor U or L as computed by ZSYTRF. !> !> On exit, if INFO = 0, the (symmetric) inverse of the original !> matrix. If UPLO = 'U', the upper triangular part of the !> inverse is formed and the part of A below the diagonal is not !> referenced; if UPLO = 'L' the lower triangular part of the !> inverse is formed and the part of A above the diagonal is !> not referenced. !>
[in]	LDA	!> LDA is INTEGER !> The leading dimension of the array A. LDA >= max(1,N). !>
[in]	I1	!> I1 is INTEGER !> Index of the first row to swap !>
[in]	I2	!> I2 is INTEGER !> Index of the second row to swap !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 101 of file zsyswapr.f.

*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER        UPLO
      INTEGER          I1, I2, LDA, N
*     ..
*     .. Array Arguments ..
      COMPLEX*16       A( LDA, N )
*
*  =====================================================================
*
*     ..
*     .. Local Scalars ..
      LOGICAL            UPPER
      INTEGER            I
      COMPLEX*16         TMP
*
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL           zswap
*     ..
*     .. Executable Statements ..
*
      upper = lsame( uplo, 'U' )
      IF (upper) THEN
*
*         UPPER
*         first swap
*          - swap column I1 and I2 from I1 to I1-1
         CALL zswap( i1-1, a(1,i1), 1, a(1,i2), 1 )
*
*          second swap :
*          - swap A(I1,I1) and A(I2,I2)
*          - swap row I1 from I1+1 to I2-1 with col I2 from I1+1 to I2-1
         tmp=a(i1,i1)
         a(i1,i1)=a(i2,i2)
         a(i2,i2)=tmp
*
         DO i=1,i2-i1-1
            tmp=a(i1,i1+i)
            a(i1,i1+i)=a(i1+i,i2)
            a(i1+i,i2)=tmp
         END DO
*
*          third swap
*          - swap row I1 and I2 from I2+1 to N
         DO i=i2+1,n
            tmp=a(i1,i)
            a(i1,i)=a(i2,i)
            a(i2,i)=tmp
         END DO
*
        ELSE
*
*         LOWER
*         first swap
*          - swap row I1 and I2 from I1 to I1-1
         CALL zswap( i1-1, a(i1,1), lda, a(i2,1), lda )
*
*         second swap :
*          - swap A(I1,I1) and A(I2,I2)
*          - swap col I1 from I1+1 to I2-1 with row I2 from I1+1 to I2-1
          tmp=a(i1,i1)
          a(i1,i1)=a(i2,i2)
          a(i2,i2)=tmp
*
          DO i=1,i2-i1-1
             tmp=a(i1+i,i1)
             a(i1+i,i1)=a(i2,i1+i)
             a(i2,i1+i)=tmp
          END DO
*
*         third swap
*          - swap col I1 and I2 from I2+1 to N
          DO i=i2+1,n
             tmp=a(i,i1)
             a(i,i1)=a(i,i2)
             a(i,i2)=tmp
          END DO
*
      ENDIF

ztgsy2()

subroutine ztgsy2	(	character	trans,
		integer	ijob,
		integer	m,
		integer	n,
		complex*16, dimension( lda, * )	a,
		integer	lda,
		complex*16, dimension( ldb, * )	b,
		integer	ldb,
		complex*16, dimension( ldc, * )	c,
		integer	ldc,
		complex*16, dimension( ldd, * )	d,
		integer	ldd,
		complex*16, dimension( lde, * )	e,
		integer	lde,
		complex*16, dimension( ldf, * )	f,
		integer	ldf,
		double precision	scale,
		double precision	rdsum,
		double precision	rdscal,
		integer	info )

ZTGSY2 solves the generalized Sylvester equation (unblocked algorithm).

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

Purpose:

!>
!> ZTGSY2 solves the generalized Sylvester equation
!>
!>             A * R - L * B = scale * C               (1)
!>             D * R - L * E = scale * F
!>
!> using Level 1 and 2 BLAS, 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. A, B, D and E are upper triangular
!> (i.e., (A,D) and (B,E) in generalized Schur form).
!>
!> The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output
!> scaling factor chosen to avoid overflow.
!>
!> In matrix notation solving equation (1) corresponds to solve
!> Zx = scale * b, where Z is defined as
!>
!>        Z = [ kron(In, A)  -kron(B**H, Im) ]             (2)
!>            [ kron(In, D)  -kron(E**H, Im) ],
!>
!> Ik is the identity matrix of size k and X**H is the conjuguate transpose of X.
!> kron(X, Y) is the Kronecker product between the matrices X and Y.
!>
!> If TRANS = 'C', y in the conjugate transposed system Z**H*y = scale*b
!> is solved for, which is equivalent to solve for R and L in
!>
!>             A**H * R  + D**H * L   = scale * C           (3)
!>             R  * B**H + L  * E**H  = scale * -F
!>
!> This case is used to compute an estimate of Dif[(A, D), (B, E)] =
!> = sigma_min(Z) using reverse communication with ZLACON.
!>
!> ZTGSY2 also (IJOB >= 1) contributes to the computation in ZTGSYL
!> of an upper bound on the separation between to matrix pairs. Then
!> the input (A, D), (B, E) are sub-pencils of two matrix pairs in
!> ZTGSYL.
!> 

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: A contribution from this subsystem to a Frobenius !> norm-based estimate of the separation between two matrix !> pairs is computed. (look ahead strategy is used). !> =2: A contribution from this subsystem to a Frobenius !> norm-based estimate of the separation between two matrix !> pairs is computed. (DGECON on sub-systems is used.) !> Not referenced if TRANS = 'T'. !>
[in]	M	!> M is INTEGER !> On entry, M specifies the order of A and D, and the row !> dimension of C, F, R and L. !>
[in]	N	!> N is INTEGER !> On entry, N specifies the order of B and E, and the column !> dimension of C, F, R and L. !>
[in]	A	!> A is COMPLEX*16 array, dimension (LDA, M) !> On entry, A contains an upper triangular matrix. !>
[in]	LDA	!> LDA is INTEGER !> The leading dimension of the matrix A. LDA >= max(1, M). !>
[in]	B	!> B is COMPLEX*16 array, dimension (LDB, N) !> On entry, B contains an upper triangular matrix. !>
[in]	LDB	!> LDB is INTEGER !> The leading dimension of the matrix B. LDB >= max(1, N). !>
[in,out]	C	!> C is COMPLEX*16 array, dimension (LDC, N) !> On entry, C contains the right-hand-side of the first matrix !> equation in (1). !> On exit, if IJOB = 0, C has been overwritten by the solution !> R. !>
[in]	LDC	!> LDC is INTEGER !> The leading dimension of the matrix C. LDC >= max(1, M). !>
[in]	D	!> D is COMPLEX*16 array, dimension (LDD, M) !> On entry, D contains an upper triangular matrix. !>
[in]	LDD	!> LDD is INTEGER !> The leading dimension of the matrix D. LDD >= max(1, M). !>
[in]	E	!> E is COMPLEX*16 array, dimension (LDE, N) !> On entry, E contains an upper triangular matrix. !>
[in]	LDE	!> LDE is INTEGER !> The leading dimension of the matrix E. LDE >= max(1, N). !>
[in,out]	F	!> F is COMPLEX*16 array, dimension (LDF, N) !> On entry, F contains the right-hand-side of the second matrix !> equation in (1). !> On exit, if IJOB = 0, F has been overwritten by the solution !> L. !>
[in]	LDF	!> LDF is INTEGER !> The leading dimension of the matrix F. LDF >= max(1, M). !>
[out]	SCALE	!> SCALE is DOUBLE PRECISION !> On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions !> R and L (C and F on entry) will hold the solutions to a !> slightly perturbed system but the input matrices A, B, D and !> E have not been changed. If SCALE = 0, R and L will hold the !> solutions to the homogeneous system with C = F = 0. !> Normally, SCALE = 1. !>
[in,out]	RDSUM	!> RDSUM is DOUBLE PRECISION !> On entry, the sum of squares of computed contributions to !> the Dif-estimate under computation by ZTGSYL, where the !> scaling factor RDSCAL (see below) has been factored out. !> On exit, the corresponding sum of squares updated with the !> contributions from the current sub-system. !> If TRANS = 'T' RDSUM is not touched. !> NOTE: RDSUM only makes sense when ZTGSY2 is called by !> ZTGSYL. !>
[in,out]	RDSCAL	!> RDSCAL is DOUBLE PRECISION !> On entry, scaling factor used to prevent overflow in RDSUM. !> On exit, RDSCAL is updated w.r.t. the current contributions !> in RDSUM. !> If TRANS = 'T', RDSCAL is not touched. !> NOTE: RDSCAL only makes sense when ZTGSY2 is called by !> ZTGSYL. !>
[out]	INFO	!> INFO is INTEGER !> On exit, if INFO is set to !> =0: Successful exit !> <0: If INFO = -i, input argument number i is illegal. !> >0: The matrix pairs (A, D) and (B, E) have common or very !> 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.

Definition at line 256 of file ztgsy2.f.

*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          TRANS
      INTEGER            IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N
      DOUBLE PRECISION   RDSCAL, RDSUM, SCALE
*     ..
*     .. Array Arguments ..
      COMPLEX*16         A( LDA, * ), B( LDB, * ), C( LDC, * ),
     $                   D( LDD, * ), E( LDE, * ), F( LDF, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      INTEGER            LDZ
      parameter( zero = 0.0d+0, one = 1.0d+0, ldz = 2 )
*     ..
*     .. Local Scalars ..
      LOGICAL            NOTRAN
      INTEGER            I, IERR, J, K
      DOUBLE PRECISION   SCALOC
      COMPLEX*16         ALPHA
*     ..
*     .. Local Arrays ..
      INTEGER            IPIV( LDZ ), JPIV( LDZ )
      COMPLEX*16         RHS( LDZ ), Z( LDZ, LDZ )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL           xerbla, zaxpy, zgesc2, zgetc2, zlatdf, zscal
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          dcmplx, dconjg, max
*     ..
*     .. Executable Statements ..
*
*     Decode and test input parameters
*
      info = 0
      ierr = 0
      notran = lsame( trans, 'N' )
      IF( .NOT.notran .AND. .NOT.lsame( trans, 'C' ) ) THEN
         info = -1
      ELSE IF( notran ) THEN
         IF( ( ijob.LT.0 ) .OR. ( ijob.GT.2 ) ) THEN
            info = -2
         END IF
      END IF
      IF( info.EQ.0 ) THEN
         IF( m.LE.0 ) THEN
            info = -3
         ELSE IF( n.LE.0 ) THEN
            info = -4
         ELSE IF( lda.LT.max( 1, m ) ) THEN
            info = -6
         ELSE IF( ldb.LT.max( 1, n ) ) THEN
            info = -8
         ELSE IF( ldc.LT.max( 1, m ) ) THEN
            info = -10
         ELSE IF( ldd.LT.max( 1, m ) ) THEN
            info = -12
         ELSE IF( lde.LT.max( 1, n ) ) THEN
            info = -14
         ELSE IF( ldf.LT.max( 1, m ) ) THEN
            info = -16
         END IF
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'ZTGSY2', -info )
         RETURN
      END IF
*
      IF( notran ) THEN
*
*        Solve (I, J) - system
*           A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
*           D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
*        for I = M, M - 1, ..., 1; J = 1, 2, ..., N
*
         scale = one
         scaloc = one
         DO 30 j = 1, n
            DO 20 i = m, 1, -1
*
*              Build 2 by 2 system
*
               z( 1, 1 ) = a( i, i )
               z( 2, 1 ) = d( i, i )
               z( 1, 2 ) = -b( j, j )
               z( 2, 2 ) = -e( j, j )
*
*              Set up right hand side(s)
*
               rhs( 1 ) = c( i, j )
               rhs( 2 ) = f( i, j )
*
*              Solve Z * x = RHS
*
               CALL zgetc2( ldz, z, ldz, ipiv, jpiv, ierr )
               IF( ierr.GT.0 )
     $            info = ierr
               IF( ijob.EQ.0 ) THEN
                  CALL zgesc2( ldz, z, ldz, rhs, ipiv, jpiv, scaloc )
                  IF( scaloc.NE.one ) THEN
                     DO 10 k = 1, n
                        CALL zscal( m, dcmplx( scaloc, zero ),
     $                              c( 1, k ), 1 )
                        CALL zscal( m, dcmplx( scaloc, zero ),
     $                              f( 1, k ), 1 )
   10                CONTINUE
                     scale = scale*scaloc
                  END IF
               ELSE
                  CALL zlatdf( ijob, ldz, z, ldz, rhs, rdsum, rdscal,
     $                         ipiv, jpiv )
               END IF
*
*              Unpack solution vector(s)
*
               c( i, j ) = rhs( 1 )
               f( i, j ) = rhs( 2 )
*
*              Substitute R(I, J) and L(I, J) into remaining equation.
*
               IF( i.GT.1 ) THEN
                  alpha = -rhs( 1 )
                  CALL zaxpy( i-1, alpha, a( 1, i ), 1, c( 1, j ), 1 )
                  CALL zaxpy( i-1, alpha, d( 1, i ), 1, f( 1, j ), 1 )
               END IF
               IF( j.LT.n ) THEN
                  CALL zaxpy( n-j, rhs( 2 ), b( j, j+1 ), ldb,
     $                        c( i, j+1 ), ldc )
                  CALL zaxpy( n-j, rhs( 2 ), e( j, j+1 ), lde,
     $                        f( i, j+1 ), ldf )
               END IF
*
   20       CONTINUE
   30    CONTINUE
      ELSE
*
*        Solve transposed (I, J) - system:
*           A(I, I)**H * R(I, J) + D(I, I)**H * L(J, J) = C(I, J)
*           R(I, I) * B(J, J) + L(I, J) * E(J, J)   = -F(I, J)
*        for I = 1, 2, ..., M, J = N, N - 1, ..., 1
*
         scale = one
         scaloc = one
         DO 80 i = 1, m
            DO 70 j = n, 1, -1
*
*              Build 2 by 2 system Z**H
*
               z( 1, 1 ) = dconjg( a( i, i ) )
               z( 2, 1 ) = -dconjg( b( j, j ) )
               z( 1, 2 ) = dconjg( d( i, i ) )
               z( 2, 2 ) = -dconjg( e( j, j ) )
*
*
*              Set up right hand side(s)
*
               rhs( 1 ) = c( i, j )
               rhs( 2 ) = f( i, j )
*
*              Solve Z**H * x = RHS
*
               CALL zgetc2( ldz, z, ldz, ipiv, jpiv, ierr )
               IF( ierr.GT.0 )
     $            info = ierr
               CALL zgesc2( ldz, z, ldz, rhs, ipiv, jpiv, scaloc )
               IF( scaloc.NE.one ) THEN
                  DO 40 k = 1, n
                     CALL zscal( m, dcmplx( scaloc, zero ), c( 1, k ),
     $                           1 )
                     CALL zscal( m, dcmplx( scaloc, zero ), f( 1, k ),
     $                           1 )
   40             CONTINUE
                  scale = scale*scaloc
               END IF
*
*              Unpack solution vector(s)
*
               c( i, j ) = rhs( 1 )
               f( i, j ) = rhs( 2 )
*
*              Substitute R(I, J) and L(I, J) into remaining equation.
*
               DO 50 k = 1, j - 1
                  f( i, k ) = f( i, k ) + rhs( 1 )*dconjg( b( k, j ) ) +
     $                        rhs( 2 )*dconjg( e( k, j ) )
   50          CONTINUE
               DO 60 k = i + 1, m
                  c( k, j ) = c( k, j ) - dconjg( a( i, k ) )*rhs( 1 ) -
     $                        dconjg( d( i, k ) )*rhs( 2 )
   60          CONTINUE
*
   70       CONTINUE
   80    CONTINUE
      END IF
      RETURN
*
*     End of ZTGSY2
*