![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CPST01 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 CPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM, 00012 * PIV, RWORK, RESID, RANK ) 00013 * 00014 * .. Scalar Arguments .. 00015 * REAL RESID 00016 * INTEGER LDA, LDAFAC, LDPERM, N, RANK 00017 * CHARACTER UPLO 00018 * .. 00019 * .. Array Arguments .. 00020 * COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), 00021 * $ PERM( LDPERM, * ) 00022 * REAL RWORK( * ) 00023 * INTEGER PIV( * ) 00024 * .. 00025 * 00026 * 00027 *> \par Purpose: 00028 * ============= 00029 *> 00030 *> \verbatim 00031 *> 00032 *> CPST01 reconstructs an Hermitian positive semidefinite matrix A 00033 *> from its L or U factors and the permutation matrix P and computes 00034 *> the residual 00035 *> norm( P*L*L'*P' - A ) / ( N * norm(A) * EPS ) or 00036 *> norm( P*U'*U*P' - A ) / ( N * norm(A) * EPS ), 00037 *> where EPS is the machine epsilon, L' is the conjugate transpose of L, 00038 *> and U' is the conjugate transpose of U. 00039 *> \endverbatim 00040 * 00041 * Arguments: 00042 * ========== 00043 * 00044 *> \param[in] UPLO 00045 *> \verbatim 00046 *> UPLO is CHARACTER*1 00047 *> Specifies whether the upper or lower triangular part of the 00048 *> Hermitian matrix A is stored: 00049 *> = 'U': Upper triangular 00050 *> = 'L': Lower triangular 00051 *> \endverbatim 00052 *> 00053 *> \param[in] N 00054 *> \verbatim 00055 *> N is INTEGER 00056 *> The number of rows and columns of the matrix A. N >= 0. 00057 *> \endverbatim 00058 *> 00059 *> \param[in] A 00060 *> \verbatim 00061 *> A is COMPLEX array, dimension (LDA,N) 00062 *> The original Hermitian matrix A. 00063 *> \endverbatim 00064 *> 00065 *> \param[in] LDA 00066 *> \verbatim 00067 *> LDA is INTEGER 00068 *> The leading dimension of the array A. LDA >= max(1,N) 00069 *> \endverbatim 00070 *> 00071 *> \param[in] AFAC 00072 *> \verbatim 00073 *> AFAC is COMPLEX array, dimension (LDAFAC,N) 00074 *> The factor L or U from the L*L' or U'*U 00075 *> factorization of A. 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[out] PERM 00085 *> \verbatim 00086 *> PERM is COMPLEX array, dimension (LDPERM,N) 00087 *> Overwritten with the reconstructed matrix, and then with the 00088 *> difference P*L*L'*P' - A (or P*U'*U*P' - A) 00089 *> \endverbatim 00090 *> 00091 *> \param[in] LDPERM 00092 *> \verbatim 00093 *> LDPERM is INTEGER 00094 *> The leading dimension of the array PERM. 00095 *> LDAPERM >= max(1,N). 00096 *> \endverbatim 00097 *> 00098 *> \param[in] PIV 00099 *> \verbatim 00100 *> PIV is INTEGER array, dimension (N) 00101 *> PIV is such that the nonzero entries are 00102 *> P( PIV( K ), K ) = 1. 00103 *> \endverbatim 00104 *> 00105 *> \param[out] RWORK 00106 *> \verbatim 00107 *> RWORK is REAL array, dimension (N) 00108 *> \endverbatim 00109 *> 00110 *> \param[out] RESID 00111 *> \verbatim 00112 *> RESID is REAL 00113 *> If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS ) 00114 *> If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS ) 00115 *> \endverbatim 00116 *> 00117 *> \param[in] RANK 00118 *> \verbatim 00119 *> RANK is INTEGER 00120 *> number of nonzero singular values of A. 00121 *> \endverbatim 00122 * 00123 * Authors: 00124 * ======== 00125 * 00126 *> \author Univ. of Tennessee 00127 *> \author Univ. of California Berkeley 00128 *> \author Univ. of Colorado Denver 00129 *> \author NAG Ltd. 00130 * 00131 *> \date November 2011 00132 * 00133 *> \ingroup complex_lin 00134 * 00135 * ===================================================================== 00136 SUBROUTINE CPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM, 00137 $ PIV, RWORK, RESID, RANK ) 00138 * 00139 * -- LAPACK test routine (version 3.4.0) -- 00140 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00141 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00142 * November 2011 00143 * 00144 * .. Scalar Arguments .. 00145 REAL RESID 00146 INTEGER LDA, LDAFAC, LDPERM, N, RANK 00147 CHARACTER UPLO 00148 * .. 00149 * .. Array Arguments .. 00150 COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), 00151 $ PERM( LDPERM, * ) 00152 REAL RWORK( * ) 00153 INTEGER PIV( * ) 00154 * .. 00155 * 00156 * ===================================================================== 00157 * 00158 * .. Parameters .. 00159 REAL ZERO, ONE 00160 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00161 COMPLEX CZERO 00162 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) ) 00163 * .. 00164 * .. Local Scalars .. 00165 COMPLEX TC 00166 REAL ANORM, EPS, TR 00167 INTEGER I, J, K 00168 * .. 00169 * .. External Functions .. 00170 COMPLEX CDOTC 00171 REAL CLANHE, SLAMCH 00172 LOGICAL LSAME 00173 EXTERNAL CDOTC, CLANHE, SLAMCH, LSAME 00174 * .. 00175 * .. External Subroutines .. 00176 EXTERNAL CHER, CSCAL, CTRMV 00177 * .. 00178 * .. Intrinsic Functions .. 00179 INTRINSIC AIMAG, CONJG, REAL 00180 * .. 00181 * .. Executable Statements .. 00182 * 00183 * Quick exit if N = 0. 00184 * 00185 IF( N.LE.0 ) THEN 00186 RESID = ZERO 00187 RETURN 00188 END IF 00189 * 00190 * Exit with RESID = 1/EPS if ANORM = 0. 00191 * 00192 EPS = SLAMCH( 'Epsilon' ) 00193 ANORM = CLANHE( '1', UPLO, N, A, LDA, RWORK ) 00194 IF( ANORM.LE.ZERO ) THEN 00195 RESID = ONE / EPS 00196 RETURN 00197 END IF 00198 * 00199 * Check the imaginary parts of the diagonal elements and return with 00200 * an error code if any are nonzero. 00201 * 00202 DO 100 J = 1, N 00203 IF( AIMAG( AFAC( J, J ) ).NE.ZERO ) THEN 00204 RESID = ONE / EPS 00205 RETURN 00206 END IF 00207 100 CONTINUE 00208 * 00209 * Compute the product U'*U, overwriting U. 00210 * 00211 IF( LSAME( UPLO, 'U' ) ) THEN 00212 * 00213 IF( RANK.LT.N ) THEN 00214 DO 120 J = RANK + 1, N 00215 DO 110 I = RANK + 1, J 00216 AFAC( I, J ) = CZERO 00217 110 CONTINUE 00218 120 CONTINUE 00219 END IF 00220 * 00221 DO 130 K = N, 1, -1 00222 * 00223 * Compute the (K,K) element of the result. 00224 * 00225 TR = CDOTC( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 ) 00226 AFAC( K, K ) = TR 00227 * 00228 * Compute the rest of column K. 00229 * 00230 CALL CTRMV( 'Upper', 'Conjugate', 'Non-unit', K-1, AFAC, 00231 $ LDAFAC, AFAC( 1, K ), 1 ) 00232 * 00233 130 CONTINUE 00234 * 00235 * Compute the product L*L', overwriting L. 00236 * 00237 ELSE 00238 * 00239 IF( RANK.LT.N ) THEN 00240 DO 150 J = RANK + 1, N 00241 DO 140 I = J, N 00242 AFAC( I, J ) = CZERO 00243 140 CONTINUE 00244 150 CONTINUE 00245 END IF 00246 * 00247 DO 160 K = N, 1, -1 00248 * Add a multiple of column K of the factor L to each of 00249 * columns K+1 through N. 00250 * 00251 IF( K+1.LE.N ) 00252 $ CALL CHER( 'Lower', N-K, ONE, AFAC( K+1, K ), 1, 00253 $ AFAC( K+1, K+1 ), LDAFAC ) 00254 * 00255 * Scale column K by the diagonal element. 00256 * 00257 TC = AFAC( K, K ) 00258 CALL CSCAL( N-K+1, TC, AFAC( K, K ), 1 ) 00259 160 CONTINUE 00260 * 00261 END IF 00262 * 00263 * Form P*L*L'*P' or P*U'*U*P' 00264 * 00265 IF( LSAME( UPLO, 'U' ) ) THEN 00266 * 00267 DO 180 J = 1, N 00268 DO 170 I = 1, N 00269 IF( PIV( I ).LE.PIV( J ) ) THEN 00270 IF( I.LE.J ) THEN 00271 PERM( PIV( I ), PIV( J ) ) = AFAC( I, J ) 00272 ELSE 00273 PERM( PIV( I ), PIV( J ) ) = CONJG( AFAC( J, I ) ) 00274 END IF 00275 END IF 00276 170 CONTINUE 00277 180 CONTINUE 00278 * 00279 * 00280 ELSE 00281 * 00282 DO 200 J = 1, N 00283 DO 190 I = 1, N 00284 IF( PIV( I ).GE.PIV( J ) ) THEN 00285 IF( I.GE.J ) THEN 00286 PERM( PIV( I ), PIV( J ) ) = AFAC( I, J ) 00287 ELSE 00288 PERM( PIV( I ), PIV( J ) ) = CONJG( AFAC( J, I ) ) 00289 END IF 00290 END IF 00291 190 CONTINUE 00292 200 CONTINUE 00293 * 00294 END IF 00295 * 00296 * Compute the difference P*L*L'*P' - A (or P*U'*U*P' - A). 00297 * 00298 IF( LSAME( UPLO, 'U' ) ) THEN 00299 DO 220 J = 1, N 00300 DO 210 I = 1, J - 1 00301 PERM( I, J ) = PERM( I, J ) - A( I, J ) 00302 210 CONTINUE 00303 PERM( J, J ) = PERM( J, J ) - REAL( A( J, J ) ) 00304 220 CONTINUE 00305 ELSE 00306 DO 240 J = 1, N 00307 PERM( J, J ) = PERM( J, J ) - REAL( A( J, J ) ) 00308 DO 230 I = J + 1, N 00309 PERM( I, J ) = PERM( I, J ) - A( I, J ) 00310 230 CONTINUE 00311 240 CONTINUE 00312 END IF 00313 * 00314 * Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or 00315 * ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ). 00316 * 00317 RESID = CLANHE( '1', UPLO, N, PERM, LDAFAC, RWORK ) 00318 * 00319 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS 00320 * 00321 RETURN 00322 * 00323 * End of CPST01 00324 * 00325 END