167
  168
  169
  170
  171
  172
  173      LOGICAL            TSTERR
  174      INTEGER            NMAX, NN, NOUT, NRHS
  175      DOUBLE PRECISION   THRESH
  176
  177
  178      LOGICAL            DOTYPE( * )
  179      INTEGER            IWORK( * ), NVAL( * )
  180      DOUBLE PRECISION   A( * ), AFAC( * ), ASAV( * ), B( * ),
  181     $                   BSAV( * ), RWORK( * ), S( * ), WORK( * ),
  182     $                   X( * ), XACT( * )
  183
  184
  185
  186
  187
  188      DOUBLE PRECISION   ONE, ZERO
  189      parameter( one = 1.0d+0, zero = 0.0d+0 )
  190      INTEGER            NTYPES
  191      parameter( ntypes = 9 )
  192      INTEGER            NTESTS
  193      parameter( ntests = 6 )
  194
  195
  196      LOGICAL            EQUIL, NOFACT, PREFAC, ZEROT
  197      CHARACTER          DIST, EQUED, FACT, TYPE, UPLO, XTYPE
  198      CHARACTER*3        PATH
  199      INTEGER            I, IEQUED, IFACT, IMAT, IN, INFO, IOFF, IUPLO,
  200     $                   IZERO, K, K1, KL, KU, LDA, MODE, N, NB, NBMIN,
  201     $                   NERRS, NFACT, NFAIL, NIMAT, NRUN, NT,
  202     $                   N_ERR_BNDS
  203      DOUBLE PRECISION   AINVNM, AMAX, ANORM, CNDNUM, RCOND, RCONDC,
  204     $                   ROLDC, SCOND, RPVGRW_SVXX
  205
  206
  207      CHARACTER          EQUEDS( 2 ), FACTS( 3 ), UPLOS( 2 )
  208      INTEGER            ISEED( 4 ), ISEEDY( 4 )
  209      DOUBLE PRECISION   RESULT( NTESTS ), BERR( NRHS ),
  210     $                   ERRBNDS_N( NRHS, 3 ), ERRBNDS_C( NRHS, 3 )
  211
  212
  213      LOGICAL            LSAME
  214      DOUBLE PRECISION   DGET06, DLANSY
  216
  217
  222
  223
  224      INTRINSIC          max
  225
  226
  227      LOGICAL            LERR, OK
  228      CHARACTER*32       SRNAMT
  229      INTEGER            INFOT, NUNIT
  230
  231
  232      COMMON             / infoc / infot, nunit, ok, lerr
  233      COMMON             / srnamc / srnamt
  234
  235
  236      DATA               iseedy / 1988, 1989, 1990, 1991 /
  237      DATA               uplos / 'U', 'L' /
  238      DATA               facts / 'F', 'N', 'E' /
  239      DATA               equeds / 'N', 'Y' /
  240
  241
  242
  243
  244
  245      path( 1: 1 ) = 'Double precision'
  246      path( 2: 3 ) = 'PO'
  247      nrun = 0
  248      nfail = 0
  249      nerrs = 0
  250      DO 10 i = 1, 4
  251         iseed( i ) = iseedy( i )
  252   10 CONTINUE
  253
  254
  255
  256      IF( tsterr )
  257     $   
