LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
dckcsd.f
Go to the documentation of this file.
1*> \brief \b DCKCSD
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 DCKCSD( NM, MVAL, PVAL, QVAL, NMATS, ISEED, THRESH,
12* MMAX, X, XF, U1, U2, V1T, V2T, THETA, IWORK,
13* WORK, RWORK, NIN, NOUT, INFO )
14*
15* .. Scalar Arguments ..
16* INTEGER INFO, NIN, NM, NMATS, MMAX, NOUT
17* DOUBLE PRECISION THRESH
18* ..
19* .. Array Arguments ..
20* INTEGER ISEED( 4 ), IWORK( * ), MVAL( * ), PVAL( * ),
21* \$ QVAL( * )
22* DOUBLE PRECISION RWORK( * ), THETA( * )
23* DOUBLE PRECISION U1( * ), U2( * ), V1T( * ), V2T( * ),
24* \$ WORK( * ), X( * ), XF( * )
25* ..
26*
27*
28*> \par Purpose:
29* =============
30*>
31*> \verbatim
32*>
33*> DCKCSD tests DORCSD:
34*> the CSD for an M-by-M orthogonal matrix X partitioned as
35*> [ X11 X12; X21 X22 ]. X11 is P-by-Q.
36*> \endverbatim
37*
38* Arguments:
39* ==========
40*
41*> \param[in] NM
42*> \verbatim
43*> NM is INTEGER
44*> The number of values of M contained in the vector MVAL.
45*> \endverbatim
46*>
47*> \param[in] MVAL
48*> \verbatim
49*> MVAL is INTEGER array, dimension (NM)
50*> The values of the matrix row dimension M.
51*> \endverbatim
52*>
53*> \param[in] PVAL
54*> \verbatim
55*> PVAL is INTEGER array, dimension (NM)
56*> The values of the matrix row dimension P.
57*> \endverbatim
58*>
59*> \param[in] QVAL
60*> \verbatim
61*> QVAL is INTEGER array, dimension (NM)
62*> The values of the matrix column dimension Q.
63*> \endverbatim
64*>
65*> \param[in] NMATS
66*> \verbatim
67*> NMATS is INTEGER
68*> The number of matrix types to be tested for each combination
69*> of matrix dimensions. If NMATS >= NTYPES (the maximum
70*> number of matrix types), then all the different types are
71*> generated for testing. If NMATS < NTYPES, another input line
72*> is read to get the numbers of the matrix types to be used.
73*> \endverbatim
74*>
75*> \param[in,out] ISEED
76*> \verbatim
77*> ISEED is INTEGER array, dimension (4)
78*> On entry, the seed of the random number generator. The array
79*> elements should be between 0 and 4095, otherwise they will be
80*> reduced mod 4096, and ISEED(4) must be odd.
81*> On exit, the next seed in the random number sequence after
82*> all the test matrices have been generated.
83*> \endverbatim
84*>
85*> \param[in] THRESH
86*> \verbatim
87*> THRESH is DOUBLE PRECISION
88*> The threshold value for the test ratios. A result is
89*> included in the output file if RESULT >= THRESH. To have
90*> every test ratio printed, use THRESH = 0.
91*> \endverbatim
92*>
93*> \param[in] MMAX
94*> \verbatim
95*> MMAX is INTEGER
96*> The maximum value permitted for M, used in dimensioning the
97*> work arrays.
98*> \endverbatim
99*>
100*> \param[out] X
101*> \verbatim
102*> X is DOUBLE PRECISION array, dimension (MMAX*MMAX)
103*> \endverbatim
104*>
105*> \param[out] XF
106*> \verbatim
107*> XF is DOUBLE PRECISION array, dimension (MMAX*MMAX)
108*> \endverbatim
109*>
110*> \param[out] U1
111*> \verbatim
112*> U1 is DOUBLE PRECISION array, dimension (MMAX*MMAX)
113*> \endverbatim
114*>
115*> \param[out] U2
116*> \verbatim
117*> U2 is DOUBLE PRECISION array, dimension (MMAX*MMAX)
118*> \endverbatim
119*>
120*> \param[out] V1T
121*> \verbatim
122*> V1T is DOUBLE PRECISION array, dimension (MMAX*MMAX)
123*> \endverbatim
124*>
125*> \param[out] V2T
126*> \verbatim
127*> V2T is DOUBLE PRECISION array, dimension (MMAX*MMAX)
128*> \endverbatim
129*>
130*> \param[out] THETA
131*> \verbatim
132*> THETA is DOUBLE PRECISION array, dimension (MMAX)
133*> \endverbatim
134*>
135*> \param[out] IWORK
136*> \verbatim
137*> IWORK is INTEGER array, dimension (MMAX)
138*> \endverbatim
139*>
140*> \param[out] WORK
141*> \verbatim
142*> WORK is DOUBLE PRECISION array
143*> \endverbatim
144*>
145*> \param[out] RWORK
146*> \verbatim
147*> RWORK is DOUBLE PRECISION array
148*> \endverbatim
149*>
150*> \param[in] NIN
151*> \verbatim
152*> NIN is INTEGER
153*> The unit number for input.
154*> \endverbatim
155*>
156*> \param[in] NOUT
157*> \verbatim
158*> NOUT is INTEGER
159*> The unit number for output.
160*> \endverbatim
161*>
162*> \param[out] INFO
163*> \verbatim
164*> INFO is INTEGER
165*> = 0 : successful exit
166*> > 0 : If DLAROR returns an error code, the absolute value
167*> of it is returned.
168*> \endverbatim
169*
170* Authors:
171* ========
172*
173*> \author Univ. of Tennessee
174*> \author Univ. of California Berkeley
175*> \author Univ. of Colorado Denver
176*> \author NAG Ltd.
177*
178*> \ingroup double_eig
179*
180* =====================================================================
181 SUBROUTINE dckcsd( NM, MVAL, PVAL, QVAL, NMATS, ISEED, THRESH,
182 \$ MMAX, X, XF, U1, U2, V1T, V2T, THETA, IWORK,
183 \$ WORK, RWORK, NIN, NOUT, INFO )
184*
185* -- LAPACK test routine --
186* -- LAPACK is a software package provided by Univ. of Tennessee, --
187* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
188*
189* .. Scalar Arguments ..
190 INTEGER INFO, NIN, NM, NMATS, MMAX, NOUT
191 DOUBLE PRECISION THRESH
192* ..
193* .. Array Arguments ..
194 INTEGER ISEED( 4 ), IWORK( * ), MVAL( * ), PVAL( * ),
195 \$ QVAL( * )
196 DOUBLE PRECISION RWORK( * ), THETA( * )
197 DOUBLE PRECISION U1( * ), U2( * ), V1T( * ), V2T( * ),
198 \$ work( * ), x( * ), xf( * )
199* ..
200*
201* =====================================================================
202*
203* .. Parameters ..
204 INTEGER NTESTS
205 PARAMETER ( NTESTS = 15 )
206 INTEGER NTYPES
207 parameter( ntypes = 4 )
208 DOUBLE PRECISION GAPDIGIT, ONE, ORTH, TEN, ZERO
209 parameter( gapdigit = 18.0d0, one = 1.0d0,
210 \$ orth = 1.0d-12,
211 \$ ten = 10.0d0, zero = 0.0d0 )
212 DOUBLE PRECISION PIOVER2
213 PARAMETER ( PIOVER2 = 1.57079632679489661923132169163975144210d0 )
214* ..
215* .. Local Scalars ..
216 LOGICAL FIRSTT
217 CHARACTER*3 PATH
218 INTEGER I, IINFO, IM, IMAT, J, LDU1, LDU2, LDV1T,
219 \$ ldv2t, ldx, lwork, m, nfail, nrun, nt, p, q, r
220* ..
221* .. Local Arrays ..
222 LOGICAL DOTYPE( NTYPES )
223 DOUBLE PRECISION RESULT( NTESTS )
224* ..
225* .. External Subroutines ..
226 EXTERNAL alahdg, alareq, alasum, dcsdts, dlacsg, dlaror,
227 \$ dlaset, drot
228* ..
229* .. Intrinsic Functions ..
230 INTRINSIC abs, min
231* ..
232* .. External Functions ..
233 DOUBLE PRECISION DLARAN, DLARND
234 EXTERNAL DLARAN, DLARND
235* ..
236* .. Executable Statements ..
237*
238* Initialize constants and the random number seed.
239*
240 path( 1: 3 ) = 'CSD'
241 info = 0
242 nrun = 0
243 nfail = 0
244 firstt = .true.
245 CALL alareq( path, nmats, dotype, ntypes, nin, nout )
246 ldx = mmax
247 ldu1 = mmax
248 ldu2 = mmax
249 ldv1t = mmax
250 ldv2t = mmax
251 lwork = mmax*mmax
252*
253* Do for each value of M in MVAL.
254*
255 DO 30 im = 1, nm
256 m = mval( im )
257 p = pval( im )
258 q = qval( im )
259*
260 DO 20 imat = 1, ntypes
261*
262* Do the tests only if DOTYPE( IMAT ) is true.
263*
264 IF( .NOT.dotype( imat ) )
265 \$ GO TO 20
266*
267* Generate X
268*
269 IF( imat.EQ.1 ) THEN
270 CALL dlaror( 'L', 'I', m, m, x, ldx, iseed, work, iinfo )
271 IF( m .NE. 0 .AND. iinfo .NE. 0 ) THEN
272 WRITE( nout, fmt = 9999 ) m, iinfo
273 info = abs( iinfo )
274 GO TO 20
275 END IF
276 ELSE IF( imat.EQ.2 ) THEN
277 r = min( p, m-p, q, m-q )
278 DO i = 1, r
279 theta(i) = piover2 * dlarnd( 1, iseed )
280 END DO
281 CALL dlacsg( m, p, q, theta, iseed, x, ldx, work )
282 DO i = 1, m
283 DO j = 1, m
284 x(i+(j-1)*ldx) = x(i+(j-1)*ldx) +
285 \$ orth*dlarnd(2,iseed)
286 END DO
287 END DO
288 ELSE IF( imat.EQ.3 ) THEN
289 r = min( p, m-p, q, m-q )
290 DO i = 1, r+1
291 theta(i) = ten**(-dlarnd(1,iseed)*gapdigit)
292 END DO
293 DO i = 2, r+1
294 theta(i) = theta(i-1) + theta(i)
295 END DO
296 DO i = 1, r
297 theta(i) = piover2 * theta(i) / theta(r+1)
298 END DO
299 CALL dlacsg( m, p, q, theta, iseed, x, ldx, work )
300 ELSE
301 CALL dlaset( 'F', m, m, zero, one, x, ldx )
302 DO i = 1, m
303 j = int( dlaran( iseed ) * m ) + 1
304 IF( j .NE. i ) THEN
305 CALL drot( m, x(1+(i-1)*ldx), 1, x(1+(j-1)*ldx), 1,
306 \$ zero, one )
307 END IF
308 END DO
309 END IF
310*
311 nt = 15
312*
313 CALL dcsdts( m, p, q, x, xf, ldx, u1, ldu1, u2, ldu2, v1t,
314 \$ ldv1t, v2t, ldv2t, theta, iwork, work, lwork,
315 \$ rwork, result )
316*
317* Print information about the tests that did not
318* pass the threshold.
319*
320 DO 10 i = 1, nt
321 IF( result( i ).GE.thresh ) THEN
322 IF( nfail.EQ.0 .AND. firstt ) THEN
323 firstt = .false.
324 CALL alahdg( nout, path )
325 END IF
326 WRITE( nout, fmt = 9998 )m, p, q, imat, i,
327 \$ result( i )
328 nfail = nfail + 1
329 END IF
330 10 CONTINUE
331 nrun = nrun + nt
332 20 CONTINUE
333 30 CONTINUE
334*
335* Print a summary of the results.
336*
337 CALL alasum( path, nout, nfail, nrun, 0 )
338*
339 9999 FORMAT( ' DLAROR in DCKCSD: M = ', i5, ', INFO = ', i15 )
340 9998 FORMAT( ' M=', i4, ' P=', i4, ', Q=', i4, ', type ', i2,
341 \$ ', test ', i2, ', ratio=', g13.6 )
342 RETURN
343*
344* End of DCKCSD
345*
346 END
347*
348*
349*
350 SUBROUTINE dlacsg( M, P, Q, THETA, ISEED, X, LDX, WORK )
351 IMPLICIT NONE
352*
353 INTEGER LDX, M, P, Q
354 INTEGER ISEED( 4 )
355 DOUBLE PRECISION THETA( * )
356 DOUBLE PRECISION WORK( * ), X( LDX, * )
357*
358 DOUBLE PRECISION ONE, ZERO
359 PARAMETER ( ONE = 1.0d0, zero = 0.0d0 )
360*
361 INTEGER I, INFO, R
362*
363 r = min( p, m-p, q, m-q )
364*
365 CALL dlaset( 'Full', m, m, zero, zero, x, ldx )
366*
367 DO i = 1, min(p,q)-r
368 x(i,i) = one
369 END DO
370 DO i = 1, r
371 x(min(p,q)-r+i,min(p,q)-r+i) = cos(theta(i))
372 END DO
373 DO i = 1, min(p,m-q)-r
374 x(p-i+1,m-i+1) = -one
375 END DO
376 DO i = 1, r
377 x(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) =
378 \$ -sin(theta(r-i+1))
379 END DO
380 DO i = 1, min(m-p,q)-r
381 x(m-i+1,q-i+1) = one
382 END DO
383 DO i = 1, r
384 x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
385 \$ sin(theta(r-i+1))
386 END DO
387 DO i = 1, min(m-p,m-q)-r
388 x(p+i,q+i) = one
389 END DO
390 DO i = 1, r
391 x(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) =
392 \$ cos(theta(i))
393 END DO
394 CALL dlaror( 'Left', 'No init', p, m, x, ldx, iseed, work, info )
395 CALL dlaror( 'Left', 'No init', m-p, m, x(p+1,1), ldx,
396 \$ iseed, work, info )
397 CALL dlaror( 'Right', 'No init', m, q, x, ldx, iseed,
398 \$ work, info )
399 CALL dlaror( 'Right', 'No init', m, m-q,
400 \$ x(1,q+1), ldx, iseed, work, info )
401*
402 END
403
subroutine dlacsg(M, P, Q, THETA, ISEED, X, LDX, WORK)
Definition: dckcsd.f:351
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine alareq(PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT)
ALAREQ
Definition: alareq.f:90
subroutine alahdg(IOUNIT, PATH)
ALAHDG
Definition: alahdg.f:62
subroutine alasum(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASUM
Definition: alasum.f:73
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:92
subroutine dcsdts(M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK, RWORK, RESULT)
DCSDTS
Definition: dcsdts.f:229
subroutine dckcsd(NM, MVAL, PVAL, QVAL, NMATS, ISEED, THRESH, MMAX, X, XF, U1, U2, V1T, V2T, THETA, IWORK, WORK, RWORK, NIN, NOUT, INFO)
DCKCSD
Definition: dckcsd.f:184
subroutine dlaror(SIDE, INIT, M, N, A, LDA, ISEED, X, INFO)
DLAROR
Definition: dlaror.f:146