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