TESTING/LIN/zhpt01.f

   1 *> \brief \b ZHPT01
   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 ZHPT01( UPLO, N, A, AFAC, IPIV, C, LDC, RWORK, RESID )
  12 *
  13 *       .. Scalar Arguments ..
  14 *       CHARACTER          UPLO
  15 *       INTEGER            LDC, N
  16 *       DOUBLE PRECISION   RESID
  17 *       ..
  18 *       .. Array Arguments ..
  19 *       INTEGER            IPIV( * )
  20 *       DOUBLE PRECISION   RWORK( * )
  21 *       COMPLEX*16         A( * ), AFAC( * ), C( LDC, * )
  22 *       ..
  23 *
  24 *
  25 *> \par Purpose:
  26 *  =============
  27 *>
  28 *> \verbatim
  29 *>
  30 *> ZHPT01 reconstructs a Hermitian indefinite packed matrix A from its
  31 *> block L*D*L' or U*D*U' factorization and computes the residual
  32 *>    norm( C - A ) / ( N * norm(A) * EPS ),
  33 *> where C is the reconstructed matrix, EPS is the machine epsilon,
  34 *> L' is the conjugate transpose of L, and U' is the conjugate transpose
  35 *> of U.
  36 *> \endverbatim
  37 *
  38 *  Arguments:
  39 *  ==========
  40 *
  41 *> \param[in] UPLO
  42 *> \verbatim
  43 *>          UPLO is CHARACTER*1
  44 *>          Specifies whether the upper or lower triangular part of the
  45 *>          Hermitian matrix A is stored:
  46 *>          = 'U':  Upper triangular
  47 *>          = 'L':  Lower triangular
  48 *> \endverbatim
  49 *>
  50 *> \param[in] N
  51 *> \verbatim
  52 *>          N is INTEGER
  53 *>          The number of rows and columns of the matrix A.  N >= 0.
  54 *> \endverbatim
  55 *>
  56 *> \param[in] A
  57 *> \verbatim
  58 *>          A is COMPLEX*16 array, dimension (N*(N+1)/2)
  59 *>          The original Hermitian matrix A, stored as a packed
  60 *>          triangular matrix.
  61 *> \endverbatim
  62 *>
  63 *> \param[in] AFAC
  64 *> \verbatim
  65 *>          AFAC is COMPLEX*16 array, dimension (N*(N+1)/2)
  66 *>          The factored form of the matrix A, stored as a packed
  67 *>          triangular matrix.  AFAC contains the block diagonal matrix D
  68 *>          and the multipliers used to obtain the factor L or U from the
  69 *>          block L*D*L' or U*D*U' factorization as computed by ZHPTRF.
  70 *> \endverbatim
  71 *>
  72 *> \param[in] IPIV
  73 *> \verbatim
  74 *>          IPIV is INTEGER array, dimension (N)
  75 *>          The pivot indices from ZHPTRF.
  76 *> \endverbatim
  77 *>
  78 *> \param[out] C
  79 *> \verbatim
  80 *>          C is COMPLEX*16 array, dimension (LDC,N)
  81 *> \endverbatim
  82 *>
  83 *> \param[in] LDC
  84 *> \verbatim
  85 *>          LDC is INTEGER
  86 *>          The leading dimension of the array C.  LDC >= max(1,N).
  87 *> \endverbatim
  88 *>
  89 *> \param[out] RWORK
  90 *> \verbatim
  91 *>          RWORK is DOUBLE PRECISION array, dimension (N)
  92 *> \endverbatim
  93 *>
  94 *> \param[out] RESID
  95 *> \verbatim
  96 *>          RESID is DOUBLE PRECISION
  97 *>          If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS )
  98 *>          If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS )
  99 *> \endverbatim
 100 *
 101 *  Authors:
 102 *  ========
 103 *
 104 *> \author Univ. of Tennessee
 105 *> \author Univ. of California Berkeley
 106 *> \author Univ. of Colorado Denver
 107 *> \author NAG Ltd.
 108 *
 109 *> \date November 2011
 110 *
 111 *> \ingroup complex16_lin
 112 *
 113 *  =====================================================================
 114       SUBROUTINE ZHPT01( UPLO, N, A, AFAC, IPIV, C, LDC, RWORK, RESID )
 115 *
 116 *  -- LAPACK test routine (version 3.4.0) --
 117 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 118 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 119 *     November 2011
 120 *
 121 *     .. Scalar Arguments ..
 122       CHARACTER          UPLO
 123       INTEGER            LDC, N
 124       DOUBLE PRECISION   RESID
 125 *     ..
 126 *     .. Array Arguments ..
 127       INTEGER            IPIV( * )
 128       DOUBLE PRECISION   RWORK( * )
 129       COMPLEX*16         A( * ), AFAC( * ), C( LDC, * )
 130 *     ..
 131 *
 132 *  =====================================================================
 133 *
 134 *     .. Parameters ..
 135       DOUBLE PRECISION   ZERO, ONE
 136       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 137       COMPLEX*16         CZERO, CONE
 138       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
 139      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
 140 *     ..
 141 *     .. Local Scalars ..
 142       INTEGER            I, INFO, J, JC
 143       DOUBLE PRECISION   ANORM, EPS
 144 *     ..
 145 *     .. External Functions ..
 146       LOGICAL            LSAME
 147       DOUBLE PRECISION   DLAMCH, ZLANHE, ZLANHP
 148       EXTERNAL           LSAME, DLAMCH, ZLANHE, ZLANHP
 149 *     ..
 150 *     .. External Subroutines ..
 151       EXTERNAL           ZLASET, ZLAVHP
 152 *     ..
 153 *     .. Intrinsic Functions ..
 154       INTRINSIC          DBLE, DIMAG
 155 *     ..
 156 *     .. Executable Statements ..
 157 *
 158 *     Quick exit if N = 0.
 159 *
 160       IF( N.LE.0 ) THEN
 161          RESID = ZERO
 162          RETURN
 163       END IF
 164 *
 165 *     Determine EPS and the norm of A.
 166 *
 167       EPS = DLAMCH( 'Epsilon' )
 168       ANORM = ZLANHP( '1', UPLO, N, A, RWORK )
 169 *
 170 *     Check the imaginary parts of the diagonal elements and return with
 171 *     an error code if any are nonzero.
 172 *
 173       JC = 1
 174       IF( LSAME( UPLO, 'U' ) ) THEN
 175          DO 10 J = 1, N
 176             IF( DIMAG( AFAC( JC ) ).NE.ZERO ) THEN
 177                RESID = ONE / EPS
 178                RETURN
 179             END IF
 180             JC = JC + J + 1
 181    10    CONTINUE
 182       ELSE
 183          DO 20 J = 1, N
 184             IF( DIMAG( AFAC( JC ) ).NE.ZERO ) THEN
 185                RESID = ONE / EPS
 186                RETURN
 187             END IF
 188             JC = JC + N - J + 1
 189    20    CONTINUE
 190       END IF
 191 *
 192 *     Initialize C to the identity matrix.
 193 *
 194       CALL ZLASET( 'Full', N, N, CZERO, CONE, C, LDC )
 195 *
 196 *     Call ZLAVHP to form the product D * U' (or D * L' ).
 197 *
 198       CALL ZLAVHP( UPLO, 'Conjugate', 'Non-unit', N, N, AFAC, IPIV, C,
 199      $             LDC, INFO )
 200 *
 201 *     Call ZLAVHP again to multiply by U ( or L ).
 202 *
 203       CALL ZLAVHP( UPLO, 'No transpose', 'Unit', N, N, AFAC, IPIV, C,
 204      $             LDC, INFO )
 205 *
 206 *     Compute the difference  C - A .
 207 *
 208       IF( LSAME( UPLO, 'U' ) ) THEN
 209          JC = 0
 210          DO 40 J = 1, N
 211             DO 30 I = 1, J - 1
 212                C( I, J ) = C( I, J ) - A( JC+I )
 213    30       CONTINUE
 214             C( J, J ) = C( J, J ) - DBLE( A( JC+J ) )
 215             JC = JC + J
 216    40    CONTINUE
 217       ELSE
 218          JC = 1
 219          DO 60 J = 1, N
 220             C( J, J ) = C( J, J ) - DBLE( A( JC ) )
 221             DO 50 I = J + 1, N
 222                C( I, J ) = C( I, J ) - A( JC+I-J )
 223    50       CONTINUE
 224             JC = JC + N - J + 1
 225    60    CONTINUE
 226       END IF
 227 *
 228 *     Compute norm( C - A ) / ( N * norm(A) * EPS )
 229 *
 230       RESID = ZLANHE( '1', UPLO, N, C, LDC, RWORK )
 231 *
 232       IF( ANORM.LE.ZERO ) THEN
 233          IF( RESID.NE.ZERO )
 234      $      RESID = ONE / EPS
 235       ELSE
 236          RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
 237       END IF
 238 *
 239       RETURN
 240 *
 241 *     End of ZHPT01
 242 *
 243       END