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