3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
11 * SUBROUTINE SORT01( ROWCOL, M, N, U, LDU, WORK, LWORK, RESID )
13 * .. Scalar Arguments ..
15 * INTEGER LDU, LWORK, M, N
18 * .. Array Arguments ..
19 * REAL U( LDU, * ), WORK( * )
28 *> SORT01 checks that the matrix U is orthogonal by computing the ratio
30 *> RESID = norm( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
32 *> RESID = norm( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
34 *> Alternatively, if there isn't sufficient workspace to form
35 *> I - U*U' or I - U'*U, the ratio is computed as
37 *> RESID = abs( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
39 *> RESID = abs( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
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
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
61 *> The number of rows of the matrix U.
67 *> The number of columns of the matrix U.
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'.
81 *> The leading dimension of the array U. LDU >= max(1,M).
86 *> WORK is REAL array, dimension (LWORK)
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.
100 *> RESID = norm( I - U * U' ) / ( n * EPS ), if ROWCOL = 'R', or
101 *> RESID = norm( I - U' * U ) / ( m * EPS ), if ROWCOL = 'C'.
107 *> \author Univ. of Tennessee
108 *> \author Univ. of California Berkeley
109 *> \author Univ. of Colorado Denver
112 *> \date November 2011
114 *> \ingroup single_eig
116 * =====================================================================
117 SUBROUTINE SORT01( ROWCOL, M, N, U, LDU, WORK, LWORK, RESID )
119 * -- LAPACK test routine (version 3.4.0) --
120 * -- LAPACK is a software package provided by Univ. of Tennessee, --
121 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
124 * .. Scalar Arguments ..
126 INTEGER LDU, LWORK, M, N
129 * .. Array Arguments ..
130 REAL U( LDU, * ), WORK( * )
133 * =====================================================================
137 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
139 * .. Local Scalars ..
141 INTEGER I, J, K, LDWORK, MNMIN
144 * .. External Functions ..
146 REAL SDOT, SLAMCH, SLANSY
147 EXTERNAL LSAME, SDOT, SLAMCH, SLANSY
149 * .. External Subroutines ..
150 EXTERNAL SLASET, SSYRK
152 * .. Intrinsic Functions ..
153 INTRINSIC MAX, MIN, REAL
155 * .. Executable Statements ..
159 * Quick return if possible
161 IF( M.LE.0 .OR. N.LE.0 )
164 EPS = SLAMCH( 'Precision' )
165 IF( M.LT.N .OR. ( M.EQ.N .AND. LSAME( ROWCOL, 'R' ) ) ) THEN
174 IF( ( MNMIN+1 )*MNMIN.LE.LWORK ) THEN
179 IF( LDWORK.GT.0 ) THEN
181 * Compute I - U*U' or I - U'*U.
183 CALL SLASET( 'Upper', MNMIN, MNMIN, ZERO, ONE, WORK, LDWORK )
184 CALL SSYRK( 'Upper', TRANSU, MNMIN, K, -ONE, U, LDU, ONE, WORK,
187 * Compute norm( I - U*U' ) / ( K * EPS ) .
189 RESID = SLANSY( '1', 'Upper', MNMIN, WORK, LDWORK,
190 $ WORK( LDWORK*MNMIN+1 ) )
191 RESID = ( RESID / REAL( K ) ) / EPS
192 ELSE IF( TRANSU.EQ.'T' ) THEN
194 * Find the maximum element in abs( I - U'*U ) / ( m * EPS )
203 TMP = TMP - SDOT( M, U( 1, I ), 1, U( 1, J ), 1 )
204 RESID = MAX( RESID, ABS( TMP ) )
207 RESID = ( RESID / REAL( M ) ) / EPS
210 * Find the maximum element in abs( I - U*U' ) / ( n * EPS )
219 TMP = TMP - SDOT( N, U( J, 1 ), LDU, U( I, 1 ), LDU )
220 RESID = MAX( RESID, ABS( TMP ) )
223 RESID = ( RESID / REAL( N ) ) / EPS