![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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