cqrt01.c
Go to the documentation of this file.
00001 /* cqrt01.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_b10 = {0.f,0.f};
00028 static complex c_b15 = {-1.f,0.f};
00029 static complex c_b16 = {1.f,0.f};
00030 static real c_b24 = -1.f;
00031 static real c_b25 = 1.f;
00032 
00033 /* Subroutine */ int cqrt01_(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;
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 cgeqrf_(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 cungqr_(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 /*  CQRT01 tests CGEQRF, which computes the QR factorization of an m-by-n */
00080 /*  matrix A, and partially tests CUNGQR which forms the m-by-m */
00081 /*  orthogonal matrix Q. */
00082 
00083 /*  CQRT01 compares R with Q'*A, 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 QR factorization of A, as returned by CGEQRF. */
00099 /*          See CGEQRF for further details. */
00100 
00101 /*  Q       (output) COMPLEX array, dimension (LDA,M) */
00102 /*          The m-by-m 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 R. */
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 CGEQRF. */
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 (M) */
00120 
00121 /*  RESULT  (output) REAL array, dimension (2) */
00122 /*          The test ratios: */
00123 /*          RESULT(1) = norm( R - Q'*A ) / ( M * norm(A) * EPS ) */
00124 /*          RESULT(2) = norm( I - Q'*Q ) / ( M * 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, "CGEQRF", (ftnlen)32, (ftnlen)6);
00173     cgeqrf_(m, n, &af[af_offset], lda, &tau[1], &work[1], lwork, &info);
00174 
00175 /*     Copy details of Q */
00176 
00177     claset_("Full", m, m, &c_b1, &c_b1, &q[q_offset], lda);
00178     i__1 = *m - 1;
00179     clacpy_("Lower", &i__1, n, &af[af_dim1 + 2], lda, &q[q_dim1 + 2], lda);
00180 
00181 /*     Generate the m-by-m matrix Q */
00182 
00183     s_copy(srnamc_1.srnamt, "CUNGQR", (ftnlen)32, (ftnlen)6);
00184     cungqr_(m, m, &minmn, &q[q_offset], lda, &tau[1], &work[1], lwork, &info);
00185 
00186 /*     Copy R */
00187 
00188     claset_("Full", m, n, &c_b10, &c_b10, &r__[r_offset], lda);
00189     clacpy_("Upper", m, n, &af[af_offset], lda, &r__[r_offset], lda);
00190 
00191 /*     Compute R - Q'*A */
00192 
00193     cgemm_("Conjugate transpose", "No transpose", m, n, m, &c_b15, &q[
00194             q_offset], lda, &a[a_offset], lda, &c_b16, &r__[r_offset], lda);
00195 
00196 /*     Compute norm( R - Q'*A ) / ( M * norm(A) * EPS ) . */
00197 
00198     anorm = clange_("1", m, n, &a[a_offset], lda, &rwork[1]);
00199     resid = clange_("1", m, n, &r__[r_offset], lda, &rwork[1]);
00200     if (anorm > 0.f) {
00201         result[1] = resid / (real) max(1,*m) / anorm / eps;
00202     } else {
00203         result[1] = 0.f;
00204     }
00205 
00206 /*     Compute I - Q'*Q */
00207 
00208     claset_("Full", m, m, &c_b10, &c_b16, &r__[r_offset], lda);
00209     cherk_("Upper", "Conjugate transpose", m, m, &c_b24, &q[q_offset], lda, &
00210             c_b25, &r__[r_offset], lda);
00211 
00212 /*     Compute norm( I - Q'*Q ) / ( M * EPS ) . */
00213 
00214     resid = clansy_("1", "Upper", m, &r__[r_offset], lda, &rwork[1]);
00215 
00216     result[2] = resid / (real) max(1,*m) / eps;
00217 
00218     return 0;
00219 
00220 /*     End of CQRT01 */
00221 
00222 } /* cqrt01_ */


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