cgelsy.c
Go to the documentation of this file.
00001 /* cgelsy.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__1 = 1;
00021 static integer c_n1 = -1;
00022 static integer c__0 = 0;
00023 static integer c__2 = 2;
00024 
00025 /* Subroutine */ int cgelsy_(integer *m, integer *n, integer *nrhs, complex *
00026         a, integer *lda, complex *b, integer *ldb, integer *jpvt, real *rcond, 
00027          integer *rank, complex *work, integer *lwork, real *rwork, integer *
00028         info)
00029 {
00030     /* System generated locals */
00031     integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4;
00032     real r__1, r__2;
00033     complex q__1;
00034 
00035     /* Builtin functions */
00036     double c_abs(complex *);
00037 
00038     /* Local variables */
00039     integer i__, j;
00040     complex c1, c2, s1, s2;
00041     integer nb, mn, nb1, nb2, nb3, nb4;
00042     real anrm, bnrm, smin, smax;
00043     integer iascl, ibscl;
00044     extern /* Subroutine */ int ccopy_(integer *, complex *, integer *, 
00045             complex *, integer *);
00046     integer ismin, ismax;
00047     extern /* Subroutine */ int ctrsm_(char *, char *, char *, char *, 
00048             integer *, integer *, complex *, complex *, integer *, complex *, 
00049             integer *), claic1_(integer *, 
00050             integer *, complex *, real *, complex *, complex *, real *, 
00051             complex *, complex *);
00052     real wsize;
00053     extern /* Subroutine */ int cgeqp3_(integer *, integer *, complex *, 
00054             integer *, integer *, complex *, complex *, integer *, real *, 
00055             integer *), slabad_(real *, real *);
00056     extern doublereal clange_(char *, integer *, integer *, complex *, 
00057             integer *, real *);
00058     extern /* Subroutine */ int clascl_(char *, integer *, integer *, real *, 
00059             real *, integer *, integer *, complex *, integer *, integer *);
00060     extern doublereal slamch_(char *);
00061     extern /* Subroutine */ int claset_(char *, integer *, integer *, complex 
00062             *, complex *, complex *, integer *), xerbla_(char *, 
00063             integer *);
00064     extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 
00065             integer *, integer *);
00066     real bignum;
00067     extern /* Subroutine */ int cunmqr_(char *, char *, integer *, integer *, 
00068             integer *, complex *, integer *, complex *, complex *, integer *, 
00069             complex *, integer *, integer *);
00070     real sminpr, smaxpr, smlnum;
00071     extern /* Subroutine */ int cunmrz_(char *, char *, integer *, integer *, 
00072             integer *, integer *, complex *, integer *, complex *, complex *, 
00073             integer *, complex *, integer *, integer *);
00074     integer lwkopt;
00075     logical lquery;
00076     extern /* Subroutine */ int ctzrzf_(integer *, integer *, complex *, 
00077             integer *, complex *, complex *, integer *, integer *);
00078 
00079 
00080 /*  -- LAPACK driver routine (version 3.2) -- */
00081 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00082 /*     November 2006 */
00083 
00084 /*     .. Scalar Arguments .. */
00085 /*     .. */
00086 /*     .. Array Arguments .. */
00087 /*     .. */
00088 
00089 /*  Purpose */
00090 /*  ======= */
00091 
00092 /*  CGELSY computes the minimum-norm solution to a complex linear least */
00093 /*  squares problem: */
00094 /*      minimize || A * X - B || */
00095 /*  using a complete orthogonal factorization of A.  A is an M-by-N */
00096 /*  matrix which may be rank-deficient. */
00097 
00098 /*  Several right hand side vectors b and solution vectors x can be */
00099 /*  handled in a single call; they are stored as the columns of the */
00100 /*  M-by-NRHS right hand side matrix B and the N-by-NRHS solution */
00101 /*  matrix X. */
00102 
00103 /*  The routine first computes a QR factorization with column pivoting: */
00104 /*      A * P = Q * [ R11 R12 ] */
00105 /*                  [  0  R22 ] */
00106 /*  with R11 defined as the largest leading submatrix whose estimated */
00107 /*  condition number is less than 1/RCOND.  The order of R11, RANK, */
00108 /*  is the effective rank of A. */
00109 
00110 /*  Then, R22 is considered to be negligible, and R12 is annihilated */
00111 /*  by unitary transformations from the right, arriving at the */
00112 /*  complete orthogonal factorization: */
00113 /*     A * P = Q * [ T11 0 ] * Z */
00114 /*                 [  0  0 ] */
00115 /*  The minimum-norm solution is then */
00116 /*     X = P * Z' [ inv(T11)*Q1'*B ] */
00117 /*                [        0       ] */
00118 /*  where Q1 consists of the first RANK columns of Q. */
00119 
00120 /*  This routine is basically identical to the original xGELSX except */
00121 /*  three differences: */
00122 /*    o The permutation of matrix B (the right hand side) is faster and */
00123 /*      more simple. */
00124 /*    o The call to the subroutine xGEQPF has been substituted by the */
00125 /*      the call to the subroutine xGEQP3. This subroutine is a Blas-3 */
00126 /*      version of the QR factorization with column pivoting. */
00127 /*    o Matrix B (the right hand side) is updated with Blas-3. */
00128 
00129 /*  Arguments */
00130 /*  ========= */
00131 
00132 /*  M       (input) INTEGER */
00133 /*          The number of rows of the matrix A.  M >= 0. */
00134 
00135 /*  N       (input) INTEGER */
00136 /*          The number of columns of the matrix A.  N >= 0. */
00137 
00138 /*  NRHS    (input) INTEGER */
00139 /*          The number of right hand sides, i.e., the number of */
00140 /*          columns of matrices B and X. NRHS >= 0. */
00141 
00142 /*  A       (input/output) COMPLEX array, dimension (LDA,N) */
00143 /*          On entry, the M-by-N matrix A. */
00144 /*          On exit, A has been overwritten by details of its */
00145 /*          complete orthogonal factorization. */
00146 
00147 /*  LDA     (input) INTEGER */
00148 /*          The leading dimension of the array A.  LDA >= max(1,M). */
00149 
00150 /*  B       (input/output) COMPLEX array, dimension (LDB,NRHS) */
00151 /*          On entry, the M-by-NRHS right hand side matrix B. */
00152 /*          On exit, the N-by-NRHS solution matrix X. */
00153 
00154 /*  LDB     (input) INTEGER */
00155 /*          The leading dimension of the array B. LDB >= max(1,M,N). */
00156 
00157 /*  JPVT    (input/output) INTEGER array, dimension (N) */
00158 /*          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted */
00159 /*          to the front of AP, otherwise column i is a free column. */
00160 /*          On exit, if JPVT(i) = k, then the i-th column of A*P */
00161 /*          was the k-th column of A. */
00162 
00163 /*  RCOND   (input) REAL */
00164 /*          RCOND is used to determine the effective rank of A, which */
00165 /*          is defined as the order of the largest leading triangular */
00166 /*          submatrix R11 in the QR factorization with pivoting of A, */
00167 /*          whose estimated condition number < 1/RCOND. */
00168 
00169 /*  RANK    (output) INTEGER */
00170 /*          The effective rank of A, i.e., the order of the submatrix */
00171 /*          R11.  This is the same as the order of the submatrix T11 */
00172 /*          in the complete orthogonal factorization of A. */
00173 
00174 /*  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) */
00175 /*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
00176 
00177 /*  LWORK   (input) INTEGER */
00178 /*          The dimension of the array WORK. */
00179 /*          The unblocked strategy requires that: */
00180 /*            LWORK >= MN + MAX( 2*MN, N+1, MN+NRHS ) */
00181 /*          where MN = min(M,N). */
00182 /*          The block algorithm requires that: */
00183 /*            LWORK >= MN + MAX( 2*MN, NB*(N+1), MN+MN*NB, MN+NB*NRHS ) */
00184 /*          where NB is an upper bound on the blocksize returned */
00185 /*          by ILAENV for the routines CGEQP3, CTZRZF, CTZRQF, CUNMQR, */
00186 /*          and CUNMRZ. */
00187 
00188 /*          If LWORK = -1, then a workspace query is assumed; the routine */
00189 /*          only calculates the optimal size of the WORK array, returns */
00190 /*          this value as the first entry of the WORK array, and no error */
00191 /*          message related to LWORK is issued by XERBLA. */
00192 
00193 /*  RWORK   (workspace) REAL array, dimension (2*N) */
00194 
00195 /*  INFO    (output) INTEGER */
00196 /*          = 0: successful exit */
00197 /*          < 0: if INFO = -i, the i-th argument had an illegal value */
00198 
00199 /*  Further Details */
00200 /*  =============== */
00201 
00202 /*  Based on contributions by */
00203 /*    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA */
00204 /*    E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain */
00205 /*    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain */
00206 
00207 /*  ===================================================================== */
00208 
00209 /*     .. Parameters .. */
00210 /*     .. */
00211 /*     .. Local Scalars .. */
00212 /*     .. */
00213 /*     .. External Subroutines .. */
00214 /*     .. */
00215 /*     .. External Functions .. */
00216 /*     .. */
00217 /*     .. Intrinsic Functions .. */
00218 /*     .. */
00219 /*     .. Executable Statements .. */
00220 
00221     /* Parameter adjustments */
00222     a_dim1 = *lda;
00223     a_offset = 1 + a_dim1;
00224     a -= a_offset;
00225     b_dim1 = *ldb;
00226     b_offset = 1 + b_dim1;
00227     b -= b_offset;
00228     --jpvt;
00229     --work;
00230     --rwork;
00231 
00232     /* Function Body */
00233     mn = min(*m,*n);
00234     ismin = mn + 1;
00235     ismax = (mn << 1) + 1;
00236 
00237 /*     Test the input arguments. */
00238 
00239     *info = 0;
00240     nb1 = ilaenv_(&c__1, "CGEQRF", " ", m, n, &c_n1, &c_n1);
00241     nb2 = ilaenv_(&c__1, "CGERQF", " ", m, n, &c_n1, &c_n1);
00242     nb3 = ilaenv_(&c__1, "CUNMQR", " ", m, n, nrhs, &c_n1);
00243     nb4 = ilaenv_(&c__1, "CUNMRQ", " ", m, n, nrhs, &c_n1);
00244 /* Computing MAX */
00245     i__1 = max(nb1,nb2), i__1 = max(i__1,nb3);
00246     nb = max(i__1,nb4);
00247 /* Computing MAX */
00248     i__1 = 1, i__2 = mn + (*n << 1) + nb * (*n + 1), i__1 = max(i__1,i__2), 
00249             i__2 = (mn << 1) + nb * *nrhs;
00250     lwkopt = max(i__1,i__2);
00251     q__1.r = (real) lwkopt, q__1.i = 0.f;
00252     work[1].r = q__1.r, work[1].i = q__1.i;
00253     lquery = *lwork == -1;
00254     if (*m < 0) {
00255         *info = -1;
00256     } else if (*n < 0) {
00257         *info = -2;
00258     } else if (*nrhs < 0) {
00259         *info = -3;
00260     } else if (*lda < max(1,*m)) {
00261         *info = -5;
00262     } else /* if(complicated condition) */ {
00263 /* Computing MAX */
00264         i__1 = max(1,*m);
00265         if (*ldb < max(i__1,*n)) {
00266             *info = -7;
00267         } else /* if(complicated condition) */ {
00268 /* Computing MAX */
00269             i__1 = mn << 1, i__2 = *n + 1, i__1 = max(i__1,i__2), i__2 = mn + 
00270                     *nrhs;
00271             if (*lwork < mn + max(i__1,i__2) && ! lquery) {
00272                 *info = -12;
00273             }
00274         }
00275     }
00276 
00277     if (*info != 0) {
00278         i__1 = -(*info);
00279         xerbla_("CGELSY", &i__1);
00280         return 0;
00281     } else if (lquery) {
00282         return 0;
00283     }
00284 
00285 /*     Quick return if possible */
00286 
00287 /* Computing MIN */
00288     i__1 = min(*m,*n);
00289     if (min(i__1,*nrhs) == 0) {
00290         *rank = 0;
00291         return 0;
00292     }
00293 
00294 /*     Get machine parameters */
00295 
00296     smlnum = slamch_("S") / slamch_("P");
00297     bignum = 1.f / smlnum;
00298     slabad_(&smlnum, &bignum);
00299 
00300 /*     Scale A, B if max entries outside range [SMLNUM,BIGNUM] */
00301 
00302     anrm = clange_("M", m, n, &a[a_offset], lda, &rwork[1]);
00303     iascl = 0;
00304     if (anrm > 0.f && anrm < smlnum) {
00305 
00306 /*        Scale matrix norm up to SMLNUM */
00307 
00308         clascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda, 
00309                 info);
00310         iascl = 1;
00311     } else if (anrm > bignum) {
00312 
00313 /*        Scale matrix norm down to BIGNUM */
00314 
00315         clascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, 
00316                 info);
00317         iascl = 2;
00318     } else if (anrm == 0.f) {
00319 
00320 /*        Matrix all zero. Return zero solution. */
00321 
00322         i__1 = max(*m,*n);
00323         claset_("F", &i__1, nrhs, &c_b1, &c_b1, &b[b_offset], ldb);
00324         *rank = 0;
00325         goto L70;
00326     }
00327 
00328     bnrm = clange_("M", m, nrhs, &b[b_offset], ldb, &rwork[1]);
00329     ibscl = 0;
00330     if (bnrm > 0.f && bnrm < smlnum) {
00331 
00332 /*        Scale matrix norm up to SMLNUM */
00333 
00334         clascl_("G", &c__0, &c__0, &bnrm, &smlnum, m, nrhs, &b[b_offset], ldb, 
00335                  info);
00336         ibscl = 1;
00337     } else if (bnrm > bignum) {
00338 
00339 /*        Scale matrix norm down to BIGNUM */
00340 
00341         clascl_("G", &c__0, &c__0, &bnrm, &bignum, m, nrhs, &b[b_offset], ldb, 
00342                  info);
00343         ibscl = 2;
00344     }
00345 
00346 /*     Compute QR factorization with column pivoting of A: */
00347 /*        A * P = Q * R */
00348 
00349     i__1 = *lwork - mn;
00350     cgeqp3_(m, n, &a[a_offset], lda, &jpvt[1], &work[1], &work[mn + 1], &i__1, 
00351              &rwork[1], info);
00352     i__1 = mn + 1;
00353     wsize = mn + work[i__1].r;
00354 
00355 /*     complex workspace: MN+NB*(N+1). real workspace 2*N. */
00356 /*     Details of Householder rotations stored in WORK(1:MN). */
00357 
00358 /*     Determine RANK using incremental condition estimation */
00359 
00360     i__1 = ismin;
00361     work[i__1].r = 1.f, work[i__1].i = 0.f;
00362     i__1 = ismax;
00363     work[i__1].r = 1.f, work[i__1].i = 0.f;
00364     smax = c_abs(&a[a_dim1 + 1]);
00365     smin = smax;
00366     if (c_abs(&a[a_dim1 + 1]) == 0.f) {
00367         *rank = 0;
00368         i__1 = max(*m,*n);
00369         claset_("F", &i__1, nrhs, &c_b1, &c_b1, &b[b_offset], ldb);
00370         goto L70;
00371     } else {
00372         *rank = 1;
00373     }
00374 
00375 L10:
00376     if (*rank < mn) {
00377         i__ = *rank + 1;
00378         claic1_(&c__2, rank, &work[ismin], &smin, &a[i__ * a_dim1 + 1], &a[
00379                 i__ + i__ * a_dim1], &sminpr, &s1, &c1);
00380         claic1_(&c__1, rank, &work[ismax], &smax, &a[i__ * a_dim1 + 1], &a[
00381                 i__ + i__ * a_dim1], &smaxpr, &s2, &c2);
00382 
00383         if (smaxpr * *rcond <= sminpr) {
00384             i__1 = *rank;
00385             for (i__ = 1; i__ <= i__1; ++i__) {
00386                 i__2 = ismin + i__ - 1;
00387                 i__3 = ismin + i__ - 1;
00388                 q__1.r = s1.r * work[i__3].r - s1.i * work[i__3].i, q__1.i = 
00389                         s1.r * work[i__3].i + s1.i * work[i__3].r;
00390                 work[i__2].r = q__1.r, work[i__2].i = q__1.i;
00391                 i__2 = ismax + i__ - 1;
00392                 i__3 = ismax + i__ - 1;
00393                 q__1.r = s2.r * work[i__3].r - s2.i * work[i__3].i, q__1.i = 
00394                         s2.r * work[i__3].i + s2.i * work[i__3].r;
00395                 work[i__2].r = q__1.r, work[i__2].i = q__1.i;
00396 /* L20: */
00397             }
00398             i__1 = ismin + *rank;
00399             work[i__1].r = c1.r, work[i__1].i = c1.i;
00400             i__1 = ismax + *rank;
00401             work[i__1].r = c2.r, work[i__1].i = c2.i;
00402             smin = sminpr;
00403             smax = smaxpr;
00404             ++(*rank);
00405             goto L10;
00406         }
00407     }
00408 
00409 /*     complex workspace: 3*MN. */
00410 
00411 /*     Logically partition R = [ R11 R12 ] */
00412 /*                             [  0  R22 ] */
00413 /*     where R11 = R(1:RANK,1:RANK) */
00414 
00415 /*     [R11,R12] = [ T11, 0 ] * Y */
00416 
00417     if (*rank < *n) {
00418         i__1 = *lwork - (mn << 1);
00419         ctzrzf_(rank, n, &a[a_offset], lda, &work[mn + 1], &work[(mn << 1) + 
00420                 1], &i__1, info);
00421     }
00422 
00423 /*     complex workspace: 2*MN. */
00424 /*     Details of Householder rotations stored in WORK(MN+1:2*MN) */
00425 
00426 /*     B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS) */
00427 
00428     i__1 = *lwork - (mn << 1);
00429     cunmqr_("Left", "Conjugate transpose", m, nrhs, &mn, &a[a_offset], lda, &
00430             work[1], &b[b_offset], ldb, &work[(mn << 1) + 1], &i__1, info);
00431 /* Computing MAX */
00432     i__1 = (mn << 1) + 1;
00433     r__1 = wsize, r__2 = (mn << 1) + work[i__1].r;
00434     wsize = dmax(r__1,r__2);
00435 
00436 /*     complex workspace: 2*MN+NB*NRHS. */
00437 
00438 /*     B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS) */
00439 
00440     ctrsm_("Left", "Upper", "No transpose", "Non-unit", rank, nrhs, &c_b2, &a[
00441             a_offset], lda, &b[b_offset], ldb);
00442 
00443     i__1 = *nrhs;
00444     for (j = 1; j <= i__1; ++j) {
00445         i__2 = *n;
00446         for (i__ = *rank + 1; i__ <= i__2; ++i__) {
00447             i__3 = i__ + j * b_dim1;
00448             b[i__3].r = 0.f, b[i__3].i = 0.f;
00449 /* L30: */
00450         }
00451 /* L40: */
00452     }
00453 
00454 /*     B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS) */
00455 
00456     if (*rank < *n) {
00457         i__1 = *n - *rank;
00458         i__2 = *lwork - (mn << 1);
00459         cunmrz_("Left", "Conjugate transpose", n, nrhs, rank, &i__1, &a[
00460                 a_offset], lda, &work[mn + 1], &b[b_offset], ldb, &work[(mn <<
00461                  1) + 1], &i__2, info);
00462     }
00463 
00464 /*     complex workspace: 2*MN+NRHS. */
00465 
00466 /*     B(1:N,1:NRHS) := P * B(1:N,1:NRHS) */
00467 
00468     i__1 = *nrhs;
00469     for (j = 1; j <= i__1; ++j) {
00470         i__2 = *n;
00471         for (i__ = 1; i__ <= i__2; ++i__) {
00472             i__3 = jpvt[i__];
00473             i__4 = i__ + j * b_dim1;
00474             work[i__3].r = b[i__4].r, work[i__3].i = b[i__4].i;
00475 /* L50: */
00476         }
00477         ccopy_(n, &work[1], &c__1, &b[j * b_dim1 + 1], &c__1);
00478 /* L60: */
00479     }
00480 
00481 /*     complex workspace: N. */
00482 
00483 /*     Undo scaling */
00484 
00485     if (iascl == 1) {
00486         clascl_("G", &c__0, &c__0, &anrm, &smlnum, n, nrhs, &b[b_offset], ldb, 
00487                  info);
00488         clascl_("U", &c__0, &c__0, &smlnum, &anrm, rank, rank, &a[a_offset], 
00489                 lda, info);
00490     } else if (iascl == 2) {
00491         clascl_("G", &c__0, &c__0, &anrm, &bignum, n, nrhs, &b[b_offset], ldb, 
00492                  info);
00493         clascl_("U", &c__0, &c__0, &bignum, &anrm, rank, rank, &a[a_offset], 
00494                 lda, info);
00495     }
00496     if (ibscl == 1) {
00497         clascl_("G", &c__0, &c__0, &smlnum, &bnrm, n, nrhs, &b[b_offset], ldb, 
00498                  info);
00499     } else if (ibscl == 2) {
00500         clascl_("G", &c__0, &c__0, &bignum, &bnrm, n, nrhs, &b[b_offset], ldb, 
00501                  info);
00502     }
00503 
00504 L70:
00505     q__1.r = (real) lwkopt, q__1.i = 0.f;
00506     work[1].r = q__1.r, work[1].i = q__1.i;
00507 
00508     return 0;
00509 
00510 /*     End of CGELSY */
00511 
00512 } /* cgelsy_ */


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