LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zhet01.f
Go to the documentation of this file.
00001 *> \brief \b ZHET01
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 ZHET01( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, LDC,
00012 *                          RWORK, RESID )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          UPLO
00016 *       INTEGER            LDA, LDAFAC, LDC, N
00017 *       DOUBLE PRECISION   RESID
00018 *       ..
00019 *       .. Array Arguments ..
00020 *       INTEGER            IPIV( * )
00021 *       DOUBLE PRECISION   RWORK( * )
00022 *       COMPLEX*16         A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * )
00023 *       ..
00024 *  
00025 *
00026 *> \par Purpose:
00027 *  =============
00028 *>
00029 *> \verbatim
00030 *>
00031 *> ZHET01 reconstructs a Hermitian indefinite matrix A from its
00032 *> block L*D*L' or U*D*U' factorization and computes the residual
00033 *>    norm( C - A ) / ( N * norm(A) * EPS ),
00034 *> where C is the reconstructed matrix, EPS is the machine epsilon,
00035 *> L' is the conjugate transpose of L, and U' is the conjugate transpose
00036 *> of U.
00037 *> \endverbatim
00038 *
00039 *  Arguments:
00040 *  ==========
00041 *
00042 *> \param[in] UPLO
00043 *> \verbatim
00044 *>          UPLO is CHARACTER*1
00045 *>          Specifies whether the upper or lower triangular part of the
00046 *>          Hermitian matrix A is stored:
00047 *>          = 'U':  Upper triangular
00048 *>          = 'L':  Lower triangular
00049 *> \endverbatim
00050 *>
00051 *> \param[in] N
00052 *> \verbatim
00053 *>          N is INTEGER
00054 *>          The number of rows and columns of the matrix A.  N >= 0.
00055 *> \endverbatim
00056 *>
00057 *> \param[in] A
00058 *> \verbatim
00059 *>          A is COMPLEX*16 array, dimension (LDA,N)
00060 *>          The original Hermitian matrix A.
00061 *> \endverbatim
00062 *>
00063 *> \param[in] LDA
00064 *> \verbatim
00065 *>          LDA is INTEGER
00066 *>          The leading dimension of the array A.  LDA >= max(1,N)
00067 *> \endverbatim
00068 *>
00069 *> \param[in] AFAC
00070 *> \verbatim
00071 *>          AFAC is COMPLEX*16 array, dimension (LDAFAC,N)
00072 *>          The factored form of the matrix A.  AFAC contains the block
00073 *>          diagonal matrix D and the multipliers used to obtain the
00074 *>          factor L or U from the block L*D*L' or U*D*U' factorization
00075 *>          as computed by ZHETRF.
00076 *> \endverbatim
00077 *>
00078 *> \param[in] LDAFAC
00079 *> \verbatim
00080 *>          LDAFAC is INTEGER
00081 *>          The leading dimension of the array AFAC.  LDAFAC >= max(1,N).
00082 *> \endverbatim
00083 *>
00084 *> \param[in] IPIV
00085 *> \verbatim
00086 *>          IPIV is INTEGER array, dimension (N)
00087 *>          The pivot indices from ZHETRF.
00088 *> \endverbatim
00089 *>
00090 *> \param[out] C
00091 *> \verbatim
00092 *>          C is COMPLEX*16 array, dimension (LDC,N)
00093 *> \endverbatim
00094 *>
00095 *> \param[in] LDC
00096 *> \verbatim
00097 *>          LDC is INTEGER
00098 *>          The leading dimension of the array C.  LDC >= max(1,N).
00099 *> \endverbatim
00100 *>
00101 *> \param[out] RWORK
00102 *> \verbatim
00103 *>          RWORK is DOUBLE PRECISION array, dimension (N)
00104 *> \endverbatim
00105 *>
00106 *> \param[out] RESID
00107 *> \verbatim
00108 *>          RESID is DOUBLE PRECISION
00109 *>          If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS )
00110 *>          If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS )
00111 *> \endverbatim
00112 *
00113 *  Authors:
00114 *  ========
00115 *
00116 *> \author Univ. of Tennessee 
00117 *> \author Univ. of California Berkeley 
00118 *> \author Univ. of Colorado Denver 
00119 *> \author NAG Ltd. 
00120 *
00121 *> \date November 2011
00122 *
00123 *> \ingroup complex16_lin
00124 *
00125 *  =====================================================================
00126       SUBROUTINE ZHET01( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, LDC,
00127      $                   RWORK, RESID )
00128 *
00129 *  -- LAPACK test routine (version 3.4.0) --
00130 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00131 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00132 *     November 2011
00133 *
00134 *     .. Scalar Arguments ..
00135       CHARACTER          UPLO
00136       INTEGER            LDA, LDAFAC, LDC, N
00137       DOUBLE PRECISION   RESID
00138 *     ..
00139 *     .. Array Arguments ..
00140       INTEGER            IPIV( * )
00141       DOUBLE PRECISION   RWORK( * )
00142       COMPLEX*16         A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * )
00143 *     ..
00144 *
00145 *  =====================================================================
00146 *
00147 *     .. Parameters ..
00148       DOUBLE PRECISION   ZERO, ONE
00149       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00150       COMPLEX*16         CZERO, CONE
00151       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00152      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00153 *     ..
00154 *     .. Local Scalars ..
00155       INTEGER            I, INFO, J
00156       DOUBLE PRECISION   ANORM, EPS
00157 *     ..
00158 *     .. External Functions ..
00159       LOGICAL            LSAME
00160       DOUBLE PRECISION   DLAMCH, ZLANHE
00161       EXTERNAL           LSAME, DLAMCH, ZLANHE
00162 *     ..
00163 *     .. External Subroutines ..
00164       EXTERNAL           ZLASET, ZLAVHE
00165 *     ..
00166 *     .. Intrinsic Functions ..
00167       INTRINSIC          DBLE, DIMAG
00168 *     ..
00169 *     .. Executable Statements ..
00170 *
00171 *     Quick exit if N = 0.
00172 *
00173       IF( N.LE.0 ) THEN
00174          RESID = ZERO
00175          RETURN
00176       END IF
00177 *
00178 *     Determine EPS and the norm of A.
00179 *
00180       EPS = DLAMCH( 'Epsilon' )
00181       ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK )
00182 *
00183 *     Check the imaginary parts of the diagonal elements and return with
00184 *     an error code if any are nonzero.
00185 *
00186       DO 10 J = 1, N
00187          IF( DIMAG( AFAC( J, J ) ).NE.ZERO ) THEN
00188             RESID = ONE / EPS
00189             RETURN
00190          END IF
00191    10 CONTINUE
00192 *
00193 *     Initialize C to the identity matrix.
00194 *
00195       CALL ZLASET( 'Full', N, N, CZERO, CONE, C, LDC )
00196 *
00197 *     Call ZLAVHE to form the product D * U' (or D * L' ).
00198 *
00199       CALL ZLAVHE( UPLO, 'Conjugate', 'Non-unit', N, N, AFAC, LDAFAC,
00200      $             IPIV, C, LDC, INFO )
00201 *
00202 *     Call ZLAVHE again to multiply by U (or L ).
00203 *
00204       CALL ZLAVHE( UPLO, 'No transpose', 'Unit', N, N, AFAC, LDAFAC,
00205      $             IPIV, C, LDC, INFO )
00206 *
00207 *     Compute the difference  C - A .
00208 *
00209       IF( LSAME( UPLO, 'U' ) ) THEN
00210          DO 30 J = 1, N
00211             DO 20 I = 1, J - 1
00212                C( I, J ) = C( I, J ) - A( I, J )
00213    20       CONTINUE
00214             C( J, J ) = C( J, J ) - DBLE( A( J, J ) )
00215    30    CONTINUE
00216       ELSE
00217          DO 50 J = 1, N
00218             C( J, J ) = C( J, J ) - DBLE( A( J, J ) )
00219             DO 40 I = J + 1, N
00220                C( I, J ) = C( I, J ) - A( I, J )
00221    40       CONTINUE
00222    50    CONTINUE
00223       END IF
00224 *
00225 *     Compute norm( C - A ) / ( N * norm(A) * EPS )
00226 *
00227       RESID = ZLANHE( '1', UPLO, N, C, LDC, RWORK )
00228 *
00229       IF( ANORM.LE.ZERO ) THEN
00230          IF( RESID.NE.ZERO )
00231      $      RESID = ONE / EPS
00232       ELSE
00233          RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
00234       END IF
00235 *
00236       RETURN
00237 *
00238 *     End of ZHET01
00239 *
00240       END
 All Files Functions