00001 SUBROUTINE CDRVES( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
00002 $ NOUNIT, A, LDA, H, HT, W, WT, VS, LDVS, RESULT,
00003 $ WORK, NWORK, RWORK, IWORK, BWORK, INFO )
00004
00005
00006
00007
00008
00009
00010 INTEGER INFO, LDA, LDVS, NOUNIT, NSIZES, NTYPES, NWORK
00011 REAL THRESH
00012
00013
00014 LOGICAL BWORK( * ), DOTYPE( * )
00015 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
00016 REAL RESULT( 13 ), RWORK( * )
00017 COMPLEX A( LDA, * ), H( LDA, * ), HT( LDA, * ),
00018 $ VS( LDVS, * ), W( * ), WORK( * ), WT( * )
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289 COMPLEX CZERO
00290 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) )
00291 COMPLEX CONE
00292 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) )
00293 REAL ZERO, ONE
00294 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00295 INTEGER MAXTYP
00296 PARAMETER ( MAXTYP = 21 )
00297
00298
00299 LOGICAL BADNN
00300 CHARACTER SORT
00301 CHARACTER*3 PATH
00302 INTEGER I, IINFO, IMODE, ISORT, ITYPE, IWK, J, JCOL,
00303 $ JSIZE, JTYPE, KNTEIG, LWORK, MTYPES, N,
00304 $ NERRS, NFAIL, NMAX, NNWORK, NTEST, NTESTF,
00305 $ NTESTT, RSUB, SDIM
00306 REAL ANORM, COND, CONDS, OVFL, RTULP, RTULPI, ULP,
00307 $ ULPINV, UNFL
00308
00309
00310 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ),
00311 $ KMAGN( MAXTYP ), KMODE( MAXTYP ),
00312 $ KTYPE( MAXTYP )
00313 REAL RES( 2 )
00314
00315
00316 LOGICAL SELVAL( 20 )
00317 REAL SELWI( 20 ), SELWR( 20 )
00318
00319
00320 INTEGER SELDIM, SELOPT
00321
00322
00323 COMMON / SSLCT / SELOPT, SELDIM, SELVAL, SELWR, SELWI
00324
00325
00326 LOGICAL CSLECT
00327 REAL SLAMCH
00328 EXTERNAL CSLECT, SLAMCH
00329
00330
00331 EXTERNAL CGEES, CHST01, CLACPY, CLATME, CLATMR, CLATMS,
00332 $ CLASET, SLABAD, SLASUM, XERBLA
00333
00334
00335 INTRINSIC ABS, CMPLX, MAX, MIN, SQRT
00336
00337
00338 DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
00339 DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
00340 $ 3, 1, 2, 3 /
00341 DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
00342 $ 1, 5, 5, 5, 4, 3, 1 /
00343 DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 /
00344
00345
00346
00347 PATH( 1: 1 ) = 'Complex precision'
00348 PATH( 2: 3 ) = 'ES'
00349
00350
00351
00352 NTESTT = 0
00353 NTESTF = 0
00354 INFO = 0
00355 SELOPT = 0
00356
00357
00358
00359 BADNN = .FALSE.
00360 NMAX = 0
00361 DO 10 J = 1, NSIZES
00362 NMAX = MAX( NMAX, NN( J ) )
00363 IF( NN( J ).LT.0 )
00364 $ BADNN = .TRUE.
00365 10 CONTINUE
00366
00367
00368
00369 IF( NSIZES.LT.0 ) THEN
00370 INFO = -1
00371 ELSE IF( BADNN ) THEN
00372 INFO = -2
00373 ELSE IF( NTYPES.LT.0 ) THEN
00374 INFO = -3
00375 ELSE IF( THRESH.LT.ZERO ) THEN
00376 INFO = -6
00377 ELSE IF( NOUNIT.LE.0 ) THEN
00378 INFO = -7
00379 ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN
00380 INFO = -9
00381 ELSE IF( LDVS.LT.1 .OR. LDVS.LT.NMAX ) THEN
00382 INFO = -15
00383 ELSE IF( 5*NMAX+2*NMAX**2.GT.NWORK ) THEN
00384 INFO = -18
00385 END IF
00386
00387 IF( INFO.NE.0 ) THEN
00388 CALL XERBLA( 'CDRVES', -INFO )
00389 RETURN
00390 END IF
00391
00392
00393
00394 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
00395 $ RETURN
00396
00397
00398
00399 UNFL = SLAMCH( 'Safe minimum' )
00400 OVFL = ONE / UNFL
00401 CALL SLABAD( UNFL, OVFL )
00402 ULP = SLAMCH( 'Precision' )
00403 ULPINV = ONE / ULP
00404 RTULP = SQRT( ULP )
00405 RTULPI = ONE / RTULP
00406
00407
00408
00409 NERRS = 0
00410
00411 DO 240 JSIZE = 1, NSIZES
00412 N = NN( JSIZE )
00413 IF( NSIZES.NE.1 ) THEN
00414 MTYPES = MIN( MAXTYP, NTYPES )
00415 ELSE
00416 MTYPES = MIN( MAXTYP+1, NTYPES )
00417 END IF
00418
00419 DO 230 JTYPE = 1, MTYPES
00420 IF( .NOT.DOTYPE( JTYPE ) )
00421 $ GO TO 230
00422
00423
00424
00425 DO 20 J = 1, 4
00426 IOLDSD( J ) = ISEED( J )
00427 20 CONTINUE
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445 IF( MTYPES.GT.MAXTYP )
00446 $ GO TO 90
00447
00448 ITYPE = KTYPE( JTYPE )
00449 IMODE = KMODE( JTYPE )
00450
00451
00452
00453 GO TO ( 30, 40, 50 )KMAGN( JTYPE )
00454
00455 30 CONTINUE
00456 ANORM = ONE
00457 GO TO 60
00458
00459 40 CONTINUE
00460 ANORM = OVFL*ULP
00461 GO TO 60
00462
00463 50 CONTINUE
00464 ANORM = UNFL*ULPINV
00465 GO TO 60
00466
00467 60 CONTINUE
00468
00469 CALL CLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA )
00470 IINFO = 0
00471 COND = ULPINV
00472
00473
00474
00475 IF( ITYPE.EQ.1 ) THEN
00476
00477
00478
00479 IINFO = 0
00480
00481 ELSE IF( ITYPE.EQ.2 ) THEN
00482
00483
00484
00485 DO 70 JCOL = 1, N
00486 A( JCOL, JCOL ) = CMPLX( ANORM )
00487 70 CONTINUE
00488
00489 ELSE IF( ITYPE.EQ.3 ) THEN
00490
00491
00492
00493 DO 80 JCOL = 1, N
00494 A( JCOL, JCOL ) = CMPLX( ANORM )
00495 IF( JCOL.GT.1 )
00496 $ A( JCOL, JCOL-1 ) = CONE
00497 80 CONTINUE
00498
00499 ELSE IF( ITYPE.EQ.4 ) THEN
00500
00501
00502
00503 CALL CLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND,
00504 $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ),
00505 $ IINFO )
00506
00507 ELSE IF( ITYPE.EQ.5 ) THEN
00508
00509
00510
00511 CALL CLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND,
00512 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ),
00513 $ IINFO )
00514
00515 ELSE IF( ITYPE.EQ.6 ) THEN
00516
00517
00518
00519 IF( KCONDS( JTYPE ).EQ.1 ) THEN
00520 CONDS = ONE
00521 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN
00522 CONDS = RTULPI
00523 ELSE
00524 CONDS = ZERO
00525 END IF
00526
00527 CALL CLATME( N, 'D', ISEED, WORK, IMODE, COND, CONE, ' ',
00528 $ 'T', 'T', 'T', RWORK, 4, CONDS, N, N, ANORM,
00529 $ A, LDA, WORK( 2*N+1 ), IINFO )
00530
00531 ELSE IF( ITYPE.EQ.7 ) THEN
00532
00533
00534
00535 CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
00536 $ 'T', 'N', WORK( N+1 ), 1, ONE,
00537 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
00538 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
00539
00540 ELSE IF( ITYPE.EQ.8 ) THEN
00541
00542
00543
00544 CALL CLATMR( N, N, 'D', ISEED, 'H', WORK, 6, ONE, CONE,
00545 $ 'T', 'N', WORK( N+1 ), 1, ONE,
00546 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
00547 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
00548
00549 ELSE IF( ITYPE.EQ.9 ) THEN
00550
00551
00552
00553 CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
00554 $ 'T', 'N', WORK( N+1 ), 1, ONE,
00555 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
00556 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
00557 IF( N.GE.4 ) THEN
00558 CALL CLASET( 'Full', 2, N, CZERO, CZERO, A, LDA )
00559 CALL CLASET( 'Full', N-3, 1, CZERO, CZERO, A( 3, 1 ),
00560 $ LDA )
00561 CALL CLASET( 'Full', N-3, 2, CZERO, CZERO,
00562 $ A( 3, N-1 ), LDA )
00563 CALL CLASET( 'Full', 1, N, CZERO, CZERO, A( N, 1 ),
00564 $ LDA )
00565 END IF
00566
00567 ELSE IF( ITYPE.EQ.10 ) THEN
00568
00569
00570
00571 CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
00572 $ 'T', 'N', WORK( N+1 ), 1, ONE,
00573 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0,
00574 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
00575
00576 ELSE
00577
00578 IINFO = 1
00579 END IF
00580
00581 IF( IINFO.NE.0 ) THEN
00582 WRITE( NOUNIT, FMT = 9992 )'Generator', IINFO, N, JTYPE,
00583 $ IOLDSD
00584 INFO = ABS( IINFO )
00585 RETURN
00586 END IF
00587
00588 90 CONTINUE
00589
00590
00591
00592 DO 220 IWK = 1, 2
00593 IF( IWK.EQ.1 ) THEN
00594 NNWORK = 3*N
00595 ELSE
00596 NNWORK = 5*N + 2*N**2
00597 END IF
00598 NNWORK = MAX( NNWORK, 1 )
00599
00600
00601
00602 DO 100 J = 1, 13
00603 RESULT( J ) = -ONE
00604 100 CONTINUE
00605
00606
00607
00608 DO 180 ISORT = 0, 1
00609 IF( ISORT.EQ.0 ) THEN
00610 SORT = 'N'
00611 RSUB = 0
00612 ELSE
00613 SORT = 'S'
00614 RSUB = 6
00615 END IF
00616
00617
00618
00619 CALL CLACPY( 'F', N, N, A, LDA, H, LDA )
00620 CALL CGEES( 'V', SORT, CSLECT, N, H, LDA, SDIM, W, VS,
00621 $ LDVS, WORK, NNWORK, RWORK, BWORK, IINFO )
00622 IF( IINFO.NE.0 ) THEN
00623 RESULT( 1+RSUB ) = ULPINV
00624 WRITE( NOUNIT, FMT = 9992 )'CGEES1', IINFO, N,
00625 $ JTYPE, IOLDSD
00626 INFO = ABS( IINFO )
00627 GO TO 190
00628 END IF
00629
00630
00631
00632 RESULT( 1+RSUB ) = ZERO
00633 DO 120 J = 1, N - 1
00634 DO 110 I = J + 1, N
00635 IF( H( I, J ).NE.ZERO )
00636 $ RESULT( 1+RSUB ) = ULPINV
00637 110 CONTINUE
00638 120 CONTINUE
00639
00640
00641
00642 LWORK = MAX( 1, 2*N*N )
00643 CALL CHST01( N, 1, N, A, LDA, H, LDA, VS, LDVS, WORK,
00644 $ LWORK, RWORK, RES )
00645 RESULT( 2+RSUB ) = RES( 1 )
00646 RESULT( 3+RSUB ) = RES( 2 )
00647
00648
00649
00650 RESULT( 4+RSUB ) = ZERO
00651 DO 130 I = 1, N
00652 IF( H( I, I ).NE.W( I ) )
00653 $ RESULT( 4+RSUB ) = ULPINV
00654 130 CONTINUE
00655
00656
00657
00658 CALL CLACPY( 'F', N, N, A, LDA, HT, LDA )
00659 CALL CGEES( 'N', SORT, CSLECT, N, HT, LDA, SDIM, WT,
00660 $ VS, LDVS, WORK, NNWORK, RWORK, BWORK,
00661 $ IINFO )
00662 IF( IINFO.NE.0 ) THEN
00663 RESULT( 5+RSUB ) = ULPINV
00664 WRITE( NOUNIT, FMT = 9992 )'CGEES2', IINFO, N,
00665 $ JTYPE, IOLDSD
00666 INFO = ABS( IINFO )
00667 GO TO 190
00668 END IF
00669
00670 RESULT( 5+RSUB ) = ZERO
00671 DO 150 J = 1, N
00672 DO 140 I = 1, N
00673 IF( H( I, J ).NE.HT( I, J ) )
00674 $ RESULT( 5+RSUB ) = ULPINV
00675 140 CONTINUE
00676 150 CONTINUE
00677
00678
00679
00680 RESULT( 6+RSUB ) = ZERO
00681 DO 160 I = 1, N
00682 IF( W( I ).NE.WT( I ) )
00683 $ RESULT( 6+RSUB ) = ULPINV
00684 160 CONTINUE
00685
00686
00687
00688 IF( ISORT.EQ.1 ) THEN
00689 RESULT( 13 ) = ZERO
00690 KNTEIG = 0
00691 DO 170 I = 1, N
00692 IF( CSLECT( W( I ) ) )
00693 $ KNTEIG = KNTEIG + 1
00694 IF( I.LT.N ) THEN
00695 IF( CSLECT( W( I+1 ) ) .AND.
00696 $ ( .NOT.CSLECT( W( I ) ) ) )RESULT( 13 )
00697 $ = ULPINV
00698 END IF
00699 170 CONTINUE
00700 IF( SDIM.NE.KNTEIG )
00701 $ RESULT( 13 ) = ULPINV
00702 END IF
00703
00704 180 CONTINUE
00705
00706
00707
00708 190 CONTINUE
00709
00710 NTEST = 0
00711 NFAIL = 0
00712 DO 200 J = 1, 13
00713 IF( RESULT( J ).GE.ZERO )
00714 $ NTEST = NTEST + 1
00715 IF( RESULT( J ).GE.THRESH )
00716 $ NFAIL = NFAIL + 1
00717 200 CONTINUE
00718
00719 IF( NFAIL.GT.0 )
00720 $ NTESTF = NTESTF + 1
00721 IF( NTESTF.EQ.1 ) THEN
00722 WRITE( NOUNIT, FMT = 9999 )PATH
00723 WRITE( NOUNIT, FMT = 9998 )
00724 WRITE( NOUNIT, FMT = 9997 )
00725 WRITE( NOUNIT, FMT = 9996 )
00726 WRITE( NOUNIT, FMT = 9995 )THRESH
00727 WRITE( NOUNIT, FMT = 9994 )
00728 NTESTF = 2
00729 END IF
00730
00731 DO 210 J = 1, 13
00732 IF( RESULT( J ).GE.THRESH ) THEN
00733 WRITE( NOUNIT, FMT = 9993 )N, IWK, IOLDSD, JTYPE,
00734 $ J, RESULT( J )
00735 END IF
00736 210 CONTINUE
00737
00738 NERRS = NERRS + NFAIL
00739 NTESTT = NTESTT + NTEST
00740
00741 220 CONTINUE
00742 230 CONTINUE
00743 240 CONTINUE
00744
00745
00746
00747 CALL SLASUM( PATH, NOUNIT, NERRS, NTESTT )
00748
00749 9999 FORMAT( / 1X, A3, ' -- Complex Schur Form Decomposition Driver',
00750 $ / ' Matrix types (see CDRVES for details): ' )
00751
00752 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ',
00753 $ ' ', ' 5=Diagonal: geometr. spaced entries.',
00754 $ / ' 2=Identity matrix. ', ' 6=Diagona',
00755 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
00756 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
00757 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
00758 $ 'mall, evenly spaced.' )
00759 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
00760 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
00761 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
00762 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
00763 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
00764 $ 'lex ', A6, / ' 12=Well-cond., random complex ', A6, ' ',
00765 $ ' 17=Ill-cond., large rand. complx ', A4, / ' 13=Ill-condi',
00766 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
00767 $ ' complx ', A4 )
00768 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
00769 $ 'with small random entries.', / ' 20=Matrix with large ran',
00770 $ 'dom entries. ', / )
00771 9995 FORMAT( ' Tests performed with test threshold =', F8.2,
00772 $ / ' ( A denotes A on input and T denotes A on output)',
00773 $ / / ' 1 = 0 if T in Schur form (no sort), ',
00774 $ ' 1/ulp otherwise', /
00775 $ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
00776 $ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ',
00777 $ / ' 4 = 0 if W are eigenvalues of T (no sort),',
00778 $ ' 1/ulp otherwise', /
00779 $ ' 5 = 0 if T same no matter if VS computed (no sort),',
00780 $ ' 1/ulp otherwise', /
00781 $ ' 6 = 0 if W same no matter if VS computed (no sort)',
00782 $ ', 1/ulp otherwise' )
00783 9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise',
00784 $ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
00785 $ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
00786 $ / ' 10 = 0 if W are eigenvalues of T (sort),',
00787 $ ' 1/ulp otherwise', /
00788 $ ' 11 = 0 if T same no matter if VS computed (sort),',
00789 $ ' 1/ulp otherwise', /
00790 $ ' 12 = 0 if W same no matter if VS computed (sort),',
00791 $ ' 1/ulp otherwise', /
00792 $ ' 13 = 0 if sorting succesful, 1/ulp otherwise', / )
00793 9993 FORMAT( ' N=', I5, ', IWK=', I2, ', seed=', 4( I4, ',' ),
00794 $ ' type ', I2, ', test(', I2, ')=', G10.3 )
00795 9992 FORMAT( ' CDRVES: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00796 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
00797
00798 RETURN
00799
00800
00801
00802 END