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