dgelss.c

提供矩阵类的函数库

	    t2 = dsecnd_();
	    lstime_1.timng[gelqf - 1] += t2 - t1;
	    il = iwork;

/*        Copy L to WORK(IL), zeroing out above it */

	    dlacpy_("L", m, m, &a[a_offset], lda, &work[il], &ldwork);
	    i__2 = *m - 1;
	    i__1 = *m - 1;
	    dlaset_("U", &i__2, &i__1, &c_b74, &c_b74, &work[il + ldwork], &
		    ldwork);
	    ie = il + ldwork * *m;
	    itauq = ie + *m;
	    itaup = itauq + *m;
	    iwork = itaup + *m;

/*        Bidiagonalize L in WORK(IL)   
          (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB) */

	    nb = ilaenv_(&c__1, "DGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)6, 
		    (ftnlen)1);
	    lstime_1.opcnt[gebrd - 1] += dopla_("DGEBRD", m, m, &c__0, &c__0, 
		    &nb);
	    t1 = dsecnd_();
	    i__2 = *lwork - iwork + 1;
	    dgebrd_(m, m, &work[il], &ldwork, &s[1], &work[ie], &work[itauq], 
		    &work[itaup], &work[iwork], &i__2, info);
	    t2 = dsecnd_();
	    lstime_1.timng[gebrd - 1] += t2 - t1;

/*        Multiply B by transpose of left bidiagonalizing vectors of L   
          (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB) */

	    nb = ilaenv_(&c__1, "DORMBR", "QLT", m, nrhs, m, &c_n1, (ftnlen)6,
		     (ftnlen)3);
	    lstime_1.opcnt[ormbr - 1] += dopla2_("DORMBR", "QLT", m, nrhs, m, 
		    &c__0, &nb);
	    t1 = dsecnd_();
	    i__2 = *lwork - iwork + 1;
	    dormbr_("Q", "L", "T", m, nrhs, m, &work[il], &ldwork, &work[
		    itauq], &b[b_offset], ldb, &work[iwork], &i__2, info);
	    t2 = dsecnd_();
	    lstime_1.timng[ormbr - 1] += t2 - t1;

/*        Generate right bidiagonalizing vectors of R in WORK(IL)   
          (Workspace: need M*M+5*M-1, prefer M*M+4*M+(M-1)*NB) */

	    nb = ilaenv_(&c__1, "DORGBR", "P", m, m, m, &c_n1, (ftnlen)6, (
		    ftnlen)1);
	    lstime_1.opcnt[orgbr - 1] += dopla2_("DORGBR", "P", m, m, m, &
		    c__0, &nb);
	    t1 = dsecnd_();
	    i__2 = *lwork - iwork + 1;
	    dorgbr_("P", m, m, m, &work[il], &ldwork, &work[itaup], &work[
		    iwork], &i__2, info);
	    t2 = dsecnd_();
	    lstime_1.timng[orgbr - 1] += t2 - t1;
	    iwork = ie + *m;

/*        Perform bidiagonal QR iteration,   
             computing right singular vectors of L in WORK(IL) and   
             multiplying B by transpose of left singular vectors   
          (Workspace: need M*M+M+BDSPAC) */

	    latime_1.ops = 0.;
	    t1 = dsecnd_();
	    dbdsqr_("U", m, m, &c__0, nrhs, &s[1], &work[ie], &work[il], &
		    ldwork, &a[a_offset], lda, &b[b_offset], ldb, &work[iwork]
		    , info);
	    t2 = dsecnd_();
	    lstime_1.timng[bdsqr - 1] += t2 - t1;
	    lstime_1.opcnt[bdsqr - 1] += latime_1.ops;
	    if (*info != 0) {
		goto L70;
	    }

/*        Multiply B by reciprocals of singular values */

	    lstime_1.opcnt[gelss - 1] += 1.;
/* Computing MAX */
	    d__1 = *rcond * s[1];
	    thr = max(d__1,sfmin);
	    if (*rcond < 0.) {
		lstime_1.opcnt[gelss - 1] += 1.;
/* Computing MAX */
		d__1 = eps * s[1];
		thr = max(d__1,sfmin);
	    }
	    *rank = 0;
	    i__2 = *m;
	    for (i__ = 1; i__ <= i__2; ++i__) {
		if (s[i__] > thr) {
		    lstime_1.opcnt[gelss - 1] += (doublereal) (*nrhs + 3);
		    drscl_(nrhs, &s[i__], &b_ref(i__, 1), ldb);
		    ++(*rank);
		} else {
		    dlaset_("F", &c__1, nrhs, &c_b74, &c_b74, &b_ref(i__, 1), 
			    ldb);
		}
/* L30: */
	    }
	    iwork = ie;

