LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dtpt02.f
Go to the documentation of this file.
00001 *> \brief \b DTPT02
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE DTPT02( UPLO, TRANS, DIAG, N, NRHS, AP, X, LDX, B, LDB,
00012 *                          WORK, RESID )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          DIAG, TRANS, UPLO
00016 *       INTEGER            LDB, LDX, N, NRHS
00017 *       DOUBLE PRECISION   RESID
00018 *       ..
00019 *       .. Array Arguments ..
00020 *       DOUBLE PRECISION   AP( * ), B( LDB, * ), WORK( * ), X( LDX, * )
00021 *       ..
00022 *  
00023 *
00024 *> \par Purpose:
00025 *  =============
00026 *>
00027 *> \verbatim
00028 *>
00029 *> DTPT02 computes the residual for the computed solution to a
00030 *> triangular system of linear equations  A*x = b  or  A'*x = b  when
00031 *> the triangular matrix A is stored in packed format.  Here A' is the
00032 *> transpose of A and x and b are N by NRHS matrices.  The test ratio is
00033 *> the maximum over the number of right hand sides of
00034 *>    norm(b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
00035 *> where op(A) denotes A or A' and EPS is the machine epsilon.
00036 *> \endverbatim
00037 *
00038 *  Arguments:
00039 *  ==========
00040 *
00041 *> \param[in] UPLO
00042 *> \verbatim
00043 *>          UPLO is CHARACTER*1
00044 *>          Specifies whether the matrix A is upper or lower triangular.
00045 *>          = 'U':  Upper triangular
00046 *>          = 'L':  Lower triangular
00047 *> \endverbatim
00048 *>
00049 *> \param[in] TRANS
00050 *> \verbatim
00051 *>          TRANS is CHARACTER*1
00052 *>          Specifies the operation applied to A.
00053 *>          = 'N':  A *x = b  (No transpose)
00054 *>          = 'T':  A'*x = b  (Transpose)
00055 *>          = 'C':  A'*x = b  (Conjugate transpose = Transpose)
00056 *> \endverbatim
00057 *>
00058 *> \param[in] DIAG
00059 *> \verbatim
00060 *>          DIAG is CHARACTER*1
00061 *>          Specifies whether or not the matrix A is unit triangular.
00062 *>          = 'N':  Non-unit triangular
00063 *>          = 'U':  Unit triangular
00064 *> \endverbatim
00065 *>
00066 *> \param[in] N
00067 *> \verbatim
00068 *>          N is INTEGER
00069 *>          The order of the matrix A.  N >= 0.
00070 *> \endverbatim
00071 *>
00072 *> \param[in] NRHS
00073 *> \verbatim
00074 *>          NRHS is INTEGER
00075 *>          The number of right hand sides, i.e., the number of columns
00076 *>          of the matrices X and B.  NRHS >= 0.
00077 *> \endverbatim
00078 *>
00079 *> \param[in] AP
00080 *> \verbatim
00081 *>          AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
00082 *>          The upper or lower triangular matrix A, packed columnwise in
00083 *>          a linear array.  The j-th column of A is stored in the array
00084 *>          AP as follows:
00085 *>          if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j;
00086 *>          if UPLO = 'L',
00087 *>             AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n.
00088 *> \endverbatim
00089 *>
00090 *> \param[in] X
00091 *> \verbatim
00092 *>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
00093 *>          The computed solution vectors for the system of linear
00094 *>          equations.
00095 *> \endverbatim
00096 *>
00097 *> \param[in] LDX
00098 *> \verbatim
00099 *>          LDX is INTEGER
00100 *>          The leading dimension of the array X.  LDX >= max(1,N).
00101 *> \endverbatim
00102 *>
00103 *> \param[in] B
00104 *> \verbatim
00105 *>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
00106 *>          The right hand side vectors for the system of linear
00107 *>          equations.
00108 *> \endverbatim
00109 *>
00110 *> \param[in] LDB
00111 *> \verbatim
00112 *>          LDB is INTEGER
00113 *>          The leading dimension of the array B.  LDB >= max(1,N).
00114 *> \endverbatim
00115 *>
00116 *> \param[out] WORK
00117 *> \verbatim
00118 *>          WORK is DOUBLE PRECISION array, dimension (N)
00119 *> \endverbatim
00120 *>
00121 *> \param[out] RESID
00122 *> \verbatim
00123 *>          RESID is DOUBLE PRECISION
00124 *>          The maximum over the number of right hand sides of
00125 *>          norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
00126 *> \endverbatim
00127 *
00128 *  Authors:
00129 *  ========
00130 *
00131 *> \author Univ. of Tennessee 
00132 *> \author Univ. of California Berkeley 
00133 *> \author Univ. of Colorado Denver 
00134 *> \author NAG Ltd. 
00135 *
00136 *> \date November 2011
00137 *
00138 *> \ingroup double_lin
00139 *
00140 *  =====================================================================
00141       SUBROUTINE DTPT02( UPLO, TRANS, DIAG, N, NRHS, AP, X, LDX, B, LDB,
00142      $                   WORK, RESID )
00143 *
00144 *  -- LAPACK test routine (version 3.4.0) --
00145 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00146 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00147 *     November 2011
00148 *
00149 *     .. Scalar Arguments ..
00150       CHARACTER          DIAG, TRANS, UPLO
00151       INTEGER            LDB, LDX, N, NRHS
00152       DOUBLE PRECISION   RESID
00153 *     ..
00154 *     .. Array Arguments ..
00155       DOUBLE PRECISION   AP( * ), B( LDB, * ), WORK( * ), X( LDX, * )
00156 *     ..
00157 *
00158 *  =====================================================================
00159 *
00160 *     .. Parameters ..
00161       DOUBLE PRECISION   ZERO, ONE
00162       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00163 *     ..
00164 *     .. Local Scalars ..
00165       INTEGER            J
00166       DOUBLE PRECISION   ANORM, BNORM, EPS, XNORM
00167 *     ..
00168 *     .. External Functions ..
00169       LOGICAL            LSAME
00170       DOUBLE PRECISION   DASUM, DLAMCH, DLANTP
00171       EXTERNAL           LSAME, DASUM, DLAMCH, DLANTP
00172 *     ..
00173 *     .. External Subroutines ..
00174       EXTERNAL           DAXPY, DCOPY, DTPMV
00175 *     ..
00176 *     .. Intrinsic Functions ..
00177       INTRINSIC          MAX
00178 *     ..
00179 *     .. Executable Statements ..
00180 *
00181 *     Quick exit if N = 0 or NRHS = 0
00182 *
00183       IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
00184          RESID = ZERO
00185          RETURN
00186       END IF
00187 *
00188 *     Compute the 1-norm of A or A'.
00189 *
00190       IF( LSAME( TRANS, 'N' ) ) THEN
00191          ANORM = DLANTP( '1', UPLO, DIAG, N, AP, WORK )
00192       ELSE
00193          ANORM = DLANTP( 'I', UPLO, DIAG, N, AP, WORK )
00194       END IF
00195 *
00196 *     Exit with RESID = 1/EPS if ANORM = 0.
00197 *
00198       EPS = DLAMCH( 'Epsilon' )
00199       IF( ANORM.LE.ZERO ) THEN
00200          RESID = ONE / EPS
00201          RETURN
00202       END IF
00203 *
00204 *     Compute the maximum over the number of right hand sides of
00205 *        norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
00206 *
00207       RESID = ZERO
00208       DO 10 J = 1, NRHS
00209          CALL DCOPY( N, X( 1, J ), 1, WORK, 1 )
00210          CALL DTPMV( UPLO, TRANS, DIAG, N, AP, WORK, 1 )
00211          CALL DAXPY( N, -ONE, B( 1, J ), 1, WORK, 1 )
00212          BNORM = DASUM( N, WORK, 1 )
00213          XNORM = DASUM( N, X( 1, J ), 1 )
00214          IF( XNORM.LE.ZERO ) THEN
00215             RESID = ONE / EPS
00216          ELSE
00217             RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
00218          END IF
00219    10 CONTINUE
00220 *
00221       RETURN
00222 *
00223 *     End of DTPT02
00224 *
00225       END
 All Files Functions