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 integer c__1 = 1;
00020 static integer c__0 = 0;
00021 static real c_b10 = 1.f;
00022 static real c_b35 = 0.f;
00023
00024 int clalsd_(char *uplo, integer *smlsiz, integer *n, integer
00025 *nrhs, real *d__, real *e, complex *b, integer *ldb, real *rcond,
00026 integer *rank, complex *work, real *rwork, integer *iwork, integer *
00027 info)
00028 {
00029
00030 integer b_dim1, b_offset, i__1, i__2, i__3, i__4, i__5, i__6;
00031 real r__1;
00032 complex q__1;
00033
00034
00035 double r_imag(complex *), log(doublereal), r_sign(real *, real *);
00036
00037
00038 integer c__, i__, j, k;
00039 real r__;
00040 integer s, u, z__;
00041 real cs;
00042 integer bx;
00043 real sn;
00044 integer st, vt, nm1, st1;
00045 real eps;
00046 integer iwk;
00047 real tol;
00048 integer difl, difr;
00049 real rcnd;
00050 integer jcol, irwb, perm, nsub, nlvl, sqre, bxst, jrow, irwu, jimag,
00051 jreal;
00052 extern int sgemm_(char *, char *, integer *, integer *,
00053 integer *, real *, real *, integer *, real *, integer *, real *,
00054 real *, integer *);
00055 integer irwib;
00056 extern int ccopy_(integer *, complex *, integer *,
00057 complex *, integer *);
00058 integer poles, sizei, irwrb, nsize;
00059 extern int csrot_(integer *, complex *, integer *,
00060 complex *, integer *, real *, real *);
00061 integer irwvt, icmpq1, icmpq2;
00062 extern int clalsa_(integer *, integer *, integer *,
00063 integer *, complex *, integer *, complex *, integer *, real *,
00064 integer *, real *, integer *, real *, real *, real *, real *,
00065 integer *, integer *, integer *, integer *, real *, real *, real *
00066 , real *, integer *, integer *), clascl_(char *, integer *,
00067 integer *, real *, real *, integer *, integer *, complex *,
00068 integer *, integer *);
00069 extern doublereal slamch_(char *);
00070 extern int slasda_(integer *, integer *, integer *,
00071 integer *, real *, real *, real *, integer *, real *, integer *,
00072 real *, real *, real *, real *, integer *, integer *, integer *,
00073 integer *, real *, real *, real *, real *, integer *, integer *),
00074 clacpy_(char *, integer *, integer *, complex *, integer *,
00075 complex *, integer *), claset_(char *, integer *, integer
00076 *, complex *, complex *, complex *, integer *), xerbla_(
00077 char *, integer *), slascl_(char *, integer *, integer *,
00078 real *, real *, integer *, integer *, real *, integer *, integer *
00079 );
00080 extern integer isamax_(integer *, real *, integer *);
00081 integer givcol;
00082 extern int slasdq_(char *, integer *, integer *, integer
00083 *, integer *, integer *, real *, real *, real *, integer *, real *
00084 , integer *, real *, integer *, real *, integer *),
00085 slaset_(char *, integer *, integer *, real *, real *, real *,
00086 integer *), slartg_(real *, real *, real *, real *, real *
00087 );
00088 real orgnrm;
00089 integer givnum;
00090 extern doublereal slanst_(char *, integer *, real *, real *);
00091 extern int slasrt_(char *, integer *, real *, integer *);
00092 integer givptr, nrwork, irwwrk, smlszp;
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
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213 --d__;
00214 --e;
00215 b_dim1 = *ldb;
00216 b_offset = 1 + b_dim1;
00217 b -= b_offset;
00218 --work;
00219 --rwork;
00220 --iwork;
00221
00222
00223 *info = 0;
00224
00225 if (*n < 0) {
00226 *info = -3;
00227 } else if (*nrhs < 1) {
00228 *info = -4;
00229 } else if (*ldb < 1 || *ldb < *n) {
00230 *info = -8;
00231 }
00232 if (*info != 0) {
00233 i__1 = -(*info);
00234 xerbla_("CLALSD", &i__1);
00235 return 0;
00236 }
00237
00238 eps = slamch_("Epsilon");
00239
00240
00241
00242 if (*rcond <= 0.f || *rcond >= 1.f) {
00243 rcnd = eps;
00244 } else {
00245 rcnd = *rcond;
00246 }
00247
00248 *rank = 0;
00249
00250
00251
00252 if (*n == 0) {
00253 return 0;
00254 } else if (*n == 1) {
00255 if (d__[1] == 0.f) {
00256 claset_("A", &c__1, nrhs, &c_b1, &c_b1, &b[b_offset], ldb);
00257 } else {
00258 *rank = 1;
00259 clascl_("G", &c__0, &c__0, &d__[1], &c_b10, &c__1, nrhs, &b[
00260 b_offset], ldb, info);
00261 d__[1] = dabs(d__[1]);
00262 }
00263 return 0;
00264 }
00265
00266
00267
00268 if (*(unsigned char *)uplo == 'L') {
00269 i__1 = *n - 1;
00270 for (i__ = 1; i__ <= i__1; ++i__) {
00271 slartg_(&d__[i__], &e[i__], &cs, &sn, &r__);
00272 d__[i__] = r__;
00273 e[i__] = sn * d__[i__ + 1];
00274 d__[i__ + 1] = cs * d__[i__ + 1];
00275 if (*nrhs == 1) {
00276 csrot_(&c__1, &b[i__ + b_dim1], &c__1, &b[i__ + 1 + b_dim1], &
00277 c__1, &cs, &sn);
00278 } else {
00279 rwork[(i__ << 1) - 1] = cs;
00280 rwork[i__ * 2] = sn;
00281 }
00282
00283 }
00284 if (*nrhs > 1) {
00285 i__1 = *nrhs;
00286 for (i__ = 1; i__ <= i__1; ++i__) {
00287 i__2 = *n - 1;
00288 for (j = 1; j <= i__2; ++j) {
00289 cs = rwork[(j << 1) - 1];
00290 sn = rwork[j * 2];
00291 csrot_(&c__1, &b[j + i__ * b_dim1], &c__1, &b[j + 1 + i__
00292 * b_dim1], &c__1, &cs, &sn);
00293
00294 }
00295
00296 }
00297 }
00298 }
00299
00300
00301
00302 nm1 = *n - 1;
00303 orgnrm = slanst_("M", n, &d__[1], &e[1]);
00304 if (orgnrm == 0.f) {
00305 claset_("A", n, nrhs, &c_b1, &c_b1, &b[b_offset], ldb);
00306 return 0;
00307 }
00308
00309 slascl_("G", &c__0, &c__0, &orgnrm, &c_b10, n, &c__1, &d__[1], n, info);
00310 slascl_("G", &c__0, &c__0, &orgnrm, &c_b10, &nm1, &c__1, &e[1], &nm1,
00311 info);
00312
00313
00314
00315
00316 if (*n <= *smlsiz) {
00317 irwu = 1;
00318 irwvt = irwu + *n * *n;
00319 irwwrk = irwvt + *n * *n;
00320 irwrb = irwwrk;
00321 irwib = irwrb + *n * *nrhs;
00322 irwb = irwib + *n * *nrhs;
00323 slaset_("A", n, n, &c_b35, &c_b10, &rwork[irwu], n);
00324 slaset_("A", n, n, &c_b35, &c_b10, &rwork[irwvt], n);
00325 slasdq_("U", &c__0, n, n, n, &c__0, &d__[1], &e[1], &rwork[irwvt], n,
00326 &rwork[irwu], n, &rwork[irwwrk], &c__1, &rwork[irwwrk], info);
00327 if (*info != 0) {
00328 return 0;
00329 }
00330
00331
00332
00333
00334
00335 j = irwb - 1;
00336 i__1 = *nrhs;
00337 for (jcol = 1; jcol <= i__1; ++jcol) {
00338 i__2 = *n;
00339 for (jrow = 1; jrow <= i__2; ++jrow) {
00340 ++j;
00341 i__3 = jrow + jcol * b_dim1;
00342 rwork[j] = b[i__3].r;
00343
00344 }
00345
00346 }
00347 sgemm_("T", "N", n, nrhs, n, &c_b10, &rwork[irwu], n, &rwork[irwb], n,
00348 &c_b35, &rwork[irwrb], n);
00349 j = irwb - 1;
00350 i__1 = *nrhs;
00351 for (jcol = 1; jcol <= i__1; ++jcol) {
00352 i__2 = *n;
00353 for (jrow = 1; jrow <= i__2; ++jrow) {
00354 ++j;
00355 rwork[j] = r_imag(&b[jrow + jcol * b_dim1]);
00356
00357 }
00358
00359 }
00360 sgemm_("T", "N", n, nrhs, n, &c_b10, &rwork[irwu], n, &rwork[irwb], n,
00361 &c_b35, &rwork[irwib], n);
00362 jreal = irwrb - 1;
00363 jimag = irwib - 1;
00364 i__1 = *nrhs;
00365 for (jcol = 1; jcol <= i__1; ++jcol) {
00366 i__2 = *n;
00367 for (jrow = 1; jrow <= i__2; ++jrow) {
00368 ++jreal;
00369 ++jimag;
00370 i__3 = jrow + jcol * b_dim1;
00371 i__4 = jreal;
00372 i__5 = jimag;
00373 q__1.r = rwork[i__4], q__1.i = rwork[i__5];
00374 b[i__3].r = q__1.r, b[i__3].i = q__1.i;
00375
00376 }
00377
00378 }
00379
00380 tol = rcnd * (r__1 = d__[isamax_(n, &d__[1], &c__1)], dabs(r__1));
00381 i__1 = *n;
00382 for (i__ = 1; i__ <= i__1; ++i__) {
00383 if (d__[i__] <= tol) {
00384 claset_("A", &c__1, nrhs, &c_b1, &c_b1, &b[i__ + b_dim1], ldb);
00385 } else {
00386 clascl_("G", &c__0, &c__0, &d__[i__], &c_b10, &c__1, nrhs, &b[
00387 i__ + b_dim1], ldb, info);
00388 ++(*rank);
00389 }
00390
00391 }
00392
00393
00394
00395
00396
00397
00398
00399
00400 j = irwb - 1;
00401 i__1 = *nrhs;
00402 for (jcol = 1; jcol <= i__1; ++jcol) {
00403 i__2 = *n;
00404 for (jrow = 1; jrow <= i__2; ++jrow) {
00405 ++j;
00406 i__3 = jrow + jcol * b_dim1;
00407 rwork[j] = b[i__3].r;
00408
00409 }
00410
00411 }
00412 sgemm_("T", "N", n, nrhs, n, &c_b10, &rwork[irwvt], n, &rwork[irwb],
00413 n, &c_b35, &rwork[irwrb], n);
00414 j = irwb - 1;
00415 i__1 = *nrhs;
00416 for (jcol = 1; jcol <= i__1; ++jcol) {
00417 i__2 = *n;
00418 for (jrow = 1; jrow <= i__2; ++jrow) {
00419 ++j;
00420 rwork[j] = r_imag(&b[jrow + jcol * b_dim1]);
00421
00422 }
00423
00424 }
00425 sgemm_("T", "N", n, nrhs, n, &c_b10, &rwork[irwvt], n, &rwork[irwb],
00426 n, &c_b35, &rwork[irwib], n);
00427 jreal = irwrb - 1;
00428 jimag = irwib - 1;
00429 i__1 = *nrhs;
00430 for (jcol = 1; jcol <= i__1; ++jcol) {
00431 i__2 = *n;
00432 for (jrow = 1; jrow <= i__2; ++jrow) {
00433 ++jreal;
00434 ++jimag;
00435 i__3 = jrow + jcol * b_dim1;
00436 i__4 = jreal;
00437 i__5 = jimag;
00438 q__1.r = rwork[i__4], q__1.i = rwork[i__5];
00439 b[i__3].r = q__1.r, b[i__3].i = q__1.i;
00440
00441 }
00442
00443 }
00444
00445
00446
00447 slascl_("G", &c__0, &c__0, &c_b10, &orgnrm, n, &c__1, &d__[1], n,
00448 info);
00449 slasrt_("D", n, &d__[1], info);
00450 clascl_("G", &c__0, &c__0, &orgnrm, &c_b10, n, nrhs, &b[b_offset],
00451 ldb, info);
00452
00453 return 0;
00454 }
00455
00456
00457
00458 nlvl = (integer) (log((real) (*n) / (real) (*smlsiz + 1)) / log(2.f)) + 1;
00459
00460 smlszp = *smlsiz + 1;
00461
00462 u = 1;
00463 vt = *smlsiz * *n + 1;
00464 difl = vt + smlszp * *n;
00465 difr = difl + nlvl * *n;
00466 z__ = difr + (nlvl * *n << 1);
00467 c__ = z__ + nlvl * *n;
00468 s = c__ + *n;
00469 poles = s + *n;
00470 givnum = poles + (nlvl << 1) * *n;
00471 nrwork = givnum + (nlvl << 1) * *n;
00472 bx = 1;
00473
00474 irwrb = nrwork;
00475 irwib = irwrb + *smlsiz * *nrhs;
00476 irwb = irwib + *smlsiz * *nrhs;
00477
00478 sizei = *n + 1;
00479 k = sizei + *n;
00480 givptr = k + *n;
00481 perm = givptr + *n;
00482 givcol = perm + nlvl * *n;
00483 iwk = givcol + (nlvl * *n << 1);
00484
00485 st = 1;
00486 sqre = 0;
00487 icmpq1 = 1;
00488 icmpq2 = 0;
00489 nsub = 0;
00490
00491 i__1 = *n;
00492 for (i__ = 1; i__ <= i__1; ++i__) {
00493 if ((r__1 = d__[i__], dabs(r__1)) < eps) {
00494 d__[i__] = r_sign(&eps, &d__[i__]);
00495 }
00496
00497 }
00498
00499 i__1 = nm1;
00500 for (i__ = 1; i__ <= i__1; ++i__) {
00501 if ((r__1 = e[i__], dabs(r__1)) < eps || i__ == nm1) {
00502 ++nsub;
00503 iwork[nsub] = st;
00504
00505
00506
00507
00508 if (i__ < nm1) {
00509
00510
00511
00512 nsize = i__ - st + 1;
00513 iwork[sizei + nsub - 1] = nsize;
00514 } else if ((r__1 = e[i__], dabs(r__1)) >= eps) {
00515
00516
00517
00518 nsize = *n - st + 1;
00519 iwork[sizei + nsub - 1] = nsize;
00520 } else {
00521
00522
00523
00524
00525
00526 nsize = i__ - st + 1;
00527 iwork[sizei + nsub - 1] = nsize;
00528 ++nsub;
00529 iwork[nsub] = *n;
00530 iwork[sizei + nsub - 1] = 1;
00531 ccopy_(nrhs, &b[*n + b_dim1], ldb, &work[bx + nm1], n);
00532 }
00533 st1 = st - 1;
00534 if (nsize == 1) {
00535
00536
00537
00538
00539 ccopy_(nrhs, &b[st + b_dim1], ldb, &work[bx + st1], n);
00540 } else if (nsize <= *smlsiz) {
00541
00542
00543
00544 slaset_("A", &nsize, &nsize, &c_b35, &c_b10, &rwork[vt + st1],
00545 n);
00546 slaset_("A", &nsize, &nsize, &c_b35, &c_b10, &rwork[u + st1],
00547 n);
00548 slasdq_("U", &c__0, &nsize, &nsize, &nsize, &c__0, &d__[st], &
00549 e[st], &rwork[vt + st1], n, &rwork[u + st1], n, &
00550 rwork[nrwork], &c__1, &rwork[nrwork], info)
00551 ;
00552 if (*info != 0) {
00553 return 0;
00554 }
00555
00556
00557
00558
00559
00560 j = irwb - 1;
00561 i__2 = *nrhs;
00562 for (jcol = 1; jcol <= i__2; ++jcol) {
00563 i__3 = st + nsize - 1;
00564 for (jrow = st; jrow <= i__3; ++jrow) {
00565 ++j;
00566 i__4 = jrow + jcol * b_dim1;
00567 rwork[j] = b[i__4].r;
00568
00569 }
00570
00571 }
00572 sgemm_("T", "N", &nsize, nrhs, &nsize, &c_b10, &rwork[u + st1]
00573 , n, &rwork[irwb], &nsize, &c_b35, &rwork[irwrb], &
00574 nsize);
00575 j = irwb - 1;
00576 i__2 = *nrhs;
00577 for (jcol = 1; jcol <= i__2; ++jcol) {
00578 i__3 = st + nsize - 1;
00579 for (jrow = st; jrow <= i__3; ++jrow) {
00580 ++j;
00581 rwork[j] = r_imag(&b[jrow + jcol * b_dim1]);
00582
00583 }
00584
00585 }
00586 sgemm_("T", "N", &nsize, nrhs, &nsize, &c_b10, &rwork[u + st1]
00587 , n, &rwork[irwb], &nsize, &c_b35, &rwork[irwib], &
00588 nsize);
00589 jreal = irwrb - 1;
00590 jimag = irwib - 1;
00591 i__2 = *nrhs;
00592 for (jcol = 1; jcol <= i__2; ++jcol) {
00593 i__3 = st + nsize - 1;
00594 for (jrow = st; jrow <= i__3; ++jrow) {
00595 ++jreal;
00596 ++jimag;
00597 i__4 = jrow + jcol * b_dim1;
00598 i__5 = jreal;
00599 i__6 = jimag;
00600 q__1.r = rwork[i__5], q__1.i = rwork[i__6];
00601 b[i__4].r = q__1.r, b[i__4].i = q__1.i;
00602
00603 }
00604
00605 }
00606
00607 clacpy_("A", &nsize, nrhs, &b[st + b_dim1], ldb, &work[bx +
00608 st1], n);
00609 } else {
00610
00611
00612
00613 slasda_(&icmpq1, smlsiz, &nsize, &sqre, &d__[st], &e[st], &
00614 rwork[u + st1], n, &rwork[vt + st1], &iwork[k + st1],
00615 &rwork[difl + st1], &rwork[difr + st1], &rwork[z__ +
00616 st1], &rwork[poles + st1], &iwork[givptr + st1], &
00617 iwork[givcol + st1], n, &iwork[perm + st1], &rwork[
00618 givnum + st1], &rwork[c__ + st1], &rwork[s + st1], &
00619 rwork[nrwork], &iwork[iwk], info);
00620 if (*info != 0) {
00621 return 0;
00622 }
00623 bxst = bx + st1;
00624 clalsa_(&icmpq2, smlsiz, &nsize, nrhs, &b[st + b_dim1], ldb, &
00625 work[bxst], n, &rwork[u + st1], n, &rwork[vt + st1], &
00626 iwork[k + st1], &rwork[difl + st1], &rwork[difr + st1]
00627 , &rwork[z__ + st1], &rwork[poles + st1], &iwork[
00628 givptr + st1], &iwork[givcol + st1], n, &iwork[perm +
00629 st1], &rwork[givnum + st1], &rwork[c__ + st1], &rwork[
00630 s + st1], &rwork[nrwork], &iwork[iwk], info);
00631 if (*info != 0) {
00632 return 0;
00633 }
00634 }
00635 st = i__ + 1;
00636 }
00637
00638 }
00639
00640
00641
00642 tol = rcnd * (r__1 = d__[isamax_(n, &d__[1], &c__1)], dabs(r__1));
00643
00644 i__1 = *n;
00645 for (i__ = 1; i__ <= i__1; ++i__) {
00646
00647
00648
00649
00650 if ((r__1 = d__[i__], dabs(r__1)) <= tol) {
00651 claset_("A", &c__1, nrhs, &c_b1, &c_b1, &work[bx + i__ - 1], n);
00652 } else {
00653 ++(*rank);
00654 clascl_("G", &c__0, &c__0, &d__[i__], &c_b10, &c__1, nrhs, &work[
00655 bx + i__ - 1], n, info);
00656 }
00657 d__[i__] = (r__1 = d__[i__], dabs(r__1));
00658
00659 }
00660
00661
00662
00663 icmpq2 = 1;
00664 i__1 = nsub;
00665 for (i__ = 1; i__ <= i__1; ++i__) {
00666 st = iwork[i__];
00667 st1 = st - 1;
00668 nsize = iwork[sizei + i__ - 1];
00669 bxst = bx + st1;
00670 if (nsize == 1) {
00671 ccopy_(nrhs, &work[bxst], n, &b[st + b_dim1], ldb);
00672 } else if (nsize <= *smlsiz) {
00673
00674
00675
00676
00677
00678
00679
00680
00681 j = bxst - *n - 1;
00682 jreal = irwb - 1;
00683 i__2 = *nrhs;
00684 for (jcol = 1; jcol <= i__2; ++jcol) {
00685 j += *n;
00686 i__3 = nsize;
00687 for (jrow = 1; jrow <= i__3; ++jrow) {
00688 ++jreal;
00689 i__4 = j + jrow;
00690 rwork[jreal] = work[i__4].r;
00691
00692 }
00693
00694 }
00695 sgemm_("T", "N", &nsize, nrhs, &nsize, &c_b10, &rwork[vt + st1],
00696 n, &rwork[irwb], &nsize, &c_b35, &rwork[irwrb], &nsize);
00697 j = bxst - *n - 1;
00698 jimag = irwb - 1;
00699 i__2 = *nrhs;
00700 for (jcol = 1; jcol <= i__2; ++jcol) {
00701 j += *n;
00702 i__3 = nsize;
00703 for (jrow = 1; jrow <= i__3; ++jrow) {
00704 ++jimag;
00705 rwork[jimag] = r_imag(&work[j + jrow]);
00706
00707 }
00708
00709 }
00710 sgemm_("T", "N", &nsize, nrhs, &nsize, &c_b10, &rwork[vt + st1],
00711 n, &rwork[irwb], &nsize, &c_b35, &rwork[irwib], &nsize);
00712 jreal = irwrb - 1;
00713 jimag = irwib - 1;
00714 i__2 = *nrhs;
00715 for (jcol = 1; jcol <= i__2; ++jcol) {
00716 i__3 = st + nsize - 1;
00717 for (jrow = st; jrow <= i__3; ++jrow) {
00718 ++jreal;
00719 ++jimag;
00720 i__4 = jrow + jcol * b_dim1;
00721 i__5 = jreal;
00722 i__6 = jimag;
00723 q__1.r = rwork[i__5], q__1.i = rwork[i__6];
00724 b[i__4].r = q__1.r, b[i__4].i = q__1.i;
00725
00726 }
00727
00728 }
00729 } else {
00730 clalsa_(&icmpq2, smlsiz, &nsize, nrhs, &work[bxst], n, &b[st +
00731 b_dim1], ldb, &rwork[u + st1], n, &rwork[vt + st1], &
00732 iwork[k + st1], &rwork[difl + st1], &rwork[difr + st1], &
00733 rwork[z__ + st1], &rwork[poles + st1], &iwork[givptr +
00734 st1], &iwork[givcol + st1], n, &iwork[perm + st1], &rwork[
00735 givnum + st1], &rwork[c__ + st1], &rwork[s + st1], &rwork[
00736 nrwork], &iwork[iwk], info);
00737 if (*info != 0) {
00738 return 0;
00739 }
00740 }
00741
00742 }
00743
00744
00745
00746 slascl_("G", &c__0, &c__0, &c_b10, &orgnrm, n, &c__1, &d__[1], n, info);
00747 slasrt_("D", n, &d__[1], info);
00748 clascl_("G", &c__0, &c__0, &orgnrm, &c_b10, n, nrhs, &b[b_offset], ldb,
00749 info);
00750
00751 return 0;
00752
00753
00754
00755 }