TESTING/LIN/zpot03.f

   1 *> \brief \b ZPOT03
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE ZPOT03( UPLO, N, A, LDA, AINV, LDAINV, WORK, LDWORK,
  12 *                          RWORK, RCOND, RESID )
  13 *
  14 *       .. Scalar Arguments ..
  15 *       CHARACTER          UPLO
  16 *       INTEGER            LDA, LDAINV, LDWORK, N
  17 *       DOUBLE PRECISION   RCOND, RESID
  18 *       ..
  19 *       .. Array Arguments ..
  20 *       DOUBLE PRECISION   RWORK( * )
  21 *       COMPLEX*16         A( LDA, * ), AINV( LDAINV, * ),
  22 *      $                   WORK( LDWORK, * )
  23 *       ..
  24 *
  25 *
  26 *> \par Purpose:
  27 *  =============
  28 *>
  29 *> \verbatim
  30 *>
  31 *> ZPOT03 computes the residual for a Hermitian matrix times its
  32 *> inverse:
  33 *>    norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
  34 *> where EPS is the machine epsilon.
  35 *> \endverbatim
  36 *
  37 *  Arguments:
  38 *  ==========
  39 *
  40 *> \param[in] UPLO
  41 *> \verbatim
  42 *>          UPLO is CHARACTER*1
  43 *>          Specifies whether the upper or lower triangular part of the
  44 *>          Hermitian matrix A is stored:
  45 *>          = 'U':  Upper triangular
  46 *>          = 'L':  Lower triangular
  47 *> \endverbatim
  48 *>
  49 *> \param[in] N
  50 *> \verbatim
  51 *>          N is INTEGER
  52 *>          The number of rows and columns of the matrix A.  N >= 0.
  53 *> \endverbatim
  54 *>
  55 *> \param[in] A
  56 *> \verbatim
  57 *>          A is COMPLEX*16 array, dimension (LDA,N)
  58 *>          The original Hermitian matrix A.
  59 *> \endverbatim
  60 *>
  61 *> \param[in] LDA
  62 *> \verbatim
  63 *>          LDA is INTEGER
  64 *>          The leading dimension of the array A.  LDA >= max(1,N)
  65 *> \endverbatim
  66 *>
  67 *> \param[in,out] AINV
  68 *> \verbatim
  69 *>          AINV is COMPLEX*16 array, dimension (LDAINV,N)
  70 *>          On entry, the inverse of the matrix A, stored as a Hermitian
  71 *>          matrix in the same format as A.
  72 *>          In this version, AINV is expanded into a full matrix and
  73 *>          multiplied by A, so the opposing triangle of AINV will be
  74 *>          changed; i.e., if the upper triangular part of AINV is
  75 *>          stored, the lower triangular part will be used as work space.
  76 *> \endverbatim
  77 *>
  78 *> \param[in] LDAINV
  79 *> \verbatim
  80 *>          LDAINV is INTEGER
  81 *>          The leading dimension of the array AINV.  LDAINV >= max(1,N).
  82 *> \endverbatim
  83 *>
  84 *> \param[out] WORK
  85 *> \verbatim
  86 *>          WORK is COMPLEX*16 array, dimension (LDWORK,N)
  87 *> \endverbatim
  88 *>
  89 *> \param[in] LDWORK
  90 *> \verbatim
  91 *>          LDWORK is INTEGER
  92 *>          The leading dimension of the array WORK.  LDWORK >= max(1,N).
  93 *> \endverbatim
  94 *>
  95 *> \param[out] RWORK
  96 *> \verbatim
  97 *>          RWORK is DOUBLE PRECISION array, dimension (N)
  98 *> \endverbatim
  99 *>
 100 *> \param[out] RCOND
 101 *> \verbatim
 102 *>          RCOND is DOUBLE PRECISION
 103 *>          The reciprocal of the condition number of A, computed as
 104 *>          ( 1/norm(A) ) / norm(AINV).
 105 *> \endverbatim
 106 *>
 107 *> \param[out] RESID
 108 *> \verbatim
 109 *>          RESID is DOUBLE PRECISION
 110 *>          norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS )
 111 *> \endverbatim
 112 *
 113 *  Authors:
 114 *  ========
 115 *
 116 *> \author Univ. of Tennessee
 117 *> \author Univ. of California Berkeley
 118 *> \author Univ. of Colorado Denver
 119 *> \author NAG Ltd.
 120 *
 121 *> \date November 2011
 122 *
 123 *> \ingroup complex16_lin
 124 *
 125 *  =====================================================================
 126       SUBROUTINE ZPOT03( UPLO, N, A, LDA, AINV, LDAINV, WORK, LDWORK,
 127      $                   RWORK, RCOND, RESID )
 128 *
 129 *  -- LAPACK test routine (version 3.4.0) --
 130 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 131 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 132 *     November 2011
 133 *
 134 *     .. Scalar Arguments ..
 135       CHARACTER          UPLO
 136       INTEGER            LDA, LDAINV, LDWORK, N
 137       DOUBLE PRECISION   RCOND, RESID
 138 *     ..
 139 *     .. Array Arguments ..
 140       DOUBLE PRECISION   RWORK( * )
 141       COMPLEX*16         A( LDA, * ), AINV( LDAINV, * ),
 142      $                   WORK( LDWORK, * )
 143 *     ..
 144 *
 145 *  =====================================================================
 146 *
 147 *     .. Parameters ..
 148       DOUBLE PRECISION   ZERO, ONE
 149       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 150       COMPLEX*16         CZERO, CONE
 151       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
 152      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
 153 *     ..
 154 *     .. Local Scalars ..
 155       INTEGER            I, J
 156       DOUBLE PRECISION   AINVNM, ANORM, EPS
 157 *     ..
 158 *     .. External Functions ..
 159       LOGICAL            LSAME
 160       DOUBLE PRECISION   DLAMCH, ZLANGE, ZLANHE
 161       EXTERNAL           LSAME, DLAMCH, ZLANGE, ZLANHE
 162 *     ..
 163 *     .. External Subroutines ..
 164       EXTERNAL           ZHEMM
 165 *     ..
 166 *     .. Intrinsic Functions ..
 167       INTRINSIC          DBLE, DCONJG
 168 *     ..
 169 *     .. Executable Statements ..
 170 *
 171 *     Quick exit if N = 0.
 172 *
 173       IF( N.LE.0 ) THEN
 174          RCOND = ONE
 175          RESID = ZERO
 176          RETURN
 177       END IF
 178 *
 179 *     Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
 180 *
 181       EPS = DLAMCH( 'Epsilon' )
 182       ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK )
 183       AINVNM = ZLANHE( '1', UPLO, N, AINV, LDAINV, RWORK )
 184       IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
 185          RCOND = ZERO
 186          RESID = ONE / EPS
 187          RETURN
 188       END IF
 189       RCOND = ( ONE / ANORM ) / AINVNM
 190 *
 191 *     Expand AINV into a full matrix and call ZHEMM to multiply
 192 *     AINV on the left by A.
 193 *
 194       IF( LSAME( UPLO, 'U' ) ) THEN
 195          DO 20 J = 1, N
 196             DO 10 I = 1, J - 1
 197                AINV( J, I ) = DCONJG( AINV( I, J ) )
 198    10       CONTINUE
 199    20    CONTINUE
 200       ELSE
 201          DO 40 J = 1, N
 202             DO 30 I = J + 1, N
 203                AINV( J, I ) = DCONJG( AINV( I, J ) )
 204    30       CONTINUE
 205    40    CONTINUE
 206       END IF
 207       CALL ZHEMM( 'Left', UPLO, N, N, -CONE, A, LDA, AINV, LDAINV,
 208      $            CZERO, WORK, LDWORK )
 209 *
 210 *     Add the identity matrix to WORK .
 211 *
 212       DO 50 I = 1, N
 213          WORK( I, I ) = WORK( I, I ) + CONE
 214    50 CONTINUE
 215 *
 216 *     Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
 217 *
 218       RESID = ZLANGE( '1', N, N, WORK, LDWORK, RWORK )
 219 *
 220       RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N )
 221 *
 222       RETURN
 223 *
 224 *     End of ZPOT03
 225 *
 226       END