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 complex c_b1 = {0.f,0.f};
00019 static complex c_b2 = {1.f,0.f};
00020 static integer c__0 = 0;
00021 static integer c__2 = 2;
00022 static integer c__1 = 1;
00023
00024 int cgelsx_(integer *m, integer *n, integer *nrhs, complex *
00025 a, integer *lda, complex *b, integer *ldb, integer *jpvt, real *rcond,
00026 integer *rank, complex *work, real *rwork, integer *info)
00027 {
00028
00029 integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
00030 complex q__1;
00031
00032
00033 double c_abs(complex *);
00034 void r_cnjg(complex *, complex *);
00035
00036
00037 integer i__, j, k;
00038 complex c1, c2, s1, s2, t1, t2;
00039 integer mn;
00040 real anrm, bnrm, smin, smax;
00041 integer iascl, ibscl, ismin, ismax;
00042 extern int ctrsm_(char *, char *, char *, char *,
00043 integer *, integer *, complex *, complex *, integer *, complex *,
00044 integer *), claic1_(integer *,
00045 integer *, complex *, real *, complex *, complex *, real *,
00046 complex *, complex *), cunm2r_(char *, char *, integer *, integer
00047 *, integer *, complex *, integer *, complex *, complex *, integer
00048 *, complex *, integer *), slabad_(real *, real *);
00049 extern doublereal clange_(char *, integer *, integer *, complex *,
00050 integer *, real *);
00051 extern int clascl_(char *, integer *, integer *, real *,
00052 real *, integer *, integer *, complex *, integer *, integer *), cgeqpf_(integer *, integer *, complex *, integer *,
00053 integer *, complex *, complex *, real *, integer *);
00054 extern doublereal slamch_(char *);
00055 extern int claset_(char *, integer *, integer *, complex
00056 *, complex *, complex *, integer *), xerbla_(char *,
00057 integer *);
00058 real bignum;
00059 extern int clatzm_(char *, integer *, integer *, complex
00060 *, integer *, complex *, complex *, complex *, integer *, complex
00061 *);
00062 real sminpr;
00063 extern int ctzrqf_(integer *, integer *, complex *,
00064 integer *, complex *, integer *);
00065 real smaxpr, smlnum;
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
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186 a_dim1 = *lda;
00187 a_offset = 1 + a_dim1;
00188 a -= a_offset;
00189 b_dim1 = *ldb;
00190 b_offset = 1 + b_dim1;
00191 b -= b_offset;
00192 --jpvt;
00193 --work;
00194 --rwork;
00195
00196
00197 mn = min(*m,*n);
00198 ismin = mn + 1;
00199 ismax = (mn << 1) + 1;
00200
00201
00202
00203 *info = 0;
00204 if (*m < 0) {
00205 *info = -1;
00206 } else if (*n < 0) {
00207 *info = -2;
00208 } else if (*nrhs < 0) {
00209 *info = -3;
00210 } else if (*lda < max(1,*m)) {
00211 *info = -5;
00212 } else {
00213
00214 i__1 = max(1,*m);
00215 if (*ldb < max(i__1,*n)) {
00216 *info = -7;
00217 }
00218 }
00219
00220 if (*info != 0) {
00221 i__1 = -(*info);
00222 xerbla_("CGELSX", &i__1);
00223 return 0;
00224 }
00225
00226
00227
00228
00229 i__1 = min(*m,*n);
00230 if (min(i__1,*nrhs) == 0) {
00231 *rank = 0;
00232 return 0;
00233 }
00234
00235
00236
00237 smlnum = slamch_("S") / slamch_("P");
00238 bignum = 1.f / smlnum;
00239 slabad_(&smlnum, &bignum);
00240
00241
00242
00243 anrm = clange_("M", m, n, &a[a_offset], lda, &rwork[1]);
00244 iascl = 0;
00245 if (anrm > 0.f && anrm < smlnum) {
00246
00247
00248
00249 clascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda,
00250 info);
00251 iascl = 1;
00252 } else if (anrm > bignum) {
00253
00254
00255
00256 clascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda,
00257 info);
00258 iascl = 2;
00259 } else if (anrm == 0.f) {
00260
00261
00262
00263 i__1 = max(*m,*n);
00264 claset_("F", &i__1, nrhs, &c_b1, &c_b1, &b[b_offset], ldb);
00265 *rank = 0;
00266 goto L100;
00267 }
00268
00269 bnrm = clange_("M", m, nrhs, &b[b_offset], ldb, &rwork[1]);
00270 ibscl = 0;
00271 if (bnrm > 0.f && bnrm < smlnum) {
00272
00273
00274
00275 clascl_("G", &c__0, &c__0, &bnrm, &smlnum, m, nrhs, &b[b_offset], ldb,
00276 info);
00277 ibscl = 1;
00278 } else if (bnrm > bignum) {
00279
00280
00281
00282 clascl_("G", &c__0, &c__0, &bnrm, &bignum, m, nrhs, &b[b_offset], ldb,
00283 info);
00284 ibscl = 2;
00285 }
00286
00287
00288
00289
00290 cgeqpf_(m, n, &a[a_offset], lda, &jpvt[1], &work[1], &work[mn + 1], &
00291 rwork[1], info);
00292
00293
00294
00295
00296
00297
00298 i__1 = ismin;
00299 work[i__1].r = 1.f, work[i__1].i = 0.f;
00300 i__1 = ismax;
00301 work[i__1].r = 1.f, work[i__1].i = 0.f;
00302 smax = c_abs(&a[a_dim1 + 1]);
00303 smin = smax;
00304 if (c_abs(&a[a_dim1 + 1]) == 0.f) {
00305 *rank = 0;
00306 i__1 = max(*m,*n);
00307 claset_("F", &i__1, nrhs, &c_b1, &c_b1, &b[b_offset], ldb);
00308 goto L100;
00309 } else {
00310 *rank = 1;
00311 }
00312
00313 L10:
00314 if (*rank < mn) {
00315 i__ = *rank + 1;
00316 claic1_(&c__2, rank, &work[ismin], &smin, &a[i__ * a_dim1 + 1], &a[
00317 i__ + i__ * a_dim1], &sminpr, &s1, &c1);
00318 claic1_(&c__1, rank, &work[ismax], &smax, &a[i__ * a_dim1 + 1], &a[
00319 i__ + i__ * a_dim1], &smaxpr, &s2, &c2);
00320
00321 if (smaxpr * *rcond <= sminpr) {
00322 i__1 = *rank;
00323 for (i__ = 1; i__ <= i__1; ++i__) {
00324 i__2 = ismin + i__ - 1;
00325 i__3 = ismin + i__ - 1;
00326 q__1.r = s1.r * work[i__3].r - s1.i * work[i__3].i, q__1.i =
00327 s1.r * work[i__3].i + s1.i * work[i__3].r;
00328 work[i__2].r = q__1.r, work[i__2].i = q__1.i;
00329 i__2 = ismax + i__ - 1;
00330 i__3 = ismax + i__ - 1;
00331 q__1.r = s2.r * work[i__3].r - s2.i * work[i__3].i, q__1.i =
00332 s2.r * work[i__3].i + s2.i * work[i__3].r;
00333 work[i__2].r = q__1.r, work[i__2].i = q__1.i;
00334
00335 }
00336 i__1 = ismin + *rank;
00337 work[i__1].r = c1.r, work[i__1].i = c1.i;
00338 i__1 = ismax + *rank;
00339 work[i__1].r = c2.r, work[i__1].i = c2.i;
00340 smin = sminpr;
00341 smax = smaxpr;
00342 ++(*rank);
00343 goto L10;
00344 }
00345 }
00346
00347
00348
00349
00350
00351
00352
00353 if (*rank < *n) {
00354 ctzrqf_(rank, n, &a[a_offset], lda, &work[mn + 1], info);
00355 }
00356
00357
00358
00359
00360
00361 cunm2r_("Left", "Conjugate transpose", m, nrhs, &mn, &a[a_offset], lda, &
00362 work[1], &b[b_offset], ldb, &work[(mn << 1) + 1], info);
00363
00364
00365
00366
00367
00368 ctrsm_("Left", "Upper", "No transpose", "Non-unit", rank, nrhs, &c_b2, &a[
00369 a_offset], lda, &b[b_offset], ldb);
00370
00371 i__1 = *n;
00372 for (i__ = *rank + 1; i__ <= i__1; ++i__) {
00373 i__2 = *nrhs;
00374 for (j = 1; j <= i__2; ++j) {
00375 i__3 = i__ + j * b_dim1;
00376 b[i__3].r = 0.f, b[i__3].i = 0.f;
00377
00378 }
00379
00380 }
00381
00382
00383
00384 if (*rank < *n) {
00385 i__1 = *rank;
00386 for (i__ = 1; i__ <= i__1; ++i__) {
00387 i__2 = *n - *rank + 1;
00388 r_cnjg(&q__1, &work[mn + i__]);
00389 clatzm_("Left", &i__2, nrhs, &a[i__ + (*rank + 1) * a_dim1], lda,
00390 &q__1, &b[i__ + b_dim1], &b[*rank + 1 + b_dim1], ldb, &
00391 work[(mn << 1) + 1]);
00392
00393 }
00394 }
00395
00396
00397
00398
00399
00400 i__1 = *nrhs;
00401 for (j = 1; j <= i__1; ++j) {
00402 i__2 = *n;
00403 for (i__ = 1; i__ <= i__2; ++i__) {
00404 i__3 = (mn << 1) + i__;
00405 work[i__3].r = 1.f, work[i__3].i = 0.f;
00406
00407 }
00408 i__2 = *n;
00409 for (i__ = 1; i__ <= i__2; ++i__) {
00410 i__3 = (mn << 1) + i__;
00411 if (work[i__3].r == 1.f && work[i__3].i == 0.f) {
00412 if (jpvt[i__] != i__) {
00413 k = i__;
00414 i__3 = k + j * b_dim1;
00415 t1.r = b[i__3].r, t1.i = b[i__3].i;
00416 i__3 = jpvt[k] + j * b_dim1;
00417 t2.r = b[i__3].r, t2.i = b[i__3].i;
00418 L70:
00419 i__3 = jpvt[k] + j * b_dim1;
00420 b[i__3].r = t1.r, b[i__3].i = t1.i;
00421 i__3 = (mn << 1) + k;
00422 work[i__3].r = 0.f, work[i__3].i = 0.f;
00423 t1.r = t2.r, t1.i = t2.i;
00424 k = jpvt[k];
00425 i__3 = jpvt[k] + j * b_dim1;
00426 t2.r = b[i__3].r, t2.i = b[i__3].i;
00427 if (jpvt[k] != i__) {
00428 goto L70;
00429 }
00430 i__3 = i__ + j * b_dim1;
00431 b[i__3].r = t1.r, b[i__3].i = t1.i;
00432 i__3 = (mn << 1) + k;
00433 work[i__3].r = 0.f, work[i__3].i = 0.f;
00434 }
00435 }
00436
00437 }
00438
00439 }
00440
00441
00442
00443 if (iascl == 1) {
00444 clascl_("G", &c__0, &c__0, &anrm, &smlnum, n, nrhs, &b[b_offset], ldb,
00445 info);
00446 clascl_("U", &c__0, &c__0, &smlnum, &anrm, rank, rank, &a[a_offset],
00447 lda, info);
00448 } else if (iascl == 2) {
00449 clascl_("G", &c__0, &c__0, &anrm, &bignum, n, nrhs, &b[b_offset], ldb,
00450 info);
00451 clascl_("U", &c__0, &c__0, &bignum, &anrm, rank, rank, &a[a_offset],
00452 lda, info);
00453 }
00454 if (ibscl == 1) {
00455 clascl_("G", &c__0, &c__0, &smlnum, &bnrm, n, nrhs, &b[b_offset], ldb,
00456 info);
00457 } else if (ibscl == 2) {
00458 clascl_("G", &c__0, &c__0, &bignum, &bnrm, n, nrhs, &b[b_offset], ldb,
00459 info);
00460 }
00461
00462 L100:
00463
00464 return 0;
00465
00466
00467
00468 }