TESTING/LIN/zchksp.f

   1 *> \brief \b ZCHKSP
   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 ZCHKSP( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR,
  12 *                          NMAX, A, AFAC, AINV, B, X, XACT, WORK, RWORK,
  13 *                          IWORK, NOUT )
  14 *
  15 *       .. Scalar Arguments ..
  16 *       LOGICAL            TSTERR
  17 *       INTEGER            NMAX, NN, NNS, NOUT
  18 *       DOUBLE PRECISION   THRESH
  19 *       ..
  20 *       .. Array Arguments ..
  21 *       LOGICAL            DOTYPE( * )
  22 *       INTEGER            IWORK( * ), NSVAL( * ), NVAL( * )
  23 *       DOUBLE PRECISION   RWORK( * )
  24 *       COMPLEX*16         A( * ), AFAC( * ), AINV( * ), B( * ),
  25 *      $                   WORK( * ), X( * ), XACT( * )
  26 *       ..
  27 *
  28 *
  29 *> \par Purpose:
  30 *  =============
  31 *>
  32 *> \verbatim
  33 *>
  34 *> ZCHKSP tests ZSPTRF, -TRI, -TRS, -RFS, and -CON
  35 *> \endverbatim
  36 *
  37 *  Arguments:
  38 *  ==========
  39 *
  40 *> \param[in] DOTYPE
  41 *> \verbatim
  42 *>          DOTYPE is LOGICAL array, dimension (NTYPES)
  43 *>          The matrix types to be used for testing.  Matrices of type j
  44 *>          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
  45 *>          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
  46 *> \endverbatim
  47 *>
  48 *> \param[in] NN
  49 *> \verbatim
  50 *>          NN is INTEGER
  51 *>          The number of values of N contained in the vector NVAL.
  52 *> \endverbatim
  53 *>
  54 *> \param[in] NVAL
  55 *> \verbatim
  56 *>          NVAL is INTEGER array, dimension (NN)
  57 *>          The values of the matrix dimension N.
  58 *> \endverbatim
  59 *>
  60 *> \param[in] NNS
  61 *> \verbatim
  62 *>          NNS is INTEGER
  63 *>          The number of values of NRHS contained in the vector NSVAL.
  64 *> \endverbatim
  65 *>
  66 *> \param[in] NSVAL
  67 *> \verbatim
  68 *>          NSVAL is INTEGER array, dimension (NNS)
  69 *>          The values of the number of right hand sides NRHS.
  70 *> \endverbatim
  71 *>
  72 *> \param[in] THRESH
  73 *> \verbatim
  74 *>          THRESH is DOUBLE PRECISION
  75 *>          The threshold value for the test ratios.  A result is
  76 *>          included in the output file if RESULT >= THRESH.  To have
  77 *>          every test ratio printed, use THRESH = 0.
  78 *> \endverbatim
  79 *>
  80 *> \param[in] TSTERR
  81 *> \verbatim
  82 *>          TSTERR is LOGICAL
  83 *>          Flag that indicates whether error exits are to be tested.
  84 *> \endverbatim
  85 *>
  86 *> \param[in] NMAX
  87 *> \verbatim
  88 *>          NMAX is INTEGER
  89 *>          The maximum value permitted for N, used in dimensioning the
  90 *>          work arrays.
  91 *> \endverbatim
  92 *>
  93 *> \param[out] A
  94 *> \verbatim
  95 *>          A is COMPLEX*16 array, dimension
  96 *>                      (NMAX*(NMAX+1)/2)
  97 *> \endverbatim
  98 *>
  99 *> \param[out] AFAC
 100 *> \verbatim
 101 *>          AFAC is COMPLEX*16 array, dimension
 102 *>                      (NMAX*(NMAX+1)/2)
 103 *> \endverbatim
 104 *>
 105 *> \param[out] AINV
 106 *> \verbatim
 107 *>          AINV is COMPLEX*16 array, dimension
 108 *>                      (NMAX*(NMAX+1)/2)
 109 *> \endverbatim
 110 *>
 111 *> \param[out] B
 112 *> \verbatim
 113 *>          B is COMPLEX*16 array, dimension (NMAX*NSMAX)
 114 *>          where NSMAX is the largest entry in NSVAL.
 115 *> \endverbatim
 116 *>
 117 *> \param[out] X
 118 *> \verbatim
 119 *>          X is COMPLEX*16 array, dimension (NMAX*NSMAX)
 120 *> \endverbatim
 121 *>
 122 *> \param[out] XACT
 123 *> \verbatim
 124 *>          XACT is COMPLEX*16 array, dimension (NMAX*NSMAX)
 125 *> \endverbatim
 126 *>
 127 *> \param[out] WORK
 128 *> \verbatim
 129 *>          WORK is COMPLEX*16 array, dimension
 130 *>                      (NMAX*max(2,NSMAX))
 131 *> \endverbatim
 132 *>
 133 *> \param[out] RWORK
 134 *> \verbatim
 135 *>          RWORK is DOUBLE PRECISION array,
 136 *>                                 dimension (NMAX+2*NSMAX)
 137 *> \endverbatim
 138 *>
 139 *> \param[out] IWORK
 140 *> \verbatim
 141 *>          IWORK is INTEGER array, dimension (NMAX)
 142 *> \endverbatim
 143 *>
 144 *> \param[in] NOUT
 145 *> \verbatim
 146 *>          NOUT is INTEGER
 147 *>          The unit number for output.
 148 *> \endverbatim
 149 *
 150 *  Authors:
 151 *  ========
 152 *
 153 *> \author Univ. of Tennessee
 154 *> \author Univ. of California Berkeley
 155 *> \author Univ. of Colorado Denver
 156 *> \author NAG Ltd.
 157 *
 158 *> \date November 2011
 159 *
 160 *> \ingroup complex16_lin
 161 *
 162 *  =====================================================================
 163       SUBROUTINE ZCHKSP( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR,
 164      $                   NMAX, A, AFAC, AINV, B, X, XACT, WORK, RWORK,
 165      $                   IWORK, NOUT )
 166 *
 167 *  -- LAPACK test routine (version 3.4.0) --
 168 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 169 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 170 *     November 2011
 171 *
 172 *     .. Scalar Arguments ..
 173       LOGICAL            TSTERR
 174       INTEGER            NMAX, NN, NNS, NOUT
 175       DOUBLE PRECISION   THRESH
 176 *     ..
 177 *     .. Array Arguments ..
 178       LOGICAL            DOTYPE( * )
 179       INTEGER            IWORK( * ), NSVAL( * ), NVAL( * )
 180       DOUBLE PRECISION   RWORK( * )
 181       COMPLEX*16         A( * ), AFAC( * ), AINV( * ), B( * ),
 182      $                   WORK( * ), X( * ), XACT( * )
 183 *     ..
 184 *
 185 *  =====================================================================
 186 *
 187 *     .. Parameters ..
 188       DOUBLE PRECISION   ZERO
 189       PARAMETER          ( ZERO = 0.0D+0 )
 190       INTEGER            NTYPES
 191       PARAMETER          ( NTYPES = 11 )
 192       INTEGER            NTESTS
 193       PARAMETER          ( NTESTS = 8 )
 194 *     ..
 195 *     .. Local Scalars ..
 196       LOGICAL            TRFCON, ZEROT
 197       CHARACTER          DIST, PACKIT, TYPE, UPLO, XTYPE
 198       CHARACTER*3        PATH
 199       INTEGER            I, I1, I2, IMAT, IN, INFO, IOFF, IRHS, IUPLO,
 200      $                   IZERO, J, K, KL, KU, LDA, MODE, N, NERRS,
 201      $                   NFAIL, NIMAT, NPP, NRHS, NRUN, NT
 202       DOUBLE PRECISION   ANORM, CNDNUM, RCOND, RCONDC
 203 *     ..
 204 *     .. Local Arrays ..
 205       CHARACTER          UPLOS( 2 )
 206       INTEGER            ISEED( 4 ), ISEEDY( 4 )
 207       DOUBLE PRECISION   RESULT( NTESTS )
 208 *     ..
 209 *     .. External Functions ..
 210       LOGICAL            LSAME
 211       DOUBLE PRECISION   DGET06, ZLANSP
 212       EXTERNAL           LSAME, DGET06, ZLANSP
 213 *     ..
 214 *     .. External Subroutines ..
 215       EXTERNAL           ALAERH, ALAHD, ALASUM, ZCOPY, ZERRSY, ZGET04,
 216      $                   ZLACPY, ZLARHS, ZLATB4, ZLATMS, ZLATSP, ZPPT05,
 217      $                   ZSPCON, ZSPRFS, ZSPT01, ZSPT02, ZSPT03, ZSPTRF,
 218      $                   ZSPTRI, ZSPTRS
 219 *     ..
 220 *     .. Intrinsic Functions ..
 221       INTRINSIC          MAX, MIN
 222 *     ..
 223 *     .. Scalars in Common ..
 224       LOGICAL            LERR, OK
 225       CHARACTER*32       SRNAMT
 226       INTEGER            INFOT, NUNIT
 227 *     ..
 228 *     .. Common blocks ..
 229       COMMON             / INFOC / INFOT, NUNIT, OK, LERR
 230       COMMON             / SRNAMC / SRNAMT
 231 *     ..
 232 *     .. Data statements ..
 233       DATA               ISEEDY / 1988, 1989, 1990, 1991 /
 234       DATA               UPLOS / 'U', 'L' /
 235 *     ..
 236 *     .. Executable Statements ..
 237 *
 238 *     Initialize constants and the random number seed.
 239 *
 240       PATH( 1: 1 ) = 'Zomplex precision'
 241       PATH( 2: 3 ) = 'SP'
 242       NRUN = 0
 243       NFAIL = 0
 244       NERRS = 0
 245       DO 10 I = 1, 4
 246          ISEED( I ) = ISEEDY( I )
 247    10 CONTINUE
 248 *
 249 *     Test the error exits
 250 *
 251       IF( TSTERR )
 252      $   CALL ZERRSY( PATH, NOUT )
 253       INFOT = 0
 254 *
 255 *     Do for each value of N in NVAL
 256 *
 257       DO 170 IN = 1, NN
 258          N = NVAL( IN )
 259          LDA = MAX( N, 1 )
 260          XTYPE = 'N'
 261          NIMAT = NTYPES
 262          IF( N.LE.0 )
 263      $      NIMAT = 1
 264 *
 265          DO 160 IMAT = 1, NIMAT
 266 *
 267 *           Do the tests only if DOTYPE( IMAT ) is true.
 268 *
 269             IF( .NOT.DOTYPE( IMAT ) )
 270      $         GO TO 160
 271 *
 272 *           Skip types 3, 4, 5, or 6 if the matrix size is too small.
 273 *
 274             ZEROT = IMAT.GE.3 .AND. IMAT.LE.6
 275             IF( ZEROT .AND. N.LT.IMAT-2 )
 276      $         GO TO 160
 277 *
 278 *           Do first for UPLO = 'U', then for UPLO = 'L'
 279 *
 280             DO 150 IUPLO = 1, 2
 281                UPLO = UPLOS( IUPLO )
 282                IF( LSAME( UPLO, 'U' ) ) THEN
 283                   PACKIT = 'C'
 284                ELSE
 285                   PACKIT = 'R'
 286                END IF
 287 *
 288                IF( IMAT.NE.NTYPES ) THEN
 289 *
 290 *                 Set up parameters with ZLATB4 and generate a test
 291 *                 matrix with ZLATMS.
 292 *
 293                   CALL ZLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM,
 294      $                         MODE, CNDNUM, DIST )
 295 *
 296                   SRNAMT = 'ZLATMS'
 297                   CALL ZLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
 298      $                         CNDNUM, ANORM, KL, KU, PACKIT, A, LDA,
 299      $                         WORK, INFO )
 300 *
 301 *                 Check error code from ZLATMS.
 302 *
 303                   IF( INFO.NE.0 ) THEN
 304                      CALL ALAERH( PATH, 'ZLATMS', INFO, 0, UPLO, N, N,
 305      $                            -1, -1, -1, IMAT, NFAIL, NERRS, NOUT )
 306                      GO TO 150
 307                   END IF
 308 *
 309 *                 For types 3-6, zero one or more rows and columns of
 310 *                 the matrix to test that INFO is returned correctly.
 311 *
 312                   IF( ZEROT ) THEN
 313                      IF( IMAT.EQ.3 ) THEN
 314                         IZERO = 1
 315                      ELSE IF( IMAT.EQ.4 ) THEN
 316                         IZERO = N
 317                      ELSE
 318                         IZERO = N / 2 + 1
 319                      END IF
 320 *
 321                      IF( IMAT.LT.6 ) THEN
 322 *
 323 *                       Set row and column IZERO to zero.
 324 *
 325                         IF( IUPLO.EQ.1 ) THEN
 326                            IOFF = ( IZERO-1 )*IZERO / 2
 327                            DO 20 I = 1, IZERO - 1
 328                               A( IOFF+I ) = ZERO
 329    20                      CONTINUE
 330                            IOFF = IOFF + IZERO
 331                            DO 30 I = IZERO, N
 332                               A( IOFF ) = ZERO
 333                               IOFF = IOFF + I
 334    30                      CONTINUE
 335                         ELSE
 336                            IOFF = IZERO
 337                            DO 40 I = 1, IZERO - 1
 338                               A( IOFF ) = ZERO
 339                               IOFF = IOFF + N - I
 340    40                      CONTINUE
 341                            IOFF = IOFF - IZERO
 342                            DO 50 I = IZERO, N
 343                               A( IOFF+I ) = ZERO
 344    50                      CONTINUE
 345                         END IF
 346                      ELSE
 347                         IF( IUPLO.EQ.1 ) THEN
 348 *
 349 *                          Set the first IZERO rows and columns to zero.
 350 *
 351                            IOFF = 0
 352                            DO 70 J = 1, N
 353                               I2 = MIN( J, IZERO )
 354                               DO 60 I = 1, I2
 355                                  A( IOFF+I ) = ZERO
 356    60                         CONTINUE
 357                               IOFF = IOFF + J
 358    70                      CONTINUE
 359                         ELSE
 360 *
 361 *                          Set the last IZERO rows and columns to zero.
 362 *
 363                            IOFF = 0
 364                            DO 90 J = 1, N
 365                               I1 = MAX( J, IZERO )
 366                               DO 80 I = I1, N
 367                                  A( IOFF+I ) = ZERO
 368    80                         CONTINUE
 369                               IOFF = IOFF + N - J
 370    90                      CONTINUE
 371                         END IF
 372                      END IF
 373                   ELSE
 374                      IZERO = 0
 375                   END IF
 376                ELSE
 377 *
 378 *                 Use a special block diagonal matrix to test alternate
 379 *                 code for the 2 x 2 blocks.
 380 *
 381                   CALL ZLATSP( UPLO, N, A, ISEED )
 382                END IF
 383 *
 384 *              Compute the L*D*L' or U*D*U' factorization of the matrix.
 385 *
 386                NPP = N*( N+1 ) / 2
 387                CALL ZCOPY( NPP, A, 1, AFAC, 1 )
 388                SRNAMT = 'ZSPTRF'
 389                CALL ZSPTRF( UPLO, N, AFAC, IWORK, INFO )
 390 *
 391 *              Adjust the expected value of INFO to account for
 392 *              pivoting.
 393 *
 394                K = IZERO
 395                IF( K.GT.0 ) THEN
 396   100             CONTINUE
 397                   IF( IWORK( K ).LT.0 ) THEN
 398                      IF( IWORK( K ).NE.-K ) THEN
 399                         K = -IWORK( K )
 400                         GO TO 100
 401                      END IF
 402                   ELSE IF( IWORK( K ).NE.K ) THEN
 403                      K = IWORK( K )
 404                      GO TO 100
 405                   END IF
 406                END IF
 407 *
 408 *              Check error code from ZSPTRF.
 409 *
 410                IF( INFO.NE.K )
 411      $            CALL ALAERH( PATH, 'ZSPTRF', INFO, K, UPLO, N, N, -1,
 412      $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
 413                IF( INFO.NE.0 ) THEN
 414                   TRFCON = .TRUE.
 415                ELSE
 416                   TRFCON = .FALSE.
 417                END IF
 418 *
 419 *+    TEST 1
 420 *              Reconstruct matrix from factors and compute residual.
 421 *
 422                CALL ZSPT01( UPLO, N, A, AFAC, IWORK, AINV, LDA, RWORK,
 423      $                      RESULT( 1 ) )
 424                NT = 1
 425 *
 426 *+    TEST 2
 427 *              Form the inverse and compute the residual.
 428 *
 429                IF( .NOT.TRFCON ) THEN
 430                   CALL ZCOPY( NPP, AFAC, 1, AINV, 1 )
 431                   SRNAMT = 'ZSPTRI'
 432                   CALL ZSPTRI( UPLO, N, AINV, IWORK, WORK, INFO )
 433 *
 434 *              Check error code from ZSPTRI.
 435 *
 436                   IF( INFO.NE.0 )
 437      $               CALL ALAERH( PATH, 'ZSPTRI', INFO, 0, UPLO, N, N,
 438      $                            -1, -1, -1, IMAT, NFAIL, NERRS, NOUT )
 439 *
 440                   CALL ZSPT03( UPLO, N, A, AINV, WORK, LDA, RWORK,
 441      $                         RCONDC, RESULT( 2 ) )
 442                   NT = 2
 443                END IF
 444 *
 445 *              Print information about the tests that did not pass
 446 *              the threshold.
 447 *
 448                DO 110 K = 1, NT
 449                   IF( RESULT( K ).GE.THRESH ) THEN
 450                      IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
 451      $                  CALL ALAHD( NOUT, PATH )
 452                      WRITE( NOUT, FMT = 9999 )UPLO, N, IMAT, K,
 453      $                  RESULT( K )
 454                      NFAIL = NFAIL + 1
 455                   END IF
 456   110          CONTINUE
 457                NRUN = NRUN + NT
 458 *
 459 *              Do only the condition estimate if INFO is not 0.
 460 *
 461                IF( TRFCON ) THEN
 462                   RCONDC = ZERO
 463                   GO TO 140
 464                END IF
 465 *
 466                DO 130 IRHS = 1, NNS
 467                   NRHS = NSVAL( IRHS )
 468 *
 469 *+    TEST 3
 470 *              Solve and compute residual for  A * X = B.
 471 *
 472                   SRNAMT = 'ZLARHS'
 473                   CALL ZLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
 474      $                         NRHS, A, LDA, XACT, LDA, B, LDA, ISEED,
 475      $                         INFO )
 476                   CALL ZLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
 477 *
 478                   SRNAMT = 'ZSPTRS'
 479                   CALL ZSPTRS( UPLO, N, NRHS, AFAC, IWORK, X, LDA,
 480      $                         INFO )
 481 *
 482 *              Check error code from ZSPTRS.
 483 *
 484                   IF( INFO.NE.0 )
 485      $               CALL ALAERH( PATH, 'ZSPTRS', INFO, 0, UPLO, N, N,
 486      $                            -1, -1, NRHS, IMAT, NFAIL, NERRS,
 487      $                            NOUT )
 488 *
 489                   CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
 490                   CALL ZSPT02( UPLO, N, NRHS, A, X, LDA, WORK, LDA,
 491      $                         RWORK, RESULT( 3 ) )
 492 *
 493 *+    TEST 4
 494 *              Check solution from generated exact solution.
 495 *
 496                   CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
 497      $                         RESULT( 4 ) )
 498 *
 499 *+    TESTS 5, 6, and 7
 500 *              Use iterative refinement to improve the solution.
 501 *
 502                   SRNAMT = 'ZSPRFS'
 503                   CALL ZSPRFS( UPLO, N, NRHS, A, AFAC, IWORK, B, LDA, X,
 504      $                         LDA, RWORK, RWORK( NRHS+1 ), WORK,
 505      $                         RWORK( 2*NRHS+1 ), INFO )
 506 *
 507 *              Check error code from ZSPRFS.
 508 *
 509                   IF( INFO.NE.0 )
 510      $               CALL ALAERH( PATH, 'ZSPRFS', INFO, 0, UPLO, N, N,
 511      $                            -1, -1, NRHS, IMAT, NFAIL, NERRS,
 512      $                            NOUT )
 513 *
 514                   CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
 515      $                         RESULT( 5 ) )
 516                   CALL ZPPT05( UPLO, N, NRHS, A, B, LDA, X, LDA, XACT,
 517      $                         LDA, RWORK, RWORK( NRHS+1 ),
 518      $                         RESULT( 6 ) )
 519 *
 520 *                 Print information about the tests that did not pass
 521 *                 the threshold.
 522 *
 523                   DO 120 K = 3, 7
 524                      IF( RESULT( K ).GE.THRESH ) THEN
 525                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
 526      $                     CALL ALAHD( NOUT, PATH )
 527                         WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS, IMAT,
 528      $                     K, RESULT( K )
 529                         NFAIL = NFAIL + 1
 530                      END IF
 531   120             CONTINUE
 532                   NRUN = NRUN + 5
 533   130          CONTINUE
 534 *
 535 *+    TEST 8
 536 *              Get an estimate of RCOND = 1/CNDNUM.
 537 *
 538   140          CONTINUE
 539                ANORM = ZLANSP( '1', UPLO, N, A, RWORK )
 540                SRNAMT = 'ZSPCON'
 541                CALL ZSPCON( UPLO, N, AFAC, IWORK, ANORM, RCOND, WORK,
 542      $                      INFO )
 543 *
 544 *              Check error code from ZSPCON.
 545 *
 546                IF( INFO.NE.0 )
 547      $            CALL ALAERH( PATH, 'ZSPCON', INFO, 0, UPLO, N, N, -1,
 548      $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
 549 *
 550                RESULT( 8 ) = DGET06( RCOND, RCONDC )
 551 *
 552 *              Print the test ratio if it is .GE. THRESH.
 553 *
 554                IF( RESULT( 8 ).GE.THRESH ) THEN
 555                   IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
 556      $               CALL ALAHD( NOUT, PATH )
 557                   WRITE( NOUT, FMT = 9999 )UPLO, N, IMAT, 8,
 558      $               RESULT( 8 )
 559                   NFAIL = NFAIL + 1
 560                END IF
 561                NRUN = NRUN + 1
 562   150       CONTINUE
 563   160    CONTINUE
 564   170 CONTINUE
 565 *
 566 *     Print a summary of the results.
 567 *
 568       CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
 569 *
 570  9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', type ', I2, ', test ',
 571      $      I2, ', ratio =', G12.5 )
 572  9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ',
 573      $      I2, ', test(', I2, ') =', G12.5 )
 574       RETURN
 575 *
 576 *     End of ZCHKSP
 577 *
 578       END