dgelss.c
Go to the documentation of this file.
00001 /* dgelss.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 integer c__6 = 6;
00019 static integer c_n1 = -1;
00020 static integer c__1 = 1;
00021 static integer c__0 = 0;
00022 static doublereal c_b74 = 0.;
00023 static doublereal c_b108 = 1.;
00024 
00025 /* Subroutine */ int dgelss_(integer *m, integer *n, integer *nrhs, 
00026         doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *
00027         s, doublereal *rcond, integer *rank, doublereal *work, integer *lwork, 
00028          integer *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     doublereal d__1;
00033 
00034     /* Local variables */
00035     integer i__, bl, ie, il, mm;
00036     doublereal eps, thr, anrm, bnrm;
00037     integer itau;
00038     doublereal vdum[1];
00039     extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *, 
00040             integer *, doublereal *, doublereal *, integer *, doublereal *, 
00041             integer *, doublereal *, doublereal *, integer *);
00042     integer iascl, ibscl;
00043     extern /* Subroutine */ int dgemv_(char *, integer *, integer *, 
00044             doublereal *, doublereal *, integer *, doublereal *, integer *, 
00045             doublereal *, doublereal *, integer *), drscl_(integer *, 
00046             doublereal *, doublereal *, integer *);
00047     integer chunk;
00048     doublereal sfmin;
00049     integer minmn;
00050     extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 
00051             doublereal *, integer *);
00052     integer maxmn, itaup, itauq, mnthr, iwork;
00053     extern /* Subroutine */ int dlabad_(doublereal *, doublereal *), dgebrd_(
00054             integer *, integer *, doublereal *, integer *, doublereal *, 
00055             doublereal *, doublereal *, doublereal *, doublereal *, integer *, 
00056              integer *);
00057     extern doublereal dlamch_(char *), dlange_(char *, integer *, 
00058             integer *, doublereal *, integer *, doublereal *);
00059     integer bdspac;
00060     extern /* Subroutine */ int dgelqf_(integer *, integer *, doublereal *, 
00061             integer *, doublereal *, doublereal *, integer *, integer *), 
00062             dlascl_(char *, integer *, integer *, doublereal *, doublereal *, 
00063             integer *, integer *, doublereal *, integer *, integer *),
00064              dgeqrf_(integer *, integer *, doublereal *, integer *, 
00065             doublereal *, doublereal *, integer *, integer *), dlacpy_(char *, 
00066              integer *, integer *, doublereal *, integer *, doublereal *, 
00067             integer *), dlaset_(char *, integer *, integer *, 
00068             doublereal *, doublereal *, doublereal *, integer *), 
00069             xerbla_(char *, integer *), dbdsqr_(char *, integer *, 
00070             integer *, integer *, integer *, doublereal *, doublereal *, 
00071             doublereal *, integer *, doublereal *, integer *, doublereal *, 
00072             integer *, doublereal *, integer *), dorgbr_(char *, 
00073             integer *, integer *, integer *, doublereal *, integer *, 
00074             doublereal *, doublereal *, integer *, integer *);
00075     doublereal bignum;
00076     extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 
00077             integer *, integer *);
00078     extern /* Subroutine */ int dormbr_(char *, char *, char *, integer *, 
00079             integer *, integer *, doublereal *, integer *, doublereal *, 
00080             doublereal *, integer *, doublereal *, integer *, integer *), dormlq_(char *, char *, integer *, 
00081             integer *, integer *, doublereal *, integer *, doublereal *, 
00082             doublereal *, integer *, doublereal *, integer *, integer *);
00083     integer ldwork;
00084     extern /* Subroutine */ int dormqr_(char *, char *, integer *, integer *, 
00085             integer *, doublereal *, integer *, doublereal *, doublereal *, 
00086             integer *, doublereal *, integer *, integer *);
00087     integer minwrk, maxwrk;
00088     doublereal smlnum;
00089     logical lquery;
00090 
00091 
00092 /*  -- LAPACK driver routine (version 3.2) -- */
00093 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00094 /*     November 2006 */
00095 
00096 /*     .. Scalar Arguments .. */
00097 /*     .. */
00098 /*     .. Array Arguments .. */
00099 /*     .. */
00100 
00101 /*  Purpose */
00102 /*  ======= */
00103 
00104 /*  DGELSS computes the minimum norm solution to a real linear least */
00105 /*  squares problem: */
00106 
00107 /*  Minimize 2-norm(| b - A*x |). */
00108 
00109 /*  using the singular value decomposition (SVD) of A. A is an M-by-N */
00110 /*  matrix which may be rank-deficient. */
00111 
00112 /*  Several right hand side vectors b and solution vectors x can be */
00113 /*  handled in a single call; they are stored as the columns of the */
00114 /*  M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix */
00115 /*  X. */
00116 
00117 /*  The effective rank of A is determined by treating as zero those */
00118 /*  singular values which are less than RCOND times the largest singular */
00119 /*  value. */
00120 
00121 /*  Arguments */
00122 /*  ========= */
00123 
00124 /*  M       (input) INTEGER */
00125 /*          The number of rows of the matrix A. M >= 0. */
00126 
00127 /*  N       (input) INTEGER */
00128 /*          The number of columns of the matrix A. N >= 0. */
00129 
00130 /*  NRHS    (input) INTEGER */
00131 /*          The number of right hand sides, i.e., the number of columns */
00132 /*          of the matrices B and X. NRHS >= 0. */
00133 
00134 /*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
00135 /*          On entry, the M-by-N matrix A. */
00136 /*          On exit, the first min(m,n) rows of A are overwritten with */
00137 /*          its right singular vectors, stored rowwise. */
00138 
00139 /*  LDA     (input) INTEGER */
00140 /*          The leading dimension of the array A.  LDA >= max(1,M). */
00141 
00142 /*  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) */
00143 /*          On entry, the M-by-NRHS right hand side matrix B. */
00144 /*          On exit, B is overwritten by the N-by-NRHS solution */
00145 /*          matrix X.  If m >= n and RANK = n, the residual */
00146 /*          sum-of-squares for the solution in the i-th column is given */
00147 /*          by the sum of squares of elements n+1:m in that column. */
00148 
00149 /*  LDB     (input) INTEGER */
00150 /*          The leading dimension of the array B. LDB >= max(1,max(M,N)). */
00151 
00152 /*  S       (output) DOUBLE PRECISION array, dimension (min(M,N)) */
00153 /*          The singular values of A in decreasing order. */
00154 /*          The condition number of A in the 2-norm = S(1)/S(min(m,n)). */
00155 
00156 /*  RCOND   (input) DOUBLE PRECISION */
00157 /*          RCOND is used to determine the effective rank of A. */
00158 /*          Singular values S(i) <= RCOND*S(1) are treated as zero. */
00159 /*          If RCOND < 0, machine precision is used instead. */
00160 
00161 /*  RANK    (output) INTEGER */
00162 /*          The effective rank of A, i.e., the number of singular values */
00163 /*          which are greater than RCOND*S(1). */
00164 
00165 /*  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
00166 /*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
00167 
00168 /*  LWORK   (input) INTEGER */
00169 /*          The dimension of the array WORK. LWORK >= 1, and also: */
00170 /*          LWORK >= 3*min(M,N) + max( 2*min(M,N), max(M,N), NRHS ) */
00171 /*          For good performance, LWORK should generally be larger. */
00172 
00173 /*          If LWORK = -1, then a workspace query is assumed; the routine */
00174 /*          only calculates the optimal size of the WORK array, returns */
00175 /*          this value as the first entry of the WORK array, and no error */
00176 /*          message related to LWORK is issued by XERBLA. */
00177 
00178 /*  INFO    (output) INTEGER */
00179 /*          = 0:  successful exit */
00180 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
00181 /*          > 0:  the algorithm for computing the SVD failed to converge; */
00182 /*                if INFO = i, i off-diagonal elements of an intermediate */
00183 /*                bidiagonal form did not converge to zero. */
00184 
00185 /*  ===================================================================== */
00186 
00187 /*     .. Parameters .. */
00188 /*     .. */
00189 /*     .. Local Scalars .. */
00190 /*     .. */
00191 /*     .. Local Arrays .. */
00192 /*     .. */
00193 /*     .. External Subroutines .. */
00194 /*     .. */
00195 /*     .. External Functions .. */
00196 /*     .. */
00197 /*     .. Intrinsic Functions .. */
00198 /*     .. */
00199 /*     .. Executable Statements .. */
00200 
00201 /*     Test the input arguments */
00202 
00203     /* Parameter adjustments */
00204     a_dim1 = *lda;
00205     a_offset = 1 + a_dim1;
00206     a -= a_offset;
00207     b_dim1 = *ldb;
00208     b_offset = 1 + b_dim1;
00209     b -= b_offset;
00210     --s;
00211     --work;
00212 
00213     /* Function Body */
00214     *info = 0;
00215     minmn = min(*m,*n);
00216     maxmn = max(*m,*n);
00217     lquery = *lwork == -1;
00218     if (*m < 0) {
00219         *info = -1;
00220     } else if (*n < 0) {
00221         *info = -2;
00222     } else if (*nrhs < 0) {
00223         *info = -3;
00224     } else if (*lda < max(1,*m)) {
00225         *info = -5;
00226     } else if (*ldb < max(1,maxmn)) {
00227         *info = -7;
00228     }
00229 
00230 /*     Compute workspace */
00231 /*      (Note: Comments in the code beginning "Workspace:" describe the */
00232 /*       minimal amount of workspace needed at that point in the code, */
00233 /*       as well as the preferred amount for good performance. */
00234 /*       NB refers to the optimal block size for the immediately */
00235 /*       following subroutine, as returned by ILAENV.) */
00236 
00237     if (*info == 0) {
00238         minwrk = 1;
00239         maxwrk = 1;
00240         if (minmn > 0) {
00241             mm = *m;
00242             mnthr = ilaenv_(&c__6, "DGELSS", " ", m, n, nrhs, &c_n1);
00243             if (*m >= *n && *m >= mnthr) {
00244 
00245 /*              Path 1a - overdetermined, with many more rows than */
00246 /*                        columns */
00247 
00248                 mm = *n;
00249 /* Computing MAX */
00250                 i__1 = maxwrk, i__2 = *n + *n * ilaenv_(&c__1, "DGEQRF", 
00251                         " ", m, n, &c_n1, &c_n1);
00252                 maxwrk = max(i__1,i__2);
00253 /* Computing MAX */
00254                 i__1 = maxwrk, i__2 = *n + *nrhs * ilaenv_(&c__1, "DORMQR", 
00255                         "LT", m, nrhs, n, &c_n1);
00256                 maxwrk = max(i__1,i__2);
00257             }
00258             if (*m >= *n) {
00259 
00260 /*              Path 1 - overdetermined or exactly determined */
00261 
00262 /*              Compute workspace needed for DBDSQR */
00263 
00264 /* Computing MAX */
00265                 i__1 = 1, i__2 = *n * 5;
00266                 bdspac = max(i__1,i__2);
00267 /* Computing MAX */
00268                 i__1 = maxwrk, i__2 = *n * 3 + (mm + *n) * ilaenv_(&c__1, 
00269                         "DGEBRD", " ", &mm, n, &c_n1, &c_n1);
00270                 maxwrk = max(i__1,i__2);
00271 /* Computing MAX */
00272                 i__1 = maxwrk, i__2 = *n * 3 + *nrhs * ilaenv_(&c__1, "DORMBR"
00273 , "QLT", &mm, nrhs, n, &c_n1);
00274                 maxwrk = max(i__1,i__2);
00275 /* Computing MAX */
00276                 i__1 = maxwrk, i__2 = *n * 3 + (*n - 1) * ilaenv_(&c__1, 
00277                         "DORGBR", "P", n, n, n, &c_n1);
00278                 maxwrk = max(i__1,i__2);
00279                 maxwrk = max(maxwrk,bdspac);
00280 /* Computing MAX */
00281                 i__1 = maxwrk, i__2 = *n * *nrhs;
00282                 maxwrk = max(i__1,i__2);
00283 /* Computing MAX */
00284                 i__1 = *n * 3 + mm, i__2 = *n * 3 + *nrhs, i__1 = max(i__1,
00285                         i__2);
00286                 minwrk = max(i__1,bdspac);
00287                 maxwrk = max(minwrk,maxwrk);
00288             }
00289             if (*n > *m) {
00290 
00291 /*              Compute workspace needed for DBDSQR */
00292 
00293 /* Computing MAX */
00294                 i__1 = 1, i__2 = *m * 5;
00295                 bdspac = max(i__1,i__2);
00296 /* Computing MAX */
00297                 i__1 = *m * 3 + *nrhs, i__2 = *m * 3 + *n, i__1 = max(i__1,
00298                         i__2);
00299                 minwrk = max(i__1,bdspac);
00300                 if (*n >= mnthr) {
00301 
00302 /*                 Path 2a - underdetermined, with many more columns */
00303 /*                 than rows */
00304 
00305                     maxwrk = *m + *m * ilaenv_(&c__1, "DGELQF", " ", m, n, &
00306                             c_n1, &c_n1);
00307 /* Computing MAX */
00308                     i__1 = maxwrk, i__2 = *m * *m + (*m << 2) + (*m << 1) * 
00309                             ilaenv_(&c__1, "DGEBRD", " ", m, m, &c_n1, &c_n1);
00310                     maxwrk = max(i__1,i__2);
00311 /* Computing MAX */
00312                     i__1 = maxwrk, i__2 = *m * *m + (*m << 2) + *nrhs * 
00313                             ilaenv_(&c__1, "DORMBR", "QLT", m, nrhs, m, &c_n1);
00314                     maxwrk = max(i__1,i__2);
00315 /* Computing MAX */
00316                     i__1 = maxwrk, i__2 = *m * *m + (*m << 2) + (*m - 1) * 
00317                             ilaenv_(&c__1, "DORGBR", "P", m, m, m, &c_n1);
00318                     maxwrk = max(i__1,i__2);
00319 /* Computing MAX */
00320                     i__1 = maxwrk, i__2 = *m * *m + *m + bdspac;
00321                     maxwrk = max(i__1,i__2);
00322                     if (*nrhs > 1) {
00323 /* Computing MAX */
00324                         i__1 = maxwrk, i__2 = *m * *m + *m + *m * *nrhs;
00325                         maxwrk = max(i__1,i__2);
00326                     } else {
00327 /* Computing MAX */
00328                         i__1 = maxwrk, i__2 = *m * *m + (*m << 1);
00329                         maxwrk = max(i__1,i__2);
00330                     }
00331 /* Computing MAX */
00332                     i__1 = maxwrk, i__2 = *m + *nrhs * ilaenv_(&c__1, "DORMLQ"
00333 , "LT", n, nrhs, m, &c_n1);
00334                     maxwrk = max(i__1,i__2);
00335                 } else {
00336 
00337 /*                 Path 2 - underdetermined */
00338 
00339                     maxwrk = *m * 3 + (*n + *m) * ilaenv_(&c__1, "DGEBRD", 
00340                             " ", m, n, &c_n1, &c_n1);
00341 /* Computing MAX */
00342                     i__1 = maxwrk, i__2 = *m * 3 + *nrhs * ilaenv_(&c__1, 
00343                             "DORMBR", "QLT", m, nrhs, m, &c_n1);
00344                     maxwrk = max(i__1,i__2);
00345 /* Computing MAX */
00346                     i__1 = maxwrk, i__2 = *m * 3 + *m * ilaenv_(&c__1, "DORG"
00347                             "BR", "P", m, n, m, &c_n1);
00348                     maxwrk = max(i__1,i__2);
00349                     maxwrk = max(maxwrk,bdspac);
00350 /* Computing MAX */
00351                     i__1 = maxwrk, i__2 = *n * *nrhs;
00352                     maxwrk = max(i__1,i__2);
00353                 }
00354             }
00355             maxwrk = max(minwrk,maxwrk);
00356         }
00357         work[1] = (doublereal) maxwrk;
00358 
00359         if (*lwork < minwrk && ! lquery) {
00360             *info = -12;
00361         }
00362     }
00363 
00364     if (*info != 0) {
00365         i__1 = -(*info);
00366         xerbla_("DGELSS", &i__1);
00367         return 0;
00368     } else if (lquery) {
00369         return 0;
00370     }
00371 
00372 /*     Quick return if possible */
00373 
00374     if (*m == 0 || *n == 0) {
00375         *rank = 0;
00376         return 0;
00377     }
00378 
00379 /*     Get machine parameters */
00380 
00381     eps = dlamch_("P");
00382     sfmin = dlamch_("S");
00383     smlnum = sfmin / eps;
00384     bignum = 1. / smlnum;
00385     dlabad_(&smlnum, &bignum);
00386 
00387 /*     Scale A if max element outside range [SMLNUM,BIGNUM] */
00388 
00389     anrm = dlange_("M", m, n, &a[a_offset], lda, &work[1]);
00390     iascl = 0;
00391     if (anrm > 0. && anrm < smlnum) {
00392 
00393 /*        Scale matrix norm up to SMLNUM */
00394 
00395         dlascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda, 
00396                 info);
00397         iascl = 1;
00398     } else if (anrm > bignum) {
00399 
00400 /*        Scale matrix norm down to BIGNUM */
00401 
00402         dlascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, 
00403                 info);
00404         iascl = 2;
00405     } else if (anrm == 0.) {
00406 
00407 /*        Matrix all zero. Return zero solution. */
00408 
00409         i__1 = max(*m,*n);
00410         dlaset_("F", &i__1, nrhs, &c_b74, &c_b74, &b[b_offset], ldb);
00411         dlaset_("F", &minmn, &c__1, &c_b74, &c_b74, &s[1], &c__1);
00412         *rank = 0;
00413         goto L70;
00414     }
00415 
00416 /*     Scale B if max element outside range [SMLNUM,BIGNUM] */
00417 
00418     bnrm = dlange_("M", m, nrhs, &b[b_offset], ldb, &work[1]);
00419     ibscl = 0;
00420     if (bnrm > 0. && bnrm < smlnum) {
00421 
00422 /*        Scale matrix norm up to SMLNUM */
00423 
00424         dlascl_("G", &c__0, &c__0, &bnrm, &smlnum, m, nrhs, &b[b_offset], ldb, 
00425                  info);
00426         ibscl = 1;
00427     } else if (bnrm > bignum) {
00428 
00429 /*        Scale matrix norm down to BIGNUM */
00430 
00431         dlascl_("G", &c__0, &c__0, &bnrm, &bignum, m, nrhs, &b[b_offset], ldb, 
00432                  info);
00433         ibscl = 2;
00434     }
00435 
00436 /*     Overdetermined case */
00437 
00438     if (*m >= *n) {
00439 
00440 /*        Path 1 - overdetermined or exactly determined */
00441 
00442         mm = *m;
00443         if (*m >= mnthr) {
00444 
00445 /*           Path 1a - overdetermined, with many more rows than columns */
00446 
00447             mm = *n;
00448             itau = 1;
00449             iwork = itau + *n;
00450 
00451 /*           Compute A=Q*R */
00452 /*           (Workspace: need 2*N, prefer N+N*NB) */
00453 
00454             i__1 = *lwork - iwork + 1;
00455             dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__1, 
00456                      info);
00457 
00458 /*           Multiply B by transpose(Q) */
00459 /*           (Workspace: need N+NRHS, prefer N+NRHS*NB) */
00460 
00461             i__1 = *lwork - iwork + 1;
00462             dormqr_("L", "T", m, nrhs, n, &a[a_offset], lda, &work[itau], &b[
00463                     b_offset], ldb, &work[iwork], &i__1, info);
00464 
00465 /*           Zero out below R */
00466 
00467             if (*n > 1) {
00468                 i__1 = *n - 1;
00469                 i__2 = *n - 1;
00470                 dlaset_("L", &i__1, &i__2, &c_b74, &c_b74, &a[a_dim1 + 2], 
00471                         lda);
00472             }
00473         }
00474 
00475         ie = 1;
00476         itauq = ie + *n;
00477         itaup = itauq + *n;
00478         iwork = itaup + *n;
00479 
00480 /*        Bidiagonalize R in A */
00481 /*        (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB) */
00482 
00483         i__1 = *lwork - iwork + 1;
00484         dgebrd_(&mm, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
00485                 work[itaup], &work[iwork], &i__1, info);
00486 
00487 /*        Multiply B by transpose of left bidiagonalizing vectors of R */
00488 /*        (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB) */
00489 
00490         i__1 = *lwork - iwork + 1;
00491         dormbr_("Q", "L", "T", &mm, nrhs, n, &a[a_offset], lda, &work[itauq], 
00492                 &b[b_offset], ldb, &work[iwork], &i__1, info);
00493 
00494 /*        Generate right bidiagonalizing vectors of R in A */
00495 /*        (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) */
00496 
00497         i__1 = *lwork - iwork + 1;
00498         dorgbr_("P", n, n, n, &a[a_offset], lda, &work[itaup], &work[iwork], &
00499                 i__1, info);
00500         iwork = ie + *n;
00501 
00502 /*        Perform bidiagonal QR iteration */
00503 /*          multiply B by transpose of left singular vectors */
00504 /*          compute right singular vectors in A */
00505 /*        (Workspace: need BDSPAC) */
00506 
00507         dbdsqr_("U", n, n, &c__0, nrhs, &s[1], &work[ie], &a[a_offset], lda, 
00508                 vdum, &c__1, &b[b_offset], ldb, &work[iwork], info)
00509                 ;
00510         if (*info != 0) {
00511             goto L70;
00512         }
00513 
00514 /*        Multiply B by reciprocals of singular values */
00515 
00516 /* Computing MAX */
00517         d__1 = *rcond * s[1];
00518         thr = max(d__1,sfmin);
00519         if (*rcond < 0.) {
00520 /* Computing MAX */
00521             d__1 = eps * s[1];
00522             thr = max(d__1,sfmin);
00523         }
00524         *rank = 0;
00525         i__1 = *n;
00526         for (i__ = 1; i__ <= i__1; ++i__) {
00527             if (s[i__] > thr) {
00528                 drscl_(nrhs, &s[i__], &b[i__ + b_dim1], ldb);
00529                 ++(*rank);
00530             } else {
00531                 dlaset_("F", &c__1, nrhs, &c_b74, &c_b74, &b[i__ + b_dim1], 
00532                         ldb);
00533             }
00534 /* L10: */
00535         }
00536 
00537 /*        Multiply B by right singular vectors */
00538 /*        (Workspace: need N, prefer N*NRHS) */
00539 
00540         if (*lwork >= *ldb * *nrhs && *nrhs > 1) {
00541             dgemm_("T", "N", n, nrhs, n, &c_b108, &a[a_offset], lda, &b[
00542                     b_offset], ldb, &c_b74, &work[1], ldb);
00543             dlacpy_("G", n, nrhs, &work[1], ldb, &b[b_offset], ldb)
00544                     ;
00545         } else if (*nrhs > 1) {
00546             chunk = *lwork / *n;
00547             i__1 = *nrhs;
00548             i__2 = chunk;
00549             for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
00550 /* Computing MIN */
00551                 i__3 = *nrhs - i__ + 1;
00552                 bl = min(i__3,chunk);
00553                 dgemm_("T", "N", n, &bl, n, &c_b108, &a[a_offset], lda, &b[
00554                         i__ * b_dim1 + 1], ldb, &c_b74, &work[1], n);
00555                 dlacpy_("G", n, &bl, &work[1], n, &b[i__ * b_dim1 + 1], ldb);
00556 /* L20: */
00557             }
00558         } else {
00559             dgemv_("T", n, n, &c_b108, &a[a_offset], lda, &b[b_offset], &c__1, 
00560                      &c_b74, &work[1], &c__1);
00561             dcopy_(n, &work[1], &c__1, &b[b_offset], &c__1);
00562         }
00563 
00564     } else /* if(complicated condition) */ {
00565 /* Computing MAX */
00566         i__2 = *m, i__1 = (*m << 1) - 4, i__2 = max(i__2,i__1), i__2 = max(
00567                 i__2,*nrhs), i__1 = *n - *m * 3;
00568         if (*n >= mnthr && *lwork >= (*m << 2) + *m * *m + max(i__2,i__1)) {
00569 
00570 /*        Path 2a - underdetermined, with many more columns than rows */
00571 /*        and sufficient workspace for an efficient algorithm */
00572 
00573             ldwork = *m;
00574 /* Computing MAX */
00575 /* Computing MAX */
00576             i__3 = *m, i__4 = (*m << 1) - 4, i__3 = max(i__3,i__4), i__3 = 
00577                     max(i__3,*nrhs), i__4 = *n - *m * 3;
00578             i__2 = (*m << 2) + *m * *lda + max(i__3,i__4), i__1 = *m * *lda + 
00579                     *m + *m * *nrhs;
00580             if (*lwork >= max(i__2,i__1)) {
00581                 ldwork = *lda;
00582             }
00583             itau = 1;
00584             iwork = *m + 1;
00585 
00586 /*        Compute A=L*Q */
00587 /*        (Workspace: need 2*M, prefer M+M*NB) */
00588 
00589             i__2 = *lwork - iwork + 1;
00590             dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, 
00591                      info);
00592             il = iwork;
00593 
00594 /*        Copy L to WORK(IL), zeroing out above it */
00595 
00596             dlacpy_("L", m, m, &a[a_offset], lda, &work[il], &ldwork);
00597             i__2 = *m - 1;
00598             i__1 = *m - 1;
00599             dlaset_("U", &i__2, &i__1, &c_b74, &c_b74, &work[il + ldwork], &
00600                     ldwork);
00601             ie = il + ldwork * *m;
00602             itauq = ie + *m;
00603             itaup = itauq + *m;
00604             iwork = itaup + *m;
00605 
00606 /*        Bidiagonalize L in WORK(IL) */
00607 /*        (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB) */
00608 
00609             i__2 = *lwork - iwork + 1;
00610             dgebrd_(m, m, &work[il], &ldwork, &s[1], &work[ie], &work[itauq], 
00611                     &work[itaup], &work[iwork], &i__2, info);
00612 
00613 /*        Multiply B by transpose of left bidiagonalizing vectors of L */
00614 /*        (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB) */
00615 
00616             i__2 = *lwork - iwork + 1;
00617             dormbr_("Q", "L", "T", m, nrhs, m, &work[il], &ldwork, &work[
00618                     itauq], &b[b_offset], ldb, &work[iwork], &i__2, info);
00619 
00620 /*        Generate right bidiagonalizing vectors of R in WORK(IL) */
00621 /*        (Workspace: need M*M+5*M-1, prefer M*M+4*M+(M-1)*NB) */
00622 
00623             i__2 = *lwork - iwork + 1;
00624             dorgbr_("P", m, m, m, &work[il], &ldwork, &work[itaup], &work[
00625                     iwork], &i__2, info);
00626             iwork = ie + *m;
00627 
00628 /*        Perform bidiagonal QR iteration, */
00629 /*           computing right singular vectors of L in WORK(IL) and */
00630 /*           multiplying B by transpose of left singular vectors */
00631 /*        (Workspace: need M*M+M+BDSPAC) */
00632 
00633             dbdsqr_("U", m, m, &c__0, nrhs, &s[1], &work[ie], &work[il], &
00634                     ldwork, &a[a_offset], lda, &b[b_offset], ldb, &work[iwork]
00635 , info);
00636             if (*info != 0) {
00637                 goto L70;
00638             }
00639 
00640 /*        Multiply B by reciprocals of singular values */
00641 
00642 /* Computing MAX */
00643             d__1 = *rcond * s[1];
00644             thr = max(d__1,sfmin);
00645             if (*rcond < 0.) {
00646 /* Computing MAX */
00647                 d__1 = eps * s[1];
00648                 thr = max(d__1,sfmin);
00649             }
00650             *rank = 0;
00651             i__2 = *m;
00652             for (i__ = 1; i__ <= i__2; ++i__) {
00653                 if (s[i__] > thr) {
00654                     drscl_(nrhs, &s[i__], &b[i__ + b_dim1], ldb);
00655                     ++(*rank);
00656                 } else {
00657                     dlaset_("F", &c__1, nrhs, &c_b74, &c_b74, &b[i__ + b_dim1]
00658 , ldb);
00659                 }
00660 /* L30: */
00661             }
00662             iwork = ie;
00663 
00664 /*        Multiply B by right singular vectors of L in WORK(IL) */
00665 /*        (Workspace: need M*M+2*M, prefer M*M+M+M*NRHS) */
00666 
00667             if (*lwork >= *ldb * *nrhs + iwork - 1 && *nrhs > 1) {
00668                 dgemm_("T", "N", m, nrhs, m, &c_b108, &work[il], &ldwork, &b[
00669                         b_offset], ldb, &c_b74, &work[iwork], ldb);
00670                 dlacpy_("G", m, nrhs, &work[iwork], ldb, &b[b_offset], ldb);
00671             } else if (*nrhs > 1) {
00672                 chunk = (*lwork - iwork + 1) / *m;
00673                 i__2 = *nrhs;
00674                 i__1 = chunk;
00675                 for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += 
00676                         i__1) {
00677 /* Computing MIN */
00678                     i__3 = *nrhs - i__ + 1;
00679                     bl = min(i__3,chunk);
00680                     dgemm_("T", "N", m, &bl, m, &c_b108, &work[il], &ldwork, &
00681                             b[i__ * b_dim1 + 1], ldb, &c_b74, &work[iwork], m);
00682                     dlacpy_("G", m, &bl, &work[iwork], m, &b[i__ * b_dim1 + 1]
00683 , ldb);
00684 /* L40: */
00685                 }
00686             } else {
00687                 dgemv_("T", m, m, &c_b108, &work[il], &ldwork, &b[b_dim1 + 1], 
00688                          &c__1, &c_b74, &work[iwork], &c__1);
00689                 dcopy_(m, &work[iwork], &c__1, &b[b_dim1 + 1], &c__1);
00690             }
00691 
00692 /*        Zero out below first M rows of B */
00693 
00694             i__1 = *n - *m;
00695             dlaset_("F", &i__1, nrhs, &c_b74, &c_b74, &b[*m + 1 + b_dim1], 
00696                     ldb);
00697             iwork = itau + *m;
00698 
00699 /*        Multiply transpose(Q) by B */
00700 /*        (Workspace: need M+NRHS, prefer M+NRHS*NB) */
00701 
00702             i__1 = *lwork - iwork + 1;
00703             dormlq_("L", "T", n, nrhs, m, &a[a_offset], lda, &work[itau], &b[
00704                     b_offset], ldb, &work[iwork], &i__1, info);
00705 
00706         } else {
00707 
00708 /*        Path 2 - remaining underdetermined cases */
00709 
00710             ie = 1;
00711             itauq = ie + *m;
00712             itaup = itauq + *m;
00713             iwork = itaup + *m;
00714 
00715 /*        Bidiagonalize A */
00716 /*        (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) */
00717 
00718             i__1 = *lwork - iwork + 1;
00719             dgebrd_(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
00720                     work[itaup], &work[iwork], &i__1, info);
00721 
00722 /*        Multiply B by transpose of left bidiagonalizing vectors */
00723 /*        (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB) */
00724 
00725             i__1 = *lwork - iwork + 1;
00726             dormbr_("Q", "L", "T", m, nrhs, n, &a[a_offset], lda, &work[itauq]
00727 , &b[b_offset], ldb, &work[iwork], &i__1, info);
00728 
00729 /*        Generate right bidiagonalizing vectors in A */
00730 /*        (Workspace: need 4*M, prefer 3*M+M*NB) */
00731 
00732             i__1 = *lwork - iwork + 1;
00733             dorgbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &work[
00734                     iwork], &i__1, info);
00735             iwork = ie + *m;
00736 
00737 /*        Perform bidiagonal QR iteration, */
00738 /*           computing right singular vectors of A in A and */
00739 /*           multiplying B by transpose of left singular vectors */
00740 /*        (Workspace: need BDSPAC) */
00741 
00742             dbdsqr_("L", m, n, &c__0, nrhs, &s[1], &work[ie], &a[a_offset], 
00743                     lda, vdum, &c__1, &b[b_offset], ldb, &work[iwork], info);
00744             if (*info != 0) {
00745                 goto L70;
00746             }
00747 
00748 /*        Multiply B by reciprocals of singular values */
00749 
00750 /* Computing MAX */
00751             d__1 = *rcond * s[1];
00752             thr = max(d__1,sfmin);
00753             if (*rcond < 0.) {
00754 /* Computing MAX */
00755                 d__1 = eps * s[1];
00756                 thr = max(d__1,sfmin);
00757             }
00758             *rank = 0;
00759             i__1 = *m;
00760             for (i__ = 1; i__ <= i__1; ++i__) {
00761                 if (s[i__] > thr) {
00762                     drscl_(nrhs, &s[i__], &b[i__ + b_dim1], ldb);
00763                     ++(*rank);
00764                 } else {
00765                     dlaset_("F", &c__1, nrhs, &c_b74, &c_b74, &b[i__ + b_dim1]
00766 , ldb);
00767                 }
00768 /* L50: */
00769             }
00770 
00771 /*        Multiply B by right singular vectors of A */
00772 /*        (Workspace: need N, prefer N*NRHS) */
00773 
00774             if (*lwork >= *ldb * *nrhs && *nrhs > 1) {
00775                 dgemm_("T", "N", n, nrhs, m, &c_b108, &a[a_offset], lda, &b[
00776                         b_offset], ldb, &c_b74, &work[1], ldb);
00777                 dlacpy_("F", n, nrhs, &work[1], ldb, &b[b_offset], ldb);
00778             } else if (*nrhs > 1) {
00779                 chunk = *lwork / *n;
00780                 i__1 = *nrhs;
00781                 i__2 = chunk;
00782                 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += 
00783                         i__2) {
00784 /* Computing MIN */
00785                     i__3 = *nrhs - i__ + 1;
00786                     bl = min(i__3,chunk);
00787                     dgemm_("T", "N", n, &bl, m, &c_b108, &a[a_offset], lda, &
00788                             b[i__ * b_dim1 + 1], ldb, &c_b74, &work[1], n);
00789                     dlacpy_("F", n, &bl, &work[1], n, &b[i__ * b_dim1 + 1], 
00790                             ldb);
00791 /* L60: */
00792                 }
00793             } else {
00794                 dgemv_("T", m, n, &c_b108, &a[a_offset], lda, &b[b_offset], &
00795                         c__1, &c_b74, &work[1], &c__1);
00796                 dcopy_(n, &work[1], &c__1, &b[b_offset], &c__1);
00797             }
00798         }
00799     }
00800 
00801 /*     Undo scaling */
00802 
00803     if (iascl == 1) {
00804         dlascl_("G", &c__0, &c__0, &anrm, &smlnum, n, nrhs, &b[b_offset], ldb, 
00805                  info);
00806         dlascl_("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], &
00807                 minmn, info);
00808     } else if (iascl == 2) {
00809         dlascl_("G", &c__0, &c__0, &anrm, &bignum, n, nrhs, &b[b_offset], ldb, 
00810                  info);
00811         dlascl_("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], &
00812                 minmn, info);
00813     }
00814     if (ibscl == 1) {
00815         dlascl_("G", &c__0, &c__0, &smlnum, &bnrm, n, nrhs, &b[b_offset], ldb, 
00816                  info);
00817     } else if (ibscl == 2) {
00818         dlascl_("G", &c__0, &c__0, &bignum, &bnrm, n, nrhs, &b[b_offset], ldb, 
00819                  info);
00820     }
00821 
00822 L70:
00823     work[1] = (doublereal) maxwrk;
00824     return 0;
00825 
00826 /*     End of DGELSS */
00827 
00828 } /* dgelss_ */


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