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