zlatm6.f

Go to the documentation of this file.

*> \brief \b ZLATM6
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*  Definition:
*  ===========
*
*       SUBROUTINE ZLATM6( TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA,
*                          BETA, WX, WY, S, DIF )
* 
*       .. Scalar Arguments ..
*       INTEGER            LDA, LDX, LDY, N, TYPE
*       COMPLEX*16         ALPHA, BETA, WX, WY
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   DIF( * ), S( * )
*       COMPLEX*16         A( LDA, * ), B( LDA, * ), X( LDX, * ),
*      $                   Y( LDY, * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> ZLATM6 generates test matrices for the generalized eigenvalue
*> problem, their corresponding right and left eigenvector matrices,
*> and also reciprocal condition numbers for all eigenvalues and
*> the reciprocal condition numbers of eigenvectors corresponding to
*> the 1th and 5th eigenvalues.
*>
*> Test Matrices
*> =============
*>
*> Two kinds of test matrix pairs
*>          (A, B) = inverse(YH) * (Da, Db) * inverse(X)
*> are used in the tests:
*>
*> Type 1:
*>    Da = 1+a   0    0    0    0    Db = 1   0   0   0   0
*>          0   2+a   0    0    0         0   1   0   0   0
*>          0    0   3+a   0    0         0   0   1   0   0
*>          0    0    0   4+a   0         0   0   0   1   0
*>          0    0    0    0   5+a ,      0   0   0   0   1
*> and Type 2:
*>    Da = 1+i   0    0       0       0    Db = 1   0   0   0   0
*>          0   1-i   0       0       0         0   1   0   0   0
*>          0    0    1       0       0         0   0   1   0   0
*>          0    0    0 (1+a)+(1+b)i  0         0   0   0   1   0
*>          0    0    0       0 (1+a)-(1+b)i,   0   0   0   0   1 .
*>
*> In both cases the same inverse(YH) and inverse(X) are used to compute
*> (A, B), giving the exact eigenvectors to (A,B) as (YH, X):
*>
*> YH:  =  1    0   -y    y   -y    X =  1   0  -x  -x   x
*>         0    1   -y    y   -y         0   1   x  -x  -x
*>         0    0    1    0    0         0   0   1   0   0
*>         0    0    0    1    0         0   0   0   1   0
*>         0    0    0    0    1,        0   0   0   0   1 , where
*>
*> a, b, x and y will have all values independently of each other.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] TYPE
*> \verbatim
*>          TYPE is INTEGER
*>          Specifies the problem type (see futher details).
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          Size of the matrices A and B.
*> \endverbatim
*>
*> \param[out] A
*> \verbatim
*>          A is COMPLEX*16 array, dimension (LDA, N).
*>          On exit A N-by-N is initialized according to TYPE.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of A and of B.
*> \endverbatim
*>
*> \param[out] B
*> \verbatim
*>          B is COMPLEX*16 array, dimension (LDA, N).
*>          On exit B N-by-N is initialized according to TYPE.
*> \endverbatim
*>
*> \param[out] X
*> \verbatim
*>          X is COMPLEX*16 array, dimension (LDX, N).
*>          On exit X is the N-by-N matrix of right eigenvectors.
*> \endverbatim
*>
*> \param[in] LDX
*> \verbatim
*>          LDX is INTEGER
*>          The leading dimension of X.
*> \endverbatim
*>
*> \param[out] Y
*> \verbatim
*>          Y is COMPLEX*16 array, dimension (LDY, N).
*>          On exit Y is the N-by-N matrix of left eigenvectors.
*> \endverbatim
*>
*> \param[in] LDY
*> \verbatim
*>          LDY is INTEGER
*>          The leading dimension of Y.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*>          ALPHA is COMPLEX*16
*> \endverbatim
*>
*> \param[in] BETA
*> \verbatim
*>          BETA is COMPLEX*16
*> \verbatim
*>          Weighting constants for matrix A.
*> \endverbatim
*>
*> \param[in] WX
*> \verbatim
*>          WX is COMPLEX*16
*>          Constant for right eigenvector matrix.
*> \endverbatim
*>
*> \param[in] WY
*> \verbatim
*>          WY is COMPLEX*16
*>          Constant for left eigenvector matrix.
*> \endverbatim
*>
*> \param[out] S
*> \verbatim
*>          S is DOUBLE PRECISION array, dimension (N)
*>          S(i) is the reciprocal condition number for eigenvalue i.
*> \endverbatim
*>
*> \param[out] DIF
*> \verbatim
*>          DIF is DOUBLE PRECISION array, dimension (N)
*>          DIF(i) is the reciprocal condition number for eigenvector i.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date November 2011
*
*> \ingroup complex16_matgen
*
*  =====================================================================
      SUBROUTINE zlatm6( TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA,
     $                   beta, wx, wy, s, dif )
*
*  -- LAPACK computational routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      INTEGER            LDA, LDX, LDY, N, TYPE
      COMPLEX*16         ALPHA, BETA, WX, WY
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   DIF( * ), S( * )
      COMPLEX*16         A( lda, * ), B( lda, * ), X( ldx, * ),
     $                   y( ldy, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   RONE, TWO, THREE
      parameter                ( rone = 1.0d+0, two = 2.0d+0, three = 3.0d+0 )
      COMPLEX*16         ZERO, ONE
      parameter                ( zero = ( 0.0d+0, 0.0d+0 ),
     $                   one = ( 1.0d+0, 0.0d+0 ) )
*     ..
*     .. Local Scalars ..
      INTEGER            I, INFO, J
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   RWORK( 50 )
      COMPLEX*16         WORK( 26 ), Z( 8, 8 )
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          cdabs, dble, dcmplx, dconjg, sqrt
*     ..
*     .. External Subroutines ..
      EXTERNAL           zgesvd, zlacpy, zlakf2
*     ..
*     .. Executable Statements ..
*
*     Generate test problem ...
*     (Da, Db) ...
*
      DO 20 i = 1, n
         DO 10 j = 1, n
*
            IF( i.EQ.j ) THEN
               a( i, i ) = dcmplx( i ) + alpha
               b( i, i ) = one
            ELSE
               a( i, j ) = zero
               b( i, j ) = zero
            END IF
*
   10    CONTINUE
   20 CONTINUE
      IF( type.EQ.2 ) THEN
         a( 1, 1 ) = dcmplx( rone, rone )
         a( 2, 2 ) = dconjg( a( 1, 1 ) )
         a( 3, 3 ) = one
         a( 4, 4 ) = dcmplx( dble( one+alpha ), dble( one+beta ) )
         a( 5, 5 ) = dconjg( a( 4, 4 ) )
      END IF
*
*     Form X and Y
*
      CALL zlacpy( 'F', n, n, b, lda, y, ldy )
      y( 3, 1 ) = -dconjg( wy )
      y( 4, 1 ) = dconjg( wy )
      y( 5, 1 ) = -dconjg( wy )
      y( 3, 2 ) = -dconjg( wy )
      y( 4, 2 ) = dconjg( wy )
      y( 5, 2 ) = -dconjg( wy )
*
      CALL zlacpy( 'F', n, n, b, lda, x, ldx )
      x( 1, 3 ) = -wx
      x( 1, 4 ) = -wx
      x( 1, 5 ) = wx
      x( 2, 3 ) = wx
      x( 2, 4 ) = -wx
      x( 2, 5 ) = -wx
*
*     Form (A, B)
*
      b( 1, 3 ) = wx + wy
      b( 2, 3 ) = -wx + wy
      b( 1, 4 ) = wx - wy
      b( 2, 4 ) = wx - wy
      b( 1, 5 ) = -wx + wy
      b( 2, 5 ) = wx + wy
      a( 1, 3 ) = wx*a( 1, 1 ) + wy*a( 3, 3 )
      a( 2, 3 ) = -wx*a( 2, 2 ) + wy*a( 3, 3 )
      a( 1, 4 ) = wx*a( 1, 1 ) - wy*a( 4, 4 )
      a( 2, 4 ) = wx*a( 2, 2 ) - wy*a( 4, 4 )
      a( 1, 5 ) = -wx*a( 1, 1 ) + wy*a( 5, 5 )
      a( 2, 5 ) = wx*a( 2, 2 ) + wy*a( 5, 5 )
*
*     Compute condition numbers
*
      s( 1 ) = rone / sqrt( ( rone+three*cdabs( wy )*cdabs( wy ) ) /
     $         ( rone+cdabs( a( 1, 1 ) )*cdabs( a( 1, 1 ) ) ) )
      s( 2 ) = rone / sqrt( ( rone+three*cdabs( wy )*cdabs( wy ) ) /
     $         ( rone+cdabs( a( 2, 2 ) )*cdabs( a( 2, 2 ) ) ) )
      s( 3 ) = rone / sqrt( ( rone+two*cdabs( wx )*cdabs( wx ) ) /
     $         ( rone+cdabs( a( 3, 3 ) )*cdabs( a( 3, 3 ) ) ) )
      s( 4 ) = rone / sqrt( ( rone+two*cdabs( wx )*cdabs( wx ) ) /
     $         ( rone+cdabs( a( 4, 4 ) )*cdabs( a( 4, 4 ) ) ) )
      s( 5 ) = rone / sqrt( ( rone+two*cdabs( wx )*cdabs( wx ) ) /
     $         ( rone+cdabs( a( 5, 5 ) )*cdabs( a( 5, 5 ) ) ) )
*
      CALL zlakf2( 1, 4, a, lda, a( 2, 2 ), b, b( 2, 2 ), z, 8 )
      CALL zgesvd( 'N', 'N', 8, 8, z, 8, rwork, work, 1, work( 2 ), 1,
     $             work( 3 ), 24, rwork( 9 ), info )
      dif( 1 ) = rwork( 8 )
*
      CALL zlakf2( 4, 1, a, lda, a( 5, 5 ), b, b( 5, 5 ), z, 8 )
      CALL zgesvd( 'N', 'N', 8, 8, z, 8, rwork, work, 1, work( 2 ), 1,
     $             work( 3 ), 24, rwork( 9 ), info )
      dif( 5 ) = rwork( 8 )
*
      RETURN
*
*     End of ZLATM6
*
      END