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