![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DPST01 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 DPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM, 00012 * PIV, RWORK, RESID, RANK ) 00013 * 00014 * .. Scalar Arguments .. 00015 * DOUBLE PRECISION RESID 00016 * INTEGER LDA, LDAFAC, LDPERM, N, RANK 00017 * CHARACTER UPLO 00018 * .. 00019 * .. Array Arguments .. 00020 * DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), 00021 * $ PERM( LDPERM, * ), RWORK( * ) 00022 * INTEGER PIV( * ) 00023 * .. 00024 * 00025 * 00026 *> \par Purpose: 00027 * ============= 00028 *> 00029 *> \verbatim 00030 *> 00031 *> DPST01 reconstructs a symmetric positive semidefinite matrix A 00032 *> from its L or U factors and the permutation matrix P and computes 00033 *> the residual 00034 *> norm( P*L*L'*P' - A ) / ( N * norm(A) * EPS ) or 00035 *> norm( P*U'*U*P' - A ) / ( N * norm(A) * EPS ), 00036 *> where EPS is the machine epsilon. 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 *> symmetric 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 DOUBLE PRECISION array, dimension (LDA,N) 00060 *> The original symmetric 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 DOUBLE PRECISION array, dimension (LDAFAC,N) 00072 *> The factor L or U from the L*L' or U'*U 00073 *> factorization of A. 00074 *> \endverbatim 00075 *> 00076 *> \param[in] LDAFAC 00077 *> \verbatim 00078 *> LDAFAC is INTEGER 00079 *> The leading dimension of the array AFAC. LDAFAC >= max(1,N). 00080 *> \endverbatim 00081 *> 00082 *> \param[out] PERM 00083 *> \verbatim 00084 *> PERM is DOUBLE PRECISION array, dimension (LDPERM,N) 00085 *> Overwritten with the reconstructed matrix, and then with the 00086 *> difference P*L*L'*P' - A (or P*U'*U*P' - A) 00087 *> \endverbatim 00088 *> 00089 *> \param[in] LDPERM 00090 *> \verbatim 00091 *> LDPERM is INTEGER 00092 *> The leading dimension of the array PERM. 00093 *> LDAPERM >= max(1,N). 00094 *> \endverbatim 00095 *> 00096 *> \param[in] PIV 00097 *> \verbatim 00098 *> PIV is INTEGER array, dimension (N) 00099 *> PIV is such that the nonzero entries are 00100 *> P( PIV( K ), K ) = 1. 00101 *> \endverbatim 00102 *> 00103 *> \param[out] RWORK 00104 *> \verbatim 00105 *> RWORK is DOUBLE PRECISION array, dimension (N) 00106 *> \endverbatim 00107 *> 00108 *> \param[out] RESID 00109 *> \verbatim 00110 *> RESID is DOUBLE PRECISION 00111 *> If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS ) 00112 *> If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS ) 00113 *> \endverbatim 00114 *> 00115 *> \param[in] RANK 00116 *> \verbatim 00117 *> RANK is INTEGER 00118 *> number of nonzero singular values of A. 00119 *> \endverbatim 00120 * 00121 * Authors: 00122 * ======== 00123 * 00124 *> \author Univ. of Tennessee 00125 *> \author Univ. of California Berkeley 00126 *> \author Univ. of Colorado Denver 00127 *> \author NAG Ltd. 00128 * 00129 *> \date November 2011 00130 * 00131 *> \ingroup double_lin 00132 * 00133 * ===================================================================== 00134 SUBROUTINE DPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM, 00135 $ PIV, RWORK, RESID, RANK ) 00136 * 00137 * -- LAPACK test routine (version 3.4.0) -- 00138 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00139 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00140 * November 2011 00141 * 00142 * .. Scalar Arguments .. 00143 DOUBLE PRECISION RESID 00144 INTEGER LDA, LDAFAC, LDPERM, N, RANK 00145 CHARACTER UPLO 00146 * .. 00147 * .. Array Arguments .. 00148 DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), 00149 $ PERM( LDPERM, * ), RWORK( * ) 00150 INTEGER PIV( * ) 00151 * .. 00152 * 00153 * ===================================================================== 00154 * 00155 * .. Parameters .. 00156 DOUBLE PRECISION ZERO, ONE 00157 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00158 * .. 00159 * .. Local Scalars .. 00160 DOUBLE PRECISION ANORM, EPS, T 00161 INTEGER I, J, K 00162 * .. 00163 * .. External Functions .. 00164 DOUBLE PRECISION DDOT, DLAMCH, DLANSY 00165 LOGICAL LSAME 00166 EXTERNAL DDOT, DLAMCH, DLANSY, LSAME 00167 * .. 00168 * .. External Subroutines .. 00169 EXTERNAL DSCAL, DSYR, DTRMV 00170 * .. 00171 * .. Intrinsic Functions .. 00172 INTRINSIC DBLE 00173 * .. 00174 * .. Executable Statements .. 00175 * 00176 * Quick exit if N = 0. 00177 * 00178 IF( N.LE.0 ) THEN 00179 RESID = ZERO 00180 RETURN 00181 END IF 00182 * 00183 * Exit with RESID = 1/EPS if ANORM = 0. 00184 * 00185 EPS = DLAMCH( 'Epsilon' ) 00186 ANORM = DLANSY( '1', UPLO, N, A, LDA, RWORK ) 00187 IF( ANORM.LE.ZERO ) THEN 00188 RESID = ONE / EPS 00189 RETURN 00190 END IF 00191 * 00192 * Compute the product U'*U, overwriting U. 00193 * 00194 IF( LSAME( UPLO, 'U' ) ) THEN 00195 * 00196 IF( RANK.LT.N ) THEN 00197 DO 110 J = RANK + 1, N 00198 DO 100 I = RANK + 1, J 00199 AFAC( I, J ) = ZERO 00200 100 CONTINUE 00201 110 CONTINUE 00202 END IF 00203 * 00204 DO 120 K = N, 1, -1 00205 * 00206 * Compute the (K,K) element of the result. 00207 * 00208 T = DDOT( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 ) 00209 AFAC( K, K ) = T 00210 * 00211 * Compute the rest of column K. 00212 * 00213 CALL DTRMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC, 00214 $ LDAFAC, AFAC( 1, K ), 1 ) 00215 * 00216 120 CONTINUE 00217 * 00218 * Compute the product L*L', overwriting L. 00219 * 00220 ELSE 00221 * 00222 IF( RANK.LT.N ) THEN 00223 DO 140 J = RANK + 1, N 00224 DO 130 I = J, N 00225 AFAC( I, J ) = ZERO 00226 130 CONTINUE 00227 140 CONTINUE 00228 END IF 00229 * 00230 DO 150 K = N, 1, -1 00231 * Add a multiple of column K of the factor L to each of 00232 * columns K+1 through N. 00233 * 00234 IF( K+1.LE.N ) 00235 $ CALL DSYR( 'Lower', N-K, ONE, AFAC( K+1, K ), 1, 00236 $ AFAC( K+1, K+1 ), LDAFAC ) 00237 * 00238 * Scale column K by the diagonal element. 00239 * 00240 T = AFAC( K, K ) 00241 CALL DSCAL( N-K+1, T, AFAC( K, K ), 1 ) 00242 150 CONTINUE 00243 * 00244 END IF 00245 * 00246 * Form P*L*L'*P' or P*U'*U*P' 00247 * 00248 IF( LSAME( UPLO, 'U' ) ) THEN 00249 * 00250 DO 170 J = 1, N 00251 DO 160 I = 1, N 00252 IF( PIV( I ).LE.PIV( J ) ) THEN 00253 IF( I.LE.J ) THEN 00254 PERM( PIV( I ), PIV( J ) ) = AFAC( I, J ) 00255 ELSE 00256 PERM( PIV( I ), PIV( J ) ) = AFAC( J, I ) 00257 END IF 00258 END IF 00259 160 CONTINUE 00260 170 CONTINUE 00261 * 00262 * 00263 ELSE 00264 * 00265 DO 190 J = 1, N 00266 DO 180 I = 1, N 00267 IF( PIV( I ).GE.PIV( J ) ) THEN 00268 IF( I.GE.J ) THEN 00269 PERM( PIV( I ), PIV( J ) ) = AFAC( I, J ) 00270 ELSE 00271 PERM( PIV( I ), PIV( J ) ) = AFAC( J, I ) 00272 END IF 00273 END IF 00274 180 CONTINUE 00275 190 CONTINUE 00276 * 00277 END IF 00278 * 00279 * Compute the difference P*L*L'*P' - A (or P*U'*U*P' - A). 00280 * 00281 IF( LSAME( UPLO, 'U' ) ) THEN 00282 DO 210 J = 1, N 00283 DO 200 I = 1, J 00284 PERM( I, J ) = PERM( I, J ) - A( I, J ) 00285 200 CONTINUE 00286 210 CONTINUE 00287 ELSE 00288 DO 230 J = 1, N 00289 DO 220 I = J, N 00290 PERM( I, J ) = PERM( I, J ) - A( I, J ) 00291 220 CONTINUE 00292 230 CONTINUE 00293 END IF 00294 * 00295 * Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or 00296 * ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ). 00297 * 00298 RESID = DLANSY( '1', UPLO, N, PERM, LDAFAC, RWORK ) 00299 * 00300 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS 00301 * 00302 RETURN 00303 * 00304 * End of DPST01 00305 * 00306 END