CALL derrvx( path, nout )
 
  258      infot = 0
  259
  260
  261
  262      nb = 1
  263      nbmin = 2
  266
  267
  268
  269      DO 130 in = 1, nn
  270         n = nval( in )
  271         lda = max( n, 1 )
  272         xtype = 'N'
  273         nimat = ntypes
  274         IF( n.LE.0 )
  275     $      nimat = 1
  276
  277         DO 120 imat = 1, nimat
  278
  279
  280
  281            IF( .NOT.dotype( imat ) )
  282     $         GO TO 120
  283
  284
  285
  286            zerot = imat.GE.3 .AND. imat.LE.5
  287            IF( zerot .AND. n.LT.imat-2 )
  288     $         GO TO 120
  289
  290
  291
  292            DO 110 iuplo = 1, 2
  293               uplo = uplos( iuplo )
  294
  295
  296
  297
  298               CALL dlatb4( path, imat, n, n, 
TYPE, KL, KU, ANORM, MODE,
 
  299     $                      CNDNUM, DIST )
  300
  301               srnamt = 'DLATMS'
  302               CALL dlatms( n, n, dist, iseed, 
TYPE, RWORK, MODE,
 
  303     $                      CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK,
  304     $                      INFO )
  305
  306
  307
  308               IF( info.NE.0 ) THEN
  309                  CALL alaerh( path, 
'DLATMS', info, 0, uplo, n, n, -1,
 
  310     $                         -1, -1, imat, nfail, nerrs, nout )
  311                  GO TO 110
  312               END IF
  313
  314
  315
  316
  317               IF( zerot ) THEN
  318                  IF( imat.EQ.3 ) THEN
  319                     izero = 1
  320                  ELSE IF( imat.EQ.4 ) THEN
  321                     izero = n
  322                  ELSE
  323                     izero = n / 2 + 1
  324                  END IF
  325                  ioff = ( izero-1 )*lda
  326
  327
  328
  329                  IF( iuplo.EQ.1 ) THEN
  330                     DO 20 i = 1, izero - 1
  331                        a( ioff+i ) = zero
  332   20                CONTINUE
  333                     ioff = ioff + izero
  334                     DO 30 i = izero, n
  335                        a( ioff ) = zero
  336                        ioff = ioff + lda
  337   30                CONTINUE
  338                  ELSE
  339                     ioff = izero
  340                     DO 40 i = 1, izero - 1
  341                        a( ioff ) = zero
  342                        ioff = ioff + lda
  343   40                CONTINUE
  344                     ioff = ioff - izero
  345                     DO 50 i = izero, n
  346                        a( ioff+i ) = zero
  347   50                CONTINUE
  348                  END IF
  349               ELSE
  350                  izero = 0
  351               END IF
  352
  353
  354
  355               CALL dlacpy( uplo, n, n, a, lda, asav, lda )
 
  356
  357               DO 100 iequed = 1, 2
  358                  equed = equeds( iequed )
  359                  IF( iequed.EQ.1 ) THEN
  360                     nfact = 3
  361                  ELSE
  362                     nfact = 1
  363                  END IF
  364
  365                  DO 90 ifact = 1, nfact
  366                     fact = facts( ifact )
  367                     prefac = 
lsame( fact, 
'F' )
 
  368                     nofact = 
lsame( fact, 
'N' )
 
  369                     equil = 
lsame( fact, 
'E' )
 
  370
  371                     IF( zerot ) THEN
  372                        IF( prefac )
  373     $                     GO TO 90
  374                        rcondc = zero
  375
  376                     ELSE IF( .NOT.
lsame( fact, 
'N' ) ) 
THEN 
  377
  378
  379
  380
  381
  382
  383                        CALL dlacpy( uplo, n, n, asav, lda, afac, lda )
 
  384                        IF( equil .OR. iequed.GT.1 ) THEN
  385
  386
  387
  388
  389                           CALL dpoequ( n, afac, lda, s, scond, amax,
 
  390     $                                  info )
  391                           IF( info.EQ.0 .AND. n.GT.0 ) THEN
  392                              IF( iequed.GT.1 )
  393     $                           scond = zero
  394
  395
  396
  397                              CALL dlaqsy( uplo, n, afac, lda, s, scond,
 
  398     $                                     amax, equed )
  399                           END IF
  400                        END IF
  401
  402
  403
  404
  405                        IF( equil )
  406     $                     roldc = rcondc
  407
  408
  409
  410                        anorm = 
dlansy( 
'1', uplo, n, afac, lda, rwork )
 
  411
  412
  413
  414                        CALL dpotrf( uplo, n, afac, lda, info )
 
  415
  416
  417
  418                        CALL dlacpy( uplo, n, n, afac, lda, a, lda )
 
  419                        CALL dpotri( uplo, n, a, lda, info )
 
  420
  421
  422
  423                        ainvnm = 
dlansy( 
'1', uplo, n, a, lda, rwork )
 
  424                        IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
  425                           rcondc = one
  426                        ELSE
  427                           rcondc = ( one / anorm ) / ainvnm
  428                        END IF
  429                     END IF
  430
  431
  432
  433                     CALL dlacpy( uplo, n, n, asav, lda, a, lda )
 
  434
  435
  436
  437                     srnamt = 'DLARHS'
  438                     CALL dlarhs( path, xtype, uplo, 
' ', n, n, kl, ku,
 
  439     $                            nrhs, a, lda, xact, lda, b, lda,
  440     $                            iseed, info )
  441                     xtype = 'C'
  442                     CALL dlacpy( 
'Full', n, nrhs, b, lda, bsav, lda )
 
  443
  444                     IF( nofact ) THEN
  445
  446
  447
  448
  449
  450
  451                        CALL dlacpy( uplo, n, n, a, lda, afac, lda )
 
  452                        CALL dlacpy( 
'Full', n, nrhs, b, lda, x, lda )
 
  453
  454                        srnamt = 'DPOSV '
  455                        CALL dposv( uplo, n, nrhs, afac, lda, x, lda,
 
  456     $                              info )
  457
  458
  459
  460                        IF( info.NE.izero ) THEN
  461                           CALL alaerh( path, 
'DPOSV ', info, izero,
 
  462     $                                  uplo, n, n, -1, -1, nrhs, imat,
  463     $                                  nfail, nerrs, nout )
  464                           GO TO 70
  465                        ELSE IF( info.NE.0 ) THEN
  466                           GO TO 70
  467                        END IF
  468
  469
  470
  471
  472                        CALL dpot01( uplo, n, a, lda, afac, lda, rwork,
 
  473     $                               result( 1 ) )
  474
  475
  476
  477                        CALL dlacpy( 
'Full', n, nrhs, b, lda, work,
 
  478     $                               lda )
  479                        CALL dpot02( uplo, n, nrhs, a, lda, x, lda,
 
  480     $                               work, lda, rwork, result( 2 ) )
  481
  482
  483
  484                        CALL dget04( n, nrhs, x, lda, xact, lda, rcondc,
 
  485     $                               result( 3 ) )
  486                        nt = 3
  487
  488
  489
  490
  491                        DO 60 k = 1, nt
  492                           IF( result( k ).GE.thresh ) THEN
  493                              IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
  494     $                           
CALL aladhd( nout, path )
 
  495                              WRITE( nout, fmt = 9999 )'DPOSV ', uplo,
  496     $                           n, imat, k, result( k )
  497                              nfail = nfail + 1
  498                           END IF
  499   60                   CONTINUE
  500                        nrun = nrun + nt
  501   70                   CONTINUE
  502                     END IF
  503
  504
  505
  506                     IF( .NOT.prefac )
  507     $                  
CALL dlaset( uplo, n, n, zero, zero, afac, lda )
 
  508                     CALL dlaset( 
'Full', n, nrhs, zero, zero, x, lda )
 
  509                     IF( iequed.GT.1 .AND. n.GT.0 ) THEN
  510
  511
  512
  513
  514                        CALL dlaqsy( uplo, n, a, lda, s, scond, amax,
 
  515     $                               equed )
  516                     END IF
  517
  518
  519
  520
  521                     srnamt = 'DPOSVX'
  522                     CALL dposvx( fact, uplo, n, nrhs, a, lda, afac,
 
  523     $                            lda, equed, s, b, lda, x, lda, rcond,
  524     $                            rwork, rwork( nrhs+1 ), work, iwork,
  525     $                            info )
  526
  527
  528
  529                     IF( info.NE.izero ) THEN
  530                        CALL alaerh( path, 
'DPOSVX', info, izero,
 
  531     $                               fact // uplo, n, n, -1, -1, nrhs,
  532     $                               imat, nfail, nerrs, nout )
  533                        GO TO 90
  534                     END IF
  535
  536                     IF( info.EQ.0 ) THEN
  537                        IF( .NOT.prefac ) THEN
  538
  539
  540
  541
  542                           CALL dpot01( uplo, n, a, lda, afac, lda,
 
  543     $                                  rwork( 2*nrhs+1 ), result( 1 ) )
  544                           k1 = 1
  545                        ELSE
  546                           k1 = 2
  547                        END IF
  548
  549
  550
  551                        CALL dlacpy( 
'Full', n, nrhs, bsav, lda, work,
 
  552     $                               lda )
  553                        CALL dpot02( uplo, n, nrhs, asav, lda, x, lda,
 
  554     $                               work, lda, rwork( 2*nrhs+1 ),
  555     $                               result( 2 ) )
  556
  557
  558
  559                        IF( nofact .OR. ( prefac .AND. 
lsame( equed,
 
  560     $                      'N' ) ) ) THEN
  561                           CALL dget04( n, nrhs, x, lda, xact, lda,
 
  562     $                                  rcondc, result( 3 ) )
  563                        ELSE
  564                           CALL dget04( n, nrhs, x, lda, xact, lda,
 
  565     $                                  roldc, result( 3 ) )
  566                        END IF
  567
  568
  569
  570
  571                        CALL dpot05( uplo, n, nrhs, asav, lda, b, lda,
 
  572     $                               x, lda, xact, lda, rwork,
  573     $                               rwork( nrhs+1 ), result( 4 ) )
  574                     ELSE
  575                        k1 = 6
  576                     END IF
  577
  578
  579
  580
  581                     result( 6 ) = 
dget06( rcond, rcondc )
 
  582
  583
  584
  585
  586                     DO 80 k = k1, 6
  587                        IF( result( k ).GE.thresh ) THEN
  588                           IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
  589     $                        
CALL aladhd( nout, path )
 
  590                           IF( prefac ) THEN
  591                              WRITE( nout, fmt = 9997 )'DPOSVX', fact,
  592     $                           uplo, n, equed, imat, k, result( k )
  593                           ELSE
  594                              WRITE( nout, fmt = 9998 )'DPOSVX', fact,
  595     $                           uplo, n, imat, k, result( k )
  596                           END IF
  597                           nfail = nfail + 1
  598                        END IF
  599   80                CONTINUE
  600                     nrun = nrun + 7 - k1
  601
  602
  603
  604
  605
  606                     CALL dlacpy( 
'Full', n, n, asav, lda, a, lda )
 
  607                     CALL dlacpy( 
'Full', n, nrhs, bsav, lda, b, lda )
 
  608 
  609                     IF( .NOT.prefac )
  610     $                  
CALL dlaset( uplo, n, n, zero, zero, afac, lda )
 
  611                     CALL dlaset( 
'Full', n, nrhs, zero, zero, x, lda )
 
  612                     IF( iequed.GT.1 .AND. n.GT.0 ) THEN
  613
  614
  615
  616
  617                        CALL dlaqsy( uplo, n, a, lda, s, scond, amax,
 
  618     $                               equed )
  619                     END IF
  620
  621
  622
  623
  624                     srnamt = 'DPOSVXX'
  625                     n_err_bnds = 3
  626                     CALL dposvxx( fact, uplo, n, nrhs, a, lda, afac,
 
  627     $                    lda, equed, s, b, lda, x,
  628     $                    lda, rcond, rpvgrw_svxx, berr, n_err_bnds,
  629     $                    errbnds_n, errbnds_c, 0, zero, work,
  630     $                    iwork, info )
  631
  632
  633
  634                     IF( info.EQ.n+1 ) GOTO 90
  635                     IF( info.NE.izero ) THEN
  636                        CALL alaerh( path, 
'DPOSVXX', info, izero,
 
  637     $                               fact // uplo, n, n, -1, -1, nrhs,
  638     $                               imat, nfail, nerrs, nout )
  639                        GO TO 90
  640                     END IF
  641
  642                     IF( info.EQ.0 ) THEN
  643                        IF( .NOT.prefac ) THEN
  644
  645
  646
  647
  648                           CALL dpot01( uplo, n, a, lda, afac, lda,
 
  649     $                                  rwork( 2*nrhs+1 ), result( 1 ) )
  650                           k1 = 1
  651                        ELSE
  652                           k1 = 2
  653                        END IF
  654
  655
  656
  657                        CALL dlacpy( 
'Full', n, nrhs, bsav, lda, work,
 
  658     $                               lda )
  659                        CALL dpot02( uplo, n, nrhs, asav, lda, x, lda,
 
  660     $                               work, lda, rwork( 2*nrhs+1 ),
  661     $                               result( 2 ) )
  662
  663
  664
  665                        IF( nofact .OR. ( prefac .AND. 
lsame( equed,
 
  666     $                      'N' ) ) ) THEN
  667                           CALL dget04( n, nrhs, x, lda, xact, lda,
 
  668     $                                  rcondc, result( 3 ) )
  669                        ELSE
  670                           CALL dget04( n, nrhs, x, lda, xact, lda,
 
  671     $                                  roldc, result( 3 ) )
  672                        END IF
  673
  674
  675
  676
  677                        CALL dpot05( uplo, n, nrhs, asav, lda, b, lda,
 
  678     $                               x, lda, xact, lda, rwork,
  679     $                               rwork( nrhs+1 ), result( 4 ) )
  680                     ELSE
  681                        k1 = 6
  682                     END IF
  683
  684
  685
  686
  687                     result( 6 ) = 
dget06( rcond, rcondc )
 
  688
  689
  690
  691
  692                     DO 85 k = k1, 6
  693                        IF( result( k ).GE.thresh ) THEN
  694                           IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
  695     $                        
CALL aladhd( nout, path )
 
  696                           IF( prefac ) THEN
  697                              WRITE( nout, fmt = 9997 )'DPOSVXX', fact,
  698     $                           uplo, n, equed, imat, k, result( k )
  699                           ELSE
  700                              WRITE( nout, fmt = 9998 )'DPOSVXX', fact,
  701     $                           uplo, n, imat, k, result( k )
  702                           END IF
  703                           nfail = nfail + 1
  704                        END IF
  705   85                CONTINUE
  706                     nrun = nrun + 7 - k1
  707   90             CONTINUE
  708  100          CONTINUE
  709  110       CONTINUE
  710  120    CONTINUE
  711  130 CONTINUE
  712
  713
  714
  715      CALL alasvm( path, nout, nfail, nrun, nerrs )
 
  716
  717 
  718
  719 
  721 
  722 9999 FORMAT( 1x, a, ', UPLO=''', a1, ''', N =', i5, ', type ', i1,
  723     $      ', test(', i1, ')=', g12.5 )
  724 9998 FORMAT( 1x, a, ', FACT=''', a1, ''', UPLO=''', a1, ''', N=', i5,
  725     $      ', type ', i1, ', test(', i1, ')=', g12.5 )
  726 9997 FORMAT( 1x, a, ', FACT=''', a1, ''', UPLO=''', a1, ''', N=', i5,
  727     $      ', EQUED=''', a1, ''', type ', i1, ', test(', i1, ') =',
  728     $      g12.5 )
  729      RETURN
  730
  731
  732
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
subroutine dlarhs(path, xtype, uplo, trans, m, n, kl, ku, nrhs, a, lda, x, ldx, b, ldb, iseed, info)
DLARHS
subroutine xlaenv(ispec, nvalue)
XLAENV
subroutine aladhd(iounit, path)
ALADHD
subroutine alaerh(path, subnam, info, infoe, opts, m, n, kl, ku, n5, imat, nfail, nerrs, nout)
ALAERH
subroutine debchvxx(thresh, path)
DEBCHVXX
subroutine derrvx(path, nunit)
DERRVX
subroutine dget04(n, nrhs, x, ldx, xact, ldxact, rcond, resid)
DGET04
double precision function dget06(rcond, rcondc)
DGET06
subroutine dlatb4(path, imat, m, n, type, kl, ku, anorm, mode, cndnum, dist)
DLATB4
subroutine dlatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
DLATMS
subroutine dpot01(uplo, n, a, lda, afac, ldafac, rwork, resid)
DPOT01
subroutine dpot02(uplo, n, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
DPOT02
subroutine dpot05(uplo, n, nrhs, a, lda, b, ldb, x, ldx, xact, ldxact, ferr, berr, reslts)
DPOT05
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
double precision function dlansy(norm, uplo, n, a, lda, work)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
subroutine dlaqsy(uplo, n, a, lda, s, scond, amax, equed)
DLAQSY scales a symmetric/Hermitian matrix, using scaling factors computed by spoequ.
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
logical function lsame(ca, cb)
LSAME
subroutine dpoequ(n, a, lda, s, scond, amax, info)
DPOEQU
subroutine dposv(uplo, n, nrhs, a, lda, b, ldb, info)
DPOSV computes the solution to system of linear equations A * X = B for PO matrices
subroutine dposvx(fact, uplo, n, nrhs, a, lda, af, ldaf, equed, s, b, ldb, x, ldx, rcond, ferr, berr, work, iwork, info)
DPOSVX computes the solution to system of linear equations A * X = B for PO matrices
subroutine dposvxx(fact, uplo, n, nrhs, a, lda, af, ldaf, equed, s, b, ldb, x, ldx, rcond, rpvgrw, berr, n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params, work, iwork, info)
DPOSVXX computes the solution to system of linear equations A * X = B for PO matrices
subroutine dpotrf(uplo, n, a, lda, info)
DPOTRF
subroutine dpotri(uplo, n, a, lda, info)
DPOTRI