![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DPOT03 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 DPOT03( UPLO, N, A, LDA, AINV, LDAINV, WORK, LDWORK, 00012 * RWORK, RCOND, RESID ) 00013 * 00014 * .. Scalar Arguments .. 00015 * CHARACTER UPLO 00016 * INTEGER LDA, LDAINV, LDWORK, N 00017 * DOUBLE PRECISION RCOND, RESID 00018 * .. 00019 * .. Array Arguments .. 00020 * DOUBLE PRECISION A( LDA, * ), AINV( LDAINV, * ), RWORK( * ), 00021 * $ WORK( LDWORK, * ) 00022 * .. 00023 * 00024 * 00025 *> \par Purpose: 00026 * ============= 00027 *> 00028 *> \verbatim 00029 *> 00030 *> DPOT03 computes the residual for a symmetric matrix times its 00031 *> inverse: 00032 *> norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ), 00033 *> where EPS is the machine epsilon. 00034 *> \endverbatim 00035 * 00036 * Arguments: 00037 * ========== 00038 * 00039 *> \param[in] UPLO 00040 *> \verbatim 00041 *> UPLO is CHARACTER*1 00042 *> Specifies whether the upper or lower triangular part of the 00043 *> symmetric matrix A is stored: 00044 *> = 'U': Upper triangular 00045 *> = 'L': Lower triangular 00046 *> \endverbatim 00047 *> 00048 *> \param[in] N 00049 *> \verbatim 00050 *> N is INTEGER 00051 *> The number of rows and columns of the matrix A. N >= 0. 00052 *> \endverbatim 00053 *> 00054 *> \param[in] A 00055 *> \verbatim 00056 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00057 *> The original symmetric matrix A. 00058 *> \endverbatim 00059 *> 00060 *> \param[in] LDA 00061 *> \verbatim 00062 *> LDA is INTEGER 00063 *> The leading dimension of the array A. LDA >= max(1,N) 00064 *> \endverbatim 00065 *> 00066 *> \param[in,out] AINV 00067 *> \verbatim 00068 *> AINV is DOUBLE PRECISION array, dimension (LDAINV,N) 00069 *> On entry, the inverse of the matrix A, stored as a symmetric 00070 *> matrix in the same format as A. 00071 *> In this version, AINV is expanded into a full matrix and 00072 *> multiplied by A, so the opposing triangle of AINV will be 00073 *> changed; i.e., if the upper triangular part of AINV is 00074 *> stored, the lower triangular part will be used as work space. 00075 *> \endverbatim 00076 *> 00077 *> \param[in] LDAINV 00078 *> \verbatim 00079 *> LDAINV is INTEGER 00080 *> The leading dimension of the array AINV. LDAINV >= max(1,N). 00081 *> \endverbatim 00082 *> 00083 *> \param[out] WORK 00084 *> \verbatim 00085 *> WORK is DOUBLE PRECISION array, dimension (LDWORK,N) 00086 *> \endverbatim 00087 *> 00088 *> \param[in] LDWORK 00089 *> \verbatim 00090 *> LDWORK is INTEGER 00091 *> The leading dimension of the array WORK. LDWORK >= max(1,N). 00092 *> \endverbatim 00093 *> 00094 *> \param[out] RWORK 00095 *> \verbatim 00096 *> RWORK is DOUBLE PRECISION array, dimension (N) 00097 *> \endverbatim 00098 *> 00099 *> \param[out] RCOND 00100 *> \verbatim 00101 *> RCOND is DOUBLE PRECISION 00102 *> The reciprocal of the condition number of A, computed as 00103 *> ( 1/norm(A) ) / norm(AINV). 00104 *> \endverbatim 00105 *> 00106 *> \param[out] RESID 00107 *> \verbatim 00108 *> RESID is DOUBLE PRECISION 00109 *> norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS ) 00110 *> \endverbatim 00111 * 00112 * Authors: 00113 * ======== 00114 * 00115 *> \author Univ. of Tennessee 00116 *> \author Univ. of California Berkeley 00117 *> \author Univ. of Colorado Denver 00118 *> \author NAG Ltd. 00119 * 00120 *> \date November 2011 00121 * 00122 *> \ingroup double_lin 00123 * 00124 * ===================================================================== 00125 SUBROUTINE DPOT03( UPLO, N, A, LDA, AINV, LDAINV, WORK, LDWORK, 00126 $ RWORK, RCOND, RESID ) 00127 * 00128 * -- LAPACK test routine (version 3.4.0) -- 00129 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00130 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00131 * November 2011 00132 * 00133 * .. Scalar Arguments .. 00134 CHARACTER UPLO 00135 INTEGER LDA, LDAINV, LDWORK, N 00136 DOUBLE PRECISION RCOND, RESID 00137 * .. 00138 * .. Array Arguments .. 00139 DOUBLE PRECISION A( LDA, * ), AINV( LDAINV, * ), RWORK( * ), 00140 $ WORK( LDWORK, * ) 00141 * .. 00142 * 00143 * ===================================================================== 00144 * 00145 * .. Parameters .. 00146 DOUBLE PRECISION ZERO, ONE 00147 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00148 * .. 00149 * .. Local Scalars .. 00150 INTEGER I, J 00151 DOUBLE PRECISION AINVNM, ANORM, EPS 00152 * .. 00153 * .. External Functions .. 00154 LOGICAL LSAME 00155 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY 00156 EXTERNAL LSAME, DLAMCH, DLANGE, DLANSY 00157 * .. 00158 * .. External Subroutines .. 00159 EXTERNAL DSYMM 00160 * .. 00161 * .. Intrinsic Functions .. 00162 INTRINSIC DBLE 00163 * .. 00164 * .. Executable Statements .. 00165 * 00166 * Quick exit if N = 0. 00167 * 00168 IF( N.LE.0 ) THEN 00169 RCOND = ONE 00170 RESID = ZERO 00171 RETURN 00172 END IF 00173 * 00174 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 00175 * 00176 EPS = DLAMCH( 'Epsilon' ) 00177 ANORM = DLANSY( '1', UPLO, N, A, LDA, RWORK ) 00178 AINVNM = DLANSY( '1', UPLO, N, AINV, LDAINV, RWORK ) 00179 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00180 RCOND = ZERO 00181 RESID = ONE / EPS 00182 RETURN 00183 END IF 00184 RCOND = ( ONE / ANORM ) / AINVNM 00185 * 00186 * Expand AINV into a full matrix and call DSYMM to multiply 00187 * AINV on the left by A. 00188 * 00189 IF( LSAME( UPLO, 'U' ) ) THEN 00190 DO 20 J = 1, N 00191 DO 10 I = 1, J - 1 00192 AINV( J, I ) = AINV( I, J ) 00193 10 CONTINUE 00194 20 CONTINUE 00195 ELSE 00196 DO 40 J = 1, N 00197 DO 30 I = J + 1, N 00198 AINV( J, I ) = AINV( I, J ) 00199 30 CONTINUE 00200 40 CONTINUE 00201 END IF 00202 CALL DSYMM( 'Left', UPLO, N, N, -ONE, A, LDA, AINV, LDAINV, ZERO, 00203 $ WORK, LDWORK ) 00204 * 00205 * Add the identity matrix to WORK . 00206 * 00207 DO 50 I = 1, N 00208 WORK( I, I ) = WORK( I, I ) + ONE 00209 50 CONTINUE 00210 * 00211 * Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS) 00212 * 00213 RESID = DLANGE( '1', N, N, WORK, LDWORK, RWORK ) 00214 * 00215 RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N ) 00216 * 00217 RETURN 00218 * 00219 * End of DPOT03 00220 * 00221 END