/*        Multiply B by right singular vectors of L in WORK(IL)   
          (Workspace: need M*M+2*M, prefer M*M+M+M*NRHS) */

	    if (*lwork >= *ldb * *nrhs + iwork - 1 && *nrhs > 1) {
		lstime_1.opcnt[gemm - 1] += dopbl3_("DGEMM ", m, nrhs, m);
		t1 = dsecnd_();
		dgemm_("T", "N", m, nrhs, m, &c_b145, &work[il], &ldwork, &b[
			b_offset], ldb, &c_b74, &work[iwork], ldb);
		t2 = dsecnd_();
		lstime_1.timng[gemm - 1] += t2 - t1;
		dlacpy_("G", m, nrhs, &work[iwork], ldb, &b[b_offset], ldb);
	    } else if (*nrhs > 1) {
		chunk = (*lwork - iwork + 1) / *m;
		i__2 = *nrhs;
		i__1 = chunk;
		for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += 
			i__1) {
/* Computing MIN */
		    i__3 = *nrhs - i__ + 1;
		    bl = min(i__3,chunk);
		    lstime_1.opcnt[gemm - 1] += dopbl3_("DGEMM ", m, &bl, m);
		    t1 = dsecnd_();
		    dgemm_("T", "N", m, &bl, m, &c_b145, &work[il], &ldwork, &
			    b_ref(1, i__), ldb, &c_b74, &work[iwork], n);
		    t2 = dsecnd_();
		    lstime_1.timng[gemm - 1] += t2 - t1;
		    dlacpy_("G", m, &bl, &work[iwork], n, &b_ref(1, i__), ldb);
/* L40: */
		}
	    } else {
		lstime_1.opcnt[gemv - 1] += dopbl2_("DGEMV ", m, m, &c__0, &
			c__0);
		t1 = dsecnd_();
		dgemv_("T", m, m, &c_b145, &work[il], &ldwork, &b_ref(1, 1), &
			c__1, &c_b74, &work[iwork], &c__1);
		t2 = dsecnd_();
		lstime_1.timng[gemv - 1] += t2 - t1;
		dcopy_(m, &work[iwork], &c__1, &b_ref(1, 1), &c__1);
	    }

/*        Zero out below first M rows of B */

	    i__1 = *n - *m;
	    dlaset_("F", &i__1, nrhs, &c_b74, &c_b74, &b_ref(*m + 1, 1), ldb);
	    iwork = itau + *m;

/*        Multiply transpose(Q) by B   
          (Workspace: need M+NRHS, prefer M+NRHS*NB) */

	    nb = ilaenv_(&c__1, "DORMLQ", "LT", n, nrhs, m, &c_n1, (ftnlen)6, 
		    (ftnlen)2);
	    lstime_1.opcnt[ormlq - 1] += dopla_("DORMLQ", n, nrhs, m, &c__0, &
		    nb);
	    t1 = dsecnd_();
	    i__1 = *lwork - iwork + 1;
	    dormlq_("L", "T", n, nrhs, m, &a[a_offset], lda, &work[itau], &b[
		    b_offset], ldb, &work[iwork], &i__1, info);
	    t2 = dsecnd_();
	    lstime_1.timng[ormlq - 1] += t2 - t1;

	} else {

/*        Path 2 - remaining underdetermined cases */

	    ie = 1;
	    itauq = ie + *m;
	    itaup = itauq + *m;
	    iwork = itaup + *m;

/*        Bidiagonalize A   
          (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) */

	    nb = ilaenv_(&c__1, "DGEBRD", " ", m, n, &c_n1, &c_n1, (ftnlen)6, 
		    (ftnlen)1);
	    lstime_1.opcnt[gebrd - 1] += dopla_("DGEBRD", m, n, &c__0, &c__0, 
		    &nb);
	    t1 = dsecnd_();
	    i__1 = *lwork - iwork + 1;
	    dgebrd_(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
		    work[itaup], &work[iwork], &i__1, info);
	    t2 = dsecnd_();
	    lstime_1.timng[gebrd - 1] += t2 - t1;

/*        Multiply B by transpose of left bidiagonalizing vectors   
          (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB) */

	    nb = ilaenv_(&c__1, "DORMBR", "QLT", m, nrhs, n, &c_n1, (ftnlen)6,
		     (ftnlen)3);
	    lstime_1.opcnt[ormbr - 1] += dopla2_("DORMBR", "QLT", m, nrhs, n, 
		    &c__0, &nb);
	    t1 = dsecnd_();
	    i__1 = *lwork - iwork + 1;
	    dormbr_("Q", "L", "T", m, nrhs, n, &a[a_offset], lda, &work[itauq]
		    , &b[b_offset], ldb, &work[iwork], &i__1, info);
	    t2 = dsecnd_();
	    lstime_1.timng[ormbr - 1] += t2 - t1;

/*        Generate right bidiagonalizing vectors in A   
          (Workspace: need 4*M, prefer 3*M+M*NB) */

	    nb = ilaenv_(&c__1, "DORGBR", "P", m, n, m, &c_n1, (ftnlen)6, (
		    ftnlen)1);
	    lstime_1.opcnt[orgbr - 1] += dopla2_("DORGBR", "P", m, n, m, &
		    c__0, &nb);
	    t1 = dsecnd_();
	    i__1 = *lwork - iwork + 1;
	    dorgbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &work[
		    iwork], &i__1, info);
	    t2 = dsecnd_();
	    lstime_1.timng[orgbr - 1] += t2 - t1;
	    iwork = ie + *m;

