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


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