LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
ztpt01.f
Go to the documentation of this file.
00001 *> \brief \b ZTPT01
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 ZTPT01( UPLO, DIAG, N, AP, AINVP, RCOND, RWORK, RESID )
00012 * 
00013 *       .. Scalar Arguments ..
00014 *       CHARACTER          DIAG, UPLO
00015 *       INTEGER            N
00016 *       DOUBLE PRECISION   RCOND, RESID
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       DOUBLE PRECISION   RWORK( * )
00020 *       COMPLEX*16         AINVP( * ), AP( * )
00021 *       ..
00022 *  
00023 *
00024 *> \par Purpose:
00025 *  =============
00026 *>
00027 *> \verbatim
00028 *>
00029 *> ZTPT01 computes the residual for a triangular matrix A times its
00030 *> inverse when A is stored in packed format:
00031 *>    RESID = norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS ),
00032 *> where EPS is the machine epsilon.
00033 *> \endverbatim
00034 *
00035 *  Arguments:
00036 *  ==========
00037 *
00038 *> \param[in] UPLO
00039 *> \verbatim
00040 *>          UPLO is CHARACTER*1
00041 *>          Specifies whether the matrix A is upper or lower triangular.
00042 *>          = 'U':  Upper triangular
00043 *>          = 'L':  Lower triangular
00044 *> \endverbatim
00045 *>
00046 *> \param[in] DIAG
00047 *> \verbatim
00048 *>          DIAG is CHARACTER*1
00049 *>          Specifies whether or not the matrix A is unit triangular.
00050 *>          = 'N':  Non-unit triangular
00051 *>          = 'U':  Unit triangular
00052 *> \endverbatim
00053 *>
00054 *> \param[in] N
00055 *> \verbatim
00056 *>          N is INTEGER
00057 *>          The order of the matrix A.  N >= 0.
00058 *> \endverbatim
00059 *>
00060 *> \param[in] AP
00061 *> \verbatim
00062 *>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
00063 *>          The original upper or lower triangular matrix A, packed
00064 *>          columnwise in a linear array.  The j-th column of A is stored
00065 *>          in the array AP as follows:
00066 *>          if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j;
00067 *>          if UPLO = 'L',
00068 *>             AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n.
00069 *> \endverbatim
00070 *>
00071 *> \param[in] AINVP
00072 *> \verbatim
00073 *>          AINVP is COMPLEX*16 array, dimension (N*(N+1)/2)
00074 *>          On entry, the (triangular) inverse of the matrix A, packed
00075 *>          columnwise in a linear array as in AP.
00076 *>          On exit, the contents of AINVP are destroyed.
00077 *> \endverbatim
00078 *>
00079 *> \param[out] RCOND
00080 *> \verbatim
00081 *>          RCOND is DOUBLE PRECISION
00082 *>          The reciprocal condition number of A, computed as
00083 *>          1/(norm(A) * norm(AINV)).
00084 *> \endverbatim
00085 *>
00086 *> \param[out] RWORK
00087 *> \verbatim
00088 *>          RWORK is DOUBLE PRECISION array, dimension (N)
00089 *> \endverbatim
00090 *>
00091 *> \param[out] RESID
00092 *> \verbatim
00093 *>          RESID is DOUBLE PRECISION
00094 *>          norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS )
00095 *> \endverbatim
00096 *
00097 *  Authors:
00098 *  ========
00099 *
00100 *> \author Univ. of Tennessee 
00101 *> \author Univ. of California Berkeley 
00102 *> \author Univ. of Colorado Denver 
00103 *> \author NAG Ltd. 
00104 *
00105 *> \date November 2011
00106 *
00107 *> \ingroup complex16_lin
00108 *
00109 *  =====================================================================
00110       SUBROUTINE ZTPT01( UPLO, DIAG, N, AP, AINVP, RCOND, RWORK, RESID )
00111 *
00112 *  -- LAPACK test routine (version 3.4.0) --
00113 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00114 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00115 *     November 2011
00116 *
00117 *     .. Scalar Arguments ..
00118       CHARACTER          DIAG, UPLO
00119       INTEGER            N
00120       DOUBLE PRECISION   RCOND, RESID
00121 *     ..
00122 *     .. Array Arguments ..
00123       DOUBLE PRECISION   RWORK( * )
00124       COMPLEX*16         AINVP( * ), AP( * )
00125 *     ..
00126 *
00127 *  =====================================================================
00128 *
00129 *     .. Parameters ..
00130       DOUBLE PRECISION   ZERO, ONE
00131       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00132 *     ..
00133 *     .. Local Scalars ..
00134       LOGICAL            UNITD
00135       INTEGER            J, JC
00136       DOUBLE PRECISION   AINVNM, ANORM, EPS
00137 *     ..
00138 *     .. External Functions ..
00139       LOGICAL            LSAME
00140       DOUBLE PRECISION   DLAMCH, ZLANTP
00141       EXTERNAL           LSAME, DLAMCH, ZLANTP
00142 *     ..
00143 *     .. External Subroutines ..
00144       EXTERNAL           ZTPMV
00145 *     ..
00146 *     .. Intrinsic Functions ..
00147       INTRINSIC          DBLE
00148 *     ..
00149 *     .. Executable Statements ..
00150 *
00151 *     Quick exit if N = 0.
00152 *
00153       IF( N.LE.0 ) THEN
00154          RCOND = ONE
00155          RESID = ZERO
00156          RETURN
00157       END IF
00158 *
00159 *     Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
00160 *
00161       EPS = DLAMCH( 'Epsilon' )
00162       ANORM = ZLANTP( '1', UPLO, DIAG, N, AP, RWORK )
00163       AINVNM = ZLANTP( '1', UPLO, DIAG, N, AINVP, RWORK )
00164       IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
00165          RCOND = ZERO
00166          RESID = ONE / EPS
00167          RETURN
00168       END IF
00169       RCOND = ( ONE / ANORM ) / AINVNM
00170 *
00171 *     Compute A * AINV, overwriting AINV.
00172 *
00173       UNITD = LSAME( DIAG, 'U' )
00174       IF( LSAME( UPLO, 'U' ) ) THEN
00175          JC = 1
00176          DO 10 J = 1, N
00177             IF( UNITD )
00178      $         AINVP( JC+J-1 ) = ONE
00179 *
00180 *           Form the j-th column of A*AINV.
00181 *
00182             CALL ZTPMV( 'Upper', 'No transpose', DIAG, J, AP,
00183      $                  AINVP( JC ), 1 )
00184 *
00185 *           Subtract 1 from the diagonal to form A*AINV - I.
00186 *
00187             AINVP( JC+J-1 ) = AINVP( JC+J-1 ) - ONE
00188             JC = JC + J
00189    10    CONTINUE
00190       ELSE
00191          JC = 1
00192          DO 20 J = 1, N
00193             IF( UNITD )
00194      $         AINVP( JC ) = ONE
00195 *
00196 *           Form the j-th column of A*AINV.
00197 *
00198             CALL ZTPMV( 'Lower', 'No transpose', DIAG, N-J+1, AP( JC ),
00199      $                  AINVP( JC ), 1 )
00200 *
00201 *           Subtract 1 from the diagonal to form A*AINV - I.
00202 *
00203             AINVP( JC ) = AINVP( JC ) - ONE
00204             JC = JC + N - J + 1
00205    20    CONTINUE
00206       END IF
00207 *
00208 *     Compute norm(A*AINV - I) / (N * norm(A) * norm(AINV) * EPS)
00209 *
00210       RESID = ZLANTP( '1', UPLO, 'Non-unit', N, AINVP, RWORK )
00211 *
00212       RESID = ( ( RESID*RCOND ) / DBLE( N ) ) / EPS
00213 *
00214       RETURN
00215 *
00216 *     End of ZTPT01
00217 *
00218       END
 All Files Functions