 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZSPT03( UPLO, N, A, AINV, WORK, LDW, RWORK, RCOND, 00002 $ RESID ) 00003 * 00004 * -- LAPACK test routine (version 3.1) -- 00005 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER UPLO 00010 INTEGER LDW, N 00011 DOUBLE PRECISION RCOND, RESID 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION RWORK( * ) 00015 COMPLEX*16 A( * ), AINV( * ), WORK( LDW, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * ZSPT03 computes the residual for a complex symmetric packed matrix 00022 * times its inverse: 00023 * norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ), 00024 * where EPS is the machine epsilon. 00025 * 00026 * Arguments 00027 * ========== 00028 * 00029 * UPLO (input) CHARACTER*1 00030 * Specifies whether the upper or lower triangular part of the 00031 * complex symmetric matrix A is stored: 00032 * = 'U': Upper triangular 00033 * = 'L': Lower triangular 00034 * 00035 * N (input) INTEGER 00036 * The number of rows and columns of the matrix A. N >= 0. 00037 * 00038 * A (input) COMPLEX*16 array, dimension (N*(N+1)/2) 00039 * The original complex symmetric matrix A, stored as a packed 00040 * triangular matrix. 00041 * 00042 * AINV (input) COMPLEX*16 array, dimension (N*(N+1)/2) 00043 * The (symmetric) inverse of the matrix A, stored as a packed 00044 * triangular matrix. 00045 * 00046 * WORK (workspace) COMPLEX*16 array, dimension (LDWORK,N) 00047 * 00048 * LDWORK (input) INTEGER 00049 * The leading dimension of the array WORK. LDWORK >= max(1,N). 00050 * 00051 * RWORK (workspace) DOUBLE PRECISION array, dimension (N) 00052 * 00053 * RCOND (output) DOUBLE PRECISION 00054 * The reciprocal of the condition number of A, computed as 00055 * ( 1/norm(A) ) / norm(AINV). 00056 * 00057 * RESID (output) DOUBLE PRECISION 00058 * norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS ) 00059 * 00060 * ===================================================================== 00061 * 00062 * .. Parameters .. 00063 DOUBLE PRECISION ZERO, ONE 00064 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00065 * .. 00066 * .. Local Scalars .. 00067 INTEGER I, ICOL, J, JCOL, K, KCOL, NALL 00068 DOUBLE PRECISION AINVNM, ANORM, EPS 00069 COMPLEX*16 T 00070 * .. 00071 * .. External Functions .. 00072 LOGICAL LSAME 00073 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANSP 00074 COMPLEX*16 ZDOTU 00075 EXTERNAL LSAME, DLAMCH, ZLANGE, ZLANSP, ZDOTU 00076 * .. 00077 * .. Intrinsic Functions .. 00078 INTRINSIC DBLE 00079 * .. 00080 * .. Executable Statements .. 00081 * 00082 * Quick exit if N = 0. 00083 * 00084 IF( N.LE.0 ) THEN 00085 RCOND = ONE 00086 RESID = ZERO 00087 RETURN 00088 END IF 00089 * 00090 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 00091 * 00092 EPS = DLAMCH( 'Epsilon' ) 00093 ANORM = ZLANSP( '1', UPLO, N, A, RWORK ) 00094 AINVNM = ZLANSP( '1', UPLO, N, AINV, RWORK ) 00095 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00096 RCOND = ZERO 00097 RESID = ONE / EPS 00098 RETURN 00099 END IF 00100 RCOND = ( ONE / ANORM ) / AINVNM 00101 * 00102 * Case where both A and AINV are upper triangular: 00103 * Each element of - A * AINV is computed by taking the dot product 00104 * of a row of A with a column of AINV. 00105 * 00106 IF( LSAME( UPLO, 'U' ) ) THEN 00107 DO 70 I = 1, N 00108 ICOL = ( ( I-1 )*I ) / 2 + 1 00109 * 00110 * Code when J <= I 00111 * 00112 DO 30 J = 1, I 00113 JCOL = ( ( J-1 )*J ) / 2 + 1 00114 T = ZDOTU( J, A( ICOL ), 1, AINV( JCOL ), 1 ) 00115 JCOL = JCOL + 2*J - 1 00116 KCOL = ICOL - 1 00117 DO 10 K = J + 1, I 00118 T = T + A( KCOL+K )*AINV( JCOL ) 00119 JCOL = JCOL + K 00120 10 CONTINUE 00121 KCOL = KCOL + 2*I 00122 DO 20 K = I + 1, N 00123 T = T + A( KCOL )*AINV( JCOL ) 00124 KCOL = KCOL + K 00125 JCOL = JCOL + K 00126 20 CONTINUE 00127 WORK( I, J ) = -T 00128 30 CONTINUE 00129 * 00130 * Code when J > I 00131 * 00132 DO 60 J = I + 1, N 00133 JCOL = ( ( J-1 )*J ) / 2 + 1 00134 T = ZDOTU( I, A( ICOL ), 1, AINV( JCOL ), 1 ) 00135 JCOL = JCOL - 1 00136 KCOL = ICOL + 2*I - 1 00137 DO 40 K = I + 1, J 00138 T = T + A( KCOL )*AINV( JCOL+K ) 00139 KCOL = KCOL + K 00140 40 CONTINUE 00141 JCOL = JCOL + 2*J 00142 DO 50 K = J + 1, N 00143 T = T + A( KCOL )*AINV( JCOL ) 00144 KCOL = KCOL + K 00145 JCOL = JCOL + K 00146 50 CONTINUE 00147 WORK( I, J ) = -T 00148 60 CONTINUE 00149 70 CONTINUE 00150 ELSE 00151 * 00152 * Case where both A and AINV are lower triangular 00153 * 00154 NALL = ( N*( N+1 ) ) / 2 00155 DO 140 I = 1, N 00156 * 00157 * Code when J <= I 00158 * 00159 ICOL = NALL - ( ( N-I+1 )*( N-I+2 ) ) / 2 + 1 00160 DO 100 J = 1, I 00161 JCOL = NALL - ( ( N-J )*( N-J+1 ) ) / 2 - ( N-I ) 00162 T = ZDOTU( N-I+1, A( ICOL ), 1, AINV( JCOL ), 1 ) 00163 KCOL = I 00164 JCOL = J 00165 DO 80 K = 1, J - 1 00166 T = T + A( KCOL )*AINV( JCOL ) 00167 JCOL = JCOL + N - K 00168 KCOL = KCOL + N - K 00169 80 CONTINUE 00170 JCOL = JCOL - J 00171 DO 90 K = J, I - 1 00172 T = T + A( KCOL )*AINV( JCOL+K ) 00173 KCOL = KCOL + N - K 00174 90 CONTINUE 00175 WORK( I, J ) = -T 00176 100 CONTINUE 00177 * 00178 * Code when J > I 00179 * 00180 ICOL = NALL - ( ( N-I )*( N-I+1 ) ) / 2 00181 DO 130 J = I + 1, N 00182 JCOL = NALL - ( ( N-J+1 )*( N-J+2 ) ) / 2 + 1 00183 T = ZDOTU( N-J+1, A( ICOL-N+J ), 1, AINV( JCOL ), 1 ) 00184 KCOL = I 00185 JCOL = J 00186 DO 110 K = 1, I - 1 00187 T = T + A( KCOL )*AINV( JCOL ) 00188 JCOL = JCOL + N - K 00189 KCOL = KCOL + N - K 00190 110 CONTINUE 00191 KCOL = KCOL - I 00192 DO 120 K = I, J - 1 00193 T = T + A( KCOL+K )*AINV( JCOL ) 00194 JCOL = JCOL + N - K 00195 120 CONTINUE 00196 WORK( I, J ) = -T 00197 130 CONTINUE 00198 140 CONTINUE 00199 END IF 00200 * 00201 * Add the identity matrix to WORK . 00202 * 00203 DO 150 I = 1, N 00204 WORK( I, I ) = WORK( I, I ) + ONE 00205 150 CONTINUE 00206 * 00207 * Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS) 00208 * 00209 RESID = ZLANGE( '1', N, N, WORK, LDW, RWORK ) 00210 * 00211 RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N ) 00212 * 00213 RETURN 00214 * 00215 * End of ZSPT03 00216 * 00217 END
1.7.3 
