LAPACK: TESTING/EIG/dstt21.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DSTT21( N, KBAND, AD, AE, SD, SE, U, LDU, WORK,
00002      $                   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       INTEGER            KBAND, LDU, N
00010 *     ..
00011 *     .. Array Arguments ..
00012       DOUBLE PRECISION   AD( * ), AE( * ), RESULT( 2 ), SD( * ),
00013      $                   SE( * ), U( LDU, * ), WORK( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DSTT21 checks a decomposition of the form
00020 *
00021 *     A = U S U'
00022 *
00023 *  where ' means transpose, A is symmetric tridiagonal, U is orthogonal,
00024 *  and S is diagonal (if KBAND=0) or symmetric tridiagonal (if KBAND=1).
00025 *  Two tests are performed:
00026 *
00027 *     RESULT(1) = | A - U S U' | / ( |A| n ulp )
00028 *
00029 *     RESULT(2) = | I - UU' | / ( n ulp )
00030 *
00031 *  Arguments
00032 *  =========
00033 *
00034 *  N       (input) INTEGER
00035 *          The size of the matrix.  If it is zero, DSTT21 does nothing.
00036 *          It must be at least zero.
00037 *
00038 *  KBAND   (input) INTEGER
00039 *          The bandwidth of the matrix S.  It may only be zero or one.
00040 *          If zero, then S is diagonal, and SE is not referenced.  If
00041 *          one, then S is symmetric tri-diagonal.
00042 *
00043 *  AD      (input) DOUBLE PRECISION array, dimension (N)
00044 *          The diagonal of the original (unfactored) matrix A.  A is
00045 *          assumed to be symmetric tridiagonal.
00046 *
00047 *  AE      (input) DOUBLE PRECISION array, dimension (N-1)
00048 *          The off-diagonal of the original (unfactored) matrix A.  A
00049 *          is assumed to be symmetric tridiagonal.  AE(1) is the (1,2)
00050 *          and (2,1) element, AE(2) is the (2,3) and (3,2) element, etc.
00051 *
00052 *  SD      (input) DOUBLE PRECISION array, dimension (N)
00053 *          The diagonal of the (symmetric tri-) diagonal matrix S.
00054 *
00055 *  SE      (input) DOUBLE PRECISION array, dimension (N-1)
00056 *          The off-diagonal of the (symmetric tri-) diagonal matrix S.
00057 *          Not referenced if KBSND=0.  If KBAND=1, then AE(1) is the
00058 *          (1,2) and (2,1) element, SE(2) is the (2,3) and (3,2)
00059 *          element, etc.
00060 *
00061 *  U       (input) DOUBLE PRECISION array, dimension (LDU, N)
00062 *          The orthogonal matrix in the decomposition.
00063 *
00064 *  LDU     (input) INTEGER
00065 *          The leading dimension of U.  LDU must be at least N.
00066 *
00067 *  WORK    (workspace) DOUBLE PRECISION array, dimension (N*(N+1))
00068 *
00069 *  RESULT  (output) DOUBLE PRECISION array, dimension (2)
00070 *          The values computed by the two tests described above.  The
00071 *          values are currently limited to 1/ulp, to avoid overflow.
00072 *          RESULT(1) is always modified.
00073 *
00074 *  =====================================================================
00075 *
00076 *     .. Parameters ..
00077       DOUBLE PRECISION   ZERO, ONE
00078       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00079 *     ..
00080 *     .. Local Scalars ..
00081       INTEGER            J
00082       DOUBLE PRECISION   ANORM, TEMP1, TEMP2, ULP, UNFL, WNORM
00083 *     ..
00084 *     .. External Functions ..
00085       DOUBLE PRECISION   DLAMCH, DLANGE, DLANSY
00086       EXTERNAL           DLAMCH, DLANGE, DLANSY
00087 *     ..
00088 *     .. External Subroutines ..
00089       EXTERNAL           DGEMM, DLASET, DSYR, DSYR2
00090 *     ..
00091 *     .. Intrinsic Functions ..
00092       INTRINSIC          ABS, DBLE, MAX, MIN
00093 *     ..
00094 *     .. Executable Statements ..
00095 *
00096 *     1)      Constants
00097 *
00098       RESULT( 1 ) = ZERO
00099       RESULT( 2 ) = ZERO
00100       IF( N.LE.0 )
00101      $   RETURN
00102 *
00103       UNFL = DLAMCH( 'Safe minimum' )
00104       ULP = DLAMCH( 'Precision' )
00105 *
00106 *     Do Test 1
00107 *
00108 *     Copy A & Compute its 1-Norm:
00109 *
00110       CALL DLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
00111 *
00112       ANORM = ZERO
00113       TEMP1 = ZERO
00114 *
00115       DO 10 J = 1, N - 1
00116          WORK( ( N+1 )*( J-1 )+1 ) = AD( J )
00117          WORK( ( N+1 )*( J-1 )+2 ) = AE( J )
00118          TEMP2 = ABS( AE( J ) )
00119          ANORM = MAX( ANORM, ABS( AD( J ) )+TEMP1+TEMP2 )
00120          TEMP1 = TEMP2
00121    10 CONTINUE
00122 *
00123       WORK( N**2 ) = AD( N )
00124       ANORM = MAX( ANORM, ABS( AD( N ) )+TEMP1, UNFL )
00125 *
00126 *     Norm of A - USU'
00127 *
00128       DO 20 J = 1, N
00129          CALL DSYR( 'L', N, -SD( J ), U( 1, J ), 1, WORK, N )
00130    20 CONTINUE
00131 *
00132       IF( N.GT.1 .AND. KBAND.EQ.1 ) THEN
00133          DO 30 J = 1, N - 1
00134             CALL DSYR2( 'L', N, -SE( J ), U( 1, J ), 1, U( 1, J+1 ), 1,
00135      $                  WORK, N )
00136    30    CONTINUE
00137       END IF
00138 *
00139       WNORM = DLANSY( '1', 'L', N, WORK, N, WORK( N**2+1 ) )
00140 *
00141       IF( ANORM.GT.WNORM ) THEN
00142          RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP )
00143       ELSE
00144          IF( ANORM.LT.ONE ) THEN
00145             RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
00146          ELSE
00147             RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP )
00148          END IF
00149       END IF
00150 *
00151 *     Do Test 2
00152 *
00153 *     Compute  UU' - I
00154 *
00155       CALL DGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK,
00156      $            N )
00157 *
00158       DO 40 J = 1, N
00159          WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - ONE
00160    40 CONTINUE
00161 *
00162       RESULT( 2 ) = MIN( DBLE( N ), DLANGE( '1', N, N, WORK, N,
00163      $              WORK( N**2+1 ) ) ) / ( N*ULP )
00164 *
00165       RETURN
00166 *
00167 *     End of DSTT21
00168 *
00169       END