00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015
00016
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
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
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 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 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
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
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
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
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
00170
00171 if (tpsd) {
00172
00173
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
00186
00187 i__1 = *n + *nrhs;
00188
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
00194
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
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
00206 }
00207
00208 }
00209
00210 } else {
00211
00212
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
00220 }
00221
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
00232
00233 dgelq2_(&ldwork, n, &work[1], &ldwork, &work[ldwork * *n + 1], &work[
00234 ldwork * (*n + 1) + 1], &info);
00235
00236
00237
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
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
00249 }
00250
00251 }
00252
00253 }
00254
00255
00256 i__1 = max(*m,*n);
00257 ret_val = err / ((doublereal) max(i__1,*nrhs) * dlamch_("Epsilon"));
00258
00259 return ret_val;
00260
00261
00262
00263 }