OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
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*> \ingroup single_eig
224*
225* =====================================================================
226 SUBROUTINE scsdts( 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 REAL RESULT( 15 ), RWORK( * ), THETA( * )
240 REAL U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
241 $ v2t( ldv2t, * ), work( lwork ), x( ldx, * ),
242 $ xf( ldx, * )
243* ..
244*
245* =====================================================================
246*
247* .. Parameters ..
248 REAL REALONE, REALZERO
249 PARAMETER ( REALONE = 1.0e0, realzero = 0.0e0 )
250 REAL ZERO, ONE
251 parameter( zero = 0.0e0, one = 1.0e0 )
252 REAL PIOVER2
253 parameter( piover2 = 1.57079632679489661923132169163975144210e0 )
254* ..
255* .. Local Scalars ..
256 INTEGER I, INFO, R
257 REAL EPS2, RESID, ULP, ULPINV
258* ..
259* .. External Functions ..
260 REAL SLAMCH, SLANGE, SLANSY
261 EXTERNAL SLAMCH, SLANGE, SLANSY
262* ..
263* .. External Subroutines ..
264 EXTERNAL sgemm, slacpy, slaset, sorcsd, sorcsd2by1,
265 $ ssyrk
266* ..
267* .. Intrinsic Functions ..
268 INTRINSIC cos, max, min, real, sin
269* ..
270* .. Executable Statements ..
271*
272 ulp = slamch( 'Precision' )
273 ulpinv = realone / ulp
274*
275* The first half of the routine checks the 2-by-2 CSD
276*
277 CALL slaset( 'Full', m, m, zero, one, work, ldx )
278 CALL ssyrk( 'Upper', 'Conjugate transpose', m, m, -one, x, ldx,
279 $ one, work, ldx )
280 IF (m.GT.0) THEN
281 eps2 = max( ulp,
282 $ slange( '1', m, m, work, ldx, rwork ) / real( 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 slacpy( 'Full', m, m, x, ldx, xf, ldx )
291*
292* Compute the CSD
293*
294 CALL sorcsd( '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, iwork, info )
298*
299* Compute XF := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
300*
301 CALL slacpy( 'Full', m, m, x, ldx, xf, ldx )
302*
303 CALL sgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
304 $ xf, ldx, v1t, ldv1t, zero, work, ldx )
305*
306 CALL sgemm( '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) - cos(theta(i))
315 END DO
316*
317 CALL sgemm( 'no transpose', 'conjugate transpose', P, M-Q, M-Q,
318 $ ONE, XF(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
319*
320 CALL SGEMM( 'conjugate transpose', 'no transpose', P, M-Q, P,
321 $ ONE, U1, LDU1, WORK, LDX, ZERO, XF(1,Q+1), LDX )
322*
323 DO I = 1, MIN(P,M-Q)-R
324 XF(P-I+1,M-I+1) = XF(P-I+1,M-I+1) + ONE
325 END DO
326 DO I = 1, R
327 XF(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
328 $ XF(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
329 $ SIN(THETA(R-I+1))
330 END DO
331*
332 CALL SGEMM( 'no transpose', 'conjugate transpose', M-P, Q, Q, ONE,
333 $ XF(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
334*
335 CALL SGEMM( 'conjugate transpose', 'no transpose', M-P, Q, M-P,
336 $ ONE, U2, LDU2, WORK, LDX, ZERO, XF(P+1,1), LDX )
337*
338 DO I = 1, MIN(M-P,Q)-R
339 XF(M-I+1,Q-I+1) = XF(M-I+1,Q-I+1) - ONE
340 END DO
341 DO I = 1, R
342 XF(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
343 $ XF(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
344 $ SIN(THETA(R-I+1))
345 END DO
346*
347 CALL SGEMM( 'no transpose', 'conjugate transpose', M-P, M-Q, M-Q,
348 $ ONE, XF(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
349*
350 CALL SGEMM( 'conjugate transpose', 'no transpose', M-P, M-Q, M-P,
351 $ ONE, U2, LDU2, WORK, LDX, ZERO, XF(P+1,Q+1), LDX )
352*
353 DO I = 1, MIN(M-P,M-Q)-R
354 XF(P+I,Q+I) = XF(P+I,Q+I) - ONE
355 END DO
356 DO I = 1, R
357 XF(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
358 $ XF(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
359 $ COS(THETA(I))
360 END DO
361*
362* Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
363*
364 RESID = SLANGE( '1', p, q, xf, ldx, rwork )
365 result( 1 ) = ( resid / real(max(1,p,q)) ) / eps2
366*
367* Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
368*
369 resid = slange( '1', p, m-q, xf(1,q+1), ldx, rwork )
370 result( 2 ) = ( resid / real(max(1,p,m-q)) ) / eps2
371*
372* Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
373*
374 resid = slange( '1', m-p, q, xf(p+1,1), ldx, rwork )
375 result( 3 ) = ( resid / real(max(1,m-p,q)) ) / eps2
376*
377* Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
378*
379 resid = slange( '1', m-p, m-q, xf(p+1,q+1), ldx, rwork )
380 result( 4 ) = ( resid / real(max(1,m-p,m-q)) ) / eps2
381*
382* Compute I - U1'*U1
383*
384 CALL slaset( 'Full', p, p, zero, one, work, ldu1 )
385 CALL ssyrk( 'Upper', 'Conjugate transpose', p, p, -one, u1, ldu1,
386 $ one, work, ldu1 )
387*
388* Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
389*
390 resid = slansy( '1', 'Upper', p, work, ldu1, rwork )
391 result( 5 ) = ( resid / real(max(1,p)) ) / ulp
392*
393* Compute I - U2'*U2
394*
395 CALL slaset( 'Full', m-p, m-p, zero, one, work, ldu2 )
396 CALL ssyrk( 'Upper', 'Conjugate transpose', m-p, m-p, -one, u2,
397 $ ldu2, one, work, ldu2 )
398*
399* Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
400*
401 resid = slansy( '1', 'Upper', m-p, work, ldu2, rwork )
402 result( 6 ) = ( resid / real(max(1,m-p)) ) / ulp
403*
404* Compute I - V1T*V1T'
405*
406 CALL slaset( 'Full', q, q, zero, one, work, ldv1t )
407 CALL ssyrk( 'Upper', 'No transpose', q, q, -one, v1t, ldv1t, one,
408 $ work, ldv1t )
409*
410* Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
411*
412 resid = slansy( '1', 'Upper', q, work, ldv1t, rwork )
413 result( 7 ) = ( resid / real(max(1,q)) ) / ulp
414*
415* Compute I - V2T*V2T'
416*
417 CALL slaset( 'Full', m-q, m-q, zero, one, work, ldv2t )
418 CALL ssyrk( 'Upper', 'No transpose', m-q, m-q, -one, v2t, ldv2t,
419 $ one, work, ldv2t )
420*
421* Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
422*
423 resid = slansy( '1', 'Upper', m-q, work, ldv2t, rwork )
424 result( 8 ) = ( resid / real(max(1,m-q)) ) / ulp
425*
426* Check sorting
427*
428 result( 9 ) = realzero
429 DO i = 1, r
430 IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
431 result( 9 ) = ulpinv
432 END IF
433 IF( i.GT.1 ) THEN
434 IF ( theta(i).LT.theta(i-1) ) THEN
435 result( 9 ) = ulpinv
436 END IF
437 END IF
438 END DO
439*
440* The second half of the routine checks the 2-by-1 CSD
441*
442 CALL slaset( 'Full', q, q, zero, one, work, ldx )
443 CALL ssyrk( 'Upper', 'Conjugate transpose', q, m, -one, x, ldx,
444 $ one, work, ldx )
445 IF (m.GT.0) THEN
446 eps2 = max( ulp,
447 $ slange( '1', q, q, work, ldx, rwork ) / real( m ) )
448 ELSE
449 eps2 = ulp
450 END IF
451 r = min( p, m-p, q, m-q )
452*
453* Copy the matrix [X11;X21] to the array XF.
454*
455 CALL slacpy( 'Full', m, q, x, ldx, xf, ldx )
456*
457* Compute the CSD
458*
459 CALL sorcsd2by1( 'Y', 'Y', 'Y', m, p, q, xf(1,1), ldx, xf(p+1,1),
460 $ ldx, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, work,
461 $ lwork, iwork, info )
462*
463* Compute [X11;X21] := diag(U1,U2)'*[X11;X21]*V1 - [D11;D21]
464*
465 CALL sgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
466 $ x, ldx, v1t, ldv1t, zero, work, ldx )
467*
468 CALL sgemm( 'Conjugate transpose', 'No transpose', p, q, p, one,
469 $ u1, ldu1, work, ldx, zero, x, ldx )
470*
471 DO i = 1, min(p,q)-r
472 x(i,i) = x(i,i) - one
473 END DO
474 DO i = 1, r
475 x(min(p,q)-r+i,min(p,q)-r+i) =
476 $ x(min(p,q)-r+i,min(p,q)-r+i) - cos(theta(i))
477 END DO
478*
479 CALL sgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
480 $ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
481*
482 CALL sgemm( 'Conjugate transpose', 'No transpose', m-p, q, m-p,
483 $ one, u2, ldu2, work, ldx, zero, x(p+1,1), ldx )
484*
485 DO i = 1, min(m-p,q)-r
486 x(m-i+1,q-i+1) = x(m-i+1,q-i+1) - one
487 END DO
488 DO i = 1, r
489 x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
490 $ x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
491 $ sin(theta(r-i+1))
492 END DO
493*
494* Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
495*
496 resid = slange( '1', p, q, x, ldx, rwork )
497 result( 10 ) = ( resid / real(max(1,p,q)) ) / eps2
498*
499* Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
500*
501 resid = slange( '1', m-p, q, x(p+1,1), ldx, rwork )
502 result( 11 ) = ( resid / real(max(1,m-p,q)) ) / eps2
503*
504* Compute I - U1'*U1
505*
506 CALL slaset( 'Full', p, p, zero, one, work, ldu1 )
507 CALL ssyrk( 'Upper', 'Conjugate transpose', p, p, -one, u1, ldu1,
508 $ one, work, ldu1 )
509*
510* Compute norm( I - U1'*U1 ) / ( MAX(1,P) * ULP ) .
511*
512 resid = slansy( '1', 'Upper', p, work, ldu1, rwork )
513 result( 12 ) = ( resid / real(max(1,p)) ) / ulp
514*
515* Compute I - U2'*U2
516*
517 CALL slaset( 'Full', m-p, m-p, zero, one, work, ldu2 )
518 CALL ssyrk( 'Upper', 'Conjugate transpose', m-p, m-p, -one, u2,
519 $ ldu2, one, work, ldu2 )
520*
521* Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
522*
523 resid = slansy( '1', 'Upper', m-p, work, ldu2, rwork )
524 result( 13 ) = ( resid / real(max(1,m-p)) ) / ulp
525*
526* Compute I - V1T*V1T'
527*
528 CALL slaset( 'Full', q, q, zero, one, work, ldv1t )
529 CALL ssyrk( 'Upper', 'No transpose', q, q, -one, v1t, ldv1t, one,
530 $ work, ldv1t )
531*
532* Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
533*
534 resid = slansy( '1', 'Upper', q, work, ldv1t, rwork )
535 result( 14 ) = ( resid / real(max(1,q)) ) / ulp
536*
537* Check sorting
538*
539 result( 15 ) = realzero
540 DO i = 1, r
541 IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
542 result( 15 ) = ulpinv
543 END IF
544 IF( i.GT.1 ) THEN
545 IF ( theta(i).LT.theta(i-1) ) THEN
546 result( 15 ) = ulpinv
547 END IF
548 END IF
549 END DO
550*
551 RETURN
552*
553* End of SCSDTS
554*
555 END
556
slaset
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
slacpy
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
sorcsd2by1
subroutine sorcsd2by1(jobu1, jobu2, jobv1t, m, p, q, x11, ldx11, x21, ldx21, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, work, lwork, iwork, info)
SORCSD2BY1
Definition sorcsd2by1.f:233
sorcsd
recursive subroutine sorcsd(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, iwork, info)
SORCSD
Definition sorcsd.f:300
ssyrk
subroutine ssyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
SSYRK
Definition ssyrk.f:169
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:187
scsdts
subroutine scsdts(m, p, q, x, xf, ldx, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t, theta, iwork, work, lwork, rwork, result)
SCSDTS
Definition scsdts.f:229
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21