TESTING/LIN/zsyt03.f

   1 *> \brief \b ZSYT03
   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 ZSYT03( 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 *> ZSYT03 computes the residual for a complex symmetric matrix times
  32 *> its 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 *>          complex symmetric 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 complex symmetric 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 symmetric
  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 *>          RCOND = 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 ZSYT03( 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 *
 148 *     .. Parameters ..
 149       DOUBLE PRECISION   ZERO, ONE
 150       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 151       COMPLEX*16         CZERO, CONE
 152       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
 153      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
 154 *     ..
 155 *     .. Local Scalars ..
 156       INTEGER            I, J
 157       DOUBLE PRECISION   AINVNM, ANORM, EPS
 158 *     ..
 159 *     .. External Functions ..
 160       LOGICAL            LSAME
 161       DOUBLE PRECISION   DLAMCH, ZLANGE, ZLANSY
 162       EXTERNAL           LSAME, DLAMCH, ZLANGE, ZLANSY
 163 *     ..
 164 *     .. External Subroutines ..
 165       EXTERNAL           ZSYMM
 166 *     ..
 167 *     .. Intrinsic Functions ..
 168       INTRINSIC          DBLE
 169 *     ..
 170 *     .. Executable Statements ..
 171 *
 172 *     Quick exit if N = 0
 173 *
 174       IF( N.LE.0 ) THEN
 175          RCOND = ONE
 176          RESID = ZERO
 177          RETURN
 178       END IF
 179 *
 180 *     Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
 181 *
 182       EPS = DLAMCH( 'Epsilon' )
 183       ANORM = ZLANSY( '1', UPLO, N, A, LDA, RWORK )
 184       AINVNM = ZLANSY( '1', UPLO, N, AINV, LDAINV, RWORK )
 185       IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
 186          RCOND = ZERO
 187          RESID = ONE / EPS
 188          RETURN
 189       END IF
 190       RCOND = ( ONE / ANORM ) / AINVNM
 191 *
 192 *     Expand AINV into a full matrix and call ZSYMM to multiply
 193 *     AINV on the left by A (store the result in WORK).
 194 *
 195       IF( LSAME( UPLO, 'U' ) ) THEN
 196          DO 20 J = 1, N
 197             DO 10 I = 1, J - 1
 198                AINV( J, I ) = AINV( I, J )
 199    10       CONTINUE
 200    20    CONTINUE
 201       ELSE
 202          DO 40 J = 1, N
 203             DO 30 I = J + 1, N
 204                AINV( J, I ) = AINV( I, J )
 205    30       CONTINUE
 206    40    CONTINUE
 207       END IF
 208       CALL ZSYMM( 'Left', UPLO, N, N, -CONE, A, LDA, AINV, LDAINV,
 209      $            CZERO, WORK, LDWORK )
 210 *
 211 *     Add the identity matrix to WORK .
 212 *
 213       DO 50 I = 1, N
 214          WORK( I, I ) = WORK( I, I ) + CONE
 215    50 CONTINUE
 216 *
 217 *     Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
 218 *
 219       RESID = ZLANGE( '1', N, N, WORK, LDWORK, RWORK )
 220 *
 221       RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N )
 222 *
 223       RETURN
 224 *
 225 *     End of ZSYT03
 226 *
 227       END