00001 /* zhet22.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static doublecomplex c_b1 = {0.,0.}; 00019 static doublecomplex c_b2 = {1.,0.}; 00020 00021 /* Subroutine */ int zhet22_(integer *itype, char *uplo, integer *n, integer * 00022 m, integer *kband, doublecomplex *a, integer *lda, doublereal *d__, 00023 doublereal *e, doublecomplex *u, integer *ldu, doublecomplex *v, 00024 integer *ldv, doublecomplex *tau, doublecomplex *work, doublereal * 00025 rwork, doublereal *result) 00026 { 00027 /* System generated locals */ 00028 integer a_dim1, a_offset, u_dim1, u_offset, v_dim1, v_offset, i__1, i__2, 00029 i__3, i__4; 00030 doublereal d__1, d__2; 00031 doublecomplex z__1; 00032 00033 /* Local variables */ 00034 integer j, jj, nn, jj1, jj2; 00035 doublereal ulp; 00036 integer nnp1; 00037 doublereal unfl, anorm; 00038 extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 00039 integer *, doublecomplex *, doublecomplex *, integer *, 00040 doublecomplex *, integer *, doublecomplex *, doublecomplex *, 00041 integer *), zhemm_(char *, char *, integer *, 00042 integer *, doublecomplex *, doublecomplex *, integer *, 00043 doublecomplex *, integer *, doublecomplex *, doublecomplex *, 00044 integer *), zunt01_(char *, integer *, integer *, 00045 doublecomplex *, integer *, doublecomplex *, integer *, 00046 doublereal *, doublereal *); 00047 doublereal wnorm; 00048 extern doublereal dlamch_(char *), zlanhe_(char *, char *, 00049 integer *, doublecomplex *, integer *, doublereal *); 00050 00051 00052 /* -- LAPACK test routine (version 3.1) -- */ 00053 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00054 /* November 2006 */ 00055 00056 /* .. Scalar Arguments .. */ 00057 /* .. */ 00058 /* .. Array Arguments .. */ 00059 /* .. */ 00060 00061 /* Purpose */ 00062 /* ======= */ 00063 00064 /* ZHET22 generally checks a decomposition of the form */ 00065 00066 /* A U = U S */ 00067 00068 /* where A is complex Hermitian, the columns of U are orthonormal, */ 00069 /* and S is diagonal (if KBAND=0) or symmetric tridiagonal (if */ 00070 /* KBAND=1). If ITYPE=1, then U is represented as a dense matrix, */ 00071 /* otherwise the U is expressed as a product of Householder */ 00072 /* transformations, whose vectors are stored in the array "V" and */ 00073 /* whose scaling constants are in "TAU"; we shall use the letter */ 00074 /* "V" to refer to the product of Householder transformations */ 00075 /* (which should be equal to U). */ 00076 00077 /* Specifically, if ITYPE=1, then: */ 00078 00079 /* RESULT(1) = | U' A U - S | / ( |A| m ulp ) *and* */ 00080 /* RESULT(2) = | I - U'U | / ( m ulp ) */ 00081 00082 /* Arguments */ 00083 /* ========= */ 00084 00085 /* ITYPE INTEGER */ 00086 /* Specifies the type of tests to be performed. */ 00087 /* 1: U expressed as a dense orthogonal matrix: */ 00088 /* RESULT(1) = | A - U S U' | / ( |A| n ulp ) *and* */ 00089 /* RESULT(2) = | I - UU' | / ( n ulp ) */ 00090 00091 /* UPLO CHARACTER */ 00092 /* If UPLO='U', the upper triangle of A will be used and the */ 00093 /* (strictly) lower triangle will not be referenced. If */ 00094 /* UPLO='L', the lower triangle of A will be used and the */ 00095 /* (strictly) upper triangle will not be referenced. */ 00096 /* Not modified. */ 00097 00098 /* N INTEGER */ 00099 /* The size of the matrix. If it is zero, ZHET22 does nothing. */ 00100 /* It must be at least zero. */ 00101 /* Not modified. */ 00102 00103 /* M INTEGER */ 00104 /* The number of columns of U. If it is zero, ZHET22 does */ 00105 /* nothing. It must be at least zero. */ 00106 /* Not modified. */ 00107 00108 /* KBAND INTEGER */ 00109 /* The bandwidth of the matrix. It may only be zero or one. */ 00110 /* If zero, then S is diagonal, and E is not referenced. If */ 00111 /* one, then S is symmetric tri-diagonal. */ 00112 /* Not modified. */ 00113 00114 /* A COMPLEX*16 array, dimension (LDA , N) */ 00115 /* The original (unfactored) matrix. It is assumed to be */ 00116 /* symmetric, and only the upper (UPLO='U') or only the lower */ 00117 /* (UPLO='L') will be referenced. */ 00118 /* Not modified. */ 00119 00120 /* LDA INTEGER */ 00121 /* The leading dimension of A. It must be at least 1 */ 00122 /* and at least N. */ 00123 /* Not modified. */ 00124 00125 /* D DOUBLE PRECISION array, dimension (N) */ 00126 /* The diagonal of the (symmetric tri-) diagonal matrix. */ 00127 /* Not modified. */ 00128 00129 /* E DOUBLE PRECISION array, dimension (N) */ 00130 /* The off-diagonal of the (symmetric tri-) diagonal matrix. */ 00131 /* E(1) is ignored, E(2) is the (1,2) and (2,1) element, etc. */ 00132 /* Not referenced if KBAND=0. */ 00133 /* Not modified. */ 00134 00135 /* U COMPLEX*16 array, dimension (LDU, N) */ 00136 /* If ITYPE=1, this contains the orthogonal matrix in */ 00137 /* the decomposition, expressed as a dense matrix. */ 00138 /* Not modified. */ 00139 00140 /* LDU INTEGER */ 00141 /* The leading dimension of U. LDU must be at least N and */ 00142 /* at least 1. */ 00143 /* Not modified. */ 00144 00145 /* V COMPLEX*16 array, dimension (LDV, N) */ 00146 /* If ITYPE=2 or 3, the lower triangle of this array contains */ 00147 /* the Householder vectors used to describe the orthogonal */ 00148 /* matrix in the decomposition. If ITYPE=1, then it is not */ 00149 /* referenced. */ 00150 /* Not modified. */ 00151 00152 /* LDV INTEGER */ 00153 /* The leading dimension of V. LDV must be at least N and */ 00154 /* at least 1. */ 00155 /* Not modified. */ 00156 00157 /* TAU COMPLEX*16 array, dimension (N) */ 00158 /* If ITYPE >= 2, then TAU(j) is the scalar factor of */ 00159 /* v(j) v(j)' in the Householder transformation H(j) of */ 00160 /* the product U = H(1)...H(n-2) */ 00161 /* If ITYPE < 2, then TAU is not referenced. */ 00162 /* Not modified. */ 00163 00164 /* WORK COMPLEX*16 array, dimension (2*N**2) */ 00165 /* Workspace. */ 00166 /* Modified. */ 00167 00168 /* RWORK DOUBLE PRECISION array, dimension (N) */ 00169 /* Workspace. */ 00170 /* Modified. */ 00171 00172 /* RESULT DOUBLE PRECISION array, dimension (2) */ 00173 /* The values computed by the two tests described above. The */ 00174 /* values are currently limited to 1/ulp, to avoid overflow. */ 00175 /* RESULT(1) is always modified. RESULT(2) is modified only */ 00176 /* if LDU is at least N. */ 00177 /* Modified. */ 00178 00179 /* ===================================================================== */ 00180 00181 /* .. Parameters .. */ 00182 /* .. */ 00183 /* .. Local Scalars .. */ 00184 /* .. */ 00185 /* .. External Functions .. */ 00186 /* .. */ 00187 /* .. External Subroutines .. */ 00188 /* .. */ 00189 /* .. Intrinsic Functions .. */ 00190 /* .. */ 00191 /* .. Executable Statements .. */ 00192 00193 /* Parameter adjustments */ 00194 a_dim1 = *lda; 00195 a_offset = 1 + a_dim1; 00196 a -= a_offset; 00197 --d__; 00198 --e; 00199 u_dim1 = *ldu; 00200 u_offset = 1 + u_dim1; 00201 u -= u_offset; 00202 v_dim1 = *ldv; 00203 v_offset = 1 + v_dim1; 00204 v -= v_offset; 00205 --tau; 00206 --work; 00207 --rwork; 00208 --result; 00209 00210 /* Function Body */ 00211 result[1] = 0.; 00212 result[2] = 0.; 00213 if (*n <= 0 || *m <= 0) { 00214 return 0; 00215 } 00216 00217 unfl = dlamch_("Safe minimum"); 00218 ulp = dlamch_("Precision"); 00219 00220 /* Do Test 1 */ 00221 00222 /* Norm of A: */ 00223 00224 /* Computing MAX */ 00225 d__1 = zlanhe_("1", uplo, n, &a[a_offset], lda, &rwork[1]); 00226 anorm = max(d__1,unfl); 00227 00228 /* Compute error matrix: */ 00229 00230 /* ITYPE=1: error = U' A U - S */ 00231 00232 zhemm_("L", uplo, n, m, &c_b2, &a[a_offset], lda, &u[u_offset], ldu, & 00233 c_b1, &work[1], n); 00234 nn = *n * *n; 00235 nnp1 = nn + 1; 00236 zgemm_("C", "N", m, m, n, &c_b2, &u[u_offset], ldu, &work[1], n, &c_b1, & 00237 work[nnp1], n); 00238 i__1 = *m; 00239 for (j = 1; j <= i__1; ++j) { 00240 jj = nn + (j - 1) * *n + j; 00241 i__2 = jj; 00242 i__3 = jj; 00243 i__4 = j; 00244 z__1.r = work[i__3].r - d__[i__4], z__1.i = work[i__3].i; 00245 work[i__2].r = z__1.r, work[i__2].i = z__1.i; 00246 /* L10: */ 00247 } 00248 if (*kband == 1 && *n > 1) { 00249 i__1 = *m; 00250 for (j = 2; j <= i__1; ++j) { 00251 jj1 = nn + (j - 1) * *n + j - 1; 00252 jj2 = nn + (j - 2) * *n + j; 00253 i__2 = jj1; 00254 i__3 = jj1; 00255 i__4 = j - 1; 00256 z__1.r = work[i__3].r - e[i__4], z__1.i = work[i__3].i; 00257 work[i__2].r = z__1.r, work[i__2].i = z__1.i; 00258 i__2 = jj2; 00259 i__3 = jj2; 00260 i__4 = j - 1; 00261 z__1.r = work[i__3].r - e[i__4], z__1.i = work[i__3].i; 00262 work[i__2].r = z__1.r, work[i__2].i = z__1.i; 00263 /* L20: */ 00264 } 00265 } 00266 wnorm = zlanhe_("1", uplo, m, &work[nnp1], n, &rwork[1]); 00267 00268 if (anorm > wnorm) { 00269 result[1] = wnorm / anorm / (*m * ulp); 00270 } else { 00271 if (anorm < 1.) { 00272 /* Computing MIN */ 00273 d__1 = wnorm, d__2 = *m * anorm; 00274 result[1] = min(d__1,d__2) / anorm / (*m * ulp); 00275 } else { 00276 /* Computing MIN */ 00277 d__1 = wnorm / anorm, d__2 = (doublereal) (*m); 00278 result[1] = min(d__1,d__2) / (*m * ulp); 00279 } 00280 } 00281 00282 /* Do Test 2 */ 00283 00284 /* Compute U'U - I */ 00285 00286 if (*itype == 1) { 00287 i__1 = (*n << 1) * *n; 00288 zunt01_("Columns", n, m, &u[u_offset], ldu, &work[1], &i__1, &rwork[1] 00289 , &result[2]); 00290 } 00291 00292 return 0; 00293 00294 /* End of ZHET22 */ 00295 00296 } /* zhet22_ */