LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ dstt21()

subroutine dstt21 ( integer  N,
integer  KBAND,
double precision, dimension( * )  AD,
double precision, dimension( * )  AE,
double precision, dimension( * )  SD,
double precision, dimension( * )  SE,
double precision, dimension( ldu, * )  U,
integer  LDU,
double precision, dimension( * )  WORK,
double precision, dimension( 2 )  RESULT 
)

DSTT21

Purpose:
 DSTT21 checks a decomposition of the form

    A = U S U'

 where ' means transpose, A is symmetric tridiagonal, U is orthogonal,
 and S is diagonal (if KBAND=0) or symmetric tridiagonal (if KBAND=1).
 Two tests are performed:

    RESULT(1) = | A - U S U' | / ( |A| n ulp )

    RESULT(2) = | I - UU' | / ( n ulp )
Parameters
[in]N
          N is INTEGER
          The size of the matrix.  If it is zero, DSTT21 does nothing.
          It must be at least zero.
[in]KBAND
          KBAND is INTEGER
          The bandwidth of the matrix S.  It may only be zero or one.
          If zero, then S is diagonal, and SE is not referenced.  If
          one, then S is symmetric tri-diagonal.
[in]AD
          AD is DOUBLE PRECISION array, dimension (N)
          The diagonal of the original (unfactored) matrix A.  A is
          assumed to be symmetric tridiagonal.
[in]AE
          AE is DOUBLE PRECISION array, dimension (N-1)
          The off-diagonal of the original (unfactored) matrix A.  A
          is assumed to be symmetric tridiagonal.  AE(1) is the (1,2)
          and (2,1) element, AE(2) is the (2,3) and (3,2) element, etc.
[in]SD
          SD is DOUBLE PRECISION array, dimension (N)
          The diagonal of the (symmetric tri-) diagonal matrix S.
[in]SE
          SE is DOUBLE PRECISION array, dimension (N-1)
          The off-diagonal of the (symmetric tri-) diagonal matrix S.
          Not referenced if KBSND=0.  If KBAND=1, then AE(1) is the
          (1,2) and (2,1) element, SE(2) is the (2,3) and (3,2)
          element, etc.
[in]U
          U is DOUBLE PRECISION array, dimension (LDU, N)
          The orthogonal matrix in the decomposition.
[in]LDU
          LDU is INTEGER
          The leading dimension of U.  LDU must be at least N.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (N*(N+1))
[out]RESULT
          RESULT is DOUBLE PRECISION array, dimension (2)
          The values computed by the two tests described above.  The
          values are currently limited to 1/ulp, to avoid overflow.
          RESULT(1) is always modified.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 129 of file dstt21.f.

129 *
130 * -- LAPACK test routine (version 3.7.0) --
131 * -- LAPACK is a software package provided by Univ. of Tennessee, --
132 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
133 * December 2016
134 *
135 * .. Scalar Arguments ..
136  INTEGER kband, ldu, n
137 * ..
138 * .. Array Arguments ..
139  DOUBLE PRECISION ad( * ), ae( * ), result( 2 ), sd( * ),
140  $ se( * ), u( ldu, * ), work( * )
141 * ..
142 *
143 * =====================================================================
144 *
145 * .. Parameters ..
146  DOUBLE PRECISION zero, one
147  parameter( zero = 0.0d0, one = 1.0d0 )
148 * ..
149 * .. Local Scalars ..
150  INTEGER j
151  DOUBLE PRECISION anorm, temp1, temp2, ulp, unfl, wnorm
152 * ..
153 * .. External Functions ..
154  DOUBLE PRECISION dlamch, dlange, dlansy
155  EXTERNAL dlamch, dlange, dlansy
156 * ..
157 * .. External Subroutines ..
158  EXTERNAL dgemm, dlaset, dsyr, dsyr2
159 * ..
160 * .. Intrinsic Functions ..
161  INTRINSIC abs, dble, max, min
162 * ..
163 * .. Executable Statements ..
164 *
165 * 1) Constants
166 *
167  result( 1 ) = zero
168  result( 2 ) = zero
169  IF( n.LE.0 )
170  $ RETURN
171 *
172  unfl = dlamch( 'Safe minimum' )
173  ulp = dlamch( 'Precision' )
174 *
175 * Do Test 1
176 *
177 * Copy A & Compute its 1-Norm:
178 *
179  CALL dlaset( 'Full', n, n, zero, zero, work, n )
180 *
181  anorm = zero
182  temp1 = zero
183 *
184  DO 10 j = 1, n - 1
185  work( ( n+1 )*( j-1 )+1 ) = ad( j )
186  work( ( n+1 )*( j-1 )+2 ) = ae( j )
187  temp2 = abs( ae( j ) )
188  anorm = max( anorm, abs( ad( j ) )+temp1+temp2 )
189  temp1 = temp2
190  10 CONTINUE
191 *
192  work( n**2 ) = ad( n )
193  anorm = max( anorm, abs( ad( n ) )+temp1, unfl )
194 *
195 * Norm of A - USU'
196 *
197  DO 20 j = 1, n
198  CALL dsyr( 'L', n, -sd( j ), u( 1, j ), 1, work, n )
199  20 CONTINUE
200 *
201  IF( n.GT.1 .AND. kband.EQ.1 ) THEN
202  DO 30 j = 1, n - 1
203  CALL dsyr2( 'L', n, -se( j ), u( 1, j ), 1, u( 1, j+1 ), 1,
204  $ work, n )
205  30 CONTINUE
206  END IF
207 *
208  wnorm = dlansy( '1', 'L', n, work, n, work( n**2+1 ) )
209 *
210  IF( anorm.GT.wnorm ) THEN
211  result( 1 ) = ( wnorm / anorm ) / ( n*ulp )
212  ELSE
213  IF( anorm.LT.one ) THEN
214  result( 1 ) = ( min( wnorm, n*anorm ) / anorm ) / ( n*ulp )
215  ELSE
216  result( 1 ) = min( wnorm / anorm, dble( n ) ) / ( n*ulp )
217  END IF
218  END IF
219 *
220 * Do Test 2
221 *
222 * Compute UU' - I
223 *
224  CALL dgemm( 'N', 'C', n, n, n, one, u, ldu, u, ldu, zero, work,
225  $ n )
226 *
227  DO 40 j = 1, n
228  work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - one
229  40 CONTINUE
230 *
231  result( 2 ) = min( dble( n ), dlange( '1', n, n, work, n,
232  $ work( n**2+1 ) ) ) / ( n*ulp )
233 *
234  RETURN
235 *
236 * End of DSTT21
237 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
double precision function dlansy(NORM, UPLO, N, A, LDA, WORK)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric matrix.
Definition: dlansy.f:124
subroutine dsyr(UPLO, N, ALPHA, X, INCX, A, LDA)
DSYR
Definition: dsyr.f:134
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: dlange.f:116
subroutine dsyr2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DSYR2
Definition: dsyr2.f:149
Here is the call graph for this function:
Here is the caller graph for this function: