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 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 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
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
00042 int s_copy(char *, char *, ftnlen, ftnlen);
00043
00044
00045 real eps;
00046 integer info;
00047 extern 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 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 int cungrq_(integer *, integer *, integer *,
00064 complex *, integer *, complex *, complex *, integer *, integer *);
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
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
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
00163 minmn = min(*m,*n);
00164 eps = slamch_("Epsilon");
00165
00166
00167
00168 clacpy_("Full", m, n, &a[a_offset], lda, &af[af_offset], lda);
00169
00170
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
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
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
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
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
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
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
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
00254
00255 }