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


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