00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015
00016
00017
00018 struct {
00019 char srnamt[32];
00020 } srnamc_;
00021
00022 #define srnamc_1 srnamc_
00023
00024
00025
00026 static doublereal c_b4 = -1e10;
00027 static integer c__2 = 2;
00028 static doublereal c_b21 = -1.;
00029 static doublereal c_b22 = 1.;
00030
00031 int dqrt03_(integer *m, integer *n, integer *k, doublereal *
00032 af, doublereal *c__, doublereal *cc, doublereal *q, integer *lda,
00033 doublereal *tau, doublereal *work, integer *lwork, doublereal *rwork,
00034 doublereal *result)
00035 {
00036
00037
00038 static integer iseed[4] = { 1988,1989,1990,1991 };
00039
00040
00041 integer af_dim1, af_offset, c_dim1, c_offset, cc_dim1, cc_offset, q_dim1,
00042 q_offset, i__1;
00043
00044
00045 int s_copy(char *, char *, ftnlen, ftnlen);
00046
00047
00048 integer j, mc, nc;
00049 doublereal eps;
00050 char side[1];
00051 integer info;
00052 extern int dgemm_(char *, char *, integer *, integer *,
00053 integer *, doublereal *, doublereal *, integer *, doublereal *,
00054 integer *, doublereal *, doublereal *, integer *);
00055 integer iside;
00056 extern logical lsame_(char *, char *);
00057 doublereal resid, cnorm;
00058 char trans[1];
00059 extern doublereal dlamch_(char *), dlange_(char *, integer *,
00060 integer *, doublereal *, integer *, doublereal *);
00061 extern int dlacpy_(char *, integer *, integer *,
00062 doublereal *, integer *, doublereal *, integer *),
00063 dlaset_(char *, integer *, integer *, doublereal *, doublereal *,
00064 doublereal *, integer *), dlarnv_(integer *, integer *,
00065 integer *, doublereal *), dorgqr_(integer *, integer *, integer *,
00066 doublereal *, integer *, doublereal *, doublereal *, integer *,
00067 integer *);
00068 integer itrans;
00069 extern int dormqr_(char *, char *, integer *, integer *,
00070 integer *, doublereal *, integer *, doublereal *, doublereal *,
00071 integer *, doublereal *, integer *, integer *);
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
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160 q_dim1 = *lda;
00161 q_offset = 1 + q_dim1;
00162 q -= q_offset;
00163 cc_dim1 = *lda;
00164 cc_offset = 1 + cc_dim1;
00165 cc -= cc_offset;
00166 c_dim1 = *lda;
00167 c_offset = 1 + c_dim1;
00168 c__ -= c_offset;
00169 af_dim1 = *lda;
00170 af_offset = 1 + af_dim1;
00171 af -= af_offset;
00172 --tau;
00173 --work;
00174 --rwork;
00175 --result;
00176
00177
00178
00179
00180
00181 eps = dlamch_("Epsilon");
00182
00183
00184
00185 dlaset_("Full", m, m, &c_b4, &c_b4, &q[q_offset], lda);
00186 i__1 = *m - 1;
00187 dlacpy_("Lower", &i__1, k, &af[af_dim1 + 2], lda, &q[q_dim1 + 2], lda);
00188
00189
00190
00191 s_copy(srnamc_1.srnamt, "DORGQR", (ftnlen)32, (ftnlen)6);
00192 dorgqr_(m, m, k, &q[q_offset], lda, &tau[1], &work[1], lwork, &info);
00193
00194 for (iside = 1; iside <= 2; ++iside) {
00195 if (iside == 1) {
00196 *(unsigned char *)side = 'L';
00197 mc = *m;
00198 nc = *n;
00199 } else {
00200 *(unsigned char *)side = 'R';
00201 mc = *n;
00202 nc = *m;
00203 }
00204
00205
00206
00207 i__1 = nc;
00208 for (j = 1; j <= i__1; ++j) {
00209 dlarnv_(&c__2, iseed, &mc, &c__[j * c_dim1 + 1]);
00210
00211 }
00212 cnorm = dlange_("1", &mc, &nc, &c__[c_offset], lda, &rwork[1]);
00213 if (cnorm == 0.) {
00214 cnorm = 1.;
00215 }
00216
00217 for (itrans = 1; itrans <= 2; ++itrans) {
00218 if (itrans == 1) {
00219 *(unsigned char *)trans = 'N';
00220 } else {
00221 *(unsigned char *)trans = 'T';
00222 }
00223
00224
00225
00226 dlacpy_("Full", &mc, &nc, &c__[c_offset], lda, &cc[cc_offset],
00227 lda);
00228
00229
00230
00231 s_copy(srnamc_1.srnamt, "DORMQR", (ftnlen)32, (ftnlen)6);
00232 dormqr_(side, trans, &mc, &nc, k, &af[af_offset], lda, &tau[1], &
00233 cc[cc_offset], lda, &work[1], lwork, &info);
00234
00235
00236
00237 if (lsame_(side, "L")) {
00238 dgemm_(trans, "No transpose", &mc, &nc, &mc, &c_b21, &q[
00239 q_offset], lda, &c__[c_offset], lda, &c_b22, &cc[
00240 cc_offset], lda);
00241 } else {
00242 dgemm_("No transpose", trans, &mc, &nc, &nc, &c_b21, &c__[
00243 c_offset], lda, &q[q_offset], lda, &c_b22, &cc[
00244 cc_offset], lda);
00245 }
00246
00247
00248
00249 resid = dlange_("1", &mc, &nc, &cc[cc_offset], lda, &rwork[1]);
00250 result[(iside - 1 << 1) + itrans] = resid / ((doublereal) max(1,*
00251 m) * cnorm * eps);
00252
00253
00254 }
00255
00256 }
00257
00258 return 0;
00259
00260
00261
00262 }