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