LAPACK 3.3.1
Linear Algebra PACKage

dsyt22.f

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