LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
zckcsd.f
Go to the documentation of this file.
1 *> \brief \b ZCKCSD
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 ZCKCSD( 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 * COMPLEX*16 U1( * ), U2( * ), V1T( * ), V2T( * ),
24 * $ WORK( * ), X( * ), XF( * )
25 * ..
26 *
27 *
28 *> \par Purpose:
29 * =============
30 *>
31 *> \verbatim
32 *>
33 *> ZCKCSD tests ZUNCSD:
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 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 COMPLEX*16 array, dimension (MMAX*MMAX)
103 *> \endverbatim
104 *>
105 *> \param[out] XF
106 *> \verbatim
107 *> XF is COMPLEX*16 array, dimension (MMAX*MMAX)
108 *> \endverbatim
109 *>
110 *> \param[out] U1
111 *> \verbatim
112 *> U1 is COMPLEX*16 array, dimension (MMAX*MMAX)
113 *> \endverbatim
114 *>
115 *> \param[out] U2
116 *> \verbatim
117 *> U2 is COMPLEX*16 array, dimension (MMAX*MMAX)
118 *> \endverbatim
119 *>
120 *> \param[out] V1T
121 *> \verbatim
122 *> V1T is COMPLEX*16 array, dimension (MMAX*MMAX)
123 *> \endverbatim
124 *>
125 *> \param[out] V2T
126 *> \verbatim
127 *> V2T is COMPLEX*16 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 COMPLEX*16 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 ZLAROR 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 December 2016
179 *
180 *> \ingroup complex16_eig
181 *
182 * =====================================================================
183  SUBROUTINE zckcsd( 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.7.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 * December 2016
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  COMPLEX*16 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, ORTH, PIOVER2, REALONE, REALZERO, TEN
212  parameter( gapdigit = 18.0d0, orth = 1.0d-12,
213  $ piover2 = 1.57079632679489662d0,
214  $ realone = 1.0d0, realzero = 0.0d0,
215  $ ten = 10.0d0 )
216  COMPLEX*16 ONE, ZERO
217  PARAMETER ( ONE = (1.0d0,0.0d0), zero = (0.0d0,0.0d0) )
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  DOUBLE PRECISION RESULT( NTESTS )
228 * ..
229 * .. External Subroutines ..
230  EXTERNAL alahdg, alareq, alasum, zcsdts, zlacsg, zlaror,
231  $ zlaset
232 * ..
233 * .. Intrinsic Functions ..
234  INTRINSIC abs, min
235 * ..
236 * .. External Functions ..
237  DOUBLE PRECISION DLARAN, DLARND
238  EXTERNAL DLARAN, DLARND
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 zlaror( '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 * dlarnd( 1, iseed )
284  END DO
285  CALL zlacsg( 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*dlarnd(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**(-dlarnd(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 zlacsg( m, p, q, theta, iseed, x, ldx, work )
304  ELSE
305  CALL zlaset( 'F', m, m, zero, one, x, ldx )
306  DO i = 1, m
307  j = int( dlaran( iseed ) * m ) + 1
308  IF( j .NE. i ) THEN
309  CALL zdrot( 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 zcsdts( 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( ' ZLAROR in ZCKCSD: 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 ZCKCSD
349 *
350  END
351 *
352 *
353 *
354  SUBROUTINE zlacsg( M, P, Q, THETA, ISEED, X, LDX, WORK )
355  IMPLICIT NONE
356 *
357  INTEGER LDX, M, P, Q
358  INTEGER ISEED( 4 )
359  DOUBLE PRECISION THETA( * )
360  COMPLEX*16 WORK( * ), X( LDX, * )
361 *
362  COMPLEX*16 ONE, ZERO
363  PARAMETER ( ONE = (1.0d0,0.0d0), zero = (0.0d0,0.0d0) )
364 *
365  INTEGER I, INFO, R
366 *
367  r = min( p, m-p, q, m-q )
368 *
369  CALL zlaset( '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) = dcmplx( cos(theta(i)), 0.0d0 )
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  $ dcmplx( -sin(theta(r-i+1)), 0.0d0 )
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  $ dcmplx( sin(theta(r-i+1)), 0.0d0 )
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  $ dcmplx( cos(theta(i)), 0.0d0 )
397  END DO
398  CALL zlaror( 'Left', 'No init', p, m, x, ldx, iseed, work, info )
399  CALL zlaror( 'Left', 'No init', m-p, m, x(p+1,1), ldx,
400  $ iseed, work, info )
401  CALL zlaror( 'Right', 'No init', m, q, x, ldx, iseed,
402  $ work, info )
403  CALL zlaror( 'Right', 'No init', m, m-q,
404  $ x(1,q+1), ldx, iseed, work, info )
405 *
406  END
407 
zckcsd
subroutine zckcsd(NM, MVAL, PVAL, QVAL, NMATS, ISEED, THRESH, MMAX, X, XF, U1, U2, V1T, V2T, THETA, IWORK, WORK, RWORK, NIN, NOUT, INFO)
ZCKCSD
Definition: zckcsd.f:186
alahdg
subroutine alahdg(IOUNIT, PATH)
ALAHDG
Definition: alahdg.f:64
zlaset
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:108
zcsdts
subroutine zcsdts(M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK, RWORK, RESULT)
ZCSDTS
Definition: zcsdts.f:231
alareq
subroutine alareq(PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT)
ALAREQ
Definition: alareq.f:92
alasum
subroutine alasum(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASUM
Definition: alasum.f:75
zlaror
subroutine zlaror(SIDE, INIT, M, N, A, LDA, ISEED, X, INFO)
ZLAROR
Definition: zlaror.f:160
zlacsg
subroutine zlacsg(M, P, Q, THETA, ISEED, X, LDX, WORK)
Definition: zckcsd.f:355
zdrot
subroutine zdrot(N, CX, INCX, CY, INCY, C, S)
ZDROT
Definition: zdrot.f:100