![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZSPT03 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 ZSPT03( UPLO, N, A, AINV, WORK, LDW, RWORK, RCOND, 00012 * RESID ) 00013 * 00014 * .. Scalar Arguments .. 00015 * CHARACTER UPLO 00016 * INTEGER LDW, N 00017 * DOUBLE PRECISION RCOND, RESID 00018 * .. 00019 * .. Array Arguments .. 00020 * DOUBLE PRECISION RWORK( * ) 00021 * COMPLEX*16 A( * ), AINV( * ), WORK( LDW, * ) 00022 * .. 00023 * 00024 * 00025 *> \par Purpose: 00026 * ============= 00027 *> 00028 *> \verbatim 00029 *> 00030 *> ZSPT03 computes the residual for a complex symmetric packed matrix 00031 *> times its 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 *> complex 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 COMPLEX*16 array, dimension (N*(N+1)/2) 00057 *> The original complex symmetric matrix A, stored as a packed 00058 *> triangular matrix. 00059 *> \endverbatim 00060 *> 00061 *> \param[in] AINV 00062 *> \verbatim 00063 *> AINV is COMPLEX*16 array, dimension (N*(N+1)/2) 00064 *> The (symmetric) inverse of the matrix A, stored as a packed 00065 *> triangular matrix. 00066 *> \endverbatim 00067 *> 00068 *> \param[out] WORK 00069 *> \verbatim 00070 *> WORK is COMPLEX*16 array, dimension (LDW,N) 00071 *> \endverbatim 00072 *> 00073 *> \param[in] LDW 00074 *> \verbatim 00075 *> LDW is INTEGER 00076 *> The leading dimension of the array WORK. LDW >= max(1,N). 00077 *> \endverbatim 00078 *> 00079 *> \param[out] RWORK 00080 *> \verbatim 00081 *> RWORK is DOUBLE PRECISION array, dimension (N) 00082 *> \endverbatim 00083 *> 00084 *> \param[out] RCOND 00085 *> \verbatim 00086 *> RCOND is DOUBLE PRECISION 00087 *> The reciprocal of the condition number of A, computed as 00088 *> ( 1/norm(A) ) / norm(AINV). 00089 *> \endverbatim 00090 *> 00091 *> \param[out] RESID 00092 *> \verbatim 00093 *> RESID is DOUBLE PRECISION 00094 *> norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS ) 00095 *> \endverbatim 00096 * 00097 * Authors: 00098 * ======== 00099 * 00100 *> \author Univ. of Tennessee 00101 *> \author Univ. of California Berkeley 00102 *> \author Univ. of Colorado Denver 00103 *> \author NAG Ltd. 00104 * 00105 *> \date November 2011 00106 * 00107 *> \ingroup complex16_lin 00108 * 00109 * ===================================================================== 00110 SUBROUTINE ZSPT03( UPLO, N, A, AINV, WORK, LDW, RWORK, RCOND, 00111 $ 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 LDW, N 00121 DOUBLE PRECISION RCOND, RESID 00122 * .. 00123 * .. Array Arguments .. 00124 DOUBLE PRECISION RWORK( * ) 00125 COMPLEX*16 A( * ), AINV( * ), WORK( LDW, * ) 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, ICOL, J, JCOL, K, KCOL, NALL 00136 DOUBLE PRECISION AINVNM, ANORM, EPS 00137 COMPLEX*16 T 00138 * .. 00139 * .. External Functions .. 00140 LOGICAL LSAME 00141 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANSP 00142 COMPLEX*16 ZDOTU 00143 EXTERNAL LSAME, DLAMCH, ZLANGE, ZLANSP, ZDOTU 00144 * .. 00145 * .. Intrinsic Functions .. 00146 INTRINSIC DBLE 00147 * .. 00148 * .. Executable Statements .. 00149 * 00150 * Quick exit if N = 0. 00151 * 00152 IF( N.LE.0 ) THEN 00153 RCOND = ONE 00154 RESID = ZERO 00155 RETURN 00156 END IF 00157 * 00158 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 00159 * 00160 EPS = DLAMCH( 'Epsilon' ) 00161 ANORM = ZLANSP( '1', UPLO, N, A, RWORK ) 00162 AINVNM = ZLANSP( '1', UPLO, N, AINV, RWORK ) 00163 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00164 RCOND = ZERO 00165 RESID = ONE / EPS 00166 RETURN 00167 END IF 00168 RCOND = ( ONE / ANORM ) / AINVNM 00169 * 00170 * Case where both A and AINV are upper triangular: 00171 * Each element of - A * AINV is computed by taking the dot product 00172 * of a row of A with a column of AINV. 00173 * 00174 IF( LSAME( UPLO, 'U' ) ) THEN 00175 DO 70 I = 1, N 00176 ICOL = ( ( I-1 )*I ) / 2 + 1 00177 * 00178 * Code when J <= I 00179 * 00180 DO 30 J = 1, I 00181 JCOL = ( ( J-1 )*J ) / 2 + 1 00182 T = ZDOTU( J, A( ICOL ), 1, AINV( JCOL ), 1 ) 00183 JCOL = JCOL + 2*J - 1 00184 KCOL = ICOL - 1 00185 DO 10 K = J + 1, I 00186 T = T + A( KCOL+K )*AINV( JCOL ) 00187 JCOL = JCOL + K 00188 10 CONTINUE 00189 KCOL = KCOL + 2*I 00190 DO 20 K = I + 1, N 00191 T = T + A( KCOL )*AINV( JCOL ) 00192 KCOL = KCOL + K 00193 JCOL = JCOL + K 00194 20 CONTINUE 00195 WORK( I, J ) = -T 00196 30 CONTINUE 00197 * 00198 * Code when J > I 00199 * 00200 DO 60 J = I + 1, N 00201 JCOL = ( ( J-1 )*J ) / 2 + 1 00202 T = ZDOTU( I, A( ICOL ), 1, AINV( JCOL ), 1 ) 00203 JCOL = JCOL - 1 00204 KCOL = ICOL + 2*I - 1 00205 DO 40 K = I + 1, J 00206 T = T + A( KCOL )*AINV( JCOL+K ) 00207 KCOL = KCOL + K 00208 40 CONTINUE 00209 JCOL = JCOL + 2*J 00210 DO 50 K = J + 1, N 00211 T = T + A( KCOL )*AINV( JCOL ) 00212 KCOL = KCOL + K 00213 JCOL = JCOL + K 00214 50 CONTINUE 00215 WORK( I, J ) = -T 00216 60 CONTINUE 00217 70 CONTINUE 00218 ELSE 00219 * 00220 * Case where both A and AINV are lower triangular 00221 * 00222 NALL = ( N*( N+1 ) ) / 2 00223 DO 140 I = 1, N 00224 * 00225 * Code when J <= I 00226 * 00227 ICOL = NALL - ( ( N-I+1 )*( N-I+2 ) ) / 2 + 1 00228 DO 100 J = 1, I 00229 JCOL = NALL - ( ( N-J )*( N-J+1 ) ) / 2 - ( N-I ) 00230 T = ZDOTU( N-I+1, A( ICOL ), 1, AINV( JCOL ), 1 ) 00231 KCOL = I 00232 JCOL = J 00233 DO 80 K = 1, J - 1 00234 T = T + A( KCOL )*AINV( JCOL ) 00235 JCOL = JCOL + N - K 00236 KCOL = KCOL + N - K 00237 80 CONTINUE 00238 JCOL = JCOL - J 00239 DO 90 K = J, I - 1 00240 T = T + A( KCOL )*AINV( JCOL+K ) 00241 KCOL = KCOL + N - K 00242 90 CONTINUE 00243 WORK( I, J ) = -T 00244 100 CONTINUE 00245 * 00246 * Code when J > I 00247 * 00248 ICOL = NALL - ( ( N-I )*( N-I+1 ) ) / 2 00249 DO 130 J = I + 1, N 00250 JCOL = NALL - ( ( N-J+1 )*( N-J+2 ) ) / 2 + 1 00251 T = ZDOTU( N-J+1, A( ICOL-N+J ), 1, AINV( JCOL ), 1 ) 00252 KCOL = I 00253 JCOL = J 00254 DO 110 K = 1, I - 1 00255 T = T + A( KCOL )*AINV( JCOL ) 00256 JCOL = JCOL + N - K 00257 KCOL = KCOL + N - K 00258 110 CONTINUE 00259 KCOL = KCOL - I 00260 DO 120 K = I, J - 1 00261 T = T + A( KCOL+K )*AINV( JCOL ) 00262 JCOL = JCOL + N - K 00263 120 CONTINUE 00264 WORK( I, J ) = -T 00265 130 CONTINUE 00266 140 CONTINUE 00267 END IF 00268 * 00269 * Add the identity matrix to WORK . 00270 * 00271 DO 150 I = 1, N 00272 WORK( I, I ) = WORK( I, I ) + ONE 00273 150 CONTINUE 00274 * 00275 * Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS) 00276 * 00277 RESID = ZLANGE( '1', N, N, WORK, LDW, RWORK ) 00278 * 00279 RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N ) 00280 * 00281 RETURN 00282 * 00283 * End of ZSPT03 00284 * 00285 END