OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
zcsdts.f
Go to the documentation of this file.
1*> \brief \b ZCSDTS
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 ZCSDTS( 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* DOUBLE PRECISION RESULT( 15 ), RWORK( * ), THETA( * )
21* COMPLEX*16 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*> ZCSDTS tests ZUNCSD, which, given an M-by-M partitioned unitary
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 COMPLEX*16 array, dimension (LDX,M)
91*> The M-by-M matrix X.
92*> \endverbatim
93*>
94*> \param[out] XF
95*> \verbatim
96*> XF is COMPLEX*16 array, dimension (LDX,M)
97*> Details of the CSD of X, as returned by ZUNCSD;
98*> see ZUNCSD 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 COMPLEX*16 array, dimension(LDU1,P)
111*> The P-by-P unitary 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 COMPLEX*16 array, dimension(LDU2,M-P)
123*> The (M-P)-by-(M-P) unitary 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 COMPLEX*16 array, dimension(LDV1T,Q)
135*> The Q-by-Q unitary 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 COMPLEX*16 array, dimension(LDV2T,M-Q)
148*> The (M-Q)-by-(M-Q) unitary 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 DOUBLE PRECISION 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 ZUNCSD 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 COMPLEX*16 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 DOUBLE PRECISION array
185*> \endverbatim
186*>
187*> \param[out] RESULT
188*> \verbatim
189*> RESULT is DOUBLE PRECISION 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*> \ingroup complex16_eig
224*
225* =====================================================================
226 SUBROUTINE zcsdts( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
227 $ LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
228 $ RWORK, RESULT )
229*
230* -- LAPACK test routine --
231* -- LAPACK is a software package provided by Univ. of Tennessee, --
232* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
233*
234* .. Scalar Arguments ..
235 INTEGER LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
236* ..
237* .. Array Arguments ..
238 INTEGER IWORK( * )
239 DOUBLE PRECISION RESULT( 15 ), RWORK( * ), THETA( * )
240 COMPLEX*16 U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
241 $ v2t( ldv2t, * ), work( lwork ), x( ldx, * ),
242 $ xf( ldx, * )
243* ..
244*
245* =====================================================================
246*
247* .. Parameters ..
248 DOUBLE PRECISION REALONE, REALZERO
249 PARAMETER ( REALONE = 1.0d0, realzero = 0.0d0 )
250 COMPLEX*16 ZERO, ONE
251 parameter( zero = (0.0d0,0.0d0), one = (1.0d0,0.0d0) )
252 DOUBLE PRECISION PIOVER2
253 parameter( piover2 = 1.57079632679489661923132169163975144210d0 )
254* ..
255* .. Local Scalars ..
256 INTEGER I, INFO, R
257 DOUBLE PRECISION EPS2, RESID, ULP, ULPINV
258* ..
259* .. External Functions ..
260 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHE
261 EXTERNAL DLAMCH, ZLANGE, ZLANHE
262* ..
263* .. External Subroutines ..
264 EXTERNAL zgemm, zherk, zlacpy, zlaset, zuncsd,
265 $ zuncsd2by1
266* ..
267* .. Intrinsic Functions ..
268 INTRINSIC cos, dble, dcmplx, max, min, real, sin
269* ..
270* .. Executable Statements ..
271*
272 ulp = dlamch( 'Precision' )
273 ulpinv = realone / ulp
274*
275* The first half of the routine checks the 2-by-2 CSD
276*
277 CALL zlaset( 'Full', m, m, zero, one, work, ldx )
278 CALL zherk( 'Upper', 'Conjugate transpose', m, m, -realone,
279 $ x, ldx, realone, work, ldx )
280 IF (m.GT.0) THEN
281 eps2 = max( ulp,
282 $ zlange( '1', m, m, work, ldx, rwork ) / dble( m ) )
283 ELSE
284 eps2 = ulp
285 END IF
286 r = min( p, m-p, q, m-q )
287*
288* Copy the matrix X to the array XF.
289*
290 CALL zlacpy( 'Full', m, m, x, ldx, xf, ldx )
291*
292* Compute the CSD
293*
294 CALL zuncsd( 'Y', 'Y', 'Y', 'Y', 'N', 'D', m, p, q, xf(1,1), ldx,
295 $ xf(1,q+1), ldx, xf(p+1,1), ldx, xf(p+1,q+1), ldx,
296 $ theta, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t,
297 $ work, lwork, rwork, 17*(r+2), iwork, info )
298*
299* Compute XF := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
300*
301 CALL zlacpy( 'Full', m, m, x, ldx, xf, ldx )
302*
303 CALL zgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
304 $ xf, ldx, v1t, ldv1t, zero, work, ldx )
305*
306 CALL zgemm( 'Conjugate transpose', 'No transpose', p, q, p, one,
307 $ u1, ldu1, work, ldx, zero, xf, ldx )
308*
309 DO i = 1, min(p,q)-r
310 xf(i,i) = xf(i,i) - one
311 END DO
312 DO i = 1, r
313 xf(min(p,q)-r+i,min(p,q)-r+i) =
314 $ xf(min(p,q)-r+i,min(p,q)-r+i) - dcmplx( cos(theta(i)),
315 $ 0.0d0 )
316 END DO
317*
318 CALL zgemm( 'No transpose', 'Conjugate transpose', p, m-q, m-q,
319 $ one, xf(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
320*
321 CALL zgemm( 'Conjugate transpose', 'No transpose', p, m-q, p,
322 $ one, u1, ldu1, work, ldx, zero, xf(1,q+1), ldx )
323*
324 DO i = 1, min(p,m-q)-r
325 xf(p-i+1,m-i+1) = xf(p-i+1,m-i+1) + one
326 END DO
327 DO i = 1, r
328 xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) =
329 $ xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) +
330 $ dcmplx( sin(theta(r-i+1)), 0.0d0 )
331 END DO
332*
333 CALL zgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
334 $ xf(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
335*
336 CALL zgemm( 'Conjugate transpose', 'No transpose', m-p, q, m-p,
337 $ one, u2, ldu2, work, ldx, zero, xf(p+1,1), ldx )
338*
339 DO i = 1, min(m-p,q)-r
340 xf(m-i+1,q-i+1) = xf(m-i+1,q-i+1) - one
341 END DO
342 DO i = 1, r
343 xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
344 $ xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
345 $ dcmplx( sin(theta(r-i+1)), 0.0d0 )
346 END DO
347*
348 CALL zgemm( 'No transpose', 'Conjugate transpose', m-p, m-q, m-q,
349 $ one, xf(p+1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
350*
351 CALL zgemm( 'Conjugate transpose', 'No transpose', m-p, m-q, m-p,
352 $ one, u2, ldu2, work, ldx, zero, xf(p+1,q+1), ldx )
353*
354 DO i = 1, min(m-p,m-q)-r
355 xf(p+i,q+i) = xf(p+i,q+i) - one
356 END DO
357 DO i = 1, r
358 xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) =
359 $ xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) -
360 $ dcmplx( cos(theta(i)), 0.0d0 )
361 END DO
362*
363* Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
364*
365 resid = zlange( '1', p, q, xf, ldx, rwork )
366 result( 1 ) = ( resid / real(max(1,p,q)) ) / eps2
367*
368* Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
369*
370 resid = zlange( '1', p, m-q, xf(1,q+1), ldx, rwork )
371 result( 2 ) = ( resid / real(max(1,p,m-q)) ) / eps2
372*
373* Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
374*
375 resid = zlange( '1', m-p, q, xf(p+1,1), ldx, rwork )
376 result( 3 ) = ( resid / real(max(1,m-p,q)) ) / eps2
377*
378* Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
379*
380 resid = zlange( '1', m-p, m-q, xf(p+1,q+1), ldx, rwork )
381 result( 4 ) = ( resid / real(max(1,m-p,m-q)) ) / eps2
382*
383* Compute I - U1'*U1
384*
385 CALL zlaset( 'Full', p, p, zero, one, work, ldu1 )
386 CALL zherk( 'Upper', 'Conjugate transpose', p, p, -realone,
387 $ u1, ldu1, realone, work, ldu1 )
388*
389* Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
390*
391 resid = zlanhe( '1', 'Upper', p, work, ldu1, rwork )
392 result( 5 ) = ( resid / real(max(1,p)) ) / ulp
393*
394* Compute I - U2'*U2
395*
396 CALL zlaset( 'Full', m-p, m-p, zero, one, work, ldu2 )
397 CALL zherk( 'Upper', 'Conjugate transpose', m-p, m-p, -realone,
398 $ u2, ldu2, realone, work, ldu2 )
399*
400* Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
401*
402 resid = zlanhe( '1', 'Upper', m-p, work, ldu2, rwork )
403 result( 6 ) = ( resid / real(max(1,m-p)) ) / ulp
404*
405* Compute I - V1T*V1T'
406*
407 CALL zlaset( 'Full', q, q, zero, one, work, ldv1t )
408 CALL zherk( 'Upper', 'No transpose', q, q, -realone,
409 $ v1t, ldv1t, realone, work, ldv1t )
410*
411* Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
412*
413 resid = zlanhe( '1', 'Upper', q, work, ldv1t, rwork )
414 result( 7 ) = ( resid / real(max(1,q)) ) / ulp
415*
416* Compute I - V2T*V2T'
417*
418 CALL zlaset( 'Full', m-q, m-q, zero, one, work, ldv2t )
419 CALL zherk( 'Upper', 'No transpose', m-q, m-q, -realone,
420 $ v2t, ldv2t, realone, work, ldv2t )
421*
422* Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
423*
424 resid = zlanhe( '1', 'Upper', m-q, work, ldv2t, rwork )
425 result( 8 ) = ( resid / real(max(1,m-q)) ) / ulp
426*
427* Check sorting
428*
429 result( 9 ) = realzero
430 DO i = 1, r
431 IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
432 result( 9 ) = ulpinv
433 END IF
434 IF( i.GT.1) THEN
435 IF ( theta(i).LT.theta(i-1) ) THEN
436 result( 9 ) = ulpinv
437 END IF
438 END IF
439 END DO
440*
441* The second half of the routine checks the 2-by-1 CSD
442*
443 CALL zlaset( 'Full', q, q, zero, one, work, ldx )
444 CALL zherk( 'Upper', 'Conjugate transpose', q, m, -realone,
445 $ x, ldx, realone, work, ldx )
446 IF (m.GT.0) THEN
447 eps2 = max( ulp,
448 $ zlange( '1', Q, Q, WORK, LDX, RWORK ) / DBLE( M ) )
449 ELSE
450 EPS2 = ULP
451 END IF
452 R = MIN( P, M-P, Q, M-Q )
453*
454* Copy the matrix X to the array XF.
455*
456 CALL ZLACPY( 'full', M, M, X, LDX, XF, LDX )
457*
458* Compute the CSD
459*
460 CALL ZUNCSD2BY1( 'y', 'y', 'y', M, P, Q, XF(1,1), LDX, XF(P+1,1),
461 $ LDX, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK,
462 $ LWORK, RWORK, 17*(R+2), IWORK, INFO )
463*
464* Compute [X11;X21] := diag(U1,U2)'*[X11;X21]*V1 - [D11;D21]
465*
466 CALL ZGEMM( 'no transpose', 'conjugate transpose', P, Q, Q, ONE,
467 $ X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
468*
469 CALL ZGEMM( 'conjugate transpose', 'no transpose', P, Q, P, ONE,
470 $ U1, LDU1, WORK, LDX, ZERO, X, LDX )
471*
472 DO I = 1, MIN(P,Q)-R
473 X(I,I) = X(I,I) - ONE
474 END DO
475 DO I = 1, R
476 X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
477 $ X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - DCMPLX( COS(THETA(I)),
478 $ 0.0D0 )
479 END DO
480*
481 CALL ZGEMM( 'no transpose', 'conjugate transpose', M-P, Q, Q, ONE,
482 $ X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
483*
484 CALL ZGEMM( '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 $ DCMPLX( SIN(THETA(R-I+1)), 0.0D0 )
494 END DO
495*
496* Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
497*
498 RESID = ZLANGE( '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 = ZLANGE( '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 ZLASET( 'full', P, P, ZERO, ONE, WORK, LDU1 )
509 CALL ZHERK( 'upper', 'conjugate transpose', P, P, -REALONE,
510 $ U1, LDU1, REALONE, WORK, LDU1 )
511*
512* Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
513*
514 RESID = ZLANHE( '1', 'upper', P, WORK, LDU1, RWORK )
515 RESULT( 12 ) = ( RESID / REAL(MAX(1,P)) ) / ULP
516*
517* Compute I - U2'*U2
518*
519 CALL ZLASET( 'full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
520 CALL ZHERK( 'upper', 'conjugate transpose', M-P, M-P, -REALONE,
521 $ U2, LDU2, REALONE, WORK, LDU2 )
522*
523* Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
524*
525 RESID = ZLANHE( '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 ZLASET( 'full', Q, Q, ZERO, ONE, WORK, LDV1T )
531 CALL ZHERK( 'upper', 'no transpose', Q, Q, -REALONE,
532 $ V1T, LDV1T, REALONE, WORK, LDV1T )
533*
534* Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
535*
536 RESID = ZLANHE( '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.LT..OR..GT. IF( THETA(I)REALZERO THETA(I)PIOVER2 ) THEN
544 RESULT( 15 ) = ULPINV
545 END IF
546.GT. IF( I1) THEN
547.LT. IF ( THETA(I)THETA(I-1) ) THEN
548 RESULT( 15 ) = ULPINV
549 END IF
550 END IF
551 END DO
552*
553 RETURN
554*
555* End of ZCSDTS
556*
557 END
558
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
zlaset
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:106
zuncsd2by1
subroutine zuncsd2by1(jobu1, jobu2, jobv1t, m, p, q, x11, ldx11, x21, ldx21, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, work, lwork, rwork, lrwork, iwork, info)
ZUNCSD2BY1
Definition zuncsd2by1.f:256
zuncsd
recursive subroutine zuncsd(jobu1, jobu2, jobv1t, jobv2t, trans, signs, m, p, q, x11, ldx11, x12, ldx12, x21, ldx21, x22, ldx22, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t, work, lwork, rwork, lrwork, iwork, info)
ZUNCSD
Definition zuncsd.f:320
zherk
subroutine zherk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
ZHERK
Definition zherk.f:173
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:187
zcsdts
subroutine zcsdts(m, p, q, x, xf, ldx, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t, theta, iwork, work, lwork, rwork, result)
ZCSDTS
Definition zcsdts.f:229
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21