LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 *> \date November 2011
179 *
180 *> \ingroup double_eig
181 *
182 * =====================================================================
183  SUBROUTINE dckcsd( NM, MVAL, PVAL, QVAL, NMATS, ISEED, THRESH,
184  $ mmax, x, xf, u1, u2, v1t, v2t, theta, iwork,
185  $ work, rwork, nin, nout, info )
186 *
187 * -- LAPACK test routine (version 3.4.0) --
188 * -- LAPACK is a software package provided by Univ. of Tennessee, --
189 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
190 * November 2011
191 *
192 * .. Scalar Arguments ..
193  INTEGER INFO, NIN, NM, NMATS, MMAX, NOUT
194  DOUBLE PRECISION THRESH
195 * ..
196 * .. Array Arguments ..
197  INTEGER ISEED( 4 ), IWORK( * ), MVAL( * ), PVAL( * ),
198  $ qval( * )
199  DOUBLE PRECISION RWORK( * ), THETA( * )
200  DOUBLE PRECISION U1( * ), U2( * ), V1T( * ), V2T( * ),
201  $ work( * ), x( * ), xf( * )
202 * ..
203 *
204 * =====================================================================
205 *
206 * .. Parameters ..
207  INTEGER NTESTS
208  parameter ( ntests = 15 )
209  INTEGER NTYPES
210  parameter ( ntypes = 4 )
211  DOUBLE PRECISION GAPDIGIT, ONE, ORTH, PIOVER2, TEN, ZERO
212  parameter ( gapdigit = 18.0d0, one = 1.0d0,
213  $ orth = 1.0d-12,
214  $ piover2 = 1.57079632679489662d0,
215  $ ten = 10.0d0, zero = 0.0d0 )
216 * ..
217 * .. Local Scalars ..
218  LOGICAL FIRSTT
219  CHARACTER*3 PATH
220  INTEGER I, IINFO, IM, IMAT, J, LDU1, LDU2, LDV1T,
221  $ ldv2t, ldx, lwork, m, nfail, nrun, nt, p, q, r
222 * ..
223 * .. Local Arrays ..
224  LOGICAL DOTYPE( ntypes )
225  DOUBLE PRECISION RESULT( ntests )
226 * ..
227 * .. External Subroutines ..
228  EXTERNAL alahdg, alareq, alasum, dcsdts, dlacsg, dlaror,
229  $ dlaset
230 * ..
231 * .. Intrinsic Functions ..
232  INTRINSIC abs, min
233 * ..
234 * .. External Functions ..
235  DOUBLE PRECISION DLARAN, DLARND
236  EXTERNAL dlaran, dlarnd
237 * ..
238 * .. Executable Statements ..
239 *
240 * Initialize constants and the random number seed.
241 *
242  path( 1: 3 ) = 'CSD'
243  info = 0
244  nrun = 0
245  nfail = 0
246  firstt = .true.
247  CALL alareq( path, nmats, dotype, ntypes, nin, nout )
248  ldx = mmax
249  ldu1 = mmax
250  ldu2 = mmax
251  ldv1t = mmax
252  ldv2t = mmax
253  lwork = mmax*mmax
254 *
255 * Do for each value of M in MVAL.
256 *
257  DO 30 im = 1, nm
258  m = mval( im )
259  p = pval( im )
260  q = qval( im )
261 *
262  DO 20 imat = 1, ntypes
263 *
264 * Do the tests only if DOTYPE( IMAT ) is true.
265 *
266  IF( .NOT.dotype( imat ) )
267  $ GO TO 20
268 *
269 * Generate X
270 *
271  IF( imat.EQ.1 ) THEN
272  CALL dlaror( 'L', 'I', m, m, x, ldx, iseed, work, iinfo )
273  IF( m .NE. 0 .AND. iinfo .NE. 0 ) THEN
274  WRITE( nout, fmt = 9999 ) m, iinfo
275  info = abs( iinfo )
276  GO TO 20
277  END IF
278  ELSE IF( imat.EQ.2 ) THEN
279  r = min( p, m-p, q, m-q )
280  DO i = 1, r
281  theta(i) = piover2 * dlarnd( 1, iseed )
282  END DO
283  CALL dlacsg( m, p, q, theta, iseed, x, ldx, work )
284  DO i = 1, m
285  DO j = 1, m
286  x(i+(j-1)*ldx) = x(i+(j-1)*ldx) +
287  $ orth*dlarnd(2,iseed)
288  END DO
289  END DO
290  ELSE IF( imat.EQ.3 ) THEN
291  r = min( p, m-p, q, m-q )
292  DO i = 1, r+1
293  theta(i) = ten**(-dlarnd(1,iseed)*gapdigit)
294  END DO
295  DO i = 2, r+1
296  theta(i) = theta(i-1) + theta(i)
297  END DO
298  DO i = 1, r
299  theta(i) = piover2 * theta(i) / theta(r+1)
300  END DO
301  CALL dlacsg( m, p, q, theta, iseed, x, ldx, work )
302  ELSE
303  CALL dlaset( 'F', m, m, zero, one, x, ldx )
304  DO i = 1, m
305  j = int( dlaran( iseed ) * m ) + 1
306  IF( j .NE. i ) THEN
307  CALL drot( m, x(1+(i-1)*ldx), 1, x(1+(j-1)*ldx), 1,
308  $ zero, one )
309  END IF
310  END DO
311  END IF
312 *
313  nt = 15
314 *
315  CALL dcsdts( m, p, q, x, xf, ldx, u1, ldu1, u2, ldu2, v1t,
316  $ ldv1t, v2t, ldv2t, theta, iwork, work, lwork,
317  $ rwork, result )
318 *
319 * Print information about the tests that did not
320 * pass the threshold.
321 *
322  DO 10 i = 1, nt
323  IF( result( i ).GE.thresh ) THEN
324  IF( nfail.EQ.0 .AND. firstt ) THEN
325  firstt = .false.
326  CALL alahdg( nout, path )
327  END IF
328  WRITE( nout, fmt = 9998 )m, p, q, imat, i,
329  $ result( i )
330  nfail = nfail + 1
331  END IF
332  10 CONTINUE
333  nrun = nrun + nt
334  20 CONTINUE
335  30 CONTINUE
336 *
337 * Print a summary of the results.
338 *
339  CALL alasum( path, nout, nfail, nrun, 0 )
340 *
341  9999 FORMAT( ' DLAROR in DCKCSD: M = ', i5, ', INFO = ', i15 )
342  9998 FORMAT( ' M=', i4, ' P=', i4, ', Q=', i4, ', type ', i2,
343  $ ', test ', i2, ', ratio=', g13.6 )
344  RETURN
345 *
346 * End of DCKCSD
347 *
348  END
349 *
350 *
351 *
352  SUBROUTINE dlacsg( M, P, Q, THETA, ISEED, X, LDX, WORK )
353  IMPLICIT NONE
354 *
355  INTEGER LDX, M, P, Q
356  INTEGER ISEED( 4 )
357  DOUBLE PRECISION THETA( * )
358  DOUBLE PRECISION WORK( * ), X( ldx, * )
359 *
360  DOUBLE PRECISION ONE, ZERO
361  parameter ( one = 1.0d0, zero = 0.0d0 )
362 *
363  INTEGER I, INFO, R
364 *
365  r = min( p, m-p, q, m-q )
366 *
367  CALL dlaset( 'Full', m, m, zero, zero, x, ldx )
368 *
369  DO i = 1, min(p,q)-r
370  x(i,i) = one
371  END DO
372  DO i = 1, r
373  x(min(p,q)-r+i,min(p,q)-r+i) = cos(theta(i))
374  END DO
375  DO i = 1, min(p,m-q)-r
376  x(p-i+1,m-i+1) = -one
377  END DO
378  DO i = 1, r
379  x(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) =
380  $ -sin(theta(r-i+1))
381  END DO
382  DO i = 1, min(m-p,q)-r
383  x(m-i+1,q-i+1) = one
384  END DO
385  DO i = 1, r
386  x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
387  $ sin(theta(r-i+1))
388  END DO
389  DO i = 1, min(m-p,m-q)-r
390  x(p+i,q+i) = one
391  END DO
392  DO i = 1, r
393  x(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) =
394  $ cos(theta(i))
395  END DO
396  CALL dlaror( 'Left', 'No init', p, m, x, ldx, iseed, work, info )
397  CALL dlaror( 'Left', 'No init', m-p, m, x(p+1,1), ldx,
398  $ iseed, work, info )
399  CALL dlaror( 'Right', 'No init', m, q, x, ldx, iseed,
400  $ work, info )
401  CALL dlaror( 'Right', 'No init', m, m-q,
402  $ x(1,q+1), ldx, iseed, work, info )
403 *
404  END
405 
dcsdts
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:231
dlaset
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:112
drot
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:53
dlacsg
subroutine dlacsg(M, P, Q, THETA, ISEED, X, LDX, WORK)
Definition: dckcsd.f:353
dlaror
subroutine dlaror(SIDE, INIT, M, N, A, LDA, ISEED, X, INFO)
DLAROR
Definition: dlaror.f:148
alareq
subroutine alareq(PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT)
ALAREQ
Definition: alareq.f:92
dckcsd
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:186
alahdg
subroutine alahdg(IOUNIT, PATH)
ALAHDG
Definition: alahdg.f:64
alasum
subroutine alasum(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASUM
Definition: alasum.f:75