LAPACK: TESTING/EIG/dsbt21.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DSBT21( UPLO, N, KA, KS, A, LDA, D, E, U, LDU, WORK,
00002      $                   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            KA, KS, LDA, LDU, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   A( LDA, * ), D( * ), E( * ), RESULT( 2 ),
00014      $                   U( LDU, * ), WORK( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  DSBT21  generally checks a decomposition of the form
00021 *
00022 *          A = U S U'
00023 *
00024 *  where ' means transpose, A is symmetric banded, U is
00025 *  orthogonal, and S is diagonal (if KS=0) or symmetric
00026 *  tridiagonal (if KS=1).
00027 *
00028 *  Specifically:
00029 *
00030 *          RESULT(1) = | A - U S U' | / ( |A| n ulp ) *and*
00031 *          RESULT(2) = | I - UU' | / ( n ulp )
00032 *
00033 *  Arguments
00034 *  =========
00035 *
00036 *  UPLO    (input) CHARACTER
00037 *          If UPLO='U', the upper triangle of A and V will be used and
00038 *          the (strictly) lower triangle will not be referenced.
00039 *          If UPLO='L', the lower triangle of A and V will be used and
00040 *          the (strictly) upper triangle will not be referenced.
00041 *
00042 *  N       (input) INTEGER
00043 *          The size of the matrix.  If it is zero, DSBT21 does nothing.
00044 *          It must be at least zero.
00045 *
00046 *  KA      (input) INTEGER
00047 *          The bandwidth of the matrix A.  It must be at least zero.  If
00048 *          it is larger than N-1, then max( 0, N-1 ) will be used.
00049 *
00050 *  KS      (input) INTEGER
00051 *          The bandwidth of the matrix S.  It may only be zero or one.
00052 *          If zero, then S is diagonal, and E is not referenced.  If
00053 *          one, then S is symmetric tri-diagonal.
00054 *
00055 *  A       (input) DOUBLE PRECISION array, dimension (LDA, N)
00056 *          The original (unfactored) matrix.  It is assumed to be
00057 *          symmetric, and only the upper (UPLO='U') or only the lower
00058 *          (UPLO='L') will be referenced.
00059 *
00060 *  LDA     (input) INTEGER
00061 *          The leading dimension of A.  It must be at least 1
00062 *          and at least min( KA, N-1 ).
00063 *
00064 *  D       (input) DOUBLE PRECISION array, dimension (N)
00065 *          The diagonal of the (symmetric tri-) diagonal matrix S.
00066 *
00067 *  E       (input) DOUBLE PRECISION array, dimension (N-1)
00068 *          The off-diagonal of the (symmetric tri-) diagonal matrix S.
00069 *          E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and
00070 *          (3,2) element, etc.
00071 *          Not referenced if KS=0.
00072 *
00073 *  U       (input) DOUBLE PRECISION array, dimension (LDU, N)
00074 *          The orthogonal matrix in the decomposition, expressed as a
00075 *          dense matrix (i.e., not as a product of Householder
00076 *          transformations, Givens transformations, etc.)
00077 *
00078 *  LDU     (input) INTEGER
00079 *          The leading dimension of U.  LDU must be at least N and
00080 *          at least 1.
00081 *
00082 *  WORK    (workspace) DOUBLE PRECISION array, dimension (N**2+N)
00083 *
00084 *  RESULT  (output) DOUBLE PRECISION array, dimension (2)
00085 *          The values computed by the two tests described above.  The
00086 *          values are currently limited to 1/ulp, to avoid overflow.
00087 *
00088 *  =====================================================================
00089 *
00090 *     .. Parameters ..
00091       DOUBLE PRECISION   ZERO, ONE
00092       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00093 *     ..
00094 *     .. Local Scalars ..
00095       LOGICAL            LOWER
00096       CHARACTER          CUPLO
00097       INTEGER            IKA, J, JC, JR, LW
00098       DOUBLE PRECISION   ANORM, ULP, UNFL, WNORM
00099 *     ..
00100 *     .. External Functions ..
00101       LOGICAL            LSAME
00102       DOUBLE PRECISION   DLAMCH, DLANGE, DLANSB, DLANSP
00103       EXTERNAL           LSAME, DLAMCH, DLANGE, DLANSB, DLANSP
00104 *     ..
00105 *     .. External Subroutines ..
00106       EXTERNAL           DGEMM, DSPR, DSPR2
00107 *     ..
00108 *     .. Intrinsic Functions ..
00109       INTRINSIC          DBLE, MAX, MIN
00110 *     ..
00111 *     .. Executable Statements ..
00112 *
00113 *     Constants
00114 *
00115       RESULT( 1 ) = ZERO
00116       RESULT( 2 ) = ZERO
00117       IF( N.LE.0 )
00118      $   RETURN
00119 *
00120       IKA = MAX( 0, MIN( N-1, KA ) )
00121       LW = ( N*( N+1 ) ) / 2
00122 *
00123       IF( LSAME( UPLO, 'U' ) ) THEN
00124          LOWER = .FALSE.
00125          CUPLO = 'U'
00126       ELSE
00127          LOWER = .TRUE.
00128          CUPLO = 'L'
00129       END IF
00130 *
00131       UNFL = DLAMCH( 'Safe minimum' )
00132       ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
00133 *
00134 *     Some Error Checks
00135 *
00136 *     Do Test 1
00137 *
00138 *     Norm of A:
00139 *
00140       ANORM = MAX( DLANSB( '1', CUPLO, N, IKA, A, LDA, WORK ), UNFL )
00141 *
00142 *     Compute error matrix:    Error = A - U S U'
00143 *
00144 *     Copy A from SB to SP storage format.
00145 *
00146       J = 0
00147       DO 50 JC = 1, N
00148          IF( LOWER ) THEN
00149             DO 10 JR = 1, MIN( IKA+1, N+1-JC )
00150                J = J + 1
00151                WORK( J ) = A( JR, JC )
00152    10       CONTINUE
00153             DO 20 JR = IKA + 2, N + 1 - JC
00154                J = J + 1
00155                WORK( J ) = ZERO
00156    20       CONTINUE
00157          ELSE
00158             DO 30 JR = IKA + 2, JC
00159                J = J + 1
00160                WORK( J ) = ZERO
00161    30       CONTINUE
00162             DO 40 JR = MIN( IKA, JC-1 ), 0, -1
00163                J = J + 1
00164                WORK( J ) = A( IKA+1-JR, JC )
00165    40       CONTINUE
00166          END IF
00167    50 CONTINUE
00168 *
00169       DO 60 J = 1, N
00170          CALL DSPR( CUPLO, N, -D( J ), U( 1, J ), 1, WORK )
00171    60 CONTINUE
00172 *
00173       IF( N.GT.1 .AND. KS.EQ.1 ) THEN
00174          DO 70 J = 1, N - 1
00175             CALL DSPR2( CUPLO, N, -E( J ), U( 1, J ), 1, U( 1, J+1 ), 1,
00176      $                  WORK )
00177    70    CONTINUE
00178       END IF
00179       WNORM = DLANSP( '1', CUPLO, N, WORK, WORK( LW+1 ) )
00180 *
00181       IF( ANORM.GT.WNORM ) THEN
00182          RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP )
00183       ELSE
00184          IF( ANORM.LT.ONE ) THEN
00185             RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
00186          ELSE
00187             RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP )
00188          END IF
00189       END IF
00190 *
00191 *     Do Test 2
00192 *
00193 *     Compute  UU' - I
00194 *
00195       CALL DGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK,
00196      $            N )
00197 *
00198       DO 80 J = 1, N
00199          WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - ONE
00200    80 CONTINUE
00201 *
00202       RESULT( 2 ) = MIN( DLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) ),
00203      $              DBLE( N ) ) / ( N*ULP )
00204 *
00205       RETURN
00206 *
00207 *     End of DSBT21
00208 *
00209       END