cppt05.f

Go to the documentation of this file.

      SUBROUTINE cppt05( UPLO, N, NRHS, AP, 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            ldb, ldx, ldxact, n, nrhs
*     ..
*     .. Array Arguments ..
      REAL               berr( * ), ferr( * ), reslts( * )
      COMPLEX            ap( * ), b( ldb, * ), x( ldx, * ),
     $                   xact( ldxact, * )
*     ..
*
*  Purpose
*  =======
*
*  CPPT05 tests the error bounds from iterative refinement for the
*  computed solution to a system of equations A*X = B, where A is a
*  Hermitian matrix in packed storage format.
*
*  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
*          Hermitian 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.
*
*  AP      (input) COMPLEX array, dimension (N*(N+1)/2)
*          The upper or lower triangle of the Hermitian matrix A, packed
*          columnwise in a linear array.  The j-th column of A is stored
*          in the array AP as follows:
*          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
*          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
*
*  B       (input) COMPLEX 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) COMPLEX 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) COMPLEX 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) REAL 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) REAL 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) REAL 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 ..
      REAL               zero, one
      parameter( zero = 0.0e+0, one = 1.0e+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            upper
      INTEGER            i, imax, j, jc, k
      REAL               axbi, diff, eps, errbnd, ovfl, tmp, unfl, xnorm
      COMPLEX            zdum
*     ..
*     .. External Functions ..
      LOGICAL            lsame
      INTEGER            icamax
      REAL               slamch
      EXTERNAL           lsame, icamax, slamch
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, aimag, max, min, real
*     ..
*     .. Statement Functions ..
      REAL               cabs1
*     ..
*     .. Statement Function definitions ..
      cabs1( zdum ) = abs( REAL( ZDUM ) ) + abs( aimag( zdum ) )
*     ..
*     .. 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 = slamch( 'Epsilon' )
      unfl = slamch( '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 = icamax( n, x( 1, j ), 1 )
         xnorm = max( cabs1( x( imax, j ) ), unfl )
         diff = zero
         DO 10 i = 1, n
            diff = max( diff, cabs1( 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 = cabs1( b( i, k ) )
            IF( upper ) THEN
               jc = ( ( i-1 )*i ) / 2
               DO 40 j = 1, i - 1
                  tmp = tmp + cabs1( ap( jc+j ) )*cabs1( x( j, k ) )
   40          continue
               tmp = tmp + abs( REAL( AP( JC+I ) ) )*cabs1( x( i, k ) )
               jc = jc + i + i
               DO 50 j = i + 1, n
                  tmp = tmp + cabs1( ap( jc ) )*cabs1( x( j, k ) )
                  jc = jc + j
   50          continue
            ELSE
               jc = i
               DO 60 j = 1, i - 1
                  tmp = tmp + cabs1( ap( jc ) )*cabs1( x( j, k ) )
                  jc = jc + n - j
   60          continue
               tmp = tmp + abs( REAL( AP( JC ) ) )*cabs1( x( i, k ) )
               DO 70 j = i + 1, n
                  tmp = tmp + cabs1( ap( jc+j-i ) )*cabs1( 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 CPPT05
*
      END