LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
scsdts.f
Go to the documentation of this file.
1 *> \brief \b SCSDTS
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 SCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
12 * LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
13 * RWORK, RESULT )
14 *
15 * .. Scalar Arguments ..
16 * INTEGER LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
17 * ..
18 * .. Array Arguments ..
19 * INTEGER IWORK( * )
20 * REAL RESULT( 15 ), RWORK( * ), THETA( * )
21 * REAL U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
22 * $ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
23 * $ XF( LDX, * )
24 * ..
25 *
26 *
27 *> \par Purpose:
28 * =============
29 *>
30 *> \verbatim
31 *>
32 *> SCSDTS tests SORCSD, which, given an M-by-M partitioned orthogonal
33 *> matrix X,
34 *> Q M-Q
35 *> X = [ X11 X12 ] P ,
36 *> [ X21 X22 ] M-P
37 *>
38 *> computes the CSD
39 *>
40 *> [ U1 ]**T * [ X11 X12 ] * [ V1 ]
41 *> [ U2 ] [ X21 X22 ] [ V2 ]
42 *>
43 *> [ I 0 0 | 0 0 0 ]
44 *> [ 0 C 0 | 0 -S 0 ]
45 *> [ 0 0 0 | 0 0 -I ]
46 *> = [---------------------] = [ D11 D12 ] .
47 *> [ 0 0 0 | I 0 0 ] [ D21 D22 ]
48 *> [ 0 S 0 | 0 C 0 ]
49 *> [ 0 0 I | 0 0 0 ]
50 *>
51 *> and also SORCSD2BY1, which, given
52 *> Q
53 *> [ X11 ] P ,
54 *> [ X21 ] M-P
55 *>
56 *> computes the 2-by-1 CSD
57 *>
58 *> [ I 0 0 ]
59 *> [ 0 C 0 ]
60 *> [ 0 0 0 ]
61 *> [ U1 ]**T * [ X11 ] * V1 = [----------] = [ D11 ] ,
62 *> [ U2 ] [ X21 ] [ 0 0 0 ] [ D21 ]
63 *> [ 0 S 0 ]
64 *> [ 0 0 I ]
65 *> \endverbatim
66 *
67 * Arguments:
68 * ==========
69 *
70 *> \param[in] M
71 *> \verbatim
72 *> M is INTEGER
73 *> The number of rows of the matrix X. M >= 0.
74 *> \endverbatim
75 *>
76 *> \param[in] P
77 *> \verbatim
78 *> P is INTEGER
79 *> The number of rows of the matrix X11. P >= 0.
80 *> \endverbatim
81 *>
82 *> \param[in] Q
83 *> \verbatim
84 *> Q is INTEGER
85 *> The number of columns of the matrix X11. Q >= 0.
86 *> \endverbatim
87 *>
88 *> \param[in] X
89 *> \verbatim
90 *> X is REAL array, dimension (LDX,M)
91 *> The M-by-M matrix X.
92 *> \endverbatim
93 *>
94 *> \param[out] XF
95 *> \verbatim
96 *> XF is REAL array, dimension (LDX,M)
97 *> Details of the CSD of X, as returned by SORCSD;
98 *> see SORCSD for further details.
99 *> \endverbatim
100 *>
101 *> \param[in] LDX
102 *> \verbatim
103 *> LDX is INTEGER
104 *> The leading dimension of the arrays X and XF.
105 *> LDX >= max( 1,M ).
106 *> \endverbatim
107 *>
108 *> \param[out] U1
109 *> \verbatim
110 *> U1 is REAL array, dimension(LDU1,P)
111 *> The P-by-P orthogonal matrix U1.
112 *> \endverbatim
113 *>
114 *> \param[in] LDU1
115 *> \verbatim
116 *> LDU1 is INTEGER
117 *> The leading dimension of the array U1. LDU >= max(1,P).
118 *> \endverbatim
119 *>
120 *> \param[out] U2
121 *> \verbatim
122 *> U2 is REAL array, dimension(LDU2,M-P)
123 *> The (M-P)-by-(M-P) orthogonal matrix U2.
124 *> \endverbatim
125 *>
126 *> \param[in] LDU2
127 *> \verbatim
128 *> LDU2 is INTEGER
129 *> The leading dimension of the array U2. LDU >= max(1,M-P).
130 *> \endverbatim
131 *>
132 *> \param[out] V1T
133 *> \verbatim
134 *> V1T is REAL array, dimension(LDV1T,Q)
135 *> The Q-by-Q orthogonal matrix V1T.
136 *> \endverbatim
137 *>
138 *> \param[in] LDV1T
139 *> \verbatim
140 *> LDV1T is INTEGER
141 *> The leading dimension of the array V1T. LDV1T >=
142 *> max(1,Q).
143 *> \endverbatim
144 *>
145 *> \param[out] V2T
146 *> \verbatim
147 *> V2T is REAL array, dimension(LDV2T,M-Q)
148 *> The (M-Q)-by-(M-Q) orthogonal matrix V2T.
149 *> \endverbatim
150 *>
151 *> \param[in] LDV2T
152 *> \verbatim
153 *> LDV2T is INTEGER
154 *> The leading dimension of the array V2T. LDV2T >=
155 *> max(1,M-Q).
156 *> \endverbatim
157 *>
158 *> \param[out] THETA
159 *> \verbatim
160 *> THETA is REAL array, dimension MIN(P,M-P,Q,M-Q)
161 *> The CS values of X; the essentially diagonal matrices C and
162 *> S are constructed from THETA; see subroutine SORCSD for
163 *> details.
164 *> \endverbatim
165 *>
166 *> \param[out] IWORK
167 *> \verbatim
168 *> IWORK is INTEGER array, dimension (M)
169 *> \endverbatim
170 *>
171 *> \param[out] WORK
172 *> \verbatim
173 *> WORK is REAL array, dimension (LWORK)
174 *> \endverbatim
175 *>
176 *> \param[in] LWORK
177 *> \verbatim
178 *> LWORK is INTEGER
179 *> The dimension of the array WORK
180 *> \endverbatim
181 *>
182 *> \param[out] RWORK
183 *> \verbatim
184 *> RWORK is REAL array
185 *> \endverbatim
186 *>
187 *> \param[out] RESULT
188 *> \verbatim
189 *> RESULT is REAL array, dimension (15)
190 *> The test ratios:
191 *> First, the 2-by-2 CSD:
192 *> RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
193 *> RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 )
194 *> RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
195 *> RESULT(4) = norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 )
196 *> RESULT(5) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
197 *> RESULT(6) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
198 *> RESULT(7) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
199 *> RESULT(8) = norm( I - V2T'*V2T ) / ( MAX(1,M-Q)*ULP )
200 *> RESULT(9) = 0 if THETA is in increasing order and
201 *> all angles are in [0,pi/2];
202 *> = ULPINV otherwise.
203 *> Then, the 2-by-1 CSD:
204 *> RESULT(10) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
205 *> RESULT(11) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
206 *> RESULT(12) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
207 *> RESULT(13) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
208 *> RESULT(14) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
209 *> RESULT(15) = 0 if THETA is in increasing order and
210 *> all angles are in [0,pi/2];
211 *> = ULPINV otherwise.
212 *> ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
213 *> \endverbatim
214 *
215 * Authors:
216 * ========
217 *
218 *> \author Univ. of Tennessee
219 *> \author Univ. of California Berkeley
220 *> \author Univ. of Colorado Denver
221 *> \author NAG Ltd.
222 *
223 *> \date November 2011
224 *
225 *> \ingroup single_eig
226 *
227 * =====================================================================
228  SUBROUTINE scsdts( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
229  $ ldv1t, v2t, ldv2t, theta, iwork, work, lwork,
230  $ rwork, result )
231 *
232 * -- LAPACK test routine (version 3.4.0) --
233 * -- LAPACK is a software package provided by Univ. of Tennessee, --
234 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
235 * November 2011
236 *
237 * .. Scalar Arguments ..
238  INTEGER ldx, ldu1, ldu2, ldv1t, ldv2t, lwork, m, p, q
239 * ..
240 * .. Array Arguments ..
241  INTEGER iwork( * )
242  REAL result( 15 ), rwork( * ), theta( * )
243  REAL u1( ldu1, * ), u2( ldu2, * ), v1t( ldv1t, * ),
244  $ v2t( ldv2t, * ), work( lwork ), x( ldx, * ),
245  $ xf( ldx, * )
246 * ..
247 *
248 * =====================================================================
249 *
250 * .. Parameters ..
251  REAL piover2, realone, realzero
252  parameter( piover2 = 1.57079632679489662e0,
253  $ realone = 1.0e0, realzero = 0.0e0 )
254  REAL zero, one
255  parameter( zero = 0.0e0, one = 1.0e0 )
256 * ..
257 * .. Local Scalars ..
258  INTEGER i, info, r
259  REAL eps2, resid, ulp, ulpinv
260 * ..
261 * .. External Functions ..
262  REAL slamch, slange, slansy
263  EXTERNAL slamch, slange, slansy
264 * ..
265 * .. External Subroutines ..
266  EXTERNAL sgemm, slacpy, slaset, sorcsd, sorcsd2by1,
267  $ ssyrk
268 * ..
269 * .. Intrinsic Functions ..
270  INTRINSIC cos, max, min, REAL, sin
271 * ..
272 * .. Executable Statements ..
273 *
274  ulp = slamch( 'Precision' )
275  ulpinv = realone / ulp
276 *
277 * The first half of the routine checks the 2-by-2 CSD
278 *
279  CALL slaset( 'Full', m, m, zero, one, work, ldx )
280  CALL ssyrk( 'Upper', 'Conjugate transpose', m, m, -one, x, ldx,
281  $ one, work, ldx )
282  IF (m.GT.0) THEN
283  eps2 = max( ulp,
284  $ slange( '1', m, m, work, ldx, rwork ) / REAL( M ) )
285  ELSE
286  eps2 = ulp
287  END IF
288  r = min( p, m-p, q, m-q )
289 *
290 * Copy the matrix X to the array XF.
291 *
292  CALL slacpy( 'Full', m, m, x, ldx, xf, ldx )
293 *
294 * Compute the CSD
295 *
296  CALL sorcsd( 'Y', 'Y', 'Y', 'Y', 'N', 'D', m, p, q, xf(1,1), ldx,
297  $ xf(1,q+1), ldx, xf(p+1,1), ldx, xf(p+1,q+1), ldx,
298  $ theta, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t,
299  $ work, lwork, iwork, info )
300 *
301 * Compute XF := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
302 *
303  CALL slacpy( 'Full', m, m, x, ldx, xf, ldx )
304 *
305  CALL sgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
306  $ xf, ldx, v1t, ldv1t, zero, work, ldx )
307 *
308  CALL sgemm( 'Conjugate transpose', 'No transpose', p, q, p, one,
309  $ u1, ldu1, work, ldx, zero, xf, ldx )
310 *
311  DO i = 1, min(p,q)-r
312  xf(i,i) = xf(i,i) - one
313  END DO
314  DO i = 1, r
315  xf(min(p,q)-r+i,min(p,q)-r+i) =
316  $ xf(min(p,q)-r+i,min(p,q)-r+i) - cos(theta(i))
317  END DO
318 *
319  CALL sgemm( 'No transpose', 'Conjugate transpose', p, m-q, m-q,
320  $ one, xf(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
321 *
322  CALL sgemm( 'Conjugate transpose', 'No transpose', p, m-q, p,
323  $ one, u1, ldu1, work, ldx, zero, xf(1,q+1), ldx )
324 *
325  DO i = 1, min(p,m-q)-r
326  xf(p-i+1,m-i+1) = xf(p-i+1,m-i+1) + one
327  END DO
328  DO i = 1, r
329  xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) =
330  $ xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) +
331  $ sin(theta(r-i+1))
332  END DO
333 *
334  CALL sgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
335  $ xf(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
336 *
337  CALL sgemm( 'Conjugate transpose', 'No transpose', m-p, q, m-p,
338  $ one, u2, ldu2, work, ldx, zero, xf(p+1,1), ldx )
339 *
340  DO i = 1, min(m-p,q)-r
341  xf(m-i+1,q-i+1) = xf(m-i+1,q-i+1) - one
342  END DO
343  DO i = 1, r
344  xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
345  $ xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
346  $ sin(theta(r-i+1))
347  END DO
348 *
349  CALL sgemm( 'No transpose', 'Conjugate transpose', m-p, m-q, m-q,
350  $ one, xf(p+1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
351 *
352  CALL sgemm( 'Conjugate transpose', 'No transpose', m-p, m-q, m-p,
353  $ one, u2, ldu2, work, ldx, zero, xf(p+1,q+1), ldx )
354 *
355  DO i = 1, min(m-p,m-q)-r
356  xf(p+i,q+i) = xf(p+i,q+i) - one
357  END DO
358  DO i = 1, r
359  xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) =
360  $ xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) -
361  $ cos(theta(i))
362  END DO
363 *
364 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
365 *
366  resid = slange( '1', p, q, xf, ldx, rwork )
367  result( 1 ) = ( resid / REAL(MAX(1,P,Q)) ) / eps2
368 *
369 * Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
370 *
371  resid = slange( '1', p, m-q, xf(1,q+1), ldx, rwork )
372  result( 2 ) = ( resid / REAL(MAX(1,P,M-Q)) ) / eps2
373 *
374 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
375 *
376  resid = slange( '1', m-p, q, xf(p+1,1), ldx, rwork )
377  result( 3 ) = ( resid / REAL(MAX(1,M-P,Q)) ) / eps2
378 *
379 * Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
380 *
381  resid = slange( '1', m-p, m-q, xf(p+1,q+1), ldx, rwork )
382  result( 4 ) = ( resid / REAL(MAX(1,M-P,M-Q)) ) / eps2
383 *
384 * Compute I - U1'*U1
385 *
386  CALL slaset( 'Full', p, p, zero, one, work, ldu1 )
387  CALL ssyrk( 'Upper', 'Conjugate transpose', p, p, -one, u1, ldu1,
388  $ one, work, ldu1 )
389 *
390 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
391 *
392  resid = slansy( '1', 'Upper', p, work, ldu1, rwork )
393  result( 5 ) = ( resid / REAL(MAX(1,P)) ) / ulp
394 *
395 * Compute I - U2'*U2
396 *
397  CALL slaset( 'Full', m-p, m-p, zero, one, work, ldu2 )
398  CALL ssyrk( 'Upper', 'Conjugate transpose', m-p, m-p, -one, u2,
399  $ ldu2, one, work, ldu2 )
400 *
401 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
402 *
403  resid = slansy( '1', 'Upper', m-p, work, ldu2, rwork )
404  result( 6 ) = ( resid / REAL(MAX(1,M-P)) ) / ulp
405 *
406 * Compute I - V1T*V1T'
407 *
408  CALL slaset( 'Full', q, q, zero, one, work, ldv1t )
409  CALL ssyrk( 'Upper', 'No transpose', q, q, -one, v1t, ldv1t, one,
410  $ work, ldv1t )
411 *
412 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
413 *
414  resid = slansy( '1', 'Upper', q, work, ldv1t, rwork )
415  result( 7 ) = ( resid / REAL(MAX(1,Q)) ) / ulp
416 *
417 * Compute I - V2T*V2T'
418 *
419  CALL slaset( 'Full', m-q, m-q, zero, one, work, ldv2t )
420  CALL ssyrk( 'Upper', 'No transpose', m-q, m-q, -one, v2t, ldv2t,
421  $ one, work, ldv2t )
422 *
423 * Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
424 *
425  resid = slansy( '1', 'Upper', m-q, work, ldv2t, rwork )
426  result( 8 ) = ( resid / REAL(MAX(1,M-Q)) ) / ulp
427 *
428 * Check sorting
429 *
430  result( 9 ) = realzero
431  DO i = 1, r
432  IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
433  result( 9 ) = ulpinv
434  END IF
435  IF( i.GT.1 ) THEN
436  IF ( theta(i).LT.theta(i-1) ) THEN
437  result( 9 ) = ulpinv
438  END IF
439  END IF
440  END DO
441 *
442 * The second half of the routine checks the 2-by-1 CSD
443 *
444  CALL slaset( 'Full', q, q, zero, one, work, ldx )
445  CALL ssyrk( 'Upper', 'Conjugate transpose', q, m, -one, x, ldx,
446  $ one, work, ldx )
447  IF (m.GT.0) THEN
448  eps2 = max( ulp,
449  $ slange( '1', q, q, work, ldx, rwork ) / REAL( M ) )
450  ELSE
451  eps2 = ulp
452  END IF
453  r = min( p, m-p, q, m-q )
454 *
455 * Copy the matrix [X11;X21] to the array XF.
456 *
457  CALL slacpy( 'Full', m, q, x, ldx, xf, ldx )
458 *
459 * Compute the CSD
460 *
461  CALL sorcsd2by1( 'Y', 'Y', 'Y', m, p, q, xf(1,1), ldx, xf(p+1,1),
462  $ ldx, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, work,
463  $ lwork, iwork, info )
464 *
465 * Compute [X11;X21] := diag(U1,U2)'*[X11;X21]*V1 - [D11;D21]
466 *
467  CALL sgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
468  $ x, ldx, v1t, ldv1t, zero, work, ldx )
469 *
470  CALL sgemm( 'Conjugate transpose', 'No transpose', p, q, p, one,
471  $ u1, ldu1, work, ldx, zero, x, ldx )
472 *
473  DO i = 1, min(p,q)-r
474  x(i,i) = x(i,i) - one
475  END DO
476  DO i = 1, r
477  x(min(p,q)-r+i,min(p,q)-r+i) =
478  $ x(min(p,q)-r+i,min(p,q)-r+i) - cos(theta(i))
479  END DO
480 *
481  CALL sgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
482  $ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
483 *
484  CALL sgemm( 'Conjugate transpose', 'No transpose', m-p, q, m-p,
485  $ one, u2, ldu2, work, ldx, zero, x(p+1,1), ldx )
486 *
487  DO i = 1, min(m-p,q)-r
488  x(m-i+1,q-i+1) = x(m-i+1,q-i+1) - one
489  END DO
490  DO i = 1, r
491  x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
492  $ x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
493  $ sin(theta(r-i+1))
494  END DO
495 *
496 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
497 *
498  resid = slange( '1', p, q, x, ldx, rwork )
499  result( 10 ) = ( resid / REAL(MAX(1,P,Q)) ) / eps2
500 *
501 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
502 *
503  resid = slange( '1', m-p, q, x(p+1,1), ldx, rwork )
504  result( 11 ) = ( resid / REAL(MAX(1,M-P,Q)) ) / eps2
505 *
506 * Compute I - U1'*U1
507 *
508  CALL slaset( 'Full', p, p, zero, one, work, ldu1 )
509  CALL ssyrk( 'Upper', 'Conjugate transpose', p, p, -one, u1, ldu1,
510  $ one, work, ldu1 )
511 *
512 * Compute norm( I - U1'*U1 ) / ( MAX(1,P) * ULP ) .
513 *
514  resid = slansy( '1', 'Upper', p, work, ldu1, rwork )
515  result( 12 ) = ( resid / REAL(MAX(1,P)) ) / ulp
516 *
517 * Compute I - U2'*U2
518 *
519  CALL slaset( 'Full', m-p, m-p, zero, one, work, ldu2 )
520  CALL ssyrk( 'Upper', 'Conjugate transpose', m-p, m-p, -one, u2,
521  $ ldu2, one, work, ldu2 )
522 *
523 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
524 *
525  resid = slansy( '1', 'Upper', m-p, work, ldu2, rwork )
526  result( 13 ) = ( resid / REAL(MAX(1,M-P)) ) / ulp
527 *
528 * Compute I - V1T*V1T'
529 *
530  CALL slaset( 'Full', q, q, zero, one, work, ldv1t )
531  CALL ssyrk( 'Upper', 'No transpose', q, q, -one, v1t, ldv1t, one,
532  $ work, ldv1t )
533 *
534 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
535 *
536  resid = slansy( '1', 'Upper', q, work, ldv1t, rwork )
537  result( 14 ) = ( resid / REAL(MAX(1,Q)) ) / ulp
538 *
539 * Check sorting
540 *
541  result( 15 ) = realzero
542  DO i = 1, r
543  IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
544  result( 15 ) = ulpinv
545  END IF
546  IF( i.GT.1 ) THEN
547  IF ( theta(i).LT.theta(i-1) ) THEN
548  result( 15 ) = ulpinv
549  END IF
550  END IF
551  END DO
552 *
553  RETURN
554 *
555 * End of SCSDTS
556 *
557  END
558