LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sort01.f
Go to the documentation of this file.
1*> \brief \b SORT01
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 SORT01( ROWCOL, M, N, U, LDU, WORK, LWORK, RESID )
12*
13* .. Scalar Arguments ..
14* CHARACTER ROWCOL
15* INTEGER LDU, LWORK, M, N
16* REAL RESID
17* ..
18* .. Array Arguments ..
19* REAL U( LDU, * ), WORK( * )
20* ..
21*
22*
23*> \par Purpose:
24* =============
25*>
26*> \verbatim
27*>
28*> SORT01 checks that the matrix U is orthogonal by computing the ratio
29*>
30*> RESID = norm( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
31*> or
32*> RESID = norm( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
33*>
34*> Alternatively, if there isn't sufficient workspace to form
35*> I - U*U' or I - U'*U, the ratio is computed as
36*>
37*> RESID = abs( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
38*> or
39*> RESID = abs( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
40*>
41*> where EPS is the machine precision. ROWCOL is used only if m = n;
42*> if m > n, ROWCOL is assumed to be 'C', and if m < n, ROWCOL is
43*> assumed to be 'R'.
44*> \endverbatim
45*
46* Arguments:
47* ==========
48*
49*> \param[in] ROWCOL
50*> \verbatim
51*> ROWCOL is CHARACTER
52*> Specifies whether the rows or columns of U should be checked
53*> for orthogonality. Used only if M = N.
54*> = 'R': Check for orthogonal rows of U
55*> = 'C': Check for orthogonal columns of U
56*> \endverbatim
57*>
58*> \param[in] M
59*> \verbatim
60*> M is INTEGER
61*> The number of rows of the matrix U.
62*> \endverbatim
63*>
64*> \param[in] N
65*> \verbatim
66*> N is INTEGER
67*> The number of columns of the matrix U.
68*> \endverbatim
69*>
70*> \param[in] U
71*> \verbatim
72*> U is REAL array, dimension (LDU,N)
73*> The orthogonal matrix U. U is checked for orthogonal columns
74*> if m > n or if m = n and ROWCOL = 'C'. U is checked for
75*> orthogonal rows if m < n or if m = n and ROWCOL = 'R'.
76*> \endverbatim
77*>
78*> \param[in] LDU
79*> \verbatim
80*> LDU is INTEGER
81*> The leading dimension of the array U. LDU >= max(1,M).
82*> \endverbatim
83*>
84*> \param[out] WORK
85*> \verbatim
86*> WORK is REAL array, dimension (LWORK)
87*> \endverbatim
88*>
89*> \param[in] LWORK
90*> \verbatim
91*> LWORK is INTEGER
92*> The length of the array WORK. For best performance, LWORK
93*> should be at least N*(N+1) if ROWCOL = 'C' or M*(M+1) if
94*> ROWCOL = 'R', but the test will be done even if LWORK is 0.
95*> \endverbatim
96*>
97*> \param[out] RESID
98*> \verbatim
99*> RESID is REAL
100*> RESID = norm( I - U * U' ) / ( n * EPS ), if ROWCOL = 'R', or
101*> RESID = norm( I - U' * U ) / ( m * EPS ), if ROWCOL = 'C'.
102*> \endverbatim
103*
104* Authors:
105* ========
106*
107*> \author Univ. of Tennessee
108*> \author Univ. of California Berkeley
109*> \author Univ. of Colorado Denver
110*> \author NAG Ltd.
111*
112*> \ingroup single_eig
113*
114* =====================================================================
115 SUBROUTINE sort01( ROWCOL, M, N, U, LDU, WORK, LWORK, RESID )
116*
117* -- LAPACK test routine --
118* -- LAPACK is a software package provided by Univ. of Tennessee, --
119* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
120*
121* .. Scalar Arguments ..
122 CHARACTER ROWCOL
123 INTEGER LDU, LWORK, M, N
124 REAL RESID
125* ..
126* .. Array Arguments ..
127 REAL U( LDU, * ), WORK( * )
128* ..
129*
130* =====================================================================
131*
132* .. Parameters ..
133 REAL ZERO, ONE
134 parameter( zero = 0.0e+0, one = 1.0e+0 )
135* ..
136* .. Local Scalars ..
137 CHARACTER TRANSU
138 INTEGER I, J, K, LDWORK, MNMIN
139 REAL EPS, TMP
140* ..
141* .. External Functions ..
142 LOGICAL LSAME
143 REAL SDOT, SLAMCH, SLANSY
144 EXTERNAL lsame, sdot, slamch, slansy
145* ..
146* .. External Subroutines ..
147 EXTERNAL slaset, ssyrk
148* ..
149* .. Intrinsic Functions ..
150 INTRINSIC max, min, real
151* ..
152* .. Executable Statements ..
153*
154 resid = zero
155*
156* Quick return if possible
157*
158 IF( m.LE.0 .OR. n.LE.0 )
159 $ RETURN
160*
161 eps = slamch( 'Precision' )
162 IF( m.LT.n .OR. ( m.EQ.n .AND. lsame( rowcol, 'R' ) ) ) THEN
163 transu = 'N'
164 k = n
165 ELSE
166 transu = 'T'
167 k = m
168 END IF
169 mnmin = min( m, n )
170*
171 IF( ( mnmin+1 )*mnmin.LE.lwork ) THEN
172 ldwork = mnmin
173 ELSE
174 ldwork = 0
175 END IF
176 IF( ldwork.GT.0 ) THEN
177*
178* Compute I - U*U' or I - U'*U.
179*
180 CALL slaset( 'Upper', mnmin, mnmin, zero, one, work, ldwork )
181 CALL ssyrk( 'Upper', transu, mnmin, k, -one, u, ldu, one, work,
182 $ ldwork )
183*
184* Compute norm( I - U*U' ) / ( K * EPS ) .
185*
186 resid = slansy( '1', 'Upper', mnmin, work, ldwork,
187 $ work( ldwork*mnmin+1 ) )
188 resid = ( resid / real( k ) ) / eps
189 ELSE IF( transu.EQ.'T' ) THEN
190*
191* Find the maximum element in abs( I - U'*U ) / ( m * EPS )
192*
193 DO 20 j = 1, n
194 DO 10 i = 1, j
195 IF( i.NE.j ) THEN
196 tmp = zero
197 ELSE
198 tmp = one
199 END IF
200 tmp = tmp - sdot( m, u( 1, i ), 1, u( 1, j ), 1 )
201 resid = max( resid, abs( tmp ) )
202 10 CONTINUE
203 20 CONTINUE
204 resid = ( resid / real( m ) ) / eps
205 ELSE
206*
207* Find the maximum element in abs( I - U*U' ) / ( n * EPS )
208*
209 DO 40 j = 1, m
210 DO 30 i = 1, j
211 IF( i.NE.j ) THEN
212 tmp = zero
213 ELSE
214 tmp = one
215 END IF
216 tmp = tmp - sdot( n, u( j, 1 ), ldu, u( i, 1 ), ldu )
217 resid = max( resid, abs( tmp ) )
218 30 CONTINUE
219 40 CONTINUE
220 resid = ( resid / real( n ) ) / eps
221 END IF
222 RETURN
223*
224* End of SORT01
225*
226 END
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
Definition: ssyrk.f:169
subroutine sort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
SORT01
Definition: sort01.f:116