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 real c_b6 = -1e10f;
00027 static real c_b13 = 0.f;
00028 static real c_b20 = -1.f;
00029 static real c_b21 = 1.f;
00030
00031 int sqlt01_(integer *m, integer *n, real *a, real *af, real *
00032 q, real *l, integer *lda, real *tau, real *work, integer *lwork, real
00033 *rwork, real *result)
00034 {
00035
00036 integer a_dim1, a_offset, af_dim1, af_offset, l_dim1, l_offset, q_dim1,
00037 q_offset, i__1, i__2;
00038
00039
00040 int s_copy(char *, char *, ftnlen, ftnlen);
00041
00042
00043 real eps;
00044 integer info;
00045 real resid;
00046 extern int sgemm_(char *, char *, integer *, integer *,
00047 integer *, real *, real *, integer *, real *, integer *, real *,
00048 real *, integer *);
00049 real anorm;
00050 integer minmn;
00051 extern int ssyrk_(char *, char *, integer *, integer *,
00052 real *, real *, integer *, real *, real *, integer *);
00053 extern doublereal slamch_(char *), slange_(char *, integer *,
00054 integer *, real *, integer *, real *);
00055 extern int sgeqlf_(integer *, integer *, real *, integer
00056 *, real *, real *, integer *, integer *), slacpy_(char *, integer
00057 *, integer *, real *, integer *, real *, integer *),
00058 slaset_(char *, integer *, integer *, real *, real *, real *,
00059 integer *), sorgql_(integer *, integer *, integer *, real
00060 *, integer *, real *, real *, integer *, integer *);
00061 extern doublereal slansy_(char *, char *, integer *, real *, integer *,
00062 real *);
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
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143 l_dim1 = *lda;
00144 l_offset = 1 + l_dim1;
00145 l -= l_offset;
00146 q_dim1 = *lda;
00147 q_offset = 1 + q_dim1;
00148 q -= q_offset;
00149 af_dim1 = *lda;
00150 af_offset = 1 + af_dim1;
00151 af -= af_offset;
00152 a_dim1 = *lda;
00153 a_offset = 1 + a_dim1;
00154 a -= a_offset;
00155 --tau;
00156 --work;
00157 --rwork;
00158 --result;
00159
00160
00161 minmn = min(*m,*n);
00162 eps = slamch_("Epsilon");
00163
00164
00165
00166 slacpy_("Full", m, n, &a[a_offset], lda, &af[af_offset], lda);
00167
00168
00169
00170 s_copy(srnamc_1.srnamt, "SGEQLF", (ftnlen)32, (ftnlen)6);
00171 sgeqlf_(m, n, &af[af_offset], lda, &tau[1], &work[1], lwork, &info);
00172
00173
00174
00175 slaset_("Full", m, m, &c_b6, &c_b6, &q[q_offset], lda);
00176 if (*m >= *n) {
00177 if (*n < *m && *n > 0) {
00178 i__1 = *m - *n;
00179 slacpy_("Full", &i__1, n, &af[af_offset], lda, &q[(*m - *n + 1) *
00180 q_dim1 + 1], lda);
00181 }
00182 if (*n > 1) {
00183 i__1 = *n - 1;
00184 i__2 = *n - 1;
00185 slacpy_("Upper", &i__1, &i__2, &af[*m - *n + 1 + (af_dim1 << 1)],
00186 lda, &q[*m - *n + 1 + (*m - *n + 2) * q_dim1], lda);
00187 }
00188 } else {
00189 if (*m > 1) {
00190 i__1 = *m - 1;
00191 i__2 = *m - 1;
00192 slacpy_("Upper", &i__1, &i__2, &af[(*n - *m + 2) * af_dim1 + 1],
00193 lda, &q[(q_dim1 << 1) + 1], lda);
00194 }
00195 }
00196
00197
00198
00199 s_copy(srnamc_1.srnamt, "SORGQL", (ftnlen)32, (ftnlen)6);
00200 sorgql_(m, m, &minmn, &q[q_offset], lda, &tau[1], &work[1], lwork, &info);
00201
00202
00203
00204 slaset_("Full", m, n, &c_b13, &c_b13, &l[l_offset], lda);
00205 if (*m >= *n) {
00206 if (*n > 0) {
00207 slacpy_("Lower", n, n, &af[*m - *n + 1 + af_dim1], lda, &l[*m - *
00208 n + 1 + l_dim1], lda);
00209 }
00210 } else {
00211 if (*n > *m && *m > 0) {
00212 i__1 = *n - *m;
00213 slacpy_("Full", m, &i__1, &af[af_offset], lda, &l[l_offset], lda);
00214 }
00215 if (*m > 0) {
00216 slacpy_("Lower", m, m, &af[(*n - *m + 1) * af_dim1 + 1], lda, &l[(
00217 *n - *m + 1) * l_dim1 + 1], lda);
00218 }
00219 }
00220
00221
00222
00223 sgemm_("Transpose", "No transpose", m, n, m, &c_b20, &q[q_offset], lda, &
00224 a[a_offset], lda, &c_b21, &l[l_offset], lda);
00225
00226
00227
00228 anorm = slange_("1", m, n, &a[a_offset], lda, &rwork[1]);
00229 resid = slange_("1", m, n, &l[l_offset], lda, &rwork[1]);
00230 if (anorm > 0.f) {
00231 result[1] = resid / (real) max(1,*m) / anorm / eps;
00232 } else {
00233 result[1] = 0.f;
00234 }
00235
00236
00237
00238 slaset_("Full", m, m, &c_b13, &c_b21, &l[l_offset], lda);
00239 ssyrk_("Upper", "Transpose", m, m, &c_b20, &q[q_offset], lda, &c_b21, &l[
00240 l_offset], lda);
00241
00242
00243
00244 resid = slansy_("1", "Upper", m, &l[l_offset], lda, &rwork[1]);
00245
00246 result[2] = resid / (real) max(1,*m) / eps;
00247
00248 return 0;
00249
00250
00251
00252 }