LAPACK: TESTING/EIG/dspt21.f Source File

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