LAPACK: TESTING/EIG/chet21.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CHET21( ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V,
00002      $                   LDV, TAU, WORK, RWORK, RESULT )
00003 *
00004 *  -- LAPACK test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          UPLO
00010       INTEGER            ITYPE, KBAND, LDA, LDU, LDV, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       REAL               D( * ), E( * ), RESULT( 2 ), RWORK( * )
00014       COMPLEX            A( LDA, * ), TAU( * ), U( LDU, * ),
00015      $                   V( LDV, * ), WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  CHET21 generally checks a decomposition of the form
00022 *
00023 *     A = U S U*
00024 *
00025 *  where * means conjugate transpose, A is hermitian, U is unitary, and
00026 *  S is diagonal (if KBAND=0) or (real) symmetric tridiagonal (if
00027 *  KBAND=1).
00028 *
00029 *  If ITYPE=1, then U is represented as a dense matrix; otherwise U is
00030 *  expressed as a product of Householder transformations, whose vectors
00031 *  are stored in the array "V" and whose scaling constants are in "TAU".
00032 *  We shall use the letter "V" to refer to the product of Householder
00033 *  transformations (which should be equal to U).
00034 *
00035 *  Specifically, if ITYPE=1, then:
00036 *
00037 *     RESULT(1) = | A - U S U* | / ( |A| n ulp ) *and*
00038 *     RESULT(2) = | I - UU* | / ( n ulp )
00039 *
00040 *  If ITYPE=2, then:
00041 *
00042 *     RESULT(1) = | A - V S V* | / ( |A| n ulp )
00043 *
00044 *  If ITYPE=3, then:
00045 *
00046 *     RESULT(1) = | I - UV* | / ( n ulp )
00047 *
00048 *  For ITYPE > 1, the transformation U is expressed as a product
00049 *  V = H(1)...H(n-2),  where H(j) = I  -  tau(j) v(j) v(j)*  and each
00050 *  vector v(j) has its first j elements 0 and the remaining n-j elements
00051 *  stored in V(j+1:n,j).
00052 *
00053 *  Arguments
00054 *  =========
00055 *
00056 *  ITYPE   (input) INTEGER
00057 *          Specifies the type of tests to be performed.
00058 *          1: U expressed as a dense unitary matrix:
00059 *             RESULT(1) = | A - U S U* | / ( |A| n ulp )   *and*
00060 *             RESULT(2) = | I - UU* | / ( n ulp )
00061 *
00062 *          2: U expressed as a product V of Housholder transformations:
00063 *             RESULT(1) = | A - V S V* | / ( |A| n ulp )
00064 *
00065 *          3: U expressed both as a dense unitary matrix and
00066 *             as a product of Housholder transformations:
00067 *             RESULT(1) = | I - UV* | / ( n ulp )
00068 *
00069 *  UPLO    (input) CHARACTER
00070 *          If UPLO='U', the upper triangle of A and V will be used and
00071 *          the (strictly) lower triangle will not be referenced.
00072 *          If UPLO='L', the lower triangle of A and V will be used and
00073 *          the (strictly) upper triangle will not be referenced.
00074 *
00075 *  N       (input) INTEGER
00076 *          The size of the matrix.  If it is zero, CHET21 does nothing.
00077 *          It must be at least zero.
00078 *
00079 *  KBAND   (input) INTEGER
00080 *          The bandwidth of the matrix.  It may only be zero or one.
00081 *          If zero, then S is diagonal, and E is not referenced.  If
00082 *          one, then S is symmetric tri-diagonal.
00083 *
00084 *  A       (input) COMPLEX array, dimension (LDA, N)
00085 *          The original (unfactored) matrix.  It is assumed to be
00086 *          hermitian, and only the upper (UPLO='U') or only the lower
00087 *          (UPLO='L') will be referenced.
00088 *
00089 *  LDA     (input) INTEGER
00090 *          The leading dimension of A.  It must be at least 1
00091 *          and at least N.
00092 *
00093 *  D       (input) REAL array, dimension (N)
00094 *          The diagonal of the (symmetric tri-) diagonal matrix.
00095 *
00096 *  E       (input) REAL array, dimension (N-1)
00097 *          The off-diagonal of the (symmetric tri-) diagonal matrix.
00098 *          E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and
00099 *          (3,2) element, etc.
00100 *          Not referenced if KBAND=0.
00101 *
00102 *  U       (input) COMPLEX array, dimension (LDU, N)
00103 *          If ITYPE=1 or 3, this contains the unitary matrix in
00104 *          the decomposition, expressed as a dense matrix.  If ITYPE=2,
00105 *          then it is not referenced.
00106 *
00107 *  LDU     (input) INTEGER
00108 *          The leading dimension of U.  LDU must be at least N and
00109 *          at least 1.
00110 *
00111 *  V       (input) COMPLEX array, dimension (LDV, N)
00112 *          If ITYPE=2 or 3, the columns of this array contain the
00113 *          Householder vectors used to describe the unitary matrix
00114 *          in the decomposition.  If UPLO='L', then the vectors are in
00115 *          the lower triangle, if UPLO='U', then in the upper
00116 *          triangle.
00117 *          *NOTE* If ITYPE=2 or 3, V is modified and restored.  The
00118 *          subdiagonal (if UPLO='L') or the superdiagonal (if UPLO='U')
00119 *          is set to one, and later reset to its original value, during
00120 *          the course of the calculation.
00121 *          If ITYPE=1, then it is neither referenced nor modified.
00122 *
00123 *  LDV     (input) INTEGER
00124 *          The leading dimension of V.  LDV must be at least N and
00125 *          at least 1.
00126 *
00127 *  TAU     (input) COMPLEX array, dimension (N)
00128 *          If ITYPE >= 2, then TAU(j) is the scalar factor of
00129 *          v(j) v(j)* in the Householder transformation H(j) of
00130 *          the product  U = H(1)...H(n-2)
00131 *          If ITYPE < 2, then TAU is not referenced.
00132 *
00133 *  WORK    (workspace) COMPLEX array, dimension (2*N**2)
00134 *
00135 *  RWORK   (workspace) REAL array, dimension (N)
00136 *
00137 *  RESULT  (output) REAL array, dimension (2)
00138 *          The values computed by the two tests described above.  The
00139 *          values are currently limited to 1/ulp, to avoid overflow.
00140 *          RESULT(1) is always modified.  RESULT(2) is modified only
00141 *          if ITYPE=1.
00142 *
00143 *  =====================================================================
00144 *
00145 *     .. Parameters ..
00146       REAL               ZERO, ONE, TEN
00147       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 10.0E+0 )
00148       COMPLEX            CZERO, CONE
00149       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
00150      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00151 *     ..
00152 *     .. Local Scalars ..
00153       LOGICAL            LOWER
00154       CHARACTER          CUPLO
00155       INTEGER            IINFO, J, JCOL, JR, JROW
00156       REAL               ANORM, ULP, UNFL, WNORM
00157       COMPLEX            VSAVE
00158 *     ..
00159 *     .. External Functions ..
00160       LOGICAL            LSAME
00161       REAL               CLANGE, CLANHE, SLAMCH
00162       EXTERNAL           LSAME, CLANGE, CLANHE, SLAMCH
00163 *     ..
00164 *     .. External Subroutines ..
00165       EXTERNAL           CGEMM, CHER, CHER2, CLACPY, CLARFY, CLASET,
00166      $                   CUNM2L, CUNM2R
00167 *     ..
00168 *     .. Intrinsic Functions ..
00169       INTRINSIC          CMPLX, MAX, MIN, REAL
00170 *     ..
00171 *     .. Executable Statements ..
00172 *
00173       RESULT( 1 ) = ZERO
00174       IF( ITYPE.EQ.1 )
00175      $   RESULT( 2 ) = ZERO
00176       IF( N.LE.0 )
00177      $   RETURN
00178 *
00179       IF( LSAME( UPLO, 'U' ) ) THEN
00180          LOWER = .FALSE.
00181          CUPLO = 'U'
00182       ELSE
00183          LOWER = .TRUE.
00184          CUPLO = 'L'
00185       END IF
00186 *
00187       UNFL = SLAMCH( 'Safe minimum' )
00188       ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
00189 *
00190 *     Some Error Checks
00191 *
00192       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
00193          RESULT( 1 ) = TEN / ULP
00194          RETURN
00195       END IF
00196 *
00197 *     Do Test 1
00198 *
00199 *     Norm of A:
00200 *
00201       IF( ITYPE.EQ.3 ) THEN
00202          ANORM = ONE
00203       ELSE
00204          ANORM = MAX( CLANHE( '1', CUPLO, N, A, LDA, RWORK ), UNFL )
00205       END IF
00206 *
00207 *     Compute error matrix:
00208 *
00209       IF( ITYPE.EQ.1 ) THEN
00210 *
00211 *        ITYPE=1: error = A - U S U*
00212 *
00213          CALL CLASET( 'Full', N, N, CZERO, CZERO, WORK, N )
00214          CALL CLACPY( CUPLO, N, N, A, LDA, WORK, N )
00215 *
00216          DO 10 J = 1, N
00217             CALL CHER( CUPLO, N, -D( J ), U( 1, J ), 1, WORK, N )
00218    10    CONTINUE
00219 *
00220          IF( N.GT.1 .AND. KBAND.EQ.1 ) THEN
00221             DO 20 J = 1, N - 1
00222                CALL CHER2( CUPLO, N, -CMPLX( E( J ) ), U( 1, J ), 1,
00223      $                     U( 1, J-1 ), 1, WORK, N )
00224    20       CONTINUE
00225          END IF
00226          WNORM = CLANHE( '1', CUPLO, N, WORK, N, RWORK )
00227 *
00228       ELSE IF( ITYPE.EQ.2 ) THEN
00229 *
00230 *        ITYPE=2: error = V S V* - A
00231 *
00232          CALL CLASET( 'Full', N, N, CZERO, CZERO, WORK, N )
00233 *
00234          IF( LOWER ) THEN
00235             WORK( N**2 ) = D( N )
00236             DO 40 J = N - 1, 1, -1
00237                IF( KBAND.EQ.1 ) THEN
00238                   WORK( ( N+1 )*( J-1 )+2 ) = ( CONE-TAU( J ) )*E( J )
00239                   DO 30 JR = J + 2, N
00240                      WORK( ( J-1 )*N+JR ) = -TAU( J )*E( J )*V( JR, J )
00241    30             CONTINUE
00242                END IF
00243 *
00244                VSAVE = V( J+1, J )
00245                V( J+1, J ) = ONE
00246                CALL CLARFY( 'L', N-J, V( J+1, J ), 1, TAU( J ),
00247      $                      WORK( ( N+1 )*J+1 ), N, WORK( N**2+1 ) )
00248                V( J+1, J ) = VSAVE
00249                WORK( ( N+1 )*( J-1 )+1 ) = D( J )
00250    40       CONTINUE
00251          ELSE
00252             WORK( 1 ) = D( 1 )
00253             DO 60 J = 1, N - 1
00254                IF( KBAND.EQ.1 ) THEN
00255                   WORK( ( N+1 )*J ) = ( CONE-TAU( J ) )*E( J )
00256                   DO 50 JR = 1, J - 1
00257                      WORK( J*N+JR ) = -TAU( J )*E( J )*V( JR, J+1 )
00258    50             CONTINUE
00259                END IF
00260 *
00261                VSAVE = V( J, J+1 )
00262                V( J, J+1 ) = ONE
00263                CALL CLARFY( 'U', J, V( 1, J+1 ), 1, TAU( J ), WORK, N,
00264      $                      WORK( N**2+1 ) )
00265                V( J, J+1 ) = VSAVE
00266                WORK( ( N+1 )*J+1 ) = D( J+1 )
00267    60       CONTINUE
00268          END IF
00269 *
00270          DO 90 JCOL = 1, N
00271             IF( LOWER ) THEN
00272                DO 70 JROW = JCOL, N
00273                   WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) )
00274      $                - A( JROW, JCOL )
00275    70          CONTINUE
00276             ELSE
00277                DO 80 JROW = 1, JCOL
00278                   WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) )
00279      $                - A( JROW, JCOL )
00280    80          CONTINUE
00281             END IF
00282    90    CONTINUE
00283          WNORM = CLANHE( '1', CUPLO, N, WORK, N, RWORK )
00284 *
00285       ELSE IF( ITYPE.EQ.3 ) THEN
00286 *
00287 *        ITYPE=3: error = U V* - I
00288 *
00289          IF( N.LT.2 )
00290      $      RETURN
00291          CALL CLACPY( ' ', N, N, U, LDU, WORK, N )
00292          IF( LOWER ) THEN
00293             CALL CUNM2R( 'R', 'C', N, N-1, N-1, V( 2, 1 ), LDV, TAU,
00294      $                   WORK( N+1 ), N, WORK( N**2+1 ), IINFO )
00295          ELSE
00296             CALL CUNM2L( 'R', 'C', N, N-1, N-1, V( 1, 2 ), LDV, TAU,
00297      $                   WORK, N, WORK( N**2+1 ), IINFO )
00298          END IF
00299          IF( IINFO.NE.0 ) THEN
00300             RESULT( 1 ) = TEN / ULP
00301             RETURN
00302          END IF
00303 *
00304          DO 100 J = 1, N
00305             WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - CONE
00306   100    CONTINUE
00307 *
00308          WNORM = CLANGE( '1', N, N, WORK, N, RWORK )
00309       END IF
00310 *
00311       IF( ANORM.GT.WNORM ) THEN
00312          RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP )
00313       ELSE
00314          IF( ANORM.LT.ONE ) THEN
00315             RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
00316          ELSE
00317             RESULT( 1 ) = MIN( WNORM / ANORM, REAL( N ) ) / ( N*ULP )
00318          END IF
00319       END IF
00320 *
00321 *     Do Test 2
00322 *
00323 *     Compute  UU* - I
00324 *
00325       IF( ITYPE.EQ.1 ) THEN
00326          CALL CGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO,
00327      $               WORK, N )
00328 *
00329          DO 110 J = 1, N
00330             WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - CONE
00331   110    CONTINUE
00332 *
00333          RESULT( 2 ) = MIN( CLANGE( '1', N, N, WORK, N, RWORK ),
00334      $                 REAL( N ) ) / ( N*ULP )
00335       END IF
00336 *
00337       RETURN
00338 *
00339 *     End of CHET21
00340 *
00341       END