complex

Functions
subroutine	claesy (a, b, c, rt1, rt2, evscal, cs1, sn1)
	CLAESY computes the eigenvalues and eigenvectors of a 2-by-2 complex symmetric matrix.
real function	clansy (norm, uplo, n, a, lda, work)
	CLANSY 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	claqsy (uplo, n, a, lda, s, scond, amax, equed)
	CLAQSY scales a symmetric/Hermitian matrix, using scaling factors computed by spoequ.
subroutine	csymv (uplo, n, alpha, a, lda, x, incx, beta, y, incy)
	CSYMV computes a matrix-vector product for a complex symmetric matrix.
subroutine	csyr (uplo, n, alpha, x, incx, a, lda)
	CSYR performs the symmetric rank-1 update of a complex symmetric matrix.
subroutine	csyswapr (uplo, n, a, lda, i1, i2)
	CSYSWAPR
subroutine	ctgsy2 (trans, ijob, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, rdsum, rdscal, info)
	CTGSY2 solves the generalized Sylvester equation (unblocked algorithm).

Detailed Description

This is the group of complex auxiliary functions for SY matrices

Function Documentation

claesy()

subroutine claesy	(	complex	a,
		complex	b,
		complex	c,
		complex	rt1,
		complex	rt2,
		complex	evscal,
		complex	cs1,
		complex	sn1 )

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

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

Purpose:

!>
!> CLAESY 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 !> The ( 1, 1 ) element of input matrix. !>
[in]	B	!> B is COMPLEX !> 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 !> The ( 2, 2 ) element of input matrix. !>
[out]	RT1	!> RT1 is COMPLEX !> The eigenvalue of larger modulus. !>
[out]	RT2	!> RT2 is COMPLEX !> The eigenvalue of smaller modulus. !>
[out]	EVSCAL	!> EVSCAL is COMPLEX !> 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 !>
[out]	SN1	!> SN1 is COMPLEX !> 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 claesy.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            A, B, C, CS1, EVSCAL, RT1, RT2, SN1
*     ..
*
* =====================================================================
*
*     .. Parameters ..
      REAL               ZERO
      parameter( zero = 0.0e0 )
      REAL               ONE
      parameter( one = 1.0e0 )
      COMPLEX            CONE
      parameter( cone = ( 1.0e0, 0.0e0 ) )
      REAL               HALF
      parameter( half = 0.5e0 )
      REAL               THRESH
      parameter( thresh = 0.1e0 )
*     ..
*     .. Local Scalars ..
      REAL               BABS, EVNORM, TABS, Z
      COMPLEX            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 CLAESY
*

clansy()

real function clansy	(	character	norm,
		character	uplo,
		integer	n,
		complex, dimension( lda, * )	a,
		integer	lda,
		real, dimension( * )	work )

CLANSY 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 CLANSY + dependencies [TGZ] [ZIP] [TXT]

Purpose:

!>
!> CLANSY  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

CLANSY

!>
!>    CLANSY = ( 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 CLANSY 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, CLANSY is !> set to zero. !>
[in]	A	!> A is COMPLEX 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 REAL 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 clansy.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 ..
      REAL               WORK( * )
      COMPLEX            A( LDA, * )
*     ..
*
* =====================================================================
*
*     .. Parameters ..
      REAL               ONE, ZERO
      parameter( one = 1.0e+0, zero = 0.0e+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, J
      REAL               ABSA, SCALE, SUM, VALUE
*     ..
*     .. External Functions ..
      LOGICAL            LSAME, SISNAN
      EXTERNAL           lsame, sisnan
*     ..
*     .. External Subroutines ..
      EXTERNAL           classq
*     ..
*     .. 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. sisnan( 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. sisnan( 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. sisnan( 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. sisnan( 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 classq( j-1, a( 1, j ), 1, scale, sum )
  110       CONTINUE
         ELSE
            DO 120 j = 1, n - 1
               CALL classq( n-j, a( j+1, j ), 1, scale, sum )
  120       CONTINUE
         END IF
         sum = 2*sum
         CALL classq( n, a, lda+1, scale, sum )
         VALUE = scale*sqrt( sum )
      END IF
*
      clansy = VALUE
      RETURN
*
*     End of CLANSY
*

claqsy()

subroutine claqsy	(	character	uplo,
		integer	n,
		complex, dimension( lda, * )	a,
		integer	lda,
		real, dimension( * )	s,
		real	scond,
		real	amax,
		character	equed )

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

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

Purpose:

!>
!> CLAQSY 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 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 REAL array, dimension (N) !> The scale factors for A. !>
[in]	SCOND	!> SCOND is REAL !> Ratio of the smallest S(i) to the largest S(i). !>
[in]	AMAX	!> AMAX is REAL !> 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 claqsy.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
      REAL               AMAX, SCOND
*     ..
*     .. Array Arguments ..
      REAL               S( * )
      COMPLEX            A( LDA, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ONE, THRESH
      parameter( one = 1.0e+0, thresh = 0.1e+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, J
      REAL               CJ, LARGE, SMALL
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      REAL               SLAMCH
      EXTERNAL           lsame, slamch
*     ..
*     .. Executable Statements ..
*
*     Quick return if possible
*
      IF( n.LE.0 ) THEN
         equed = 'N'
         RETURN
      END IF
*
*     Initialize LARGE and SMALL.
*
      small = slamch( 'Safe minimum' ) / slamch( '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 CLAQSY
*

csymv()

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

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

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

Purpose:

!>
!> CSYMV  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 !> On entry, ALPHA specifies the scalar alpha. !> Unchanged on exit. !>
[in]	A	!> A is COMPLEX 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 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 !> 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 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 csymv.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            ALPHA, BETA
*     ..
*     .. Array Arguments ..
      COMPLEX            A( LDA, * ), X( * ), Y( * )
*     ..
*
* =====================================================================
*
*     .. Parameters ..
      COMPLEX            ONE
      parameter( one = ( 1.0e+0, 0.0e+0 ) )
      COMPLEX            ZERO
      parameter( zero = ( 0.0e+0, 0.0e+0 ) )
*     ..
*     .. Local Scalars ..
      INTEGER            I, INFO, IX, IY, J, JX, JY, KX, KY
      COMPLEX            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( 'CSYMV ', 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 CSYMV
*

csyr()

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

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

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

Purpose:

!>
!> CSYR   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 !> On entry, ALPHA specifies the scalar alpha. !> Unchanged on exit. !>
[in]	X	!> X is COMPLEX 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 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 csyr.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            ALPHA
*     ..
*     .. Array Arguments ..
      COMPLEX            A( LDA, * ), X( * )
*     ..
*
* =====================================================================
*
*     .. Parameters ..
      COMPLEX            ZERO
      parameter( zero = ( 0.0e+0, 0.0e+0 ) )
*     ..
*     .. Local Scalars ..
      INTEGER            I, INFO, IX, J, JX, KX
      COMPLEX            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( 'CSYR  ', 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 CSYR
*

csyswapr()

subroutine csyswapr	(	character	uplo,
		integer	n,
		complex, dimension( lda, n )	a,
		integer	lda,
		integer	i1,
		integer	i2 )

CSYSWAPR

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

Purpose:

!>
!> CSYSWAPR 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 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 CSYTRF. !> !> 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 csyswapr.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          A( LDA, N )
*
*  =====================================================================
*
*     ..
*     .. Local Scalars ..
      LOGICAL            UPPER
      INTEGER            I
      COMPLEX            TMP
*
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL           cswap
*     ..
*     .. Executable Statements ..
*
      upper = lsame( uplo, 'U' )
      IF (upper) THEN
*
*         UPPER
*         first swap
*          - swap column I1 and I2 from I1 to I1-1
         CALL cswap( 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 cswap ( 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

ctgsy2()

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

CTGSY2 solves the generalized Sylvester equation (unblocked algorithm).

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

Purpose:

!>
!> CTGSY2 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 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 CLACON.
!>
!> CTGSY2 also (IJOB >= 1) contributes to the computation in CTGSYL
!> 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
!> CTGSYL.
!> 

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. (SGECON 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 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 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 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 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 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 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 REAL !> 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 REAL !> On entry, the sum of squares of computed contributions to !> the Dif-estimate under computation by CTGSYL, 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 CTGSY2 is called by !> CTGSYL. !>
[in,out]	RDSCAL	!> RDSCAL is REAL !> 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 CTGSY2 is called by !> CTGSYL. !>
[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 ctgsy2.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
      REAL               RDSCAL, RDSUM, SCALE
*     ..
*     .. Array Arguments ..
      COMPLEX            A( LDA, * ), B( LDB, * ), C( LDC, * ),
     $                   D( LDD, * ), E( LDE, * ), F( LDF, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      INTEGER            LDZ
      parameter( zero = 0.0e+0, one = 1.0e+0, ldz = 2 )
*     ..
*     .. Local Scalars ..
      LOGICAL            NOTRAN
      INTEGER            I, IERR, J, K
      REAL               SCALOC
      COMPLEX            ALPHA
*     ..
*     .. Local Arrays ..
      INTEGER            IPIV( LDZ ), JPIV( LDZ )
      COMPLEX            RHS( LDZ ), Z( LDZ, LDZ )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL           caxpy, cgesc2, cgetc2, cscal, clatdf, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          cmplx, conjg, 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( 'CTGSY2', -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 cgetc2( ldz, z, ldz, ipiv, jpiv, ierr )
               IF( ierr.GT.0 )
     $            info = ierr
               IF( ijob.EQ.0 ) THEN
                  CALL cgesc2( ldz, z, ldz, rhs, ipiv, jpiv, scaloc )
                  IF( scaloc.NE.one ) THEN
                     DO 10 k = 1, n
                        CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
     $                              1 )
                        CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
     $                              1 )
   10                CONTINUE
                     scale = scale*scaloc
                  END IF
               ELSE
                  CALL clatdf( 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 caxpy( i-1, alpha, a( 1, i ), 1, c( 1, j ), 1 )
                  CALL caxpy( i-1, alpha, d( 1, i ), 1, f( 1, j ), 1 )
               END IF
               IF( j.LT.n ) THEN
                  CALL caxpy( n-j, rhs( 2 ), b( j, j+1 ), ldb,
     $                        c( i, j+1 ), ldc )
                  CALL caxpy( 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 ) = conjg( a( i, i ) )
               z( 2, 1 ) = -conjg( b( j, j ) )
               z( 1, 2 ) = conjg( d( i, i ) )
               z( 2, 2 ) = -conjg( 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 cgetc2( ldz, z, ldz, ipiv, jpiv, ierr )
               IF( ierr.GT.0 )
     $            info = ierr
               CALL cgesc2( ldz, z, ldz, rhs, ipiv, jpiv, scaloc )
               IF( scaloc.NE.one ) THEN
                  DO 40 k = 1, n
                     CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
     $                           1 )
                     CALL cscal( m, cmplx( 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 )*conjg( b( k, j ) ) +
     $                        rhs( 2 )*conjg( e( k, j ) )
   50          CONTINUE
               DO 60 k = i + 1, m
                  c( k, j ) = c( k, j ) - conjg( a( i, k ) )*rhs( 1 ) -
     $                        conjg( d( i, k ) )*rhs( 2 )
   60          CONTINUE
*
   70       CONTINUE
   80    CONTINUE
      END IF
      RETURN
*
*     End of CTGSY2
*