sqlt01.c
Go to the documentation of this file.
00001 /* sqlt01.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_b6 = -1e10f;
00027 static real c_b13 = 0.f;
00028 static real c_b20 = -1.f;
00029 static real c_b21 = 1.f;
00030 
00031 /* Subroutine */ int sqlt01_(integer *m, integer *n, real *a, real *af, real *
00032         q, real *l, integer *lda, real *tau, real *work, integer *lwork, real 
00033         *rwork, real *result)
00034 {
00035     /* System generated locals */
00036     integer a_dim1, a_offset, af_dim1, af_offset, l_dim1, l_offset, q_dim1, 
00037             q_offset, i__1, i__2;
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     integer minmn;
00051     extern /* Subroutine */ int ssyrk_(char *, char *, integer *, integer *, 
00052             real *, real *, integer *, real *, real *, integer *);
00053     extern doublereal slamch_(char *), slange_(char *, integer *, 
00054             integer *, real *, integer *, real *);
00055     extern /* Subroutine */ int sgeqlf_(integer *, integer *, real *, integer 
00056             *, real *, real *, integer *, integer *), slacpy_(char *, integer 
00057             *, integer *, real *, integer *, real *, integer *), 
00058             slaset_(char *, integer *, integer *, real *, real *, real *, 
00059             integer *), sorgql_(integer *, integer *, integer *, real 
00060             *, integer *, real *, real *, integer *, integer *);
00061     extern doublereal slansy_(char *, char *, integer *, real *, integer *, 
00062             real *);
00063 
00064 
00065 /*  -- LAPACK test routine (version 3.1) -- */
00066 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00067 /*     November 2006 */
00068 
00069 /*     .. Scalar Arguments .. */
00070 /*     .. */
00071 /*     .. Array Arguments .. */
00072 /*     .. */
00073 
00074 /*  Purpose */
00075 /*  ======= */
00076 
00077 /*  SQLT01 tests SGEQLF, which computes the QL factorization of an m-by-n */
00078 /*  matrix A, and partially tests SORGQL which forms the m-by-m */
00079 /*  orthogonal matrix Q. */
00080 
00081 /*  SQLT01 compares L with Q'*A, and checks that Q is orthogonal. */
00082 
00083 /*  Arguments */
00084 /*  ========= */
00085 
00086 /*  M       (input) INTEGER */
00087 /*          The number of rows of the matrix A.  M >= 0. */
00088 
00089 /*  N       (input) INTEGER */
00090 /*          The number of columns of the matrix A.  N >= 0. */
00091 
00092 /*  A       (input) REAL array, dimension (LDA,N) */
00093 /*          The m-by-n matrix A. */
00094 
00095 /*  AF      (output) REAL array, dimension (LDA,N) */
00096 /*          Details of the QL factorization of A, as returned by SGEQLF. */
00097 /*          See SGEQLF for further details. */
00098 
00099 /*  Q       (output) REAL array, dimension (LDA,M) */
00100 /*          The m-by-m orthogonal matrix Q. */
00101 
00102 /*  L       (workspace) REAL array, dimension (LDA,max(M,N)) */
00103 
00104 /*  LDA     (input) INTEGER */
00105 /*          The leading dimension of the arrays A, AF, Q and R. */
00106 /*          LDA >= max(M,N). */
00107 
00108 /*  TAU     (output) REAL array, dimension (min(M,N)) */
00109 /*          The scalar factors of the elementary reflectors, as returned */
00110 /*          by SGEQLF. */
00111 
00112 /*  WORK    (workspace) REAL array, dimension (LWORK) */
00113 
00114 /*  LWORK   (input) INTEGER */
00115 /*          The dimension of the array WORK. */
00116 
00117 /*  RWORK   (workspace) REAL array, dimension (M) */
00118 
00119 /*  RESULT  (output) REAL array, dimension (2) */
00120 /*          The test ratios: */
00121 /*          RESULT(1) = norm( L - Q'*A ) / ( M * norm(A) * EPS ) */
00122 /*          RESULT(2) = norm( I - Q'*Q ) / ( M * EPS ) */
00123 
00124 /*  ===================================================================== */
00125 
00126 /*     .. Parameters .. */
00127 /*     .. */
00128 /*     .. Local Scalars .. */
00129 /*     .. */
00130 /*     .. External Functions .. */
00131 /*     .. */
00132 /*     .. External Subroutines .. */
00133 /*     .. */
00134 /*     .. Intrinsic Functions .. */
00135 /*     .. */
00136 /*     .. Scalars in Common .. */
00137 /*     .. */
00138 /*     .. Common blocks .. */
00139 /*     .. */
00140 /*     .. Executable Statements .. */
00141 
00142     /* Parameter adjustments */
00143     l_dim1 = *lda;
00144     l_offset = 1 + l_dim1;
00145     l -= l_offset;
00146     q_dim1 = *lda;
00147     q_offset = 1 + q_dim1;
00148     q -= q_offset;
00149     af_dim1 = *lda;
00150     af_offset = 1 + af_dim1;
00151     af -= af_offset;
00152     a_dim1 = *lda;
00153     a_offset = 1 + a_dim1;
00154     a -= a_offset;
00155     --tau;
00156     --work;
00157     --rwork;
00158     --result;
00159 
00160     /* Function Body */
00161     minmn = min(*m,*n);
00162     eps = slamch_("Epsilon");
00163 
00164 /*     Copy the matrix A to the array AF. */
00165 
00166     slacpy_("Full", m, n, &a[a_offset], lda, &af[af_offset], lda);
00167 
00168 /*     Factorize the matrix A in the array AF. */
00169 
00170     s_copy(srnamc_1.srnamt, "SGEQLF", (ftnlen)32, (ftnlen)6);
00171     sgeqlf_(m, n, &af[af_offset], lda, &tau[1], &work[1], lwork, &info);
00172 
00173 /*     Copy details of Q */
00174 
00175     slaset_("Full", m, m, &c_b6, &c_b6, &q[q_offset], lda);
00176     if (*m >= *n) {
00177         if (*n < *m && *n > 0) {
00178             i__1 = *m - *n;
00179             slacpy_("Full", &i__1, n, &af[af_offset], lda, &q[(*m - *n + 1) * 
00180                     q_dim1 + 1], lda);
00181         }
00182         if (*n > 1) {
00183             i__1 = *n - 1;
00184             i__2 = *n - 1;
00185             slacpy_("Upper", &i__1, &i__2, &af[*m - *n + 1 + (af_dim1 << 1)], 
00186                     lda, &q[*m - *n + 1 + (*m - *n + 2) * q_dim1], lda);
00187         }
00188     } else {
00189         if (*m > 1) {
00190             i__1 = *m - 1;
00191             i__2 = *m - 1;
00192             slacpy_("Upper", &i__1, &i__2, &af[(*n - *m + 2) * af_dim1 + 1], 
00193                     lda, &q[(q_dim1 << 1) + 1], lda);
00194         }
00195     }
00196 
00197 /*     Generate the m-by-m matrix Q */
00198 
00199     s_copy(srnamc_1.srnamt, "SORGQL", (ftnlen)32, (ftnlen)6);
00200     sorgql_(m, m, &minmn, &q[q_offset], lda, &tau[1], &work[1], lwork, &info);
00201 
00202 /*     Copy L */
00203 
00204     slaset_("Full", m, n, &c_b13, &c_b13, &l[l_offset], lda);
00205     if (*m >= *n) {
00206         if (*n > 0) {
00207             slacpy_("Lower", n, n, &af[*m - *n + 1 + af_dim1], lda, &l[*m - *
00208                     n + 1 + l_dim1], lda);
00209         }
00210     } else {
00211         if (*n > *m && *m > 0) {
00212             i__1 = *n - *m;
00213             slacpy_("Full", m, &i__1, &af[af_offset], lda, &l[l_offset], lda);
00214         }
00215         if (*m > 0) {
00216             slacpy_("Lower", m, m, &af[(*n - *m + 1) * af_dim1 + 1], lda, &l[(
00217                     *n - *m + 1) * l_dim1 + 1], lda);
00218         }
00219     }
00220 
00221 /*     Compute L - Q'*A */
00222 
00223     sgemm_("Transpose", "No transpose", m, n, m, &c_b20, &q[q_offset], lda, &
00224             a[a_offset], lda, &c_b21, &l[l_offset], lda);
00225 
00226 /*     Compute norm( L - Q'*A ) / ( M * norm(A) * EPS ) . */
00227 
00228     anorm = slange_("1", m, n, &a[a_offset], lda, &rwork[1]);
00229     resid = slange_("1", m, n, &l[l_offset], lda, &rwork[1]);
00230     if (anorm > 0.f) {
00231         result[1] = resid / (real) max(1,*m) / anorm / eps;
00232     } else {
00233         result[1] = 0.f;
00234     }
00235 
00236 /*     Compute I - Q'*Q */
00237 
00238     slaset_("Full", m, m, &c_b13, &c_b21, &l[l_offset], lda);
00239     ssyrk_("Upper", "Transpose", m, m, &c_b20, &q[q_offset], lda, &c_b21, &l[
00240             l_offset], lda);
00241 
00242 /*     Compute norm( I - Q'*Q ) / ( M * EPS ) . */
00243 
00244     resid = slansy_("1", "Upper", m, &l[l_offset], lda, &rwork[1]);
00245 
00246     result[2] = resid / (real) max(1,*m) / eps;
00247 
00248     return 0;
00249 
00250 /*     End of SQLT01 */
00251 
00252 } /* sqlt01_ */


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