LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
ssyt22.f
Go to the documentation of this file.
00001 *> \brief \b SSYT22
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 SSYT22( ITYPE, UPLO, N, M, KBAND, A, LDA, D, E, U, LDU,
00012 *                          V, LDV, TAU, WORK, RESULT )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          UPLO
00016 *       INTEGER            ITYPE, KBAND, LDA, LDU, LDV, M, N
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       REAL               A( LDA, * ), D( * ), E( * ), RESULT( 2 ),
00020 *      $                   TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
00021 *       ..
00022 *  
00023 *
00024 *> \par Purpose:
00025 *  =============
00026 *>
00027 *> \verbatim
00028 *>
00029 *>      SSYT22  generally checks a decomposition of the form
00030 *>
00031 *>              A U = U S
00032 *>
00033 *>      where A is symmetric, the columns of U are orthonormal, and S
00034 *>      is diagonal (if KBAND=0) or symmetric tridiagonal (if
00035 *>      KBAND=1).  If ITYPE=1, then U is represented as a dense matrix,
00036 *>      otherwise the U is expressed as a product of Householder
00037 *>      transformations, whose vectors are stored in the array "V" and
00038 *>      whose scaling constants are in "TAU"; we shall use the letter
00039 *>      "V" to refer to the product of Householder transformations
00040 *>      (which should be equal to U).
00041 *>
00042 *>      Specifically, if ITYPE=1, then:
00043 *>
00044 *>              RESULT(1) = | U' A U - S | / ( |A| m ulp ) *andC>              RESULT(2) = | I - U'U | / ( m ulp )
00045 *> \endverbatim
00046 *
00047 *  Arguments:
00048 *  ==========
00049 *
00050 *> \verbatim
00051 *>  ITYPE   INTEGER
00052 *>          Specifies the type of tests to be performed.
00053 *>          1: U expressed as a dense orthogonal matrix:
00054 *>             RESULT(1) = | A - U S U' | / ( |A| n ulp )   *andC>             RESULT(2) = | I - UU' | / ( n ulp )
00055 *>
00056 *>  UPLO    CHARACTER
00057 *>          If UPLO='U', the upper triangle of A will be used and the
00058 *>          (strictly) lower triangle will not be referenced.  If
00059 *>          UPLO='L', the lower triangle of A will be used and the
00060 *>          (strictly) upper triangle will not be referenced.
00061 *>          Not modified.
00062 *>
00063 *>  N       INTEGER
00064 *>          The size of the matrix.  If it is zero, SSYT22 does nothing.
00065 *>          It must be at least zero.
00066 *>          Not modified.
00067 *>
00068 *>  M       INTEGER
00069 *>          The number of columns of U.  If it is zero, SSYT22 does
00070 *>          nothing.  It must be at least zero.
00071 *>          Not modified.
00072 *>
00073 *>  KBAND   INTEGER
00074 *>          The bandwidth of the matrix.  It may only be zero or one.
00075 *>          If zero, then S is diagonal, and E is not referenced.  If
00076 *>          one, then S is symmetric tri-diagonal.
00077 *>          Not modified.
00078 *>
00079 *>  A       REAL array, dimension (LDA , N)
00080 *>          The original (unfactored) matrix.  It is assumed to be
00081 *>          symmetric, and only the upper (UPLO='U') or only the lower
00082 *>          (UPLO='L') will be referenced.
00083 *>          Not modified.
00084 *>
00085 *>  LDA     INTEGER
00086 *>          The leading dimension of A.  It must be at least 1
00087 *>          and at least N.
00088 *>          Not modified.
00089 *>
00090 *>  D       REAL array, dimension (N)
00091 *>          The diagonal of the (symmetric tri-) diagonal matrix.
00092 *>          Not modified.
00093 *>
00094 *>  E       REAL array, dimension (N)
00095 *>          The off-diagonal of the (symmetric tri-) diagonal matrix.
00096 *>          E(1) is ignored, E(2) is the (1,2) and (2,1) element, etc.
00097 *>          Not referenced if KBAND=0.
00098 *>          Not modified.
00099 *>
00100 *>  U       REAL array, dimension (LDU, N)
00101 *>          If ITYPE=1 or 3, this contains the orthogonal matrix in
00102 *>          the decomposition, expressed as a dense matrix.  If ITYPE=2,
00103 *>          then it is not referenced.
00104 *>          Not modified.
00105 *>
00106 *>  LDU     INTEGER
00107 *>          The leading dimension of U.  LDU must be at least N and
00108 *>          at least 1.
00109 *>          Not modified.
00110 *>
00111 *>  V       REAL array, dimension (LDV, N)
00112 *>          If ITYPE=2 or 3, the lower triangle of this array contains
00113 *>          the Householder vectors used to describe the orthogonal
00114 *>          matrix in the decomposition.  If ITYPE=1, then it is not
00115 *>          referenced.
00116 *>          Not modified.
00117 *>
00118 *>  LDV     INTEGER
00119 *>          The leading dimension of V.  LDV must be at least N and
00120 *>          at least 1.
00121 *>          Not modified.
00122 *>
00123 *>  TAU     REAL array, dimension (N)
00124 *>          If ITYPE >= 2, then TAU(j) is the scalar factor of
00125 *>          v(j) v(j)' in the Householder transformation H(j) of
00126 *>          the product  U = H(1)...H(n-2)
00127 *>          If ITYPE < 2, then TAU is not referenced.
00128 *>          Not modified.
00129 *>
00130 *>  WORK    REAL array, dimension (2*N**2)
00131 *>          Workspace.
00132 *>          Modified.
00133 *>
00134 *>  RESULT  REAL array, dimension (2)
00135 *>          The values computed by the two tests described above.  The
00136 *>          values are currently limited to 1/ulp, to avoid overflow.
00137 *>          RESULT(1) is always modified.  RESULT(2) is modified only
00138 *>          if LDU is at least N.
00139 *>          Modified.
00140 *> \endverbatim
00141 *
00142 *  Authors:
00143 *  ========
00144 *
00145 *> \author Univ. of Tennessee 
00146 *> \author Univ. of California Berkeley 
00147 *> \author Univ. of Colorado Denver 
00148 *> \author NAG Ltd. 
00149 *
00150 *> \date November 2011
00151 *
00152 *> \ingroup single_eig
00153 *
00154 *  =====================================================================
00155       SUBROUTINE SSYT22( ITYPE, UPLO, N, M, KBAND, A, LDA, D, E, U, LDU,
00156      $                   V, LDV, TAU, WORK, RESULT )
00157 *
00158 *  -- LAPACK test routine (version 3.4.0) --
00159 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00160 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00161 *     November 2011
00162 *
00163 *     .. Scalar Arguments ..
00164       CHARACTER          UPLO
00165       INTEGER            ITYPE, KBAND, LDA, LDU, LDV, M, N
00166 *     ..
00167 *     .. Array Arguments ..
00168       REAL               A( LDA, * ), D( * ), E( * ), RESULT( 2 ),
00169      $                   TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
00170 *     ..
00171 *
00172 *  =====================================================================
00173 *
00174 *     .. Parameters ..
00175       REAL               ZERO, ONE
00176       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00177 *     ..
00178 *     .. Local Scalars ..
00179       INTEGER            J, JJ, JJ1, JJ2, NN, NNP1
00180       REAL               ANORM, ULP, UNFL, WNORM
00181 *     ..
00182 *     .. External Functions ..
00183       REAL               SLAMCH, SLANSY
00184       EXTERNAL           SLAMCH, SLANSY
00185 *     ..
00186 *     .. External Subroutines ..
00187       EXTERNAL           SGEMM, SSYMM
00188 *     ..
00189 *     .. Intrinsic Functions ..
00190       INTRINSIC          MAX, MIN, REAL
00191 *     ..
00192 *     .. Executable Statements ..
00193 *
00194       RESULT( 1 ) = ZERO
00195       RESULT( 2 ) = ZERO
00196       IF( N.LE.0 .OR. M.LE.0 )
00197      $   RETURN
00198 *
00199       UNFL = SLAMCH( 'Safe minimum' )
00200       ULP = SLAMCH( 'Precision' )
00201 *
00202 *     Do Test 1
00203 *
00204 *     Norm of A:
00205 *
00206       ANORM = MAX( SLANSY( '1', UPLO, N, A, LDA, WORK ), UNFL )
00207 *
00208 *     Compute error matrix:
00209 *
00210 *     ITYPE=1: error = U' A U - S
00211 *
00212       CALL SSYMM( 'L', UPLO, N, M, ONE, A, LDA, U, LDU, ZERO, WORK, N )
00213       NN = N*N
00214       NNP1 = NN + 1
00215       CALL SGEMM( 'T', 'N', M, M, N, ONE, U, LDU, WORK, N, ZERO,
00216      $            WORK( NNP1 ), N )
00217       DO 10 J = 1, M
00218          JJ = NN + ( J-1 )*N + J
00219          WORK( JJ ) = WORK( JJ ) - D( J )
00220    10 CONTINUE
00221       IF( KBAND.EQ.1 .AND. N.GT.1 ) THEN
00222          DO 20 J = 2, M
00223             JJ1 = NN + ( J-1 )*N + J - 1
00224             JJ2 = NN + ( J-2 )*N + J
00225             WORK( JJ1 ) = WORK( JJ1 ) - E( J-1 )
00226             WORK( JJ2 ) = WORK( JJ2 ) - E( J-1 )
00227    20    CONTINUE
00228       END IF
00229       WNORM = SLANSY( '1', UPLO, M, WORK( NNP1 ), N, WORK( 1 ) )
00230 *
00231       IF( ANORM.GT.WNORM ) THEN
00232          RESULT( 1 ) = ( WNORM / ANORM ) / ( M*ULP )
00233       ELSE
00234          IF( ANORM.LT.ONE ) THEN
00235             RESULT( 1 ) = ( MIN( WNORM, M*ANORM ) / ANORM ) / ( M*ULP )
00236          ELSE
00237             RESULT( 1 ) = MIN( WNORM / ANORM, REAL( M ) ) / ( M*ULP )
00238          END IF
00239       END IF
00240 *
00241 *     Do Test 2
00242 *
00243 *     Compute  U'U - I
00244 *
00245       IF( ITYPE.EQ.1 )
00246      $   CALL SORT01( 'Columns', N, M, U, LDU, WORK, 2*N*N,
00247      $                RESULT( 2 ) )
00248 *
00249       RETURN
00250 *
00251 *     End of SSYT22
00252 *
00253       END
 All Files Functions