LAPACK 3.3.0

chet22.f

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