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
227