LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 = 9 )
209  INTEGER ntypes
210  parameter( ntypes = 3 )
211  DOUBLE PRECISION gapdigit, orth, piover2, ten
212  parameter( gapdigit = 18.0d0, orth = 1.0d-12,
213  $ piover2 = 1.57079632679489662d0,
214  $ ten = 10.0d0 )
215 * ..
216 * .. Local Scalars ..
217  LOGICAL firstt
218  CHARACTER*3 path
219  INTEGER i, iinfo, im, imat, j, ldu1, ldu2, ldv1t,
220  $ ldv2t, ldx, lwork, m, nfail, nrun, nt, p, q, r
221 * ..
222 * .. Local Arrays ..
223  LOGICAL dotype( ntypes )
224  DOUBLE PRECISION result( ntests )
225 * ..
226 * .. External Subroutines ..
227  EXTERNAL alahdg, alareq, alasum, dcsdts, dlacsg, dlaror,
228  $ dlaset
229 * ..
230 * .. Intrinsic Functions ..
231  INTRINSIC abs, min
232 * ..
233 * .. External Functions ..
234  DOUBLE PRECISION dlarnd
235  EXTERNAL dlarnd
236 * ..
237 * .. Executable Statements ..
238 *
239 * Initialize constants and the random number seed.
240 *
241  path( 1: 3 ) = 'CSD'
242  info = 0
243  nrun = 0
244  nfail = 0
245  firstt = .true.
246  CALL alareq( path, nmats, dotype, ntypes, nin, nout )
247  ldx = mmax
248  ldu1 = mmax
249  ldu2 = mmax
250  ldv1t = mmax
251  ldv2t = mmax
252  lwork = mmax*mmax
253 *
254 * Do for each value of M in MVAL.
255 *
256  DO 30 im = 1, nm
257  m = mval( im )
258  p = pval( im )
259  q = qval( im )
260 *
261  DO 20 imat = 1, ntypes
262 *
263 * Do the tests only if DOTYPE( IMAT ) is true.
264 *
265  IF( .NOT.dotype( imat ) )
266  $ go to 20
267 *
268 * Generate X
269 *
270  IF( imat.EQ.1 ) THEN
271  CALL dlaror( 'L', 'I', m, m, x, ldx, iseed, work, iinfo )
272  IF( m .NE. 0 .AND. iinfo .NE. 0 ) THEN
273  WRITE( nout, fmt = 9999 ) m, iinfo
274  info = abs( iinfo )
275  go to 20
276  END IF
277  ELSE IF( imat.EQ.2 ) THEN
278  r = min( p, m-p, q, m-q )
279  DO i = 1, r
280  theta(i) = piover2 * dlarnd( 1, iseed )
281  END DO
282  CALL dlacsg( m, p, q, theta, iseed, x, ldx, work )
283  DO i = 1, m
284  DO j = 1, m
285  x(i+(j-1)*ldx) = x(i+(j-1)*ldx) +
286  $ orth*dlarnd(2,iseed)
287  END DO
288  END DO
289  ELSE
290  r = min( p, m-p, q, m-q )
291  DO i = 1, r+1
292  theta(i) = ten**(-dlarnd(1,iseed)*gapdigit)
293  END DO
294  DO i = 2, r+1
295  theta(i) = theta(i-1) + theta(i)
296  END DO
297  DO i = 1, r
298  theta(i) = piover2 * theta(i) / theta(r+1)
299  END DO
300  CALL dlacsg( m, p, q, theta, iseed, x, ldx, work )
301  END IF
302 *
303  nt = 9
304 *
305  CALL dcsdts( m, p, q, x, xf, ldx, u1, ldu1, u2, ldu2, v1t,
306  $ ldv1t, v2t, ldv2t, theta, iwork, work, lwork,
307  $ rwork, result )
308 *
309 * Print information about the tests that did not
310 * pass the threshold.
311 *
312  DO 10 i = 1, nt
313  IF( result( i ).GE.thresh ) THEN
314  IF( nfail.EQ.0 .AND. firstt ) THEN
315  firstt = .false.
316  CALL alahdg( nout, path )
317  END IF
318  WRITE( nout, fmt = 9998 )m, p, q, imat, i,
319  $ result( i )
320  nfail = nfail + 1
321  END IF
322  10 continue
323  nrun = nrun + nt
324  20 continue
325  30 continue
326 *
327 * Print a summary of the results.
328 *
329  CALL alasum( path, nout, nfail, nrun, 0 )
330 *
331  9999 format( ' DLAROR in DCKCSD: M = ', i5, ', INFO = ', i15 )
332  9998 format( ' M=', i4, ' P=', i4, ', Q=', i4, ', type ', i2,
333  $ ', test ', i2, ', ratio=', g13.6 )
334  return
335 *
336 * End of DCKCSD
337 *
338  END
339 *
340 *
341 *
342  SUBROUTINE dlacsg( M, P, Q, THETA, ISEED, X, LDX, WORK )
343  IMPLICIT NONE
344 *
345  INTEGER ldx, m, p, q
346  INTEGER iseed( 4 )
347  DOUBLE PRECISION theta( * )
348  DOUBLE PRECISION work( * ), x( ldx, * )
349 *
350  DOUBLE PRECISION one, zero
351  parameter( one = 1.0d0, zero = 0.0d0 )
352 *
353  INTEGER i, info, r
354 *
355  r = min( p, m-p, q, m-q )
356 *
357  CALL dlaset( 'Full', m, m, zero, zero, x, ldx )
358 *
359  DO i = 1, min(p,q)-r
360  x(i,i) = one
361  END DO
362  DO i = 1, r
363  x(min(p,q)-r+i,min(p,q)-r+i) = cos(theta(i))
364  END DO
365  DO i = 1, min(p,m-q)-r
366  x(p-i+1,m-i+1) = -one
367  END DO
368  DO i = 1, r
369  x(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) =
370  $ -sin(theta(r-i+1))
371  END DO
372  DO i = 1, min(m-p,q)-r
373  x(m-i+1,q-i+1) = one
374  END DO
375  DO i = 1, r
376  x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
377  $ sin(theta(r-i+1))
378  END DO
379  DO i = 1, min(m-p,m-q)-r
380  x(p+i,q+i) = one
381  END DO
382  DO i = 1, r
383  x(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) =
384  $ cos(theta(i))
385  END DO
386  CALL dlaror( 'Left', 'No init', p, m, x, ldx, iseed, work, info )
387  CALL dlaror( 'Left', 'No init', m-p, m, x(p+1,1), ldx,
388  $ iseed, work, info )
389  CALL dlaror( 'Right', 'No init', m, q, x, ldx, iseed,
390  $ work, info )
391  CALL dlaror( 'Right', 'No init', m, m-q,
392  $ x(1,q+1), ldx, iseed, work, info )
393 *
394  END
395