dlqt03.c
Go to the documentation of this file.
00001 /* dlqt03.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 doublereal c_b4 = -1e10;
00027 static integer c__2 = 2;
00028 static doublereal c_b21 = -1.;
00029 static doublereal c_b22 = 1.;
00030 
00031 /* Subroutine */ int dlqt03_(integer *m, integer *n, integer *k, doublereal *
00032         af, doublereal *c__, doublereal *cc, doublereal *q, integer *lda, 
00033         doublereal *tau, doublereal *work, integer *lwork, doublereal *rwork, 
00034         doublereal *result)
00035 {
00036     /* Initialized data */
00037 
00038     static integer iseed[4] = { 1988,1989,1990,1991 };
00039 
00040     /* System generated locals */
00041     integer af_dim1, af_offset, c_dim1, c_offset, cc_dim1, cc_offset, q_dim1, 
00042             q_offset, i__1;
00043 
00044     /* Builtin functions */
00045     /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
00046 
00047     /* Local variables */
00048     integer j, mc, nc;
00049     doublereal eps;
00050     char side[1];
00051     integer info;
00052     extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *, 
00053             integer *, doublereal *, doublereal *, integer *, doublereal *, 
00054             integer *, doublereal *, doublereal *, integer *);
00055     integer iside;
00056     extern logical lsame_(char *, char *);
00057     doublereal resid, cnorm;
00058     char trans[1];
00059     extern doublereal dlamch_(char *), dlange_(char *, integer *, 
00060             integer *, doublereal *, integer *, doublereal *);
00061     extern /* Subroutine */ int dlacpy_(char *, integer *, integer *, 
00062             doublereal *, integer *, doublereal *, integer *), 
00063             dlaset_(char *, integer *, integer *, doublereal *, doublereal *, 
00064             doublereal *, integer *), dlarnv_(integer *, integer *, 
00065             integer *, doublereal *), dorglq_(integer *, integer *, integer *, 
00066              doublereal *, integer *, doublereal *, doublereal *, integer *, 
00067             integer *), dormlq_(char *, char *, integer *, integer *, integer 
00068             *, doublereal *, integer *, doublereal *, doublereal *, integer *, 
00069              doublereal *, integer *, integer *);
00070     integer itrans;
00071 
00072 
00073 /*  -- LAPACK test routine (version 3.1) -- */
00074 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00075 /*     November 2006 */
00076 
00077 /*     .. Scalar Arguments .. */
00078 /*     .. */
00079 /*     .. Array Arguments .. */
00080 /*     .. */
00081 
00082 /*  Purpose */
00083 /*  ======= */
00084 
00085 /*  DLQT03 tests DORMLQ, which computes Q*C, Q'*C, C*Q or C*Q'. */
00086 
00087 /*  DLQT03 compares the results of a call to DORMLQ with the results of */
00088 /*  forming Q explicitly by a call to DORGLQ and then performing matrix */
00089 /*  multiplication by a call to DGEMM. */
00090 
00091 /*  Arguments */
00092 /*  ========= */
00093 
00094 /*  M       (input) INTEGER */
00095 /*          The number of rows or columns of the matrix C; C is n-by-m if */
00096 /*          Q is applied from the left, or m-by-n if Q is applied from */
00097 /*          the right.  M >= 0. */
00098 
00099 /*  N       (input) INTEGER */
00100 /*          The order of the orthogonal matrix Q.  N >= 0. */
00101 
00102 /*  K       (input) INTEGER */
00103 /*          The number of elementary reflectors whose product defines the */
00104 /*          orthogonal matrix Q.  N >= K >= 0. */
00105 
00106 /*  AF      (input) DOUBLE PRECISION array, dimension (LDA,N) */
00107 /*          Details of the LQ factorization of an m-by-n matrix, as */
00108 /*          returned by DGELQF. See SGELQF for further details. */
00109 
00110 /*  C       (workspace) DOUBLE PRECISION array, dimension (LDA,N) */
00111 
00112 /*  CC      (workspace) DOUBLE PRECISION array, dimension (LDA,N) */
00113 
00114 /*  Q       (workspace) DOUBLE PRECISION array, dimension (LDA,N) */
00115 
00116 /*  LDA     (input) INTEGER */
00117 /*          The leading dimension of the arrays AF, C, CC, and Q. */
00118 
00119 /*  TAU     (input) DOUBLE PRECISION array, dimension (min(M,N)) */
00120 /*          The scalar factors of the elementary reflectors corresponding */
00121 /*          to the LQ factorization in AF. */
00122 
00123 /*  WORK    (workspace) DOUBLE PRECISION array, dimension (LWORK) */
00124 
00125 /*  LWORK   (input) INTEGER */
00126 /*          The length of WORK.  LWORK must be at least M, and should be */
00127 /*          M*NB, where NB is the blocksize for this environment. */
00128 
00129 /*  RWORK   (workspace) DOUBLE PRECISION array, dimension (M) */
00130 
00131 /*  RESULT  (output) DOUBLE PRECISION array, dimension (4) */
00132 /*          The test ratios compare two techniques for multiplying a */
00133 /*          random matrix C by an n-by-n orthogonal matrix Q. */
00134 /*          RESULT(1) = norm( Q*C - Q*C )  / ( N * norm(C) * EPS ) */
00135 /*          RESULT(2) = norm( C*Q - C*Q )  / ( N * norm(C) * EPS ) */
00136 /*          RESULT(3) = norm( Q'*C - Q'*C )/ ( N * norm(C) * EPS ) */
00137 /*          RESULT(4) = norm( C*Q' - C*Q' )/ ( N * norm(C) * EPS ) */
00138 
00139 /*  ===================================================================== */
00140 
00141 /*     .. Parameters .. */
00142 /*     .. */
00143 /*     .. Local Scalars .. */
00144 /*     .. */
00145 /*     .. External Functions .. */
00146 /*     .. */
00147 /*     .. External Subroutines .. */
00148 /*     .. */
00149 /*     .. Local Arrays .. */
00150 /*     .. */
00151 /*     .. Intrinsic Functions .. */
00152 /*     .. */
00153 /*     .. Scalars in Common .. */
00154 /*     .. */
00155 /*     .. Common blocks .. */
00156 /*     .. */
00157 /*     .. Data statements .. */
00158     /* Parameter adjustments */
00159     q_dim1 = *lda;
00160     q_offset = 1 + q_dim1;
00161     q -= q_offset;
00162     cc_dim1 = *lda;
00163     cc_offset = 1 + cc_dim1;
00164     cc -= cc_offset;
00165     c_dim1 = *lda;
00166     c_offset = 1 + c_dim1;
00167     c__ -= c_offset;
00168     af_dim1 = *lda;
00169     af_offset = 1 + af_dim1;
00170     af -= af_offset;
00171     --tau;
00172     --work;
00173     --rwork;
00174     --result;
00175 
00176     /* Function Body */
00177 /*     .. */
00178 /*     .. Executable Statements .. */
00179 
00180     eps = dlamch_("Epsilon");
00181 
00182 /*     Copy the first k rows of the factorization to the array Q */
00183 
00184     dlaset_("Full", n, n, &c_b4, &c_b4, &q[q_offset], lda);
00185     i__1 = *n - 1;
00186     dlacpy_("Upper", k, &i__1, &af[(af_dim1 << 1) + 1], lda, &q[(q_dim1 << 1) 
00187             + 1], lda);
00188 
00189 /*     Generate the n-by-n matrix Q */
00190 
00191     s_copy(srnamc_1.srnamt, "DORGLQ", (ftnlen)32, (ftnlen)6);
00192     dorglq_(n, n, k, &q[q_offset], lda, &tau[1], &work[1], lwork, &info);
00193 
00194     for (iside = 1; iside <= 2; ++iside) {
00195         if (iside == 1) {
00196             *(unsigned char *)side = 'L';
00197             mc = *n;
00198             nc = *m;
00199         } else {
00200             *(unsigned char *)side = 'R';
00201             mc = *m;
00202             nc = *n;
00203         }
00204 
00205 /*        Generate MC by NC matrix C */
00206 
00207         i__1 = nc;
00208         for (j = 1; j <= i__1; ++j) {
00209             dlarnv_(&c__2, iseed, &mc, &c__[j * c_dim1 + 1]);
00210 /* L10: */
00211         }
00212         cnorm = dlange_("1", &mc, &nc, &c__[c_offset], lda, &rwork[1]);
00213         if (cnorm == 0.) {
00214             cnorm = 1.;
00215         }
00216 
00217         for (itrans = 1; itrans <= 2; ++itrans) {
00218             if (itrans == 1) {
00219                 *(unsigned char *)trans = 'N';
00220             } else {
00221                 *(unsigned char *)trans = 'T';
00222             }
00223 
00224 /*           Copy C */
00225 
00226             dlacpy_("Full", &mc, &nc, &c__[c_offset], lda, &cc[cc_offset], 
00227                     lda);
00228 
00229 /*           Apply Q or Q' to C */
00230 
00231             s_copy(srnamc_1.srnamt, "DORMLQ", (ftnlen)32, (ftnlen)6);
00232             dormlq_(side, trans, &mc, &nc, k, &af[af_offset], lda, &tau[1], &
00233                     cc[cc_offset], lda, &work[1], lwork, &info);
00234 
00235 /*           Form explicit product and subtract */
00236 
00237             if (lsame_(side, "L")) {
00238                 dgemm_(trans, "No transpose", &mc, &nc, &mc, &c_b21, &q[
00239                         q_offset], lda, &c__[c_offset], lda, &c_b22, &cc[
00240                         cc_offset], lda);
00241             } else {
00242                 dgemm_("No transpose", trans, &mc, &nc, &nc, &c_b21, &c__[
00243                         c_offset], lda, &q[q_offset], lda, &c_b22, &cc[
00244                         cc_offset], lda);
00245             }
00246 
00247 /*           Compute error in the difference */
00248 
00249             resid = dlange_("1", &mc, &nc, &cc[cc_offset], lda, &rwork[1]);
00250             result[(iside - 1 << 1) + itrans] = resid / ((doublereal) max(1,*
00251                     n) * cnorm * eps);
00252 
00253 /* L20: */
00254         }
00255 /* L30: */
00256     }
00257 
00258     return 0;
00259 
00260 /*     End of DLQT03 */
00261 
00262 } /* dlqt03_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:47