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


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