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*> \ingroup complex_eig
179*
180*  =====================================================================
181      SUBROUTINE CCKCSD( 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      REAL               THRESH
192*     ..
193*     .. Array Arguments ..
194      INTEGER            ISEED( 4 ), IWORK( * ), MVAL( * ), PVAL( * ),
195     $                   QVAL( * )
196      REAL               RWORK( * ), THETA( * )
197      COMPLEX            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      REAL               GAPDIGIT, ORTH, REALONE, REALZERO, TEN
209      PARAMETER          ( GAPDIGIT = 10.0E0, ORTH = 1.0E-4,
210     $                     REALONE = 1.0E0, REALZERO = 0.0E0,
211     $                     TEN = 10.0E0 )
212      COMPLEX            ONE, ZERO
213      PARAMETER          ( ONE = (1.0E0,0.0E0), ZERO = (0.0E0,0.0E0) )
214      REAL               PIOVER2
215      PARAMETER ( PIOVER2 = 1.57079632679489661923132169163975144210E0 )
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      REAL               RESULT( NTESTS )
226*     ..
227*     .. External Subroutines ..
228      EXTERNAL           ALAHDG, ALAREQ, ALASUM, CCSDTS, CLACSG, CLAROR,
229     $                   CLASET, CSROT
230*     ..
231*     .. Intrinsic Functions ..
232      INTRINSIC          ABS, MIN
233*     ..
234*     .. External Functions ..
235      REAL               SLARAN, SLARND
236      EXTERNAL           SLARAN, SLARND
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 CLAROR( '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 * SLARND( 1, ISEED )
282               END DO
283               CALL CLACSG( 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*SLARND(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**(-SLARND(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 CLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
302            ELSE
303               CALL CLASET( 'F', M, M, ZERO, ONE, X, LDX )
304               DO I = 1, M
305                  J = INT( SLARAN( ISEED ) * M ) + 1
306                  IF( J .NE. I ) THEN
307                     CALL CSROT( M, X(1+(I-1)*LDX), 1, X(1+(J-1)*LDX),
308     $                 1, REALZERO, REALONE )
309                  END IF
310               END DO
311            END IF
312*
313            NT = 15
314*
315            CALL CCSDTS( 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( ' CLAROR in CCKCSD: 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 CCKCSD
347*
348      END
349*
350*
351*
352      SUBROUTINE CLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
353      IMPLICIT NONE
354*
355      INTEGER            LDX, M, P, Q
356      INTEGER            ISEED( 4 )
357      REAL               THETA( * )
358      COMPLEX            WORK( * ), X( LDX, * )
359*
360      COMPLEX            ONE, ZERO
361      PARAMETER          ( ONE = (1.0E0,0.0E0), ZERO = (0.0E0,0.0E0) )
362*
363      INTEGER            I, INFO, R
364*
365      R = MIN( P, M-P, Q, M-Q )
366*
367      CALL CLASET( '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) = CMPLX ( COS(THETA(I)), 0.0E0 )
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     $      CMPLX( -SIN(THETA(R-I+1)), 0.0E0 )
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     $      CMPLX( SIN(THETA(R-I+1)), 0.0E0 )
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     $      CMPLX( COS(THETA(I)), 0.0E0 )
395      END DO
396      CALL CLAROR( 'Left', 'No init', P, M, X, LDX, ISEED, WORK, INFO )
397      CALL CLAROR( 'Left', 'No init', M-P, M, X(P+1,1), LDX,
398     $             ISEED, WORK, INFO )
399      CALL CLAROR( 'Right', 'No init', M, Q, X, LDX, ISEED,
400     $             WORK, INFO )
401      CALL CLAROR( 'Right', 'No init', M, M-Q,
402     $             X(1,Q+1), LDX, ISEED, WORK, INFO )
403*
404      END
405
406