LAPACK 3.3.0

ctbt03.f

Go to the documentation of this file.
00001       SUBROUTINE CTBT03( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB,
00002      $                   SCALE, CNORM, TSCAL, X, LDX, B, LDB, WORK,
00003      $                   RESID )
00004 *
00005 *  -- LAPACK test routine (version 3.1) --
00006 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          DIAG, TRANS, UPLO
00011       INTEGER            KD, LDAB, LDB, LDX, N, NRHS
00012       REAL               RESID, SCALE, TSCAL
00013 *     ..
00014 *     .. Array Arguments ..
00015       REAL               CNORM( * )
00016       COMPLEX            AB( LDAB, * ), B( LDB, * ), WORK( * ),
00017      $                   X( LDX, * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  CTBT03 computes the residual for the solution to a scaled triangular
00024 *  system of equations  A*x = s*b,  A**T *x = s*b,  or  A**H *x = s*b
00025 *  when A is a triangular band matrix.  Here A**T  denotes the transpose
00026 *  of A, A**H denotes the conjugate transpose of A, s is a scalar, and
00027 *  x and b are N by NRHS matrices.  The test ratio is the maximum over
00028 *  the number of right hand sides of
00029 *     norm(s*b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
00030 *  where op(A) denotes A, A**T, or A**H, and EPS is the machine epsilon.
00031 *
00032 *  Arguments
00033 *  =========
00034 *
00035 *  UPLO    (input) CHARACTER*1
00036 *          Specifies whether the matrix A is upper or lower triangular.
00037 *          = 'U':  Upper triangular
00038 *          = 'L':  Lower triangular
00039 *
00040 *  TRANS   (input) CHARACTER*1
00041 *          Specifies the operation applied to A.
00042 *          = 'N':  A *x = s*b     (No transpose)
00043 *          = 'T':  A**T *x = s*b  (Transpose)
00044 *          = 'C':  A**H *x = s*b  (Conjugate transpose)
00045 *
00046 *  DIAG    (input) CHARACTER*1
00047 *          Specifies whether or not the matrix A is unit triangular.
00048 *          = 'N':  Non-unit triangular
00049 *          = 'U':  Unit triangular
00050 *
00051 *  N       (input) INTEGER
00052 *          The order of the matrix A.  N >= 0.
00053 *
00054 *  KD      (input) INTEGER
00055 *          The number of superdiagonals or subdiagonals of the
00056 *          triangular band matrix A.  KD >= 0.
00057 *
00058 *  NRHS    (input) INTEGER
00059 *          The number of right hand sides, i.e., the number of columns
00060 *          of the matrices X and B.  NRHS >= 0.
00061 *
00062 *  AB      (input) COMPLEX array, dimension (LDAB,N)
00063 *          The upper or lower triangular band matrix A, stored in the
00064 *          first kd+1 rows of the array. The j-th column of A is stored
00065 *          in the j-th column of the array AB as follows:
00066 *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
00067 *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
00068 *
00069 *  LDAB    (input) INTEGER
00070 *          The leading dimension of the array AB.  LDAB >= KD+1.
00071 *
00072 *  SCALE   (input) REAL
00073 *          The scaling factor s used in solving the triangular system.
00074 *
00075 *  CNORM   (input) REAL array, dimension (N)
00076 *          The 1-norms of the columns of A, not counting the diagonal.
00077 *
00078 *  TSCAL   (input) REAL
00079 *          The scaling factor used in computing the 1-norms in CNORM.
00080 *          CNORM actually contains the column norms of TSCAL*A.
00081 *
00082 *  X       (input) COMPLEX array, dimension (LDX,NRHS)
00083 *          The computed solution vectors for the system of linear
00084 *          equations.
00085 *
00086 *  LDX     (input) INTEGER
00087 *          The leading dimension of the array X.  LDX >= max(1,N).
00088 *
00089 *  B       (input) COMPLEX array, dimension (LDB,NRHS)
00090 *          The right hand side vectors for the system of linear
00091 *          equations.
00092 *
00093 *  LDB     (input) INTEGER
00094 *          The leading dimension of the array B.  LDB >= max(1,N).
00095 *
00096 *  WORK    (workspace) COMPLEX array, dimension (N)
00097 *
00098 *  RESID   (output) REAL
00099 *          The maximum over the number of right hand sides of
00100 *          norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
00101 *
00102 *  =====================================================================
00103 *
00104 *
00105 *     .. Parameters ..
00106       REAL               ONE, ZERO
00107       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00108 *     ..
00109 *     .. Local Scalars ..
00110       INTEGER            IX, J
00111       REAL               EPS, ERR, SMLNUM, TNORM, XNORM, XSCAL
00112 *     ..
00113 *     .. External Functions ..
00114       LOGICAL            LSAME
00115       INTEGER            ICAMAX
00116       REAL               SLAMCH
00117       EXTERNAL           LSAME, ICAMAX, SLAMCH
00118 *     ..
00119 *     .. External Subroutines ..
00120       EXTERNAL           CAXPY, CCOPY, CSSCAL, CTBMV
00121 *     ..
00122 *     .. Intrinsic Functions ..
00123       INTRINSIC          ABS, CMPLX, MAX, REAL
00124 *     ..
00125 *     .. Executable Statements ..
00126 *
00127 *     Quick exit if N = 0
00128 *
00129       IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
00130          RESID = ZERO
00131          RETURN
00132       END IF
00133       EPS = SLAMCH( 'Epsilon' )
00134       SMLNUM = SLAMCH( 'Safe minimum' )
00135 *
00136 *     Compute the norm of the triangular matrix A using the column
00137 *     norms already computed by CLATBS.
00138 *
00139       TNORM = ZERO
00140       IF( LSAME( DIAG, 'N' ) ) THEN
00141          IF( LSAME( UPLO, 'U' ) ) THEN
00142             DO 10 J = 1, N
00143                TNORM = MAX( TNORM, TSCAL*ABS( AB( KD+1, J ) )+
00144      $                 CNORM( J ) )
00145    10       CONTINUE
00146          ELSE
00147             DO 20 J = 1, N
00148                TNORM = MAX( TNORM, TSCAL*ABS( AB( 1, J ) )+CNORM( J ) )
00149    20       CONTINUE
00150          END IF
00151       ELSE
00152          DO 30 J = 1, N
00153             TNORM = MAX( TNORM, TSCAL+CNORM( J ) )
00154    30    CONTINUE
00155       END IF
00156 *
00157 *     Compute the maximum over the number of right hand sides of
00158 *        norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
00159 *
00160       RESID = ZERO
00161       DO 40 J = 1, NRHS
00162          CALL CCOPY( N, X( 1, J ), 1, WORK, 1 )
00163          IX = ICAMAX( N, WORK, 1 )
00164          XNORM = MAX( ONE, ABS( X( IX, J ) ) )
00165          XSCAL = ( ONE / XNORM ) / REAL( KD+1 )
00166          CALL CSSCAL( N, XSCAL, WORK, 1 )
00167          CALL CTBMV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, WORK, 1 )
00168          CALL CAXPY( N, CMPLX( -SCALE*XSCAL ), B( 1, J ), 1, WORK, 1 )
00169          IX = ICAMAX( N, WORK, 1 )
00170          ERR = TSCAL*ABS( WORK( IX ) )
00171          IX = ICAMAX( N, X( 1, J ), 1 )
00172          XNORM = ABS( X( IX, J ) )
00173          IF( ERR*SMLNUM.LE.XNORM ) THEN
00174             IF( XNORM.GT.ZERO )
00175      $         ERR = ERR / XNORM
00176          ELSE
00177             IF( ERR.GT.ZERO )
00178      $         ERR = ONE / EPS
00179          END IF
00180          IF( ERR*SMLNUM.LE.TNORM ) THEN
00181             IF( TNORM.GT.ZERO )
00182      $         ERR = ERR / TNORM
00183          ELSE
00184             IF( ERR.GT.ZERO )
00185      $         ERR = ONE / EPS
00186          END IF
00187          RESID = MAX( RESID, ERR )
00188    40 CONTINUE
00189 *
00190       RETURN
00191 *
00192 *     End of CTBT03
00193 *
00194       END
 All Files Functions