LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sstt21.f
Go to the documentation of this file.
00001 *> \brief \b SSTT21
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE SSTT21( N, KBAND, AD, AE, SD, SE, U, LDU, WORK,
00012 *                          RESULT )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       INTEGER            KBAND, LDU, N
00016 *       ..
00017 *       .. Array Arguments ..
00018 *       REAL               AD( * ), AE( * ), RESULT( 2 ), SD( * ),
00019 *      $                   SE( * ), U( LDU, * ), WORK( * )
00020 *       ..
00021 *  
00022 *
00023 *> \par Purpose:
00024 *  =============
00025 *>
00026 *> \verbatim
00027 *>
00028 *> SSTT21 checks a decomposition of the form
00029 *>
00030 *>    A = U S U'
00031 *>
00032 *> where ' means transpose, A is symmetric tridiagonal, U is orthogonal,
00033 *> and S is diagonal (if KBAND=0) or symmetric tridiagonal (if KBAND=1).
00034 *> Two tests are performed:
00035 *>
00036 *>    RESULT(1) = | A - U S U' | / ( |A| n ulp )
00037 *>
00038 *>    RESULT(2) = | I - UU' | / ( n ulp )
00039 *> \endverbatim
00040 *
00041 *  Arguments:
00042 *  ==========
00043 *
00044 *> \param[in] N
00045 *> \verbatim
00046 *>          N is INTEGER
00047 *>          The size of the matrix.  If it is zero, SSTT21 does nothing.
00048 *>          It must be at least zero.
00049 *> \endverbatim
00050 *>
00051 *> \param[in] KBAND
00052 *> \verbatim
00053 *>          KBAND is INTEGER
00054 *>          The bandwidth of the matrix S.  It may only be zero or one.
00055 *>          If zero, then S is diagonal, and SE is not referenced.  If
00056 *>          one, then S is symmetric tri-diagonal.
00057 *> \endverbatim
00058 *>
00059 *> \param[in] AD
00060 *> \verbatim
00061 *>          AD is REAL array, dimension (N)
00062 *>          The diagonal of the original (unfactored) matrix A.  A is
00063 *>          assumed to be symmetric tridiagonal.
00064 *> \endverbatim
00065 *>
00066 *> \param[in] AE
00067 *> \verbatim
00068 *>          AE is REAL array, dimension (N-1)
00069 *>          The off-diagonal of the original (unfactored) matrix A.  A
00070 *>          is assumed to be symmetric tridiagonal.  AE(1) is the (1,2)
00071 *>          and (2,1) element, AE(2) is the (2,3) and (3,2) element, etc.
00072 *> \endverbatim
00073 *>
00074 *> \param[in] SD
00075 *> \verbatim
00076 *>          SD is REAL array, dimension (N)
00077 *>          The diagonal of the (symmetric tri-) diagonal matrix S.
00078 *> \endverbatim
00079 *>
00080 *> \param[in] SE
00081 *> \verbatim
00082 *>          SE is REAL array, dimension (N-1)
00083 *>          The off-diagonal of the (symmetric tri-) diagonal matrix S.
00084 *>          Not referenced if KBSND=0.  If KBAND=1, then AE(1) is the
00085 *>          (1,2) and (2,1) element, SE(2) is the (2,3) and (3,2)
00086 *>          element, etc.
00087 *> \endverbatim
00088 *>
00089 *> \param[in] U
00090 *> \verbatim
00091 *>          U is REAL array, dimension (LDU, N)
00092 *>          The orthogonal matrix in the decomposition.
00093 *> \endverbatim
00094 *>
00095 *> \param[in] LDU
00096 *> \verbatim
00097 *>          LDU is INTEGER
00098 *>          The leading dimension of U.  LDU must be at least N.
00099 *> \endverbatim
00100 *>
00101 *> \param[out] WORK
00102 *> \verbatim
00103 *>          WORK is REAL array, dimension (N*(N+1))
00104 *> \endverbatim
00105 *>
00106 *> \param[out] RESULT
00107 *> \verbatim
00108 *>          RESULT is REAL array, dimension (2)
00109 *>          The values computed by the two tests described above.  The
00110 *>          values are currently limited to 1/ulp, to avoid overflow.
00111 *>          RESULT(1) is always modified.
00112 *> \endverbatim
00113 *
00114 *  Authors:
00115 *  ========
00116 *
00117 *> \author Univ. of Tennessee 
00118 *> \author Univ. of California Berkeley 
00119 *> \author Univ. of Colorado Denver 
00120 *> \author NAG Ltd. 
00121 *
00122 *> \date November 2011
00123 *
00124 *> \ingroup single_eig
00125 *
00126 *  =====================================================================
00127       SUBROUTINE SSTT21( N, KBAND, AD, AE, SD, SE, U, LDU, WORK,
00128      $                   RESULT )
00129 *
00130 *  -- LAPACK test routine (version 3.4.0) --
00131 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00132 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00133 *     November 2011
00134 *
00135 *     .. Scalar Arguments ..
00136       INTEGER            KBAND, LDU, N
00137 *     ..
00138 *     .. Array Arguments ..
00139       REAL               AD( * ), AE( * ), RESULT( 2 ), SD( * ),
00140      $                   SE( * ), U( LDU, * ), WORK( * )
00141 *     ..
00142 *
00143 *  =====================================================================
00144 *
00145 *     .. Parameters ..
00146       REAL               ZERO, ONE
00147       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00148 *     ..
00149 *     .. Local Scalars ..
00150       INTEGER            J
00151       REAL               ANORM, TEMP1, TEMP2, ULP, UNFL, WNORM
00152 *     ..
00153 *     .. External Functions ..
00154       REAL               SLAMCH, SLANGE, SLANSY
00155       EXTERNAL           SLAMCH, SLANGE, SLANSY
00156 *     ..
00157 *     .. External Subroutines ..
00158       EXTERNAL           SGEMM, SLASET, SSYR, SSYR2
00159 *     ..
00160 *     .. Intrinsic Functions ..
00161       INTRINSIC          ABS, MAX, MIN, REAL
00162 *     ..
00163 *     .. Executable Statements ..
00164 *
00165 *     1)      Constants
00166 *
00167       RESULT( 1 ) = ZERO
00168       RESULT( 2 ) = ZERO
00169       IF( N.LE.0 )
00170      $   RETURN
00171 *
00172       UNFL = SLAMCH( 'Safe minimum' )
00173       ULP = SLAMCH( 'Precision' )
00174 *
00175 *     Do Test 1
00176 *
00177 *     Copy A & Compute its 1-Norm:
00178 *
00179       CALL SLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
00180 *
00181       ANORM = ZERO
00182       TEMP1 = ZERO
00183 *
00184       DO 10 J = 1, N - 1
00185          WORK( ( N+1 )*( J-1 )+1 ) = AD( J )
00186          WORK( ( N+1 )*( J-1 )+2 ) = AE( J )
00187          TEMP2 = ABS( AE( J ) )
00188          ANORM = MAX( ANORM, ABS( AD( J ) )+TEMP1+TEMP2 )
00189          TEMP1 = TEMP2
00190    10 CONTINUE
00191 *
00192       WORK( N**2 ) = AD( N )
00193       ANORM = MAX( ANORM, ABS( AD( N ) )+TEMP1, UNFL )
00194 *
00195 *     Norm of A - USU'
00196 *
00197       DO 20 J = 1, N
00198          CALL SSYR( 'L', N, -SD( J ), U( 1, J ), 1, WORK, N )
00199    20 CONTINUE
00200 *
00201       IF( N.GT.1 .AND. KBAND.EQ.1 ) THEN
00202          DO 30 J = 1, N - 1
00203             CALL SSYR2( 'L', N, -SE( J ), U( 1, J ), 1, U( 1, J+1 ), 1,
00204      $                  WORK, N )
00205    30    CONTINUE
00206       END IF
00207 *
00208       WNORM = SLANSY( '1', 'L', N, WORK, N, WORK( N**2+1 ) )
00209 *
00210       IF( ANORM.GT.WNORM ) THEN
00211          RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP )
00212       ELSE
00213          IF( ANORM.LT.ONE ) THEN
00214             RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
00215          ELSE
00216             RESULT( 1 ) = MIN( WNORM / ANORM, REAL( N ) ) / ( N*ULP )
00217          END IF
00218       END IF
00219 *
00220 *     Do Test 2
00221 *
00222 *     Compute  UU' - I
00223 *
00224       CALL SGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK,
00225      $            N )
00226 *
00227       DO 40 J = 1, N
00228          WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - ONE
00229    40 CONTINUE
00230 *
00231       RESULT( 2 ) = MIN( REAL( N ), SLANGE( '1', N, N, WORK, N,
00232      $              WORK( N**2+1 ) ) ) / ( N*ULP )
00233 *
00234       RETURN
00235 *
00236 *     End of SSTT21
00237 *
00238       END
 All Files Functions