/*        Perform bidiagonal QR iteration,   
             computing right singular vectors of A in A and   
             multiplying B by transpose of left singular vectors   
          (Workspace: need BDSPAC) */

	    latime_1.ops = 0.;
	    t1 = dsecnd_();
	    dbdsqr_("L", m, n, &c__0, nrhs, &s[1], &work[ie], &a[a_offset], 
		    lda, vdum, &c__1, &b[b_offset], ldb, &work[iwork], info);
	    t2 = dsecnd_();
	    lstime_1.timng[bdsqr - 1] += t2 - t1;
	    lstime_1.opcnt[bdsqr - 1] += latime_1.ops;
	    if (*info != 0) {
		goto L70;
	    }

/*        Multiply B by reciprocals of singular values */

	    lstime_1.opcnt[gelss - 1] += 1.;
/* Computing MAX */
	    d__1 = *rcond * s[1];
	    thr = max(d__1,sfmin);
	    if (*rcond < 0.) {
		lstime_1.opcnt[gelss - 1] += 1.;
/* Computing MAX */
		d__1 = eps * s[1];
		thr = max(d__1,sfmin);
	    }
	    *rank = 0;
	    i__1 = *m;
	    for (i__ = 1; i__ <= i__1; ++i__) {
		if (s[i__] > thr) {
		    lstime_1.opcnt[gelss - 1] += (doublereal) (*nrhs + 3);
		    drscl_(nrhs, &s[i__], &b_ref(i__, 1), ldb);
		    ++(*rank);
		} else {
		    dlaset_("F", &c__1, nrhs, &c_b74, &c_b74, &b_ref(i__, 1), 
			    ldb);
		}
/* L50: */
	    }

/*        Multiply B by right singular vectors of A   
          (Workspace: need N, prefer N*NRHS) */

	    if (*lwork >= *ldb * *nrhs && *nrhs > 1) {
		lstime_1.opcnt[gemm - 1] += dopbl3_("DGEMM ", n, nrhs, m);
		t1 = dsecnd_();
		dgemm_("T", "N", n, nrhs, m, &c_b145, &a[a_offset], lda, &b[
			b_offset], ldb, &c_b74, &work[1], ldb);
		t2 = dsecnd_();
		lstime_1.timng[gemm - 1] += t2 - t1;
		dlacpy_("F", n, nrhs, &work[1], ldb, &b[b_offset], ldb);
	    } else if (*nrhs > 1) {
		chunk = *lwork / *n;
		i__1 = *nrhs;
		i__2 = chunk;
		for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += 
			i__2) {
/* Computing MIN */
		    i__3 = *nrhs - i__ + 1;
		    bl = min(i__3,chunk);
		    lstime_1.opcnt[gemm - 1] += dopbl3_("DGEMM ", n, &bl, m);
		    t1 = dsecnd_();
		    dgemm_("T", "N", n, &bl, m, &c_b145, &a[a_offset], lda, &
			    b_ref(1, i__), ldb, &c_b74, &work[1], n);
		    t2 = dsecnd_();
		    lstime_1.timng[gemm - 1] += t2 - t1;
		    dlacpy_("F", n, &bl, &work[1], n, &b_ref(1, i__), ldb);
/* L60: */
		}
	    } else {
		lstime_1.opcnt[gemv - 1] += dopbl2_("DGEMV ", m, n, &c__0, &
			c__0);
		t1 = dsecnd_();
		dgemv_("T", m, n, &c_b145, &a[a_offset], lda, &b[b_offset], &
			c__1, &c_b74, &work[1], &c__1);
		t2 = dsecnd_();
		lstime_1.timng[gemv - 1] += t2 - t1;
		dcopy_(n, &work[1], &c__1, &b[b_offset], &c__1);
	    }
	}
    }

/*     Undo scaling */

    if (iascl == 1) {
	lstime_1.opcnt[gelss - 1] += (doublereal) (*n * *nrhs + minmn);
	dlascl_("G", &c__0, &c__0, &anrm, &smlnum, n, nrhs, &b[b_offset], ldb,
		 info);
	dlascl_("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], &
		minmn, info);
    } else if (iascl == 2) {
	lstime_1.opcnt[gelss - 1] += (doublereal) (*n * *nrhs + minmn);
	dlascl_("G", &c__0, &c__0, &anrm, &bignum, n, nrhs, &b[b_offset], ldb,
		 info);
	dlascl_("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], &
		minmn, info);
    }
    if (ibscl == 1) {
	lstime_1.opcnt[gelss - 1] += (doublereal) (*n * *nrhs);
	dlascl_("G", &c__0, &c__0, &smlnum, &bnrm, n, nrhs, &b[b_offset], ldb,
		 info);
    } else if (ibscl == 2) {
	lstime_1.opcnt[gelss - 1] += (doublereal) (*n * *nrhs);
	dlascl_("G", &c__0, &c__0, &bignum, &bnrm, n, nrhs, &b[b_offset], ldb,
		 info);
    }

L70:
    work[1] = (doublereal) maxwrk;
    return 0;

/*     End of DGELSS */

} /* dgelss_ */

#undef b_ref
#undef a_ref