LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zppt03.f
Go to the documentation of this file.
00001 *> \brief \b ZPPT03
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 ZPPT03( UPLO, N, A, AINV, WORK, LDWORK, RWORK, RCOND,
00012 *                          RESID )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          UPLO
00016 *       INTEGER            LDWORK, N
00017 *       DOUBLE PRECISION   RCOND, RESID
00018 *       ..
00019 *       .. Array Arguments ..
00020 *       DOUBLE PRECISION   RWORK( * )
00021 *       COMPLEX*16         A( * ), AINV( * ), WORK( LDWORK, * )
00022 *       ..
00023 *  
00024 *
00025 *> \par Purpose:
00026 *  =============
00027 *>
00028 *> \verbatim
00029 *>
00030 *> ZPPT03 computes the residual for a Hermitian packed matrix times its
00031 *> inverse:
00032 *>    norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
00033 *> where EPS is the machine epsilon.
00034 *> \endverbatim
00035 *
00036 *  Arguments:
00037 *  ==========
00038 *
00039 *> \param[in] UPLO
00040 *> \verbatim
00041 *>          UPLO is CHARACTER*1
00042 *>          Specifies whether the upper or lower triangular part of the
00043 *>          Hermitian matrix A is stored:
00044 *>          = 'U':  Upper triangular
00045 *>          = 'L':  Lower triangular
00046 *> \endverbatim
00047 *>
00048 *> \param[in] N
00049 *> \verbatim
00050 *>          N is INTEGER
00051 *>          The number of rows and columns of the matrix A.  N >= 0.
00052 *> \endverbatim
00053 *>
00054 *> \param[in] A
00055 *> \verbatim
00056 *>          A is COMPLEX*16 array, dimension (N*(N+1)/2)
00057 *>          The original Hermitian matrix A, stored as a packed
00058 *>          triangular matrix.
00059 *> \endverbatim
00060 *>
00061 *> \param[in] AINV
00062 *> \verbatim
00063 *>          AINV is COMPLEX*16 array, dimension (N*(N+1)/2)
00064 *>          The (Hermitian) inverse of the matrix A, stored as a packed
00065 *>          triangular matrix.
00066 *> \endverbatim
00067 *>
00068 *> \param[out] WORK
00069 *> \verbatim
00070 *>          WORK is COMPLEX*16 array, dimension (LDWORK,N)
00071 *> \endverbatim
00072 *>
00073 *> \param[in] LDWORK
00074 *> \verbatim
00075 *>          LDWORK is INTEGER
00076 *>          The leading dimension of the array WORK.  LDWORK >= max(1,N).
00077 *> \endverbatim
00078 *>
00079 *> \param[out] RWORK
00080 *> \verbatim
00081 *>          RWORK is DOUBLE PRECISION array, dimension (N)
00082 *> \endverbatim
00083 *>
00084 *> \param[out] RCOND
00085 *> \verbatim
00086 *>          RCOND is DOUBLE PRECISION
00087 *>          The reciprocal of the condition number of A, computed as
00088 *>          ( 1/norm(A) ) / norm(AINV).
00089 *> \endverbatim
00090 *>
00091 *> \param[out] RESID
00092 *> \verbatim
00093 *>          RESID is DOUBLE PRECISION
00094 *>          norm(I - A*AINV) / ( 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 ZPPT03( UPLO, N, A, AINV, WORK, LDWORK, RWORK, RCOND,
00111      $                   RESID )
00112 *
00113 *  -- LAPACK test routine (version 3.4.0) --
00114 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00115 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00116 *     November 2011
00117 *
00118 *     .. Scalar Arguments ..
00119       CHARACTER          UPLO
00120       INTEGER            LDWORK, N
00121       DOUBLE PRECISION   RCOND, RESID
00122 *     ..
00123 *     .. Array Arguments ..
00124       DOUBLE PRECISION   RWORK( * )
00125       COMPLEX*16         A( * ), AINV( * ), WORK( LDWORK, * )
00126 *     ..
00127 *
00128 *  =====================================================================
00129 *
00130 *     .. Parameters ..
00131       DOUBLE PRECISION   ZERO, ONE
00132       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00133       COMPLEX*16         CZERO, CONE
00134       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00135      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00136 *     ..
00137 *     .. Local Scalars ..
00138       INTEGER            I, J, JJ
00139       DOUBLE PRECISION   AINVNM, ANORM, EPS
00140 *     ..
00141 *     .. External Functions ..
00142       LOGICAL            LSAME
00143       DOUBLE PRECISION   DLAMCH, ZLANGE, ZLANHP
00144       EXTERNAL           LSAME, DLAMCH, ZLANGE, ZLANHP
00145 *     ..
00146 *     .. Intrinsic Functions ..
00147       INTRINSIC          DBLE, DCONJG
00148 *     ..
00149 *     .. External Subroutines ..
00150       EXTERNAL           ZCOPY, ZHPMV
00151 *     ..
00152 *     .. Executable Statements ..
00153 *
00154 *     Quick exit if N = 0.
00155 *
00156       IF( N.LE.0 ) THEN
00157          RCOND = ONE
00158          RESID = ZERO
00159          RETURN
00160       END IF
00161 *
00162 *     Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
00163 *
00164       EPS = DLAMCH( 'Epsilon' )
00165       ANORM = ZLANHP( '1', UPLO, N, A, RWORK )
00166       AINVNM = ZLANHP( '1', UPLO, N, AINV, RWORK )
00167       IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
00168          RCOND = ZERO
00169          RESID = ONE / EPS
00170          RETURN
00171       END IF
00172       RCOND = ( ONE / ANORM ) / AINVNM
00173 *
00174 *     UPLO = 'U':
00175 *     Copy the leading N-1 x N-1 submatrix of AINV to WORK(1:N,2:N) and
00176 *     expand it to a full matrix, then multiply by A one column at a
00177 *     time, moving the result one column to the left.
00178 *
00179       IF( LSAME( UPLO, 'U' ) ) THEN
00180 *
00181 *        Copy AINV
00182 *
00183          JJ = 1
00184          DO 20 J = 1, N - 1
00185             CALL ZCOPY( J, AINV( JJ ), 1, WORK( 1, J+1 ), 1 )
00186             DO 10 I = 1, J - 1
00187                WORK( J, I+1 ) = DCONJG( AINV( JJ+I-1 ) )
00188    10       CONTINUE
00189             JJ = JJ + J
00190    20    CONTINUE
00191          JJ = ( ( N-1 )*N ) / 2 + 1
00192          DO 30 I = 1, N - 1
00193             WORK( N, I+1 ) = DCONJG( AINV( JJ+I-1 ) )
00194    30    CONTINUE
00195 *
00196 *        Multiply by A
00197 *
00198          DO 40 J = 1, N - 1
00199             CALL ZHPMV( 'Upper', N, -CONE, A, WORK( 1, J+1 ), 1, CZERO,
00200      $                  WORK( 1, J ), 1 )
00201    40    CONTINUE
00202          CALL ZHPMV( 'Upper', N, -CONE, A, AINV( JJ ), 1, CZERO,
00203      $               WORK( 1, N ), 1 )
00204 *
00205 *     UPLO = 'L':
00206 *     Copy the trailing N-1 x N-1 submatrix of AINV to WORK(1:N,1:N-1)
00207 *     and multiply by A, moving each column to the right.
00208 *
00209       ELSE
00210 *
00211 *        Copy AINV
00212 *
00213          DO 50 I = 1, N - 1
00214             WORK( 1, I ) = DCONJG( AINV( I+1 ) )
00215    50    CONTINUE
00216          JJ = N + 1
00217          DO 70 J = 2, N
00218             CALL ZCOPY( N-J+1, AINV( JJ ), 1, WORK( J, J-1 ), 1 )
00219             DO 60 I = 1, N - J
00220                WORK( J, J+I-1 ) = DCONJG( AINV( JJ+I ) )
00221    60       CONTINUE
00222             JJ = JJ + N - J + 1
00223    70    CONTINUE
00224 *
00225 *        Multiply by A
00226 *
00227          DO 80 J = N, 2, -1
00228             CALL ZHPMV( 'Lower', N, -CONE, A, WORK( 1, J-1 ), 1, CZERO,
00229      $                  WORK( 1, J ), 1 )
00230    80    CONTINUE
00231          CALL ZHPMV( 'Lower', N, -CONE, A, AINV( 1 ), 1, CZERO,
00232      $               WORK( 1, 1 ), 1 )
00233 *
00234       END IF
00235 *
00236 *     Add the identity matrix to WORK .
00237 *
00238       DO 90 I = 1, N
00239          WORK( I, I ) = WORK( I, I ) + CONE
00240    90 CONTINUE
00241 *
00242 *     Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
00243 *
00244       RESID = ZLANGE( '1', N, N, WORK, LDWORK, RWORK )
00245 *
00246       RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N )
00247 *
00248       RETURN
00249 *
00250 *     End of ZPPT03
00251 *
00252       END
 All Files Functions