sqrt02.c
Go to the documentation of this file.
00001 /* sqrt02.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 /* Common Block Declarations */
00017 
00018 struct {
00019     char srnamt[32];
00020 } srnamc_;
00021 
00022 #define srnamc_1 srnamc_
00023 
00024 /* Table of constant values */
00025 
00026 static real c_b4 = -1e10f;
00027 static real c_b9 = 0.f;
00028 static real c_b14 = -1.f;
00029 static real c_b15 = 1.f;
00030 
00031 /* Subroutine */ int sqrt02_(integer *m, integer *n, integer *k, real *a, 
00032         real *af, real *q, real *r__, integer *lda, real *tau, real *work, 
00033         integer *lwork, real *rwork, real *result)
00034 {
00035     /* System generated locals */
00036     integer a_dim1, a_offset, af_dim1, af_offset, q_dim1, q_offset, r_dim1, 
00037             r_offset, i__1;
00038 
00039     /* Builtin functions */
00040     /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
00041 
00042     /* Local variables */
00043     real eps;
00044     integer info;
00045     real resid;
00046     extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *, 
00047             integer *, real *, real *, integer *, real *, integer *, real *, 
00048             real *, integer *);
00049     real anorm;
00050     extern /* Subroutine */ int ssyrk_(char *, char *, integer *, integer *, 
00051             real *, real *, integer *, real *, real *, integer *);
00052     extern doublereal slamch_(char *), slange_(char *, integer *, 
00053             integer *, real *, integer *, real *);
00054     extern /* Subroutine */ int slacpy_(char *, integer *, integer *, real *, 
00055             integer *, real *, integer *), slaset_(char *, integer *, 
00056             integer *, real *, real *, real *, integer *);
00057     extern doublereal slansy_(char *, char *, integer *, real *, integer *, 
00058             real *);
00059     extern /* Subroutine */ int sorgqr_(integer *, integer *, integer *, real 
00060             *, integer *, real *, real *, integer *, integer *);
00061 
00062 
00063 /*  -- LAPACK test routine (version 3.1) -- */
00064 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00065 /*     November 2006 */
00066 
00067 /*     .. Scalar Arguments .. */
00068 /*     .. */
00069 /*     .. Array Arguments .. */
00070 /*     .. */
00071 
00072 /*  Purpose */
00073 /*  ======= */
00074 
00075 /*  SQRT02 tests SORGQR, which generates an m-by-n matrix Q with */
00076 /*  orthonornmal columns that is defined as the product of k elementary */
00077 /*  reflectors. */
00078 
00079 /*  Given the QR factorization of an m-by-n matrix A, SQRT02 generates */
00080 /*  the orthogonal matrix Q defined by the factorization of the first k */
00081 /*  columns of A; it compares R(1:n,1:k) with Q(1:m,1:n)'*A(1:m,1:k), */
00082 /*  and checks that the columns of Q are orthonormal. */
00083 
00084 /*  Arguments */
00085 /*  ========= */
00086 
00087 /*  M       (input) INTEGER */
00088 /*          The number of rows of the matrix Q to be generated.  M >= 0. */
00089 
00090 /*  N       (input) INTEGER */
00091 /*          The number of columns of the matrix Q to be generated. */
00092 /*          M >= N >= 0. */
00093 
00094 /*  K       (input) INTEGER */
00095 /*          The number of elementary reflectors whose product defines the */
00096 /*          matrix Q. N >= K >= 0. */
00097 
00098 /*  A       (input) REAL array, dimension (LDA,N) */
00099 /*          The m-by-n matrix A which was factorized by SQRT01. */
00100 
00101 /*  AF      (input) REAL array, dimension (LDA,N) */
00102 /*          Details of the QR factorization of A, as returned by SGEQRF. */
00103 /*          See SGEQRF for further details. */
00104 
00105 /*  Q       (workspace) REAL array, dimension (LDA,N) */
00106 
00107 /*  R       (workspace) REAL array, dimension (LDA,N) */
00108 
00109 /*  LDA     (input) INTEGER */
00110 /*          The leading dimension of the arrays A, AF, Q and R. LDA >= M. */
00111 
00112 /*  TAU     (input) REAL array, dimension (N) */
00113 /*          The scalar factors of the elementary reflectors corresponding */
00114 /*          to the QR factorization in AF. */
00115 
00116 /*  WORK    (workspace) REAL array, dimension (LWORK) */
00117 
00118 /*  LWORK   (input) INTEGER */
00119 /*          The dimension of the array WORK. */
00120 
00121 /*  RWORK   (workspace) REAL array, dimension (M) */
00122 
00123 /*  RESULT  (output) REAL array, dimension (2) */
00124 /*          The test ratios: */
00125 /*          RESULT(1) = norm( R - Q'*A ) / ( M * norm(A) * EPS ) */
00126 /*          RESULT(2) = norm( I - Q'*Q ) / ( M * EPS ) */
00127 
00128 /*  ===================================================================== */
00129 
00130 /*     .. Parameters .. */
00131 /*     .. */
00132 /*     .. Local Scalars .. */
00133 /*     .. */
00134 /*     .. External Functions .. */
00135 /*     .. */
00136 /*     .. External Subroutines .. */
00137 /*     .. */
00138 /*     .. Intrinsic Functions .. */
00139 /*     .. */
00140 /*     .. Scalars in Common .. */
00141 /*     .. */
00142 /*     .. Common blocks .. */
00143 /*     .. */
00144 /*     .. Executable Statements .. */
00145 
00146     /* Parameter adjustments */
00147     r_dim1 = *lda;
00148     r_offset = 1 + r_dim1;
00149     r__ -= r_offset;
00150     q_dim1 = *lda;
00151     q_offset = 1 + q_dim1;
00152     q -= q_offset;
00153     af_dim1 = *lda;
00154     af_offset = 1 + af_dim1;
00155     af -= af_offset;
00156     a_dim1 = *lda;
00157     a_offset = 1 + a_dim1;
00158     a -= a_offset;
00159     --tau;
00160     --work;
00161     --rwork;
00162     --result;
00163 
00164     /* Function Body */
00165     eps = slamch_("Epsilon");
00166 
00167 /*     Copy the first k columns of the factorization to the array Q */
00168 
00169     slaset_("Full", m, n, &c_b4, &c_b4, &q[q_offset], lda);
00170     i__1 = *m - 1;
00171     slacpy_("Lower", &i__1, k, &af[af_dim1 + 2], lda, &q[q_dim1 + 2], lda);
00172 
00173 /*     Generate the first n columns of the matrix Q */
00174 
00175     s_copy(srnamc_1.srnamt, "SORGQR", (ftnlen)32, (ftnlen)6);
00176     sorgqr_(m, n, k, &q[q_offset], lda, &tau[1], &work[1], lwork, &info);
00177 
00178 /*     Copy R(1:n,1:k) */
00179 
00180     slaset_("Full", n, k, &c_b9, &c_b9, &r__[r_offset], lda);
00181     slacpy_("Upper", n, k, &af[af_offset], lda, &r__[r_offset], lda);
00182 
00183 /*     Compute R(1:n,1:k) - Q(1:m,1:n)' * A(1:m,1:k) */
00184 
00185     sgemm_("Transpose", "No transpose", n, k, m, &c_b14, &q[q_offset], lda, &
00186             a[a_offset], lda, &c_b15, &r__[r_offset], lda);
00187 
00188 /*     Compute norm( R - Q'*A ) / ( M * norm(A) * EPS ) . */
00189 
00190     anorm = slange_("1", m, k, &a[a_offset], lda, &rwork[1]);
00191     resid = slange_("1", n, k, &r__[r_offset], lda, &rwork[1]);
00192     if (anorm > 0.f) {
00193         result[1] = resid / (real) max(1,*m) / anorm / eps;
00194     } else {
00195         result[1] = 0.f;
00196     }
00197 
00198 /*     Compute I - Q'*Q */
00199 
00200     slaset_("Full", n, n, &c_b9, &c_b15, &r__[r_offset], lda);
00201     ssyrk_("Upper", "Transpose", n, m, &c_b14, &q[q_offset], lda, &c_b15, &
00202             r__[r_offset], lda);
00203 
00204 /*     Compute norm( I - Q'*Q ) / ( M * EPS ) . */
00205 
00206     resid = slansy_("1", "Upper", n, &r__[r_offset], lda, &rwork[1]);
00207 
00208     result[2] = resid / (real) max(1,*m) / eps;
00209 
00210     return 0;
00211 
00212 /*     End of SQRT02 */
00213 
00214 } /* sqrt02_ */


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