LAPACK: TESTING/EIG/zstt21.f Source File

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