00001 /* shst01.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_b7 = 1.f; 00019 static real c_b8 = 0.f; 00020 static real c_b11 = -1.f; 00021 00022 /* Subroutine */ int shst01_(integer *n, integer *ilo, integer *ihi, real *a, 00023 integer *lda, real *h__, integer *ldh, real *q, integer *ldq, real * 00024 work, integer *lwork, real *result) 00025 { 00026 /* System generated locals */ 00027 integer a_dim1, a_offset, h_dim1, h_offset, q_dim1, q_offset; 00028 real r__1, r__2; 00029 00030 /* Local variables */ 00031 real eps, unfl, ovfl; 00032 extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *, 00033 integer *, real *, real *, integer *, real *, integer *, real *, 00034 real *, integer *); 00035 real anorm; 00036 extern /* Subroutine */ int sort01_(char *, integer *, integer *, real *, 00037 integer *, real *, integer *, real *); 00038 real wnorm; 00039 extern /* Subroutine */ int slabad_(real *, real *); 00040 extern doublereal slamch_(char *), slange_(char *, integer *, 00041 integer *, real *, integer *, real *); 00042 extern /* Subroutine */ int slacpy_(char *, integer *, integer *, real *, 00043 integer *, real *, integer *); 00044 integer ldwork; 00045 real smlnum; 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 /* SHST01 tests the reduction of a general matrix A to upper Hessenberg */ 00061 /* form: A = Q*H*Q'. Two test ratios are computed; */ 00062 00063 /* RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS ) */ 00064 /* RESULT(2) = norm( I - Q'*Q ) / ( N * EPS ) */ 00065 00066 /* The matrix Q is assumed to be given explicitly as it would be */ 00067 /* following SGEHRD + SORGHR. */ 00068 00069 /* In this version, ILO and IHI are not used and are assumed to be 1 and */ 00070 /* N, respectively. */ 00071 00072 /* Arguments */ 00073 /* ========= */ 00074 00075 /* N (input) INTEGER */ 00076 /* The order of the matrix A. N >= 0. */ 00077 00078 /* ILO (input) INTEGER */ 00079 /* IHI (input) INTEGER */ 00080 /* A is assumed to be upper triangular in rows and columns */ 00081 /* 1:ILO-1 and IHI+1:N, so Q differs from the identity only in */ 00082 /* rows and columns ILO+1:IHI. */ 00083 00084 /* A (input) REAL array, dimension (LDA,N) */ 00085 /* The original n by n matrix A. */ 00086 00087 /* LDA (input) INTEGER */ 00088 /* The leading dimension of the array A. LDA >= max(1,N). */ 00089 00090 /* H (input) REAL array, dimension (LDH,N) */ 00091 /* The upper Hessenberg matrix H from the reduction A = Q*H*Q' */ 00092 /* as computed by SGEHRD. H is assumed to be zero below the */ 00093 /* first subdiagonal. */ 00094 00095 /* LDH (input) INTEGER */ 00096 /* The leading dimension of the array H. LDH >= max(1,N). */ 00097 00098 /* Q (input) REAL array, dimension (LDQ,N) */ 00099 /* The orthogonal matrix Q from the reduction A = Q*H*Q' as */ 00100 /* computed by SGEHRD + SORGHR. */ 00101 00102 /* LDQ (input) INTEGER */ 00103 /* The leading dimension of the array Q. LDQ >= max(1,N). */ 00104 00105 /* WORK (workspace) REAL array, dimension (LWORK) */ 00106 00107 /* LWORK (input) INTEGER */ 00108 /* The length of the array WORK. LWORK >= 2*N*N. */ 00109 00110 /* RESULT (output) REAL array, dimension (2) */ 00111 /* RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS ) */ 00112 /* RESULT(2) = norm( I - Q'*Q ) / ( N * EPS ) */ 00113 00114 /* ===================================================================== */ 00115 00116 /* .. Parameters .. */ 00117 /* .. */ 00118 /* .. Local Scalars .. */ 00119 /* .. */ 00120 /* .. External Functions .. */ 00121 /* .. */ 00122 /* .. External Subroutines .. */ 00123 /* .. */ 00124 /* .. Intrinsic Functions .. */ 00125 /* .. */ 00126 /* .. Executable Statements .. */ 00127 00128 /* Quick return if possible */ 00129 00130 /* Parameter adjustments */ 00131 a_dim1 = *lda; 00132 a_offset = 1 + a_dim1; 00133 a -= a_offset; 00134 h_dim1 = *ldh; 00135 h_offset = 1 + h_dim1; 00136 h__ -= h_offset; 00137 q_dim1 = *ldq; 00138 q_offset = 1 + q_dim1; 00139 q -= q_offset; 00140 --work; 00141 --result; 00142 00143 /* Function Body */ 00144 if (*n <= 0) { 00145 result[1] = 0.f; 00146 result[2] = 0.f; 00147 return 0; 00148 } 00149 00150 unfl = slamch_("Safe minimum"); 00151 eps = slamch_("Precision"); 00152 ovfl = 1.f / unfl; 00153 slabad_(&unfl, &ovfl); 00154 smlnum = unfl * *n / eps; 00155 00156 /* Test 1: Compute norm( A - Q*H*Q' ) / ( norm(A) * N * EPS ) */ 00157 00158 /* Copy A to WORK */ 00159 00160 ldwork = max(1,*n); 00161 slacpy_(" ", n, n, &a[a_offset], lda, &work[1], &ldwork); 00162 00163 /* Compute Q*H */ 00164 00165 sgemm_("No transpose", "No transpose", n, n, n, &c_b7, &q[q_offset], ldq, 00166 &h__[h_offset], ldh, &c_b8, &work[ldwork * *n + 1], &ldwork); 00167 00168 /* Compute A - Q*H*Q' */ 00169 00170 sgemm_("No transpose", "Transpose", n, n, n, &c_b11, &work[ldwork * *n + 00171 1], &ldwork, &q[q_offset], ldq, &c_b7, &work[1], &ldwork); 00172 00173 /* Computing MAX */ 00174 r__1 = slange_("1", n, n, &a[a_offset], lda, &work[ldwork * *n + 1]); 00175 anorm = dmax(r__1,unfl); 00176 wnorm = slange_("1", n, n, &work[1], &ldwork, &work[ldwork * *n + 1]); 00177 00178 /* Note that RESULT(1) cannot overflow and is bounded by 1/(N*EPS) */ 00179 00180 /* Computing MAX */ 00181 r__1 = smlnum, r__2 = anorm * eps; 00182 result[1] = dmin(wnorm,anorm) / dmax(r__1,r__2) / *n; 00183 00184 /* Test 2: Compute norm( I - Q'*Q ) / ( N * EPS ) */ 00185 00186 sort01_("Columns", n, n, &q[q_offset], ldq, &work[1], lwork, &result[2]); 00187 00188 return 0; 00189 00190 /* End of SHST01 */ 00191 00192 } /* shst01_ */