dqrt14.c
Go to the documentation of this file.
00001 /* dqrt14.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 /* Table of constant values */
00017 
00018 static integer c__10 = 10;
00019 static integer c__1 = 1;
00020 static integer c__0 = 0;
00021 static doublereal c_b15 = 1.;
00022 
00023 doublereal dqrt14_(char *trans, integer *m, integer *n, integer *nrhs, 
00024         doublereal *a, integer *lda, doublereal *x, integer *ldx, doublereal *
00025         work, integer *lwork)
00026 {
00027     /* System generated locals */
00028     integer a_dim1, a_offset, x_dim1, x_offset, i__1, i__2, i__3;
00029     doublereal ret_val, d__1, d__2, d__3;
00030 
00031     /* Local variables */
00032     integer i__, j;
00033     doublereal err;
00034     integer info;
00035     doublereal anrm;
00036     logical tpsd;
00037     doublereal xnrm;
00038     extern logical lsame_(char *, char *);
00039     doublereal rwork[1];
00040     extern /* Subroutine */ int dgelq2_(integer *, integer *, doublereal *, 
00041             integer *, doublereal *, doublereal *, integer *), dgeqr2_(
00042             integer *, integer *, doublereal *, integer *, doublereal *, 
00043             doublereal *, integer *);
00044     extern doublereal dlamch_(char *), dlange_(char *, integer *, 
00045             integer *, doublereal *, integer *, doublereal *);
00046     extern /* Subroutine */ int dlascl_(char *, integer *, integer *, 
00047             doublereal *, doublereal *, integer *, integer *, doublereal *, 
00048             integer *, integer *), dlacpy_(char *, integer *, integer 
00049             *, doublereal *, integer *, doublereal *, integer *), 
00050             xerbla_(char *, integer *);
00051     integer ldwork;
00052 
00053 
00054 /*  -- LAPACK test routine (version 3.1) -- */
00055 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00056 /*     November 2006 */
00057 
00058 /*     .. Scalar Arguments .. */
00059 /*     .. */
00060 /*     .. Array Arguments .. */
00061 /*     .. */
00062 
00063 /*  Purpose */
00064 /*  ======= */
00065 
00066 /*  DQRT14 checks whether X is in the row space of A or A'.  It does so */
00067 /*  by scaling both X and A such that their norms are in the range */
00068 /*  [sqrt(eps), 1/sqrt(eps)], then computing a QR factorization of [A,X] */
00069 /*  (if TRANS = 'T') or an LQ factorization of [A',X]' (if TRANS = 'N'), */
00070 /*  and returning the norm of the trailing triangle, scaled by */
00071 /*  MAX(M,N,NRHS)*eps. */
00072 
00073 /*  Arguments */
00074 /*  ========= */
00075 
00076 /*  TRANS   (input) CHARACTER*1 */
00077 /*          = 'N':  No transpose, check for X in the row space of A */
00078 /*          = 'T':  Transpose, check for X in the row space of A'. */
00079 
00080 /*  M       (input) INTEGER */
00081 /*          The number of rows of the matrix A. */
00082 
00083 /*  N       (input) INTEGER */
00084 /*          The number of columns of the matrix A. */
00085 
00086 /*  NRHS    (input) INTEGER */
00087 /*          The number of right hand sides, i.e., the number of columns */
00088 /*          of X. */
00089 
00090 /*  A       (input) DOUBLE PRECISION array, dimension (LDA,N) */
00091 /*          The M-by-N matrix A. */
00092 
00093 /*  LDA     (input) INTEGER */
00094 /*          The leading dimension of the array A. */
00095 
00096 /*  X       (input) DOUBLE PRECISION array, dimension (LDX,NRHS) */
00097 /*          If TRANS = 'N', the N-by-NRHS matrix X. */
00098 /*          IF TRANS = 'T', the M-by-NRHS matrix X. */
00099 
00100 /*  LDX     (input) INTEGER */
00101 /*          The leading dimension of the array X. */
00102 
00103 /*  WORK    (workspace) DOUBLE PRECISION array dimension (LWORK) */
00104 
00105 /*  LWORK   (input) INTEGER */
00106 /*          length of workspace array required */
00107 /*          If TRANS = 'N', LWORK >= (M+NRHS)*(N+2); */
00108 /*          if TRANS = 'T', LWORK >= (N+NRHS)*(M+2). */
00109 
00110 /*  ===================================================================== */
00111 
00112 /*     .. Parameters .. */
00113 /*     .. */
00114 /*     .. Local Scalars .. */
00115 /*     .. */
00116 /*     .. Local Arrays .. */
00117 /*     .. */
00118 /*     .. External Functions .. */
00119 /*     .. */
00120 /*     .. External Subroutines .. */
00121 /*     .. */
00122 /*     .. Intrinsic Functions .. */
00123 /*     .. */
00124 /*     .. Executable Statements .. */
00125 
00126     /* Parameter adjustments */
00127     a_dim1 = *lda;
00128     a_offset = 1 + a_dim1;
00129     a -= a_offset;
00130     x_dim1 = *ldx;
00131     x_offset = 1 + x_dim1;
00132     x -= x_offset;
00133     --work;
00134 
00135     /* Function Body */
00136     ret_val = 0.;
00137     if (lsame_(trans, "N")) {
00138         ldwork = *m + *nrhs;
00139         tpsd = FALSE_;
00140         if (*lwork < (*m + *nrhs) * (*n + 2)) {
00141             xerbla_("DQRT14", &c__10);
00142             return ret_val;
00143         } else if (*n <= 0 || *nrhs <= 0) {
00144             return ret_val;
00145         }
00146     } else if (lsame_(trans, "T")) {
00147         ldwork = *m;
00148         tpsd = TRUE_;
00149         if (*lwork < (*n + *nrhs) * (*m + 2)) {
00150             xerbla_("DQRT14", &c__10);
00151             return ret_val;
00152         } else if (*m <= 0 || *nrhs <= 0) {
00153             return ret_val;
00154         }
00155     } else {
00156         xerbla_("DQRT14", &c__1);
00157         return ret_val;
00158     }
00159 
00160 /*     Copy and scale A */
00161 
00162     dlacpy_("All", m, n, &a[a_offset], lda, &work[1], &ldwork);
00163     anrm = dlange_("M", m, n, &work[1], &ldwork, rwork);
00164     if (anrm != 0.) {
00165         dlascl_("G", &c__0, &c__0, &anrm, &c_b15, m, n, &work[1], &ldwork, &
00166                 info);
00167     }
00168 
00169 /*     Copy X or X' into the right place and scale it */
00170 
00171     if (tpsd) {
00172 
00173 /*        Copy X into columns n+1:n+nrhs of work */
00174 
00175         dlacpy_("All", m, nrhs, &x[x_offset], ldx, &work[*n * ldwork + 1], &
00176                 ldwork);
00177         xnrm = dlange_("M", m, nrhs, &work[*n * ldwork + 1], &ldwork, rwork);
00178         if (xnrm != 0.) {
00179             dlascl_("G", &c__0, &c__0, &xnrm, &c_b15, m, nrhs, &work[*n * 
00180                     ldwork + 1], &ldwork, &info);
00181         }
00182         i__1 = *n + *nrhs;
00183         anrm = dlange_("One-norm", m, &i__1, &work[1], &ldwork, rwork);
00184 
00185 /*        Compute QR factorization of X */
00186 
00187         i__1 = *n + *nrhs;
00188 /* Computing MIN */
00189         i__2 = *m, i__3 = *n + *nrhs;
00190         dgeqr2_(m, &i__1, &work[1], &ldwork, &work[ldwork * (*n + *nrhs) + 1], 
00191                  &work[ldwork * (*n + *nrhs) + min(i__2, i__3)+ 1], &info);
00192 
00193 /*        Compute largest entry in upper triangle of */
00194 /*        work(n+1:m,n+1:n+nrhs) */
00195 
00196         err = 0.;
00197         i__1 = *n + *nrhs;
00198         for (j = *n + 1; j <= i__1; ++j) {
00199             i__2 = min(*m,j);
00200             for (i__ = *n + 1; i__ <= i__2; ++i__) {
00201 /* Computing MAX */
00202                 d__2 = err, d__3 = (d__1 = work[i__ + (j - 1) * *m], abs(d__1)
00203                         );
00204                 err = max(d__2,d__3);
00205 /* L10: */
00206             }
00207 /* L20: */
00208         }
00209 
00210     } else {
00211 
00212 /*        Copy X' into rows m+1:m+nrhs of work */
00213 
00214         i__1 = *n;
00215         for (i__ = 1; i__ <= i__1; ++i__) {
00216             i__2 = *nrhs;
00217             for (j = 1; j <= i__2; ++j) {
00218                 work[*m + j + (i__ - 1) * ldwork] = x[i__ + j * x_dim1];
00219 /* L30: */
00220             }
00221 /* L40: */
00222         }
00223 
00224         xnrm = dlange_("M", nrhs, n, &work[*m + 1], &ldwork, rwork)
00225                 ;
00226         if (xnrm != 0.) {
00227             dlascl_("G", &c__0, &c__0, &xnrm, &c_b15, nrhs, n, &work[*m + 1], 
00228                     &ldwork, &info);
00229         }
00230 
00231 /*        Compute LQ factorization of work */
00232 
00233         dgelq2_(&ldwork, n, &work[1], &ldwork, &work[ldwork * *n + 1], &work[
00234                 ldwork * (*n + 1) + 1], &info);
00235 
00236 /*        Compute largest entry in lower triangle in */
00237 /*        work(m+1:m+nrhs,m+1:n) */
00238 
00239         err = 0.;
00240         i__1 = *n;
00241         for (j = *m + 1; j <= i__1; ++j) {
00242             i__2 = ldwork;
00243             for (i__ = j; i__ <= i__2; ++i__) {
00244 /* Computing MAX */
00245                 d__2 = err, d__3 = (d__1 = work[i__ + (j - 1) * ldwork], abs(
00246                         d__1));
00247                 err = max(d__2,d__3);
00248 /* L50: */
00249             }
00250 /* L60: */
00251         }
00252 
00253     }
00254 
00255 /* Computing MAX */
00256     i__1 = max(*m,*n);
00257     ret_val = err / ((doublereal) max(i__1,*nrhs) * dlamch_("Epsilon"));
00258 
00259     return ret_val;
00260 
00261 /*     End of DQRT14 */
00262 
00263 } /* dqrt14_ */


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