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__1 = 1;
00019 static real c_b8 = -1.f;
00020 static real c_b9 = 1.f;
00021 static real c_b16 = 0.f;
00022
00023 int slaqps_(integer *m, integer *n, integer *offset, integer
00024 *nb, integer *kb, real *a, integer *lda, integer *jpvt, real *tau,
00025 real *vn1, real *vn2, real *auxv, real *f, integer *ldf)
00026 {
00027
00028 integer a_dim1, a_offset, f_dim1, f_offset, i__1, i__2;
00029 real r__1, r__2;
00030
00031
00032 double sqrt(doublereal);
00033 integer i_nint(real *);
00034
00035
00036 integer j, k, rk;
00037 real akk;
00038 integer pvt;
00039 real temp, temp2;
00040 extern doublereal snrm2_(integer *, real *, integer *);
00041 real tol3z;
00042 extern int sgemm_(char *, char *, integer *, integer *,
00043 integer *, real *, real *, integer *, real *, integer *, real *,
00044 real *, integer *);
00045 integer itemp;
00046 extern int sgemv_(char *, integer *, integer *, real *,
00047 real *, integer *, real *, integer *, real *, real *, integer *), sswap_(integer *, real *, integer *, real *, integer *);
00048 extern doublereal slamch_(char *);
00049 integer lsticc;
00050 extern integer isamax_(integer *, real *, integer *);
00051 extern int slarfp_(integer *, real *, real *, integer *,
00052 real *);
00053 integer lastrk;
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
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 a_dim1 = *lda;
00159 a_offset = 1 + a_dim1;
00160 a -= a_offset;
00161 --jpvt;
00162 --tau;
00163 --vn1;
00164 --vn2;
00165 --auxv;
00166 f_dim1 = *ldf;
00167 f_offset = 1 + f_dim1;
00168 f -= f_offset;
00169
00170
00171
00172 i__1 = *m, i__2 = *n + *offset;
00173 lastrk = min(i__1,i__2);
00174 lsticc = 0;
00175 k = 0;
00176 tol3z = sqrt(slamch_("Epsilon"));
00177
00178
00179
00180 L10:
00181 if (k < *nb && lsticc == 0) {
00182 ++k;
00183 rk = *offset + k;
00184
00185
00186
00187 i__1 = *n - k + 1;
00188 pvt = k - 1 + isamax_(&i__1, &vn1[k], &c__1);
00189 if (pvt != k) {
00190 sswap_(m, &a[pvt * a_dim1 + 1], &c__1, &a[k * a_dim1 + 1], &c__1);
00191 i__1 = k - 1;
00192 sswap_(&i__1, &f[pvt + f_dim1], ldf, &f[k + f_dim1], ldf);
00193 itemp = jpvt[pvt];
00194 jpvt[pvt] = jpvt[k];
00195 jpvt[k] = itemp;
00196 vn1[pvt] = vn1[k];
00197 vn2[pvt] = vn2[k];
00198 }
00199
00200
00201
00202
00203 if (k > 1) {
00204 i__1 = *m - rk + 1;
00205 i__2 = k - 1;
00206 sgemv_("No transpose", &i__1, &i__2, &c_b8, &a[rk + a_dim1], lda,
00207 &f[k + f_dim1], ldf, &c_b9, &a[rk + k * a_dim1], &c__1);
00208 }
00209
00210
00211
00212 if (rk < *m) {
00213 i__1 = *m - rk + 1;
00214 slarfp_(&i__1, &a[rk + k * a_dim1], &a[rk + 1 + k * a_dim1], &
00215 c__1, &tau[k]);
00216 } else {
00217 slarfp_(&c__1, &a[rk + k * a_dim1], &a[rk + k * a_dim1], &c__1, &
00218 tau[k]);
00219 }
00220
00221 akk = a[rk + k * a_dim1];
00222 a[rk + k * a_dim1] = 1.f;
00223
00224
00225
00226
00227
00228 if (k < *n) {
00229 i__1 = *m - rk + 1;
00230 i__2 = *n - k;
00231 sgemv_("Transpose", &i__1, &i__2, &tau[k], &a[rk + (k + 1) *
00232 a_dim1], lda, &a[rk + k * a_dim1], &c__1, &c_b16, &f[k +
00233 1 + k * f_dim1], &c__1);
00234 }
00235
00236
00237
00238 i__1 = k;
00239 for (j = 1; j <= i__1; ++j) {
00240 f[j + k * f_dim1] = 0.f;
00241
00242 }
00243
00244
00245
00246
00247
00248 if (k > 1) {
00249 i__1 = *m - rk + 1;
00250 i__2 = k - 1;
00251 r__1 = -tau[k];
00252 sgemv_("Transpose", &i__1, &i__2, &r__1, &a[rk + a_dim1], lda, &a[
00253 rk + k * a_dim1], &c__1, &c_b16, &auxv[1], &c__1);
00254
00255 i__1 = k - 1;
00256 sgemv_("No transpose", n, &i__1, &c_b9, &f[f_dim1 + 1], ldf, &
00257 auxv[1], &c__1, &c_b9, &f[k * f_dim1 + 1], &c__1);
00258 }
00259
00260
00261
00262
00263 if (k < *n) {
00264 i__1 = *n - k;
00265 sgemv_("No transpose", &i__1, &k, &c_b8, &f[k + 1 + f_dim1], ldf,
00266 &a[rk + a_dim1], lda, &c_b9, &a[rk + (k + 1) * a_dim1],
00267 lda);
00268 }
00269
00270
00271
00272 if (rk < lastrk) {
00273 i__1 = *n;
00274 for (j = k + 1; j <= i__1; ++j) {
00275 if (vn1[j] != 0.f) {
00276
00277
00278
00279
00280 temp = (r__1 = a[rk + j * a_dim1], dabs(r__1)) / vn1[j];
00281
00282 r__1 = 0.f, r__2 = (temp + 1.f) * (1.f - temp);
00283 temp = dmax(r__1,r__2);
00284
00285 r__1 = vn1[j] / vn2[j];
00286 temp2 = temp * (r__1 * r__1);
00287 if (temp2 <= tol3z) {
00288 vn2[j] = (real) lsticc;
00289 lsticc = j;
00290 } else {
00291 vn1[j] *= sqrt(temp);
00292 }
00293 }
00294
00295 }
00296 }
00297
00298 a[rk + k * a_dim1] = akk;
00299
00300
00301
00302 goto L10;
00303 }
00304 *kb = k;
00305 rk = *offset + *kb;
00306
00307
00308
00309
00310
00311
00312 i__1 = *n, i__2 = *m - *offset;
00313 if (*kb < min(i__1,i__2)) {
00314 i__1 = *m - rk;
00315 i__2 = *n - *kb;
00316 sgemm_("No transpose", "Transpose", &i__1, &i__2, kb, &c_b8, &a[rk +
00317 1 + a_dim1], lda, &f[*kb + 1 + f_dim1], ldf, &c_b9, &a[rk + 1
00318 + (*kb + 1) * a_dim1], lda);
00319 }
00320
00321
00322
00323 L40:
00324 if (lsticc > 0) {
00325 itemp = i_nint(&vn2[lsticc]);
00326 i__1 = *m - rk;
00327 vn1[lsticc] = snrm2_(&i__1, &a[rk + 1 + lsticc * a_dim1], &c__1);
00328
00329
00330
00331
00332
00333 vn2[lsticc] = vn1[lsticc];
00334 lsticc = itemp;
00335 goto L40;
00336 }
00337
00338 return 0;
00339
00340
00341
00342 }