LAPACK: TESTING/EIG/dsyt21.f Source File

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