LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cstt21.f
Go to the documentation of this file.
00001 *> \brief \b CSTT21
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 CSTT21( N, KBAND, AD, AE, SD, SE, U, LDU, WORK, RWORK,
00012 *                          RESULT )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       INTEGER            KBAND, LDU, N
00016 *       ..
00017 *       .. Array Arguments ..
00018 *       REAL               AD( * ), AE( * ), RESULT( 2 ), RWORK( * ),
00019 *      $                   SD( * ), SE( * )
00020 *       COMPLEX            U( LDU, * ), WORK( * )
00021 *       ..
00022 *  
00023 *
00024 *> \par Purpose:
00025 *  =============
00026 *>
00027 *> \verbatim
00028 *>
00029 *> CSTT21  checks a decomposition of the form
00030 *>
00031 *>    A = U S UC>
00032 *> where * means conjugate transpose, A is real symmetric tridiagonal,
00033 *> U is unitary, and S is real and diagonal (if KBAND=0) or symmetric
00034 *> tridiagonal (if KBAND=1).  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, CSTT21 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 real 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 real (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 COMPLEX array, dimension (LDU, N)
00092 *>          The unitary 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 COMPLEX array, dimension (N**2)
00104 *> \endverbatim
00105 *>
00106 *> \param[out] RWORK
00107 *> \verbatim
00108 *>          RWORK is REAL array, dimension (N)
00109 *> \endverbatim
00110 *>
00111 *> \param[out] RESULT
00112 *> \verbatim
00113 *>          RESULT is REAL array, dimension (2)
00114 *>          The values computed by the two tests described above.  The
00115 *>          values are currently limited to 1/ulp, to avoid overflow.
00116 *>          RESULT(1) is always modified.
00117 *> \endverbatim
00118 *
00119 *  Authors:
00120 *  ========
00121 *
00122 *> \author Univ. of Tennessee 
00123 *> \author Univ. of California Berkeley 
00124 *> \author Univ. of Colorado Denver 
00125 *> \author NAG Ltd. 
00126 *
00127 *> \date November 2011
00128 *
00129 *> \ingroup complex_eig
00130 *
00131 *  =====================================================================
00132       SUBROUTINE CSTT21( N, KBAND, AD, AE, SD, SE, U, LDU, WORK, RWORK,
00133      $                   RESULT )
00134 *
00135 *  -- LAPACK test routine (version 3.4.0) --
00136 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00137 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00138 *     November 2011
00139 *
00140 *     .. Scalar Arguments ..
00141       INTEGER            KBAND, LDU, N
00142 *     ..
00143 *     .. Array Arguments ..
00144       REAL               AD( * ), AE( * ), RESULT( 2 ), RWORK( * ),
00145      $                   SD( * ), SE( * )
00146       COMPLEX            U( LDU, * ), WORK( * )
00147 *     ..
00148 *
00149 *  =====================================================================
00150 *
00151 *     .. Parameters ..
00152       REAL               ZERO, ONE
00153       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00154       COMPLEX            CZERO, CONE
00155       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
00156      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00157 *     ..
00158 *     .. Local Scalars ..
00159       INTEGER            J
00160       REAL               ANORM, TEMP1, TEMP2, ULP, UNFL, WNORM
00161 *     ..
00162 *     .. External Functions ..
00163       REAL               CLANGE, CLANHE, SLAMCH
00164       EXTERNAL           CLANGE, CLANHE, SLAMCH
00165 *     ..
00166 *     .. External Subroutines ..
00167       EXTERNAL           CGEMM, CHER, CHER2, CLASET
00168 *     ..
00169 *     .. Intrinsic Functions ..
00170       INTRINSIC          ABS, CMPLX, MAX, MIN, REAL
00171 *     ..
00172 *     .. Executable Statements ..
00173 *
00174 *     1)      Constants
00175 *
00176       RESULT( 1 ) = ZERO
00177       RESULT( 2 ) = ZERO
00178       IF( N.LE.0 )
00179      $   RETURN
00180 *
00181       UNFL = SLAMCH( 'Safe minimum' )
00182       ULP = SLAMCH( 'Precision' )
00183 *
00184 *     Do Test 1
00185 *
00186 *     Copy A & Compute its 1-Norm:
00187 *
00188       CALL CLASET( 'Full', N, N, CZERO, CZERO, WORK, N )
00189 *
00190       ANORM = ZERO
00191       TEMP1 = ZERO
00192 *
00193       DO 10 J = 1, N - 1
00194          WORK( ( N+1 )*( J-1 )+1 ) = AD( J )
00195          WORK( ( N+1 )*( J-1 )+2 ) = AE( J )
00196          TEMP2 = ABS( AE( J ) )
00197          ANORM = MAX( ANORM, ABS( AD( J ) )+TEMP1+TEMP2 )
00198          TEMP1 = TEMP2
00199    10 CONTINUE
00200 *
00201       WORK( N**2 ) = AD( N )
00202       ANORM = MAX( ANORM, ABS( AD( N ) )+TEMP1, UNFL )
00203 *
00204 *     Norm of A - USU*
00205 *
00206       DO 20 J = 1, N
00207          CALL CHER( 'L', N, -SD( J ), U( 1, J ), 1, WORK, N )
00208    20 CONTINUE
00209 *
00210       IF( N.GT.1 .AND. KBAND.EQ.1 ) THEN
00211          DO 30 J = 1, N - 1
00212             CALL CHER2( 'L', N, -CMPLX( SE( J ) ), U( 1, J ), 1,
00213      $                  U( 1, J+1 ), 1, WORK, N )
00214    30    CONTINUE
00215       END IF
00216 *
00217       WNORM = CLANHE( '1', 'L', N, WORK, N, RWORK )
00218 *
00219       IF( ANORM.GT.WNORM ) THEN
00220          RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP )
00221       ELSE
00222          IF( ANORM.LT.ONE ) THEN
00223             RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
00224          ELSE
00225             RESULT( 1 ) = MIN( WNORM / ANORM, REAL( N ) ) / ( N*ULP )
00226          END IF
00227       END IF
00228 *
00229 *     Do Test 2
00230 *
00231 *     Compute  UU* - I
00232 *
00233       CALL CGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO, WORK,
00234      $            N )
00235 *
00236       DO 40 J = 1, N
00237          WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - CONE
00238    40 CONTINUE
00239 *
00240       RESULT( 2 ) = MIN( REAL( N ), CLANGE( '1', N, N, WORK, N,
00241      $              RWORK ) ) / ( N*ULP )
00242 *
00243       RETURN
00244 *
00245 *     End of CSTT21
00246 *
00247       END
 All Files Functions