LAPACK: TESTING/EIG/zhpt21.f Source File

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