LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
ctpt03.f
Go to the documentation of this file.
00001 *> \brief \b CTPT03
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 CTPT03( UPLO, TRANS, DIAG, N, NRHS, AP, SCALE, CNORM,
00012 *                          TSCAL, X, LDX, B, LDB, WORK, RESID )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          DIAG, TRANS, UPLO
00016 *       INTEGER            LDB, LDX, N, NRHS
00017 *       REAL               RESID, SCALE, TSCAL
00018 *       ..
00019 *       .. Array Arguments ..
00020 *       REAL               CNORM( * )
00021 *       COMPLEX            AP( * ), B( LDB, * ), WORK( * ), X( LDX, * )
00022 *       ..
00023 *  
00024 *
00025 *> \par Purpose:
00026 *  =============
00027 *>
00028 *> \verbatim
00029 *>
00030 *> CTPT03 computes the residual for the solution to a scaled triangular
00031 *> system of equations A*x = s*b,  A**T *x = s*b,  or  A**H *x = s*b,
00032 *> when the triangular matrix A is stored in packed format.  Here A**T
00033 *> denotes the transpose of A, A**H denotes the conjugate transpose of
00034 *> A, s is a scalar, and x and b are N by NRHS matrices.  The test ratio
00035 *> is the maximum over the number of right hand sides of
00036 *>    norm(s*b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
00037 *> where op(A) denotes A, A**T, or A**H, and EPS is the machine epsilon.
00038 *> \endverbatim
00039 *
00040 *  Arguments:
00041 *  ==========
00042 *
00043 *> \param[in] UPLO
00044 *> \verbatim
00045 *>          UPLO is CHARACTER*1
00046 *>          Specifies whether the matrix A is upper or lower triangular.
00047 *>          = 'U':  Upper triangular
00048 *>          = 'L':  Lower triangular
00049 *> \endverbatim
00050 *>
00051 *> \param[in] TRANS
00052 *> \verbatim
00053 *>          TRANS is CHARACTER*1
00054 *>          Specifies the operation applied to A.
00055 *>          = 'N':  A *x = s*b     (No transpose)
00056 *>          = 'T':  A**T *x = s*b  (Transpose)
00057 *>          = 'C':  A**H *x = s*b  (Conjugate transpose)
00058 *> \endverbatim
00059 *>
00060 *> \param[in] DIAG
00061 *> \verbatim
00062 *>          DIAG is CHARACTER*1
00063 *>          Specifies whether or not the matrix A is unit triangular.
00064 *>          = 'N':  Non-unit triangular
00065 *>          = 'U':  Unit triangular
00066 *> \endverbatim
00067 *>
00068 *> \param[in] N
00069 *> \verbatim
00070 *>          N is INTEGER
00071 *>          The order of the matrix A.  N >= 0.
00072 *> \endverbatim
00073 *>
00074 *> \param[in] NRHS
00075 *> \verbatim
00076 *>          NRHS is INTEGER
00077 *>          The number of right hand sides, i.e., the number of columns
00078 *>          of the matrices X and B.  NRHS >= 0.
00079 *> \endverbatim
00080 *>
00081 *> \param[in] AP
00082 *> \verbatim
00083 *>          AP is COMPLEX array, dimension (N*(N+1)/2)
00084 *>          The upper or lower triangular matrix A, packed columnwise in
00085 *>          a linear array.  The j-th column of A is stored in the array
00086 *>          AP as follows:
00087 *>          if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j;
00088 *>          if UPLO = 'L',
00089 *>             AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n.
00090 *> \endverbatim
00091 *>
00092 *> \param[in] SCALE
00093 *> \verbatim
00094 *>          SCALE is REAL
00095 *>          The scaling factor s used in solving the triangular system.
00096 *> \endverbatim
00097 *>
00098 *> \param[in] CNORM
00099 *> \verbatim
00100 *>          CNORM is REAL array, dimension (N)
00101 *>          The 1-norms of the columns of A, not counting the diagonal.
00102 *> \endverbatim
00103 *>
00104 *> \param[in] TSCAL
00105 *> \verbatim
00106 *>          TSCAL is REAL
00107 *>          The scaling factor used in computing the 1-norms in CNORM.
00108 *>          CNORM actually contains the column norms of TSCAL*A.
00109 *> \endverbatim
00110 *>
00111 *> \param[in] X
00112 *> \verbatim
00113 *>          X is COMPLEX array, dimension (LDX,NRHS)
00114 *>          The computed solution vectors for the system of linear
00115 *>          equations.
00116 *> \endverbatim
00117 *>
00118 *> \param[in] LDX
00119 *> \verbatim
00120 *>          LDX is INTEGER
00121 *>          The leading dimension of the array X.  LDX >= max(1,N).
00122 *> \endverbatim
00123 *>
00124 *> \param[in] B
00125 *> \verbatim
00126 *>          B is COMPLEX array, dimension (LDB,NRHS)
00127 *>          The right hand side vectors for the system of linear
00128 *>          equations.
00129 *> \endverbatim
00130 *>
00131 *> \param[in] LDB
00132 *> \verbatim
00133 *>          LDB is INTEGER
00134 *>          The leading dimension of the array B.  LDB >= max(1,N).
00135 *> \endverbatim
00136 *>
00137 *> \param[out] WORK
00138 *> \verbatim
00139 *>          WORK is COMPLEX array, dimension (N)
00140 *> \endverbatim
00141 *>
00142 *> \param[out] RESID
00143 *> \verbatim
00144 *>          RESID is REAL
00145 *>          The maximum over the number of right hand sides of
00146 *>          norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
00147 *> \endverbatim
00148 *
00149 *  Authors:
00150 *  ========
00151 *
00152 *> \author Univ. of Tennessee 
00153 *> \author Univ. of California Berkeley 
00154 *> \author Univ. of Colorado Denver 
00155 *> \author NAG Ltd. 
00156 *
00157 *> \date November 2011
00158 *
00159 *> \ingroup complex_lin
00160 *
00161 *  =====================================================================
00162       SUBROUTINE CTPT03( UPLO, TRANS, DIAG, N, NRHS, AP, SCALE, CNORM,
00163      $                   TSCAL, X, LDX, B, LDB, WORK, RESID )
00164 *
00165 *  -- LAPACK test routine (version 3.4.0) --
00166 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00167 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00168 *     November 2011
00169 *
00170 *     .. Scalar Arguments ..
00171       CHARACTER          DIAG, TRANS, UPLO
00172       INTEGER            LDB, LDX, N, NRHS
00173       REAL               RESID, SCALE, TSCAL
00174 *     ..
00175 *     .. Array Arguments ..
00176       REAL               CNORM( * )
00177       COMPLEX            AP( * ), B( LDB, * ), WORK( * ), X( LDX, * )
00178 *     ..
00179 *
00180 *  =====================================================================
00181 *
00182 *     .. Parameters ..
00183       REAL               ONE, ZERO
00184       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00185 *     ..
00186 *     .. Local Scalars ..
00187       INTEGER            IX, J, JJ
00188       REAL               EPS, ERR, SMLNUM, TNORM, XNORM, XSCAL
00189 *     ..
00190 *     .. External Functions ..
00191       LOGICAL            LSAME
00192       INTEGER            ICAMAX
00193       REAL               SLAMCH
00194       EXTERNAL           LSAME, ICAMAX, SLAMCH
00195 *     ..
00196 *     .. External Subroutines ..
00197       EXTERNAL           CAXPY, CCOPY, CSSCAL, CTPMV
00198 *     ..
00199 *     .. Intrinsic Functions ..
00200       INTRINSIC          ABS, CMPLX, MAX, REAL
00201 *     ..
00202 *     .. Executable Statements ..
00203 *
00204 *     Quick exit if N = 0.
00205 *
00206       IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
00207          RESID = ZERO
00208          RETURN
00209       END IF
00210       EPS = SLAMCH( 'Epsilon' )
00211       SMLNUM = SLAMCH( 'Safe minimum' )
00212 *
00213 *     Compute the norm of the triangular matrix A using the column
00214 *     norms already computed by CLATPS.
00215 *
00216       TNORM = 0.
00217       IF( LSAME( DIAG, 'N' ) ) THEN
00218          IF( LSAME( UPLO, 'U' ) ) THEN
00219             JJ = 1
00220             DO 10 J = 1, N
00221                TNORM = MAX( TNORM, TSCAL*ABS( AP( JJ ) )+CNORM( J ) )
00222                JJ = JJ + J + 1
00223    10       CONTINUE
00224          ELSE
00225             JJ = 1
00226             DO 20 J = 1, N
00227                TNORM = MAX( TNORM, TSCAL*ABS( AP( JJ ) )+CNORM( J ) )
00228                JJ = JJ + N - J + 1
00229    20       CONTINUE
00230          END IF
00231       ELSE
00232          DO 30 J = 1, N
00233             TNORM = MAX( TNORM, TSCAL+CNORM( J ) )
00234    30    CONTINUE
00235       END IF
00236 *
00237 *     Compute the maximum over the number of right hand sides of
00238 *        norm(op(A)*x - s*b) / ( norm(A) * norm(x) * EPS ).
00239 *
00240       RESID = ZERO
00241       DO 40 J = 1, NRHS
00242          CALL CCOPY( N, X( 1, J ), 1, WORK, 1 )
00243          IX = ICAMAX( N, WORK, 1 )
00244          XNORM = MAX( ONE, ABS( X( IX, J ) ) )
00245          XSCAL = ( ONE / XNORM ) / REAL( N )
00246          CALL CSSCAL( N, XSCAL, WORK, 1 )
00247          CALL CTPMV( UPLO, TRANS, DIAG, N, AP, WORK, 1 )
00248          CALL CAXPY( N, CMPLX( -SCALE*XSCAL ), B( 1, J ), 1, WORK, 1 )
00249          IX = ICAMAX( N, WORK, 1 )
00250          ERR = TSCAL*ABS( WORK( IX ) )
00251          IX = ICAMAX( N, X( 1, J ), 1 )
00252          XNORM = ABS( X( IX, J ) )
00253          IF( ERR*SMLNUM.LE.XNORM ) THEN
00254             IF( XNORM.GT.ZERO )
00255      $         ERR = ERR / XNORM
00256          ELSE
00257             IF( ERR.GT.ZERO )
00258      $         ERR = ONE / EPS
00259          END IF
00260          IF( ERR*SMLNUM.LE.TNORM ) THEN
00261             IF( TNORM.GT.ZERO )
00262      $         ERR = ERR / TNORM
00263          ELSE
00264             IF( ERR.GT.ZERO )
00265      $         ERR = ONE / EPS
00266          END IF
00267          RESID = MAX( RESID, ERR )
00268    40 CONTINUE
00269 *
00270       RETURN
00271 *
00272 *     End of CTPT03
00273 *
00274       END
 All Files Functions