LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
cckcsd.f
Go to the documentation of this file.
1 *> \brief \b CCKCSD
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 CCKCSD( 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 * REAL THRESH
18 * ..
19 * .. Array Arguments ..
20 * INTEGER ISEED( 4 ), IWORK( * ), MVAL( * ), PVAL( * ),
21 * $ QVAL( * )
22 * REAL RWORK( * ), THETA( * )
23 * COMPLEX U1( * ), U2( * ), V1T( * ), V2T( * ),
24 * $ WORK( * ), X( * ), XF( * )
25 * ..
26 *
27 *
28 *> \par Purpose:
29 * =============
30 *>
31 *> \verbatim
32 *>
33 *> CCKCSD tests CUNCSD:
34 *> the CSD for an M-by-M unitary 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 REAL
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 COMPLEX array, dimension (MMAX*MMAX)
103 *> \endverbatim
104 *>
105 *> \param[out] XF
106 *> \verbatim
107 *> XF is COMPLEX array, dimension (MMAX*MMAX)
108 *> \endverbatim
109 *>
110 *> \param[out] U1
111 *> \verbatim
112 *> U1 is COMPLEX array, dimension (MMAX*MMAX)
113 *> \endverbatim
114 *>
115 *> \param[out] U2
116 *> \verbatim
117 *> U2 is COMPLEX array, dimension (MMAX*MMAX)
118 *> \endverbatim
119 *>
120 *> \param[out] V1T
121 *> \verbatim
122 *> V1T is COMPLEX array, dimension (MMAX*MMAX)
123 *> \endverbatim
124 *>
125 *> \param[out] V2T
126 *> \verbatim
127 *> V2T is COMPLEX array, dimension (MMAX*MMAX)
128 *> \endverbatim
129 *>
130 *> \param[out] THETA
131 *> \verbatim
132 *> THETA is REAL 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 COMPLEX array
143 *> \endverbatim
144 *>
145 *> \param[out] RWORK
146 *> \verbatim
147 *> RWORK is REAL 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 CLAROR 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 complex_eig
181 *
182 * =====================================================================
183  SUBROUTINE cckcsd( 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  REAL thresh
195 * ..
196 * .. Array Arguments ..
197  INTEGER iseed( 4 ), iwork( * ), mval( * ), pval( * ),
198  $ qval( * )
199  REAL rwork( * ), theta( * )
200  COMPLEX 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  REAL gapdigit, orth, piover2, realone, realzero, ten
212  parameter( gapdigit = 10.0e0, orth = 1.0e-4,
213  $ piover2 = 1.57079632679489662e0,
214  $ realone = 1.0e0, realzero = 0.0e0,
215  $ ten = 10.0e0 )
216  COMPLEX one, zero
217  parameter( one = (1.0e0,0.0e0), zero = (0.0e0,0.0e0) )
218 * ..
219 * .. Local Scalars ..
220  LOGICAL firstt
221  CHARACTER*3 path
222  INTEGER i, iinfo, im, imat, j, ldu1, ldu2, ldv1t,
223  $ ldv2t, ldx, lwork, m, nfail, nrun, nt, p, q, r
224 * ..
225 * .. Local Arrays ..
226  LOGICAL dotype( ntypes )
227  REAL result( ntests )
228 * ..
229 * .. External Subroutines ..
230  EXTERNAL alahdg, alareq, alasum, ccsdts, clacsg, claror,
231  $ claset
232 * ..
233 * .. Intrinsic Functions ..
234  INTRINSIC abs, min
235 * ..
236 * .. External Functions ..
237  REAL slaran, slarnd
238  EXTERNAL slaran, slarnd
239 * ..
240 * .. Executable Statements ..
241 *
242 * Initialize constants and the random number seed.
243 *
244  path( 1: 3 ) = 'CSD'
245  info = 0
246  nrun = 0
247  nfail = 0
248  firstt = .true.
249  CALL alareq( path, nmats, dotype, ntypes, nin, nout )
250  ldx = mmax
251  ldu1 = mmax
252  ldu2 = mmax
253  ldv1t = mmax
254  ldv2t = mmax
255  lwork = mmax*mmax
256 *
257 * Do for each value of M in MVAL.
258 *
259  DO 30 im = 1, nm
260  m = mval( im )
261  p = pval( im )
262  q = qval( im )
263 *
264  DO 20 imat = 1, ntypes
265 *
266 * Do the tests only if DOTYPE( IMAT ) is true.
267 *
268  IF( .NOT.dotype( imat ) )
269  $ go to 20
270 *
271 * Generate X
272 *
273  IF( imat.EQ.1 ) THEN
274  CALL claror( 'L', 'I', m, m, x, ldx, iseed, work, iinfo )
275  IF( m .NE. 0 .AND. iinfo .NE. 0 ) THEN
276  WRITE( nout, fmt = 9999 ) m, iinfo
277  info = abs( iinfo )
278  go to 20
279  END IF
280  ELSE IF( imat.EQ.2 ) THEN
281  r = min( p, m-p, q, m-q )
282  DO i = 1, r
283  theta(i) = piover2 * slarnd( 1, iseed )
284  END DO
285  CALL clacsg( m, p, q, theta, iseed, x, ldx, work )
286  DO i = 1, m
287  DO j = 1, m
288  x(i+(j-1)*ldx) = x(i+(j-1)*ldx) +
289  $ orth*slarnd(2,iseed)
290  END DO
291  END DO
292  ELSE IF( imat.EQ.3 ) THEN
293  r = min( p, m-p, q, m-q )
294  DO i = 1, r+1
295  theta(i) = ten**(-slarnd(1,iseed)*gapdigit)
296  END DO
297  DO i = 2, r+1
298  theta(i) = theta(i-1) + theta(i)
299  END DO
300  DO i = 1, r
301  theta(i) = piover2 * theta(i) / theta(r+1)
302  END DO
303  CALL clacsg( m, p, q, theta, iseed, x, ldx, work )
304  ELSE
305  CALL claset( 'F', m, m, zero, one, x, ldx )
306  DO i = 1, m
307  j = int( slaran( iseed ) * m ) + 1
308  IF( j .NE. i ) THEN
309  CALL csrot( m, x(1+(i-1)*ldx), 1, x(1+(j-1)*ldx),
310  $ 1, realzero, realone )
311  END IF
312  END DO
313  END IF
314 *
315  nt = 15
316 *
317  CALL ccsdts( m, p, q, x, xf, ldx, u1, ldu1, u2, ldu2, v1t,
318  $ ldv1t, v2t, ldv2t, theta, iwork, work, lwork,
319  $ rwork, result )
320 *
321 * Print information about the tests that did not
322 * pass the threshold.
323 *
324  DO 10 i = 1, nt
325  IF( result( i ).GE.thresh ) THEN
326  IF( nfail.EQ.0 .AND. firstt ) THEN
327  firstt = .false.
328  CALL alahdg( nout, path )
329  END IF
330  WRITE( nout, fmt = 9998 )m, p, q, imat, i,
331  $ result( i )
332  nfail = nfail + 1
333  END IF
334  10 CONTINUE
335  nrun = nrun + nt
336  20 CONTINUE
337  30 CONTINUE
338 *
339 * Print a summary of the results.
340 *
341  CALL alasum( path, nout, nfail, nrun, 0 )
342 *
343  9999 FORMAT( ' CLAROR in CCKCSD: M = ', i5, ', INFO = ', i15 )
344  9998 FORMAT( ' M=', i4, ' P=', i4, ', Q=', i4, ', type ', i2,
345  $ ', test ', i2, ', ratio=', g13.6 )
346  RETURN
347 *
348 * End of CCKCSD
349 *
350  END
351 *
352 *
353 *
354  SUBROUTINE clacsg( M, P, Q, THETA, ISEED, X, LDX, WORK )
355  IMPLICIT NONE
356 *
357  INTEGER ldx, m, p, q
358  INTEGER iseed( 4 )
359  REAL theta( * )
360  COMPLEX work( * ), x( ldx, * )
361 *
362  COMPLEX one, zero
363  parameter( one = (1.0e0,0.0e0), zero = (0.0e0,0.0e0) )
364 *
365  INTEGER i, info, r
366 *
367  r = min( p, m-p, q, m-q )
368 *
369  CALL claset( 'Full', m, m, zero, zero, x, ldx )
370 *
371  DO i = 1, min(p,q)-r
372  x(i,i) = one
373  END DO
374  DO i = 1, r
375  x(min(p,q)-r+i,min(p,q)-r+i) = cmplx( cos(theta(i)), 0.0e0 )
376  END DO
377  DO i = 1, min(p,m-q)-r
378  x(p-i+1,m-i+1) = -one
379  END DO
380  DO i = 1, r
381  x(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) =
382  $ cmplx( -sin(theta(r-i+1)), 0.0e0 )
383  END DO
384  DO i = 1, min(m-p,q)-r
385  x(m-i+1,q-i+1) = one
386  END DO
387  DO i = 1, r
388  x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
389  $ cmplx( sin(theta(r-i+1)), 0.0e0 )
390  END DO
391  DO i = 1, min(m-p,m-q)-r
392  x(p+i,q+i) = one
393  END DO
394  DO i = 1, r
395  x(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) =
396  $ cmplx( cos(theta(i)), 0.0e0 )
397  END DO
398  CALL claror( 'Left', 'No init', p, m, x, ldx, iseed, work, info )
399  CALL claror( 'Left', 'No init', m-p, m, x(p+1,1), ldx,
400  $ iseed, work, info )
401  CALL claror( 'Right', 'No init', m, q, x, ldx, iseed,
402  $ work, info )
403  CALL claror( 'Right', 'No init', m, m-q,
404  $ x(1,q+1), ldx, iseed, work, info )
405 *
406  END
407