shst01.c
Go to the documentation of this file.
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_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:09