dpot05.f

Go to the documentation of this file.

      SUBROUTINE dpot05( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, XACT,
     $                   ldxact, ferr, berr, reslts )
*
*  -- LAPACK test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      CHARACTER          uplo
      INTEGER            lda, ldb, ldx, ldxact, n, nrhs
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   a( lda, * ), b( ldb, * ), berr( * ), ferr( * ),
     $                   reslts( * ), x( ldx, * ), xact( ldxact, * )
*     ..
*
*  Purpose
*  =======
*
*  DPOT05 tests the error bounds from iterative refinement for the
*  computed solution to a system of equations A*X = B, where A is a
*  symmetric n by n matrix.
*
*  RESLTS(1) = test of the error bound
*            = norm(X - XACT) / ( norm(X) * FERR )
*
*  A large value is returned if this ratio is not less than one.
*
*  RESLTS(2) = residual from the iterative refinement routine
*            = the maximum of BERR / ( (n+1)*EPS + (*) ), where
*              (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
*
*  Arguments
*  =========
*
*  UPLO    (input) CHARACTER*1
*          Specifies whether the upper or lower triangular part of the
*          symmetric matrix A is stored.
*          = 'U':  Upper triangular
*          = 'L':  Lower triangular
*
*  N       (input) INTEGER
*          The number of rows of the matrices X, B, and XACT, and the
*          order of the matrix A.  N >= 0.
*
*  NRHS    (input) INTEGER
*          The number of columns of the matrices X, B, and XACT.
*          NRHS >= 0.
*
*  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
*          The symmetric matrix A.  If UPLO = 'U', the leading n by n
*          upper triangular part of A contains the upper triangular part
*          of the matrix A, and the strictly lower triangular part of A
*          is not referenced.  If UPLO = 'L', the leading n by n lower
*          triangular part of A contains the lower triangular part of
*          the matrix A, and the strictly upper triangular part of A is
*          not referenced.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A.  LDA >= max(1,N).
*
*  B       (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
*          The right hand side vectors for the system of linear
*          equations.
*
*  LDB     (input) INTEGER
*          The leading dimension of the array B.  LDB >= max(1,N).
*
*  X       (input) DOUBLE PRECISION array, dimension (LDX,NRHS)
*          The computed solution vectors.  Each vector is stored as a
*          column of the matrix X.
*
*  LDX     (input) INTEGER
*          The leading dimension of the array X.  LDX >= max(1,N).
*
*  XACT    (input) DOUBLE PRECISION array, dimension (LDX,NRHS)
*          The exact solution vectors.  Each vector is stored as a
*          column of the matrix XACT.
*
*  LDXACT  (input) INTEGER
*          The leading dimension of the array XACT.  LDXACT >= max(1,N).
*
*  FERR    (input) DOUBLE PRECISION array, dimension (NRHS)
*          The estimated forward error bounds for each solution vector
*          X.  If XTRUE is the true solution, FERR bounds the magnitude
*          of the largest entry in (X - XTRUE) divided by the magnitude
*          of the largest entry in X.
*
*  BERR    (input) DOUBLE PRECISION array, dimension (NRHS)
*          The componentwise relative backward error of each solution
*          vector (i.e., the smallest relative change in any entry of A
*          or B that makes X an exact solution).
*
*  RESLTS  (output) DOUBLE PRECISION array, dimension (2)
*          The maximum over the NRHS solution vectors of the ratios:
*          RESLTS(1) = norm(X - XACT) / ( norm(X) * FERR )
*          RESLTS(2) = BERR / ( (n+1)*EPS + (*) )
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   zero, one
      parameter( zero = 0.0d+0, one = 1.0d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            upper
      INTEGER            i, imax, j, k
      DOUBLE PRECISION   axbi, diff, eps, errbnd, ovfl, tmp, unfl, xnorm
*     ..
*     .. External Functions ..
      LOGICAL            lsame
      INTEGER            idamax
      DOUBLE PRECISION   dlamch
      EXTERNAL           lsame, idamax, dlamch
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min
*     ..
*     .. Executable Statements ..
*
*     Quick exit if N = 0 or NRHS = 0.
*
      IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
         reslts( 1 ) = zero
         reslts( 2 ) = zero
         return
      END IF
*
      eps = dlamch( 'Epsilon' )
      unfl = dlamch( 'Safe minimum' )
      ovfl = one / unfl
      upper = lsame( uplo, 'U' )
*
*     Test 1:  Compute the maximum of
*        norm(X - XACT) / ( norm(X) * FERR )
*     over all the vectors X and XACT using the infinity-norm.
*
      errbnd = zero
      DO 30 j = 1, nrhs
         imax = idamax( n, x( 1, j ), 1 )
         xnorm = max( abs( x( imax, j ) ), unfl )
         diff = zero
         DO 10 i = 1, n
            diff = max( diff, abs( x( i, j )-xact( i, j ) ) )
   10    continue
*
         IF( xnorm.GT.one ) THEN
            go to 20
         ELSE IF( diff.LE.ovfl*xnorm ) THEN
            go to 20
         ELSE
            errbnd = one / eps
            go to 30
         END IF
*
   20    continue
         IF( diff / xnorm.LE.ferr( j ) ) THEN
            errbnd = max( errbnd, ( diff / xnorm ) / ferr( j ) )
         ELSE
            errbnd = one / eps
         END IF
   30 continue
      reslts( 1 ) = errbnd
*
*     Test 2:  Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where
*     (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
*
      DO 90 k = 1, nrhs
         DO 80 i = 1, n
            tmp = abs( b( i, k ) )
            IF( upper ) THEN
               DO 40 j = 1, i
                  tmp = tmp + abs( a( j, i ) )*abs( x( j, k ) )
   40          continue
               DO 50 j = i + 1, n
                  tmp = tmp + abs( a( i, j ) )*abs( x( j, k ) )
   50          continue
            ELSE
               DO 60 j = 1, i - 1
                  tmp = tmp + abs( a( i, j ) )*abs( x( j, k ) )
   60          continue
               DO 70 j = i, n
                  tmp = tmp + abs( a( j, i ) )*abs( x( j, k ) )
   70          continue
            END IF
            IF( i.EQ.1 ) THEN
               axbi = tmp
            ELSE
               axbi = min( axbi, tmp )
            END IF
   80    continue
         tmp = berr( k ) / ( ( n+1 )*eps+( n+1 )*unfl /
     $         max( axbi, ( n+1 )*unfl ) )
         IF( k.EQ.1 ) THEN
            reslts( 2 ) = tmp
         ELSE
            reslts( 2 ) = max( reslts( 2 ), tmp )
         END IF
   90 continue
*
      return
*
*     End of DPOT05
*
      END