LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ssyt22.f
Go to the documentation of this file.
1*> \brief \b SSYT22
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE SSYT22( ITYPE, UPLO, N, M, KBAND, A, LDA, D, E, U, LDU,
12* V, LDV, TAU, WORK, RESULT )
13*
14* .. Scalar Arguments ..
15* CHARACTER UPLO
16* INTEGER ITYPE, KBAND, LDA, LDU, LDV, M, N
17* ..
18* .. Array Arguments ..
19* REAL A( LDA, * ), D( * ), E( * ), RESULT( 2 ),
20* $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
21* ..
22*
23*
24*> \par Purpose:
25* =============
26*>
27*> \verbatim
28*>
29*> SSYT22 generally checks a decomposition of the form
30*>
31*> A U = U S
32*>
33*> where A is symmetric, the columns of U are orthonormal, and S
34*> is diagonal (if KBAND=0) or symmetric tridiagonal (if
35*> KBAND=1). If ITYPE=1, then U is represented as a dense matrix,
36*> otherwise the U is expressed as a product of Householder
37*> transformations, whose vectors are stored in the array "V" and
38*> whose scaling constants are in "TAU"; we shall use the letter
39*> "V" to refer to the product of Householder transformations
40*> (which should be equal to U).
41*>
42*> Specifically, if ITYPE=1, then:
43*>
44*> RESULT(1) = | U**T A U - S | / ( |A| m ulp ) and
45*> RESULT(2) = | I - U**T U | / ( m ulp )
46*> \endverbatim
47*
48* Arguments:
49* ==========
50*
51*> \verbatim
52*> ITYPE INTEGER
53*> Specifies the type of tests to be performed.
54*> 1: U expressed as a dense orthogonal matrix:
55*> RESULT(1) = | A - U S U**T | / ( |A| n ulp ) and
56*> RESULT(2) = | I - U U**T | / ( n ulp )
57*>
58*> UPLO CHARACTER
59*> If UPLO='U', the upper triangle of A will be used and the
60*> (strictly) lower triangle will not be referenced. If
61*> UPLO='L', the lower triangle of A will be used and the
62*> (strictly) upper triangle will not be referenced.
63*> Not modified.
64*>
65*> N INTEGER
66*> The size of the matrix. If it is zero, SSYT22 does nothing.
67*> It must be at least zero.
68*> Not modified.
69*>
70*> M INTEGER
71*> The number of columns of U. If it is zero, SSYT22 does
72*> nothing. It must be at least zero.
73*> Not modified.
74*>
75*> KBAND INTEGER
76*> The bandwidth of the matrix. It may only be zero or one.
77*> If zero, then S is diagonal, and E is not referenced. If
78*> one, then S is symmetric tri-diagonal.
79*> Not modified.
80*>
81*> A REAL array, dimension (LDA , N)
82*> The original (unfactored) matrix. It is assumed to be
83*> symmetric, and only the upper (UPLO='U') or only the lower
84*> (UPLO='L') will be referenced.
85*> Not modified.
86*>
87*> LDA INTEGER
88*> The leading dimension of A. It must be at least 1
89*> and at least N.
90*> Not modified.
91*>
92*> D REAL array, dimension (N)
93*> The diagonal of the (symmetric tri-) diagonal matrix.
94*> Not modified.
95*>
96*> E REAL array, dimension (N)
97*> The off-diagonal of the (symmetric tri-) diagonal matrix.
98*> E(1) is ignored, E(2) is the (1,2) and (2,1) element, etc.
99*> Not referenced if KBAND=0.
100*> Not modified.
101*>
102*> U REAL array, dimension (LDU, N)
103*> If ITYPE=1 or 3, this contains the orthogonal matrix in
104*> the decomposition, expressed as a dense matrix. If ITYPE=2,
105*> then it is not referenced.
106*> Not modified.
107*>
108*> LDU INTEGER
109*> The leading dimension of U. LDU must be at least N and
110*> at least 1.
111*> Not modified.
112*>
113*> V REAL array, dimension (LDV, N)
114*> If ITYPE=2 or 3, the lower triangle of this array contains
115*> the Householder vectors used to describe the orthogonal
116*> matrix in the decomposition. If ITYPE=1, then it is not
117*> referenced.
118*> Not modified.
119*>
120*> LDV INTEGER
121*> The leading dimension of V. LDV must be at least N and
122*> at least 1.
123*> Not modified.
124*>
125*> TAU REAL array, dimension (N)
126*> If ITYPE >= 2, then TAU(j) is the scalar factor of
127*> v(j) v(j)**T in the Householder transformation H(j) of
128*> the product U = H(1)...H(n-2)
129*> If ITYPE < 2, then TAU is not referenced.
130*> Not modified.
131*>
132*> WORK REAL array, dimension (2*N**2)
133*> Workspace.
134*> Modified.
135*>
136*> RESULT REAL array, dimension (2)
137*> The values computed by the two tests described above. The
138*> values are currently limited to 1/ulp, to avoid overflow.
139*> RESULT(1) is always modified. RESULT(2) is modified only
140*> if LDU is at least N.
141*> Modified.
142*> \endverbatim
143*
144* Authors:
145* ========
146*
147*> \author Univ. of Tennessee
148*> \author Univ. of California Berkeley
149*> \author Univ. of Colorado Denver
150*> \author NAG Ltd.
151*
152*> \ingroup single_eig
153*
154* =====================================================================
155 SUBROUTINE ssyt22( ITYPE, UPLO, N, M, KBAND, A, LDA, D, E, U, LDU,
156 $ V, LDV, TAU, WORK, RESULT )
157*
158* -- LAPACK test routine --
159* -- LAPACK is a software package provided by Univ. of Tennessee, --
160* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
161*
162* .. Scalar Arguments ..
163 CHARACTER UPLO
164 INTEGER ITYPE, KBAND, LDA, LDU, LDV, M, N
165* ..
166* .. Array Arguments ..
167 REAL A( LDA, * ), D( * ), E( * ), RESULT( 2 ),
168 $ tau( * ), u( ldu, * ), v( ldv, * ), work( * )
169* ..
170*
171* =====================================================================
172*
173* .. Parameters ..
174 REAL ZERO, ONE
175 parameter( zero = 0.0e0, one = 1.0e0 )
176* ..
177* .. Local Scalars ..
178 INTEGER J, JJ, JJ1, JJ2, NN, NNP1
179 REAL ANORM, ULP, UNFL, WNORM
180* ..
181* .. External Functions ..
182 REAL SLAMCH, SLANSY
183 EXTERNAL slamch, slansy
184* ..
185* .. External Subroutines ..
186 EXTERNAL sgemm, ssymm
187* ..
188* .. Intrinsic Functions ..
189 INTRINSIC max, min, real
190* ..
191* .. Executable Statements ..
192*
193 result( 1 ) = zero
194 result( 2 ) = zero
195 IF( n.LE.0 .OR. m.LE.0 )
196 $ RETURN
197*
198 unfl = slamch( 'Safe minimum' )
199 ulp = slamch( 'Precision' )
200*
201* Do Test 1
202*
203* Norm of A:
204*
205 anorm = max( slansy( '1', uplo, n, a, lda, work ), unfl )
206*
207* Compute error matrix:
208*
209* ITYPE=1: error = U**T A U - S
210*
211 CALL ssymm( 'L', uplo, n, m, one, a, lda, u, ldu, zero, work, n )
212 nn = n*n
213 nnp1 = nn + 1
214 CALL sgemm( 'T', 'N', m, m, n, one, u, ldu, work, n, zero,
215 $ work( nnp1 ), n )
216 DO 10 j = 1, m
217 jj = nn + ( j-1 )*n + j
218 work( jj ) = work( jj ) - d( j )
219 10 CONTINUE
220 IF( kband.EQ.1 .AND. n.GT.1 ) THEN
221 DO 20 j = 2, m
222 jj1 = nn + ( j-1 )*n + j - 1
223 jj2 = nn + ( j-2 )*n + j
224 work( jj1 ) = work( jj1 ) - e( j-1 )
225 work( jj2 ) = work( jj2 ) - e( j-1 )
226 20 CONTINUE
227 END IF
228 wnorm = slansy( '1', uplo, m, work( nnp1 ), n, work( 1 ) )
229*
230 IF( anorm.GT.wnorm ) THEN
231 result( 1 ) = ( wnorm / anorm ) / ( m*ulp )
232 ELSE
233 IF( anorm.LT.one ) THEN
234 result( 1 ) = ( min( wnorm, m*anorm ) / anorm ) / ( m*ulp )
235 ELSE
236 result( 1 ) = min( wnorm / anorm, real( m ) ) / ( m*ulp )
237 END IF
238 END IF
239*
240* Do Test 2
241*
242* Compute U**T U - I
243*
244 IF( itype.EQ.1 )
245 $ CALL sort01( 'Columns', n, m, u, ldu, work, 2*n*n,
246 $ result( 2 ) )
247*
248 RETURN
249*
250* End of SSYT22
251*
252 END
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine ssymm(side, uplo, m, n, alpha, a, lda, b, ldb, beta, c, ldc)
SSYMM
Definition ssymm.f:189
subroutine sort01(rowcol, m, n, u, ldu, work, lwork, resid)
SORT01
Definition sort01.f:116
subroutine ssyt22(itype, uplo, n, m, kband, a, lda, d, e, u, ldu, v, ldv, tau, work, result)
SSYT22
Definition ssyt22.f:157