LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zhet21.f
Go to the documentation of this file.
00001 *> \brief \b ZHET21
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 ZHET21( ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V,
00012 *                          LDV, TAU, WORK, RWORK, RESULT )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          UPLO
00016 *       INTEGER            ITYPE, KBAND, LDA, LDU, LDV, 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 *> ZHET21 generally checks a decomposition of the form
00031 *>
00032 *>    A = U S UC>
00033 *> where * means conjugate transpose, A is hermitian, U is unitary, and
00034 *> S is diagonal (if KBAND=0) or (real) symmetric tridiagonal (if
00035 *> KBAND=1).
00036 *>
00037 *> If ITYPE=1, then U is represented as a dense matrix; otherwise U is
00038 *> expressed as a product of Householder transformations, whose vectors
00039 *> are stored in the array "V" and whose scaling constants are in "TAU".
00040 *> We shall use the letter "V" to refer to the product of Householder
00041 *> transformations (which should be equal to U).
00042 *>
00043 *> Specifically, if ITYPE=1, then:
00044 *>
00045 *>    RESULT(1) = | A - U S U* | / ( |A| n ulp ) *andC>    RESULT(2) = | I - UU* | / ( n ulp )
00046 *>
00047 *> If ITYPE=2, then:
00048 *>
00049 *>    RESULT(1) = | A - V S V* | / ( |A| n ulp )
00050 *>
00051 *> If ITYPE=3, then:
00052 *>
00053 *>    RESULT(1) = | I - UV* | / ( n ulp )
00054 *>
00055 *> For ITYPE > 1, the transformation U is expressed as a product
00056 *> V = H(1)...H(n-2),  where H(j) = I  -  tau(j) v(j) v(j)C> and each
00057 *> vector v(j) has its first j elements 0 and the remaining n-j elements
00058 *> stored in V(j+1:n,j).
00059 *> \endverbatim
00060 *
00061 *  Arguments:
00062 *  ==========
00063 *
00064 *> \param[in] ITYPE
00065 *> \verbatim
00066 *>          ITYPE is INTEGER
00067 *>          Specifies the type of tests to be performed.
00068 *>          1: U expressed as a dense unitary matrix:
00069 *>             RESULT(1) = | A - U S U* | / ( |A| n ulp )   *andC>             RESULT(2) = | I - UU* | / ( n ulp )
00070 *>
00071 *>          2: U expressed as a product V of Housholder transformations:
00072 *>             RESULT(1) = | A - V S V* | / ( |A| n ulp )
00073 *>
00074 *>          3: U expressed both as a dense unitary matrix and
00075 *>             as a product of Housholder transformations:
00076 *>             RESULT(1) = | I - UV* | / ( n ulp )
00077 *> \endverbatim
00078 *>
00079 *> \param[in] UPLO
00080 *> \verbatim
00081 *>          UPLO is CHARACTER
00082 *>          If UPLO='U', the upper triangle of A and V will be used and
00083 *>          the (strictly) lower triangle will not be referenced.
00084 *>          If UPLO='L', the lower triangle of A and V will be used and
00085 *>          the (strictly) upper triangle will not be referenced.
00086 *> \endverbatim
00087 *>
00088 *> \param[in] N
00089 *> \verbatim
00090 *>          N is INTEGER
00091 *>          The size of the matrix.  If it is zero, ZHET21 does nothing.
00092 *>          It must be at least zero.
00093 *> \endverbatim
00094 *>
00095 *> \param[in] KBAND
00096 *> \verbatim
00097 *>          KBAND is INTEGER
00098 *>          The bandwidth of the matrix.  It may only be zero or one.
00099 *>          If zero, then S is diagonal, and E is not referenced.  If
00100 *>          one, then S is symmetric tri-diagonal.
00101 *> \endverbatim
00102 *>
00103 *> \param[in] A
00104 *> \verbatim
00105 *>          A is COMPLEX*16 array, dimension (LDA, N)
00106 *>          The original (unfactored) matrix.  It is assumed to be
00107 *>          hermitian, and only the upper (UPLO='U') or only the lower
00108 *>          (UPLO='L') will be referenced.
00109 *> \endverbatim
00110 *>
00111 *> \param[in] LDA
00112 *> \verbatim
00113 *>          LDA is INTEGER
00114 *>          The leading dimension of A.  It must be at least 1
00115 *>          and at least N.
00116 *> \endverbatim
00117 *>
00118 *> \param[in] D
00119 *> \verbatim
00120 *>          D is DOUBLE PRECISION array, dimension (N)
00121 *>          The diagonal of the (symmetric tri-) diagonal matrix.
00122 *> \endverbatim
00123 *>
00124 *> \param[in] E
00125 *> \verbatim
00126 *>          E is DOUBLE PRECISION array, dimension (N-1)
00127 *>          The off-diagonal of the (symmetric tri-) diagonal matrix.
00128 *>          E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and
00129 *>          (3,2) element, etc.
00130 *>          Not referenced if KBAND=0.
00131 *> \endverbatim
00132 *>
00133 *> \param[in] U
00134 *> \verbatim
00135 *>          U is COMPLEX*16 array, dimension (LDU, N)
00136 *>          If ITYPE=1 or 3, this contains the unitary matrix in
00137 *>          the decomposition, expressed as a dense matrix.  If ITYPE=2,
00138 *>          then it is not referenced.
00139 *> \endverbatim
00140 *>
00141 *> \param[in] LDU
00142 *> \verbatim
00143 *>          LDU is INTEGER
00144 *>          The leading dimension of U.  LDU must be at least N and
00145 *>          at least 1.
00146 *> \endverbatim
00147 *>
00148 *> \param[in] V
00149 *> \verbatim
00150 *>          V is COMPLEX*16 array, dimension (LDV, N)
00151 *>          If ITYPE=2 or 3, the columns of this array contain the
00152 *>          Householder vectors used to describe the unitary matrix
00153 *>          in the decomposition.  If UPLO='L', then the vectors are in
00154 *>          the lower triangle, if UPLO='U', then in the upper
00155 *>          triangle.
00156 *>          *NOTE* If ITYPE=2 or 3, V is modified and restored.  The
00157 *>          subdiagonal (if UPLO='L') or the superdiagonal (if UPLO='U')
00158 *>          is set to one, and later reset to its original value, during
00159 *>          the course of the calculation.
00160 *>          If ITYPE=1, then it is neither referenced nor modified.
00161 *> \endverbatim
00162 *>
00163 *> \param[in] LDV
00164 *> \verbatim
00165 *>          LDV is INTEGER
00166 *>          The leading dimension of V.  LDV must be at least N and
00167 *>          at least 1.
00168 *> \endverbatim
00169 *>
00170 *> \param[in] TAU
00171 *> \verbatim
00172 *>          TAU is COMPLEX*16 array, dimension (N)
00173 *>          If ITYPE >= 2, then TAU(j) is the scalar factor of
00174 *>          v(j) v(j)* in the Householder transformation H(j) of
00175 *>          the product  U = H(1)...H(n-2)
00176 *>          If ITYPE < 2, then TAU is not referenced.
00177 *> \endverbatim
00178 *>
00179 *> \param[out] WORK
00180 *> \verbatim
00181 *>          WORK is COMPLEX*16 array, dimension (2*N**2)
00182 *> \endverbatim
00183 *>
00184 *> \param[out] RWORK
00185 *> \verbatim
00186 *>          RWORK is DOUBLE PRECISION array, dimension (N)
00187 *> \endverbatim
00188 *>
00189 *> \param[out] RESULT
00190 *> \verbatim
00191 *>          RESULT is DOUBLE PRECISION array, dimension (2)
00192 *>          The values computed by the two tests described above.  The
00193 *>          values are currently limited to 1/ulp, to avoid overflow.
00194 *>          RESULT(1) is always modified.  RESULT(2) is modified only
00195 *>          if ITYPE=1.
00196 *> \endverbatim
00197 *
00198 *  Authors:
00199 *  ========
00200 *
00201 *> \author Univ. of Tennessee 
00202 *> \author Univ. of California Berkeley 
00203 *> \author Univ. of Colorado Denver 
00204 *> \author NAG Ltd. 
00205 *
00206 *> \date November 2011
00207 *
00208 *> \ingroup complex16_eig
00209 *
00210 *  =====================================================================
00211       SUBROUTINE ZHET21( ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V,
00212      $                   LDV, TAU, WORK, RWORK, RESULT )
00213 *
00214 *  -- LAPACK test routine (version 3.4.0) --
00215 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00216 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00217 *     November 2011
00218 *
00219 *     .. Scalar Arguments ..
00220       CHARACTER          UPLO
00221       INTEGER            ITYPE, KBAND, LDA, LDU, LDV, N
00222 *     ..
00223 *     .. Array Arguments ..
00224       DOUBLE PRECISION   D( * ), E( * ), RESULT( 2 ), RWORK( * )
00225       COMPLEX*16         A( LDA, * ), TAU( * ), U( LDU, * ),
00226      $                   V( LDV, * ), WORK( * )
00227 *     ..
00228 *
00229 *  =====================================================================
00230 *
00231 *     .. Parameters ..
00232       DOUBLE PRECISION   ZERO, ONE, TEN
00233       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 10.0D+0 )
00234       COMPLEX*16         CZERO, CONE
00235       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00236      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00237 *     ..
00238 *     .. Local Scalars ..
00239       LOGICAL            LOWER
00240       CHARACTER          CUPLO
00241       INTEGER            IINFO, J, JCOL, JR, JROW
00242       DOUBLE PRECISION   ANORM, ULP, UNFL, WNORM
00243       COMPLEX*16         VSAVE
00244 *     ..
00245 *     .. External Functions ..
00246       LOGICAL            LSAME
00247       DOUBLE PRECISION   DLAMCH, ZLANGE, ZLANHE
00248       EXTERNAL           LSAME, DLAMCH, ZLANGE, ZLANHE
00249 *     ..
00250 *     .. External Subroutines ..
00251       EXTERNAL           ZGEMM, ZHER, ZHER2, ZLACPY, ZLARFY, ZLASET,
00252      $                   ZUNM2L, ZUNM2R
00253 *     ..
00254 *     .. Intrinsic Functions ..
00255       INTRINSIC          DBLE, DCMPLX, MAX, MIN
00256 *     ..
00257 *     .. Executable Statements ..
00258 *
00259       RESULT( 1 ) = ZERO
00260       IF( ITYPE.EQ.1 )
00261      $   RESULT( 2 ) = ZERO
00262       IF( N.LE.0 )
00263      $   RETURN
00264 *
00265       IF( LSAME( UPLO, 'U' ) ) THEN
00266          LOWER = .FALSE.
00267          CUPLO = 'U'
00268       ELSE
00269          LOWER = .TRUE.
00270          CUPLO = 'L'
00271       END IF
00272 *
00273       UNFL = DLAMCH( 'Safe minimum' )
00274       ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
00275 *
00276 *     Some Error Checks
00277 *
00278       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
00279          RESULT( 1 ) = TEN / ULP
00280          RETURN
00281       END IF
00282 *
00283 *     Do Test 1
00284 *
00285 *     Norm of A:
00286 *
00287       IF( ITYPE.EQ.3 ) THEN
00288          ANORM = ONE
00289       ELSE
00290          ANORM = MAX( ZLANHE( '1', CUPLO, N, A, LDA, RWORK ), UNFL )
00291       END IF
00292 *
00293 *     Compute error matrix:
00294 *
00295       IF( ITYPE.EQ.1 ) THEN
00296 *
00297 *        ITYPE=1: error = A - U S U*
00298 *
00299          CALL ZLASET( 'Full', N, N, CZERO, CZERO, WORK, N )
00300          CALL ZLACPY( CUPLO, N, N, A, LDA, WORK, N )
00301 *
00302          DO 10 J = 1, N
00303             CALL ZHER( CUPLO, N, -D( J ), U( 1, J ), 1, WORK, N )
00304    10    CONTINUE
00305 *
00306          IF( N.GT.1 .AND. KBAND.EQ.1 ) THEN
00307             DO 20 J = 1, N - 1
00308                CALL ZHER2( CUPLO, N, -DCMPLX( E( J ) ), U( 1, J ), 1,
00309      $                     U( 1, J-1 ), 1, WORK, N )
00310    20       CONTINUE
00311          END IF
00312          WNORM = ZLANHE( '1', CUPLO, N, WORK, N, RWORK )
00313 *
00314       ELSE IF( ITYPE.EQ.2 ) THEN
00315 *
00316 *        ITYPE=2: error = V S V* - A
00317 *
00318          CALL ZLASET( 'Full', N, N, CZERO, CZERO, WORK, N )
00319 *
00320          IF( LOWER ) THEN
00321             WORK( N**2 ) = D( N )
00322             DO 40 J = N - 1, 1, -1
00323                IF( KBAND.EQ.1 ) THEN
00324                   WORK( ( N+1 )*( J-1 )+2 ) = ( CONE-TAU( J ) )*E( J )
00325                   DO 30 JR = J + 2, N
00326                      WORK( ( J-1 )*N+JR ) = -TAU( J )*E( J )*V( JR, J )
00327    30             CONTINUE
00328                END IF
00329 *
00330                VSAVE = V( J+1, J )
00331                V( J+1, J ) = ONE
00332                CALL ZLARFY( 'L', N-J, V( J+1, J ), 1, TAU( J ),
00333      $                      WORK( ( N+1 )*J+1 ), N, WORK( N**2+1 ) )
00334                V( J+1, J ) = VSAVE
00335                WORK( ( N+1 )*( J-1 )+1 ) = D( J )
00336    40       CONTINUE
00337          ELSE
00338             WORK( 1 ) = D( 1 )
00339             DO 60 J = 1, N - 1
00340                IF( KBAND.EQ.1 ) THEN
00341                   WORK( ( N+1 )*J ) = ( CONE-TAU( J ) )*E( J )
00342                   DO 50 JR = 1, J - 1
00343                      WORK( J*N+JR ) = -TAU( J )*E( J )*V( JR, J+1 )
00344    50             CONTINUE
00345                END IF
00346 *
00347                VSAVE = V( J, J+1 )
00348                V( J, J+1 ) = ONE
00349                CALL ZLARFY( 'U', J, V( 1, J+1 ), 1, TAU( J ), WORK, N,
00350      $                      WORK( N**2+1 ) )
00351                V( J, J+1 ) = VSAVE
00352                WORK( ( N+1 )*J+1 ) = D( J+1 )
00353    60       CONTINUE
00354          END IF
00355 *
00356          DO 90 JCOL = 1, N
00357             IF( LOWER ) THEN
00358                DO 70 JROW = JCOL, N
00359                   WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) )
00360      $                - A( JROW, JCOL )
00361    70          CONTINUE
00362             ELSE
00363                DO 80 JROW = 1, JCOL
00364                   WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) )
00365      $                - A( JROW, JCOL )
00366    80          CONTINUE
00367             END IF
00368    90    CONTINUE
00369          WNORM = ZLANHE( '1', CUPLO, N, WORK, N, RWORK )
00370 *
00371       ELSE IF( ITYPE.EQ.3 ) THEN
00372 *
00373 *        ITYPE=3: error = U V* - I
00374 *
00375          IF( N.LT.2 )
00376      $      RETURN
00377          CALL ZLACPY( ' ', N, N, U, LDU, WORK, N )
00378          IF( LOWER ) THEN
00379             CALL ZUNM2R( 'R', 'C', N, N-1, N-1, V( 2, 1 ), LDV, TAU,
00380      $                   WORK( N+1 ), N, WORK( N**2+1 ), IINFO )
00381          ELSE
00382             CALL ZUNM2L( 'R', 'C', N, N-1, N-1, V( 1, 2 ), LDV, TAU,
00383      $                   WORK, N, WORK( N**2+1 ), IINFO )
00384          END IF
00385          IF( IINFO.NE.0 ) THEN
00386             RESULT( 1 ) = TEN / ULP
00387             RETURN
00388          END IF
00389 *
00390          DO 100 J = 1, N
00391             WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - CONE
00392   100    CONTINUE
00393 *
00394          WNORM = ZLANGE( '1', N, N, WORK, N, RWORK )
00395       END IF
00396 *
00397       IF( ANORM.GT.WNORM ) THEN
00398          RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP )
00399       ELSE
00400          IF( ANORM.LT.ONE ) THEN
00401             RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
00402          ELSE
00403             RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP )
00404          END IF
00405       END IF
00406 *
00407 *     Do Test 2
00408 *
00409 *     Compute  UU* - I
00410 *
00411       IF( ITYPE.EQ.1 ) THEN
00412          CALL ZGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO,
00413      $               WORK, N )
00414 *
00415          DO 110 J = 1, N
00416             WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - CONE
00417   110    CONTINUE
00418 *
00419          RESULT( 2 ) = MIN( ZLANGE( '1', N, N, WORK, N, RWORK ),
00420      $                 DBLE( N ) ) / ( N*ULP )
00421       END IF
00422 *
00423       RETURN
00424 *
00425 *     End of ZHET21
00426 *
00427       END
 All Files Functions