LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zhet22.f
Go to the documentation of this file.
00001 *> \brief \b ZHET22
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE ZHET22( ITYPE, UPLO, N, M, KBAND, A, LDA, D, E, U, LDU,
00012 *                          V, LDV, TAU, WORK, RWORK, RESULT )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          UPLO
00016 *       INTEGER            ITYPE, KBAND, LDA, LDU, LDV, M, N
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       DOUBLE PRECISION   D( * ), E( * ), RESULT( 2 ), RWORK( * )
00020 *       COMPLEX*16         A( LDA, * ), TAU( * ), U( LDU, * ),
00021 *      $                   V( LDV, * ), WORK( * )
00022 *       ..
00023 *  
00024 *
00025 *> \par Purpose:
00026 *  =============
00027 *>
00028 *> \verbatim
00029 *>
00030 *>      ZHET22  generally checks a decomposition of the form
00031 *>
00032 *>              A U = U S
00033 *>
00034 *>      where A is complex Hermitian, the columns of U are orthonormal,
00035 *>      and S is diagonal (if KBAND=0) or symmetric tridiagonal (if
00036 *>      KBAND=1).  If ITYPE=1, then U is represented as a dense matrix,
00037 *>      otherwise the U is expressed as a product of Householder
00038 *>      transformations, whose vectors are stored in the array "V" and
00039 *>      whose scaling constants are in "TAU"; we shall use the letter
00040 *>      "V" to refer to the product of Householder transformations
00041 *>      (which should be equal to U).
00042 *>
00043 *>      Specifically, if ITYPE=1, then:
00044 *>
00045 *>              RESULT(1) = | U' A U - S | / ( |A| m ulp ) *andC>              RESULT(2) = | I - U'U | / ( m ulp )
00046 *> \endverbatim
00047 *
00048 *  Arguments:
00049 *  ==========
00050 *
00051 *> \verbatim
00052 *>  ITYPE   INTEGER
00053 *>          Specifies the type of tests to be performed.
00054 *>          1: U expressed as a dense orthogonal matrix:
00055 *>             RESULT(1) = | A - U S U' | / ( |A| n ulp )   *andC>             RESULT(2) = | I - UU' | / ( n ulp )
00056 *>
00057 *>  UPLO    CHARACTER
00058 *>          If UPLO='U', the upper triangle of A will be used and the
00059 *>          (strictly) lower triangle will not be referenced.  If
00060 *>          UPLO='L', the lower triangle of A will be used and the
00061 *>          (strictly) upper triangle will not be referenced.
00062 *>          Not modified.
00063 *>
00064 *>  N       INTEGER
00065 *>          The size of the matrix.  If it is zero, ZHET22 does nothing.
00066 *>          It must be at least zero.
00067 *>          Not modified.
00068 *>
00069 *>  M       INTEGER
00070 *>          The number of columns of U.  If it is zero, ZHET22 does
00071 *>          nothing.  It must be at least zero.
00072 *>          Not modified.
00073 *>
00074 *>  KBAND   INTEGER
00075 *>          The bandwidth of the matrix.  It may only be zero or one.
00076 *>          If zero, then S is diagonal, and E is not referenced.  If
00077 *>          one, then S is symmetric tri-diagonal.
00078 *>          Not modified.
00079 *>
00080 *>  A       COMPLEX*16 array, dimension (LDA , N)
00081 *>          The original (unfactored) matrix.  It is assumed to be
00082 *>          symmetric, and only the upper (UPLO='U') or only the lower
00083 *>          (UPLO='L') will be referenced.
00084 *>          Not modified.
00085 *>
00086 *>  LDA     INTEGER
00087 *>          The leading dimension of A.  It must be at least 1
00088 *>          and at least N.
00089 *>          Not modified.
00090 *>
00091 *>  D       DOUBLE PRECISION array, dimension (N)
00092 *>          The diagonal of the (symmetric tri-) diagonal matrix.
00093 *>          Not modified.
00094 *>
00095 *>  E       DOUBLE PRECISION array, dimension (N)
00096 *>          The off-diagonal of the (symmetric tri-) diagonal matrix.
00097 *>          E(1) is ignored, E(2) is the (1,2) and (2,1) element, etc.
00098 *>          Not referenced if KBAND=0.
00099 *>          Not modified.
00100 *>
00101 *>  U       COMPLEX*16 array, dimension (LDU, N)
00102 *>          If ITYPE=1, this contains the orthogonal matrix in
00103 *>          the decomposition, expressed as a dense matrix.
00104 *>          Not modified.
00105 *>
00106 *>  LDU     INTEGER
00107 *>          The leading dimension of U.  LDU must be at least N and
00108 *>          at least 1.
00109 *>          Not modified.
00110 *>
00111 *>  V       COMPLEX*16 array, dimension (LDV, N)
00112 *>          If ITYPE=2 or 3, the lower triangle of this array contains
00113 *>          the Householder vectors used to describe the orthogonal
00114 *>          matrix in the decomposition.  If ITYPE=1, then it is not
00115 *>          referenced.
00116 *>          Not modified.
00117 *>
00118 *>  LDV     INTEGER
00119 *>          The leading dimension of V.  LDV must be at least N and
00120 *>          at least 1.
00121 *>          Not modified.
00122 *>
00123 *>  TAU     COMPLEX*16 array, dimension (N)
00124 *>          If ITYPE >= 2, then TAU(j) is the scalar factor of
00125 *>          v(j) v(j)' in the Householder transformation H(j) of
00126 *>          the product  U = H(1)...H(n-2)
00127 *>          If ITYPE < 2, then TAU is not referenced.
00128 *>          Not modified.
00129 *>
00130 *>  WORK    COMPLEX*16 array, dimension (2*N**2)
00131 *>          Workspace.
00132 *>          Modified.
00133 *>
00134 *>  RWORK   DOUBLE PRECISION array, dimension (N)
00135 *>          Workspace.
00136 *>          Modified.
00137 *>
00138 *>  RESULT  DOUBLE PRECISION array, dimension (2)
00139 *>          The values computed by the two tests described above.  The
00140 *>          values are currently limited to 1/ulp, to avoid overflow.
00141 *>          RESULT(1) is always modified.  RESULT(2) is modified only
00142 *>          if LDU is at least N.
00143 *>          Modified.
00144 *> \endverbatim
00145 *
00146 *  Authors:
00147 *  ========
00148 *
00149 *> \author Univ. of Tennessee 
00150 *> \author Univ. of California Berkeley 
00151 *> \author Univ. of Colorado Denver 
00152 *> \author NAG Ltd. 
00153 *
00154 *> \date November 2011
00155 *
00156 *> \ingroup complex16_eig
00157 *
00158 *  =====================================================================
00159       SUBROUTINE ZHET22( ITYPE, UPLO, N, M, KBAND, A, LDA, D, E, U, LDU,
00160      $                   V, LDV, TAU, WORK, RWORK, RESULT )
00161 *
00162 *  -- LAPACK test routine (version 3.4.0) --
00163 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00164 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00165 *     November 2011
00166 *
00167 *     .. Scalar Arguments ..
00168       CHARACTER          UPLO
00169       INTEGER            ITYPE, KBAND, LDA, LDU, LDV, M, N
00170 *     ..
00171 *     .. Array Arguments ..
00172       DOUBLE PRECISION   D( * ), E( * ), RESULT( 2 ), RWORK( * )
00173       COMPLEX*16         A( LDA, * ), TAU( * ), U( LDU, * ),
00174      $                   V( LDV, * ), WORK( * )
00175 *     ..
00176 *
00177 *  =====================================================================
00178 *
00179 *     .. Parameters ..
00180       DOUBLE PRECISION   ZERO, ONE
00181       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00182       COMPLEX*16         CZERO, CONE
00183       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
00184      $                   CONE = ( 1.0D0, 0.0D0 ) )
00185 *     ..
00186 *     .. Local Scalars ..
00187       INTEGER            J, JJ, JJ1, JJ2, NN, NNP1
00188       DOUBLE PRECISION   ANORM, ULP, UNFL, WNORM
00189 *     ..
00190 *     .. External Functions ..
00191       DOUBLE PRECISION   DLAMCH, ZLANHE
00192       EXTERNAL           DLAMCH, ZLANHE
00193 *     ..
00194 *     .. External Subroutines ..
00195       EXTERNAL           ZGEMM, ZHEMM, ZUNT01
00196 *     ..
00197 *     .. Intrinsic Functions ..
00198       INTRINSIC          DBLE, MAX, MIN
00199 *     ..
00200 *     .. Executable Statements ..
00201 *
00202       RESULT( 1 ) = ZERO
00203       RESULT( 2 ) = ZERO
00204       IF( N.LE.0 .OR. M.LE.0 )
00205      $   RETURN
00206 *
00207       UNFL = DLAMCH( 'Safe minimum' )
00208       ULP = DLAMCH( 'Precision' )
00209 *
00210 *     Do Test 1
00211 *
00212 *     Norm of A:
00213 *
00214       ANORM = MAX( ZLANHE( '1', UPLO, N, A, LDA, RWORK ), UNFL )
00215 *
00216 *     Compute error matrix:
00217 *
00218 *     ITYPE=1: error = U' A U - S
00219 *
00220       CALL ZHEMM( 'L', UPLO, N, M, CONE, A, LDA, U, LDU, CZERO, WORK,
00221      $            N )
00222       NN = N*N
00223       NNP1 = NN + 1
00224       CALL ZGEMM( 'C', 'N', M, M, N, CONE, U, LDU, WORK, N, CZERO,
00225      $            WORK( NNP1 ), N )
00226       DO 10 J = 1, M
00227          JJ = NN + ( J-1 )*N + J
00228          WORK( JJ ) = WORK( JJ ) - D( J )
00229    10 CONTINUE
00230       IF( KBAND.EQ.1 .AND. N.GT.1 ) THEN
00231          DO 20 J = 2, M
00232             JJ1 = NN + ( J-1 )*N + J - 1
00233             JJ2 = NN + ( J-2 )*N + J
00234             WORK( JJ1 ) = WORK( JJ1 ) - E( J-1 )
00235             WORK( JJ2 ) = WORK( JJ2 ) - E( J-1 )
00236    20    CONTINUE
00237       END IF
00238       WNORM = ZLANHE( '1', UPLO, M, WORK( NNP1 ), N, RWORK )
00239 *
00240       IF( ANORM.GT.WNORM ) THEN
00241          RESULT( 1 ) = ( WNORM / ANORM ) / ( M*ULP )
00242       ELSE
00243          IF( ANORM.LT.ONE ) THEN
00244             RESULT( 1 ) = ( MIN( WNORM, M*ANORM ) / ANORM ) / ( M*ULP )
00245          ELSE
00246             RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( M ) ) / ( M*ULP )
00247          END IF
00248       END IF
00249 *
00250 *     Do Test 2
00251 *
00252 *     Compute  U'U - I
00253 *
00254       IF( ITYPE.EQ.1 )
00255      $   CALL ZUNT01( 'Columns', N, M, U, LDU, WORK, 2*N*N, RWORK,
00256      $                RESULT( 2 ) )
00257 *
00258       RETURN
00259 *
00260 *     End of ZHET22
00261 *
00262       END
 All Files Functions