cgelsx.c
Go to the documentation of this file.
00001 /* cgelsx.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Table of constant values */
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 /* Subroutine */ 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     /* System generated locals */
00029     integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
00030     complex q__1;
00031 
00032     /* Builtin functions */
00033     double c_abs(complex *);
00034     void r_cnjg(complex *, complex *);
00035 
00036     /* Local variables */
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 /* Subroutine */ 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 /* Subroutine */ 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 /* Subroutine */ int claset_(char *, integer *, integer *, complex 
00056             *, complex *, complex *, integer *), xerbla_(char *, 
00057             integer *);
00058     real bignum;
00059     extern /* Subroutine */ int clatzm_(char *, integer *, integer *, complex 
00060             *, integer *, complex *, complex *, complex *, integer *, complex 
00061             *);
00062     real sminpr;
00063     extern /* Subroutine */ int ctzrqf_(integer *, integer *, complex *, 
00064             integer *, complex *, integer *);
00065     real smaxpr, smlnum;
00066 
00067 
00068 /*  -- LAPACK driver routine (version 3.2) -- */
00069 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00070 /*     November 2006 */
00071 
00072 /*     .. Scalar Arguments .. */
00073 /*     .. */
00074 /*     .. Array Arguments .. */
00075 /*     .. */
00076 
00077 /*  Purpose */
00078 /*  ======= */
00079 
00080 /*  This routine is deprecated and has been replaced by routine CGELSY. */
00081 
00082 /*  CGELSX computes the minimum-norm solution to a complex linear least */
00083 /*  squares problem: */
00084 /*      minimize || A * X - B || */
00085 /*  using a complete orthogonal factorization of A.  A is an M-by-N */
00086 /*  matrix which may be rank-deficient. */
00087 
00088 /*  Several right hand side vectors b and solution vectors x can be */
00089 /*  handled in a single call; they are stored as the columns of the */
00090 /*  M-by-NRHS right hand side matrix B and the N-by-NRHS solution */
00091 /*  matrix X. */
00092 
00093 /*  The routine first computes a QR factorization with column pivoting: */
00094 /*      A * P = Q * [ R11 R12 ] */
00095 /*                  [  0  R22 ] */
00096 /*  with R11 defined as the largest leading submatrix whose estimated */
00097 /*  condition number is less than 1/RCOND.  The order of R11, RANK, */
00098 /*  is the effective rank of A. */
00099 
00100 /*  Then, R22 is considered to be negligible, and R12 is annihilated */
00101 /*  by unitary transformations from the right, arriving at the */
00102 /*  complete orthogonal factorization: */
00103 /*     A * P = Q * [ T11 0 ] * Z */
00104 /*                 [  0  0 ] */
00105 /*  The minimum-norm solution is then */
00106 /*     X = P * Z' [ inv(T11)*Q1'*B ] */
00107 /*                [        0       ] */
00108 /*  where Q1 consists of the first RANK columns of Q. */
00109 
00110 /*  Arguments */
00111 /*  ========= */
00112 
00113 /*  M       (input) INTEGER */
00114 /*          The number of rows of the matrix A.  M >= 0. */
00115 
00116 /*  N       (input) INTEGER */
00117 /*          The number of columns of the matrix A.  N >= 0. */
00118 
00119 /*  NRHS    (input) INTEGER */
00120 /*          The number of right hand sides, i.e., the number of */
00121 /*          columns of matrices B and X. NRHS >= 0. */
00122 
00123 /*  A       (input/output) COMPLEX array, dimension (LDA,N) */
00124 /*          On entry, the M-by-N matrix A. */
00125 /*          On exit, A has been overwritten by details of its */
00126 /*          complete orthogonal factorization. */
00127 
00128 /*  LDA     (input) INTEGER */
00129 /*          The leading dimension of the array A.  LDA >= max(1,M). */
00130 
00131 /*  B       (input/output) COMPLEX array, dimension (LDB,NRHS) */
00132 /*          On entry, the M-by-NRHS right hand side matrix B. */
00133 /*          On exit, the N-by-NRHS solution matrix X. */
00134 /*          If m >= n and RANK = n, the residual sum-of-squares for */
00135 /*          the solution in the i-th column is given by the sum of */
00136 /*          squares of elements N+1:M in that column. */
00137 
00138 /*  LDB     (input) INTEGER */
00139 /*          The leading dimension of the array B. LDB >= max(1,M,N). */
00140 
00141 /*  JPVT    (input/output) INTEGER array, dimension (N) */
00142 /*          On entry, if JPVT(i) .ne. 0, the i-th column of A is an */
00143 /*          initial column, otherwise it is a free column.  Before */
00144 /*          the QR factorization of A, all initial columns are */
00145 /*          permuted to the leading positions; only the remaining */
00146 /*          free columns are moved as a result of column pivoting */
00147 /*          during the factorization. */
00148 /*          On exit, if JPVT(i) = k, then the i-th column of A*P */
00149 /*          was the k-th column of A. */
00150 
00151 /*  RCOND   (input) REAL */
00152 /*          RCOND is used to determine the effective rank of A, which */
00153 /*          is defined as the order of the largest leading triangular */
00154 /*          submatrix R11 in the QR factorization with pivoting of A, */
00155 /*          whose estimated condition number < 1/RCOND. */
00156 
00157 /*  RANK    (output) INTEGER */
00158 /*          The effective rank of A, i.e., the order of the submatrix */
00159 /*          R11.  This is the same as the order of the submatrix T11 */
00160 /*          in the complete orthogonal factorization of A. */
00161 
00162 /*  WORK    (workspace) COMPLEX array, dimension */
00163 /*                      (min(M,N) + max( N, 2*min(M,N)+NRHS )), */
00164 
00165 /*  RWORK   (workspace) REAL array, dimension (2*N) */
00166 
00167 /*  INFO    (output) INTEGER */
00168 /*          = 0:  successful exit */
00169 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
00170 
00171 /*  ===================================================================== */
00172 
00173 /*     .. Parameters .. */
00174 /*     .. */
00175 /*     .. Local Scalars .. */
00176 /*     .. */
00177 /*     .. External Subroutines .. */
00178 /*     .. */
00179 /*     .. External Functions .. */
00180 /*     .. */
00181 /*     .. Intrinsic Functions .. */
00182 /*     .. */
00183 /*     .. Executable Statements .. */
00184 
00185     /* Parameter adjustments */
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     /* Function Body */
00197     mn = min(*m,*n);
00198     ismin = mn + 1;
00199     ismax = (mn << 1) + 1;
00200 
00201 /*     Test the input arguments. */
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 /* if(complicated condition) */ {
00213 /* Computing MAX */
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 /*     Quick return if possible */
00227 
00228 /* Computing MIN */
00229     i__1 = min(*m,*n);
00230     if (min(i__1,*nrhs) == 0) {
00231         *rank = 0;
00232         return 0;
00233     }
00234 
00235 /*     Get machine parameters */
00236 
00237     smlnum = slamch_("S") / slamch_("P");
00238     bignum = 1.f / smlnum;
00239     slabad_(&smlnum, &bignum);
00240 
00241 /*     Scale A, B if max elements outside range [SMLNUM,BIGNUM] */
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 /*        Scale matrix norm up to SMLNUM */
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 /*        Scale matrix norm down to BIGNUM */
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 /*        Matrix all zero. Return zero solution. */
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 /*        Scale matrix norm up to SMLNUM */
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 /*        Scale matrix norm down to BIGNUM */
00281 
00282         clascl_("G", &c__0, &c__0, &bnrm, &bignum, m, nrhs, &b[b_offset], ldb, 
00283                  info);
00284         ibscl = 2;
00285     }
00286 
00287 /*     Compute QR factorization with column pivoting of A: */
00288 /*        A * P = Q * R */
00289 
00290     cgeqpf_(m, n, &a[a_offset], lda, &jpvt[1], &work[1], &work[mn + 1], &
00291             rwork[1], info);
00292 
00293 /*     complex workspace MN+N. Real workspace 2*N. Details of Householder */
00294 /*     rotations stored in WORK(1:MN). */
00295 
00296 /*     Determine RANK using incremental condition estimation */
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 /* L20: */
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 /*     Logically partition R = [ R11 R12 ] */
00348 /*                             [  0  R22 ] */
00349 /*     where R11 = R(1:RANK,1:RANK) */
00350 
00351 /*     [R11,R12] = [ T11, 0 ] * Y */
00352 
00353     if (*rank < *n) {
00354         ctzrqf_(rank, n, &a[a_offset], lda, &work[mn + 1], info);
00355     }
00356 
00357 /*     Details of Householder rotations stored in WORK(MN+1:2*MN) */
00358 
00359 /*     B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS) */
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 /*     workspace NRHS */
00365 
00366 /*      B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS) */
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 /* L30: */
00378         }
00379 /* L40: */
00380     }
00381 
00382 /*     B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS) */
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 /* L50: */
00393         }
00394     }
00395 
00396 /*     workspace NRHS */
00397 
00398 /*     B(1:N,1:NRHS) := P * B(1:N,1:NRHS) */
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 /* L60: */
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 /* L80: */
00437         }
00438 /* L90: */
00439     }
00440 
00441 /*     Undo scaling */
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 /*     End of CGELSX */
00467 
00468 } /* cgelsx_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:26