cgelss.c

提供矩阵类的函数库

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

/*        Bidiagonalize L in WORK(IL)   
          (CWorkspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)   
          (RWorkspace: need M) */

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

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

	    nb = ilaenv_(&c__1, "CUNMBR", "QLC", m, nrhs, m, &c_n1, (ftnlen)6,
		     (ftnlen)3);
	    lstime_1.opcnt[unmbr - 1] += sopla2_("CUNMBR", "QLC", m, nrhs, m, 
		    &c__0, &nb);
	    t1 = second_();
	    i__2 = *lwork - iwork + 1;
	    cunmbr_("Q", "L", "C", m, nrhs, m, &work[il], &ldwork, &work[
		    itauq], &b[b_offset], ldb, &work[iwork], &i__2, info);
	    t2 = second_();
	    lstime_1.timng[unmbr - 1] += t2 - t1;

/*        Generate right bidiagonalizing vectors of R in WORK(IL)   
          (CWorkspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)   
          (RWorkspace: none) */

	    nb = ilaenv_(&c__1, "CUNGBR", "P", m, m, m, &c_n1, (ftnlen)6, (
		    ftnlen)1);
	    lstime_1.opcnt[ungbr - 1] += sopla2_("CUNGBR", "P", m, m, m, &
		    c__0, &nb);
	    t1 = second_();
	    i__2 = *lwork - iwork + 1;
	    cungbr_("P", m, m, m, &work[il], &ldwork, &work[itaup], &work[
		    iwork], &i__2, info);
	    t2 = second_();
	    lstime_1.timng[ungbr - 1] += t2 - t1;
	    irwork = ie + *m;

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

	    latime_1.ops = 0.f;
	    t1 = second_();
	    cbdsqr_("U", m, m, &c__0, nrhs, &s[1], &rwork[ie], &work[il], &
		    ldwork, &a[a_offset], lda, &b[b_offset], ldb, &rwork[
		    irwork], info);
	    t2 = second_();
	    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.f;
/* Computing MAX */
	    r__1 = *rcond * s[1];
	    thr = dmax(r__1,sfmin);
	    if (*rcond < 0.f) {
		lstime_1.opcnt[gelss - 1] += 1.f;
/* Computing MAX */
		r__1 = eps * s[1];
		thr = dmax(r__1,sfmin);
	    }
	    *rank = 0;
	    i__2 = *m;
	    for (i__ = 1; i__ <= i__2; ++i__) {
		if (s[i__] > thr) {
		    lstime_1.opcnt[gelss - 1] += (real) (*nrhs * 6 + 3);
		    csrscl_(nrhs, &s[i__], &b_ref(i__, 1), ldb);
		    ++(*rank);
		} else {
		    claset_("F", &c__1, nrhs, &c_b1, &c_b1, &b_ref(i__, 1), 
			    ldb);
		}
/* L30: */
	    }
	    iwork = il + *m * ldwork;

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

	    if (*lwork >= *ldb * *nrhs + iwork - 1 && *nrhs > 1) {
		lstime_1.opcnt[gemm - 1] += sopbl3_("CGEMM ", m, nrhs, m);
		cgemm_("C", "N", m, nrhs, m, &c_b2, &work[il], &ldwork, &b[
			b_offset], ldb, &c_b1, &work[iwork], ldb);
		clacpy_("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] += sopbl3_("CGEMM ", m, &bl, m);
		    t1 = second_();
		    cgemm_("C", "N", m, &bl, m, &c_b2, &work[il], &ldwork, &
			    b_ref(1, i__), ldb, &c_b1, &work[iwork], n);
		    t2 = second_();
		    lstime_1.timng[gemm - 1] += t2 - t1;
		    clacpy_("G", m, &bl, &work[iwork], n, &b_ref(1, i__), ldb);
/* L40: */
		}
	    } else {
		lstime_1.opcnt[gemv - 1] += sopbl2_("CGEMV ", m, m, &c__0, &
			c__0);
		t1 = second_();
		cgemv_("C", m, m, &c_b2, &work[il], &ldwork, &b_ref(1, 1), &
			c__1, &c_b1, &work[iwork], &c__1);
		t2 = second_();
		lstime_1.timng[gemv - 1] += t2 - t1;
		ccopy_(m, &work[iwork], &c__1, &b_ref(1, 1), &c__1);
	    }

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

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

/*        Multiply transpose(Q) by B   
          (CWorkspace: need M+NRHS, prefer M+NHRS*NB)   
          (RWorkspace: none) */

	    nb = ilaenv_(&c__1, "CUNMLQ", "LC", n, nrhs, m, &c_n1, (ftnlen)6, 
		    (ftnlen)2);
	    lstime_1.opcnt[unmlq - 1] += sopla_("CUNMLQ", n, nrhs, m, &c__0, &
		    nb);
	    t1 = second_();
	    i__1 = *lwork - iwork + 1;
	    cunmlq_("L", "C", n, nrhs, m, &a[a_offset], lda, &work[itau], &b[
		    b_offset], ldb, &work[iwork], &i__1, info);
	    t2 = second_();
	    lstime_1.timng[unmlq - 1] += t2 - t1;

	} else {

/*        Path 2 - remaining underdetermined cases */

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

/*        Bidiagonalize A   
          (CWorkspace: need 3*M, prefer 2*M+(M+N)*NB)   
          (RWorkspace: need N) */

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

/*        Multiply B by transpose of left bidiagonalizing vectors   
          (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)   
          (RWorkspace: none) */

	    nb = ilaenv_(&c__1, "CUNMBR", "QLC", m, nrhs, n, &c_n1, (ftnlen)6,
		     (ftnlen)3);
	    lstime_1.opcnt[unmbr - 1] += sopla2_("CUNMBR", "QLC", m, nrhs, n, 
		    &c__0, &nb);
	    t1 = second_();
	    i__1 = *lwork - iwork + 1;
	    cunmbr_("Q", "L", "C", m, nrhs, n, &a[a_offset], lda, &work[itauq]
		    , &b[b_offset], ldb, &work[iwork], &i__1, info);
	    t2 = second_();
	    lstime_1.timng[unmbr - 1] += t2 - t1;

/*        Generate right bidiagonalizing vectors in A   
          (CWorkspace: need 3*M, prefer 2*M+M*NB)   
          (RWorkspace: none) */

	    nb = ilaenv_(&c__1, "CUNGBR", "P", m, n, m, &c_n1, (ftnlen)6, (
		    ftnlen)1);
	    lstime_1.opcnt[ungbr - 1] += sopla2_("CUNGBR", "P", m, n, m, &
		    c__0, &nb);
	    t1 = second_();
	    i__1 = *lwork - iwork + 1;
	    cungbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &work[
		    iwork], &i__1, info);
	    t2 = second_();
	    lstime_1.timng[ungbr - 1] += t2 - t1;
	    irwork = ie + *m;

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

	    latime_1.ops = 0.f;
	    t1 = second_();
	    cbdsqr_("L", m, n, &c__0, nrhs, &s[1], &rwork[ie], &a[a_offset], 
		    lda, vdum, &c__1, &b[b_offset], ldb, &rwork[irwork], info);
	    t2 = second_();
	    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.f;
/* Computing MAX */
	    r__1 = *rcond * s[1];
	    thr = dmax(r__1,sfmin);
	    if (*rcond < 0.f) {
		lstime_1.opcnt[gelss - 1] += 1.f;
/* Computing MAX */
		r__1 = eps * s[1];
		thr = dmax(r__1,sfmin);
	    }
	    *rank = 0;
	    i__1 = *m;
	    for (i__ = 1; i__ <= i__1; ++i__) {
		if (s[i__] > thr) {
		    lstime_1.opcnt[gelss - 1] += (real) (*nrhs * 6 + 3);
		    csrscl_(nrhs, &s[i__], &b_ref(i__, 1), ldb);
		    ++(*rank);
		} else {
		    claset_("F", &c__1, nrhs, &c_b1, &c_b1, &b_ref(i__, 1), 
			    ldb);
		}
/* L50: */
	    }

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

	    if (*lwork >= *ldb * *nrhs && *nrhs > 1) {
		lstime_1.opcnt[gemm - 1] += sopbl3_("CGEMM ", n, nrhs, m);
		t1 = second_();
		cgemm_("C", "N", n, nrhs, m, &c_b2, &a[a_offset], lda, &b[
			b_offset], ldb, &c_b1, &work[1], ldb);
		t2 = second_();
		lstime_1.timng[gemm - 1] += t2 - t1;
		clacpy_("G", 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] += sopbl3_("CGEMM ", n, &bl, n);
		    t1 = second_();
		    cgemm_("C", "N", n, &bl, m, &c_b2, &a[a_offset], lda, &
			    b_ref(1, i__), ldb, &c_b1, &work[1], n);
		    t2 = second_();
		    lstime_1.timng[gemm - 1] += t2 - t1;
		    clacpy_("F", n, &bl, &work[1], n, &b_ref(1, i__), ldb);
/* L60: */
		}
	    } else {
		lstime_1.opcnt[gelss - 1] += sopbl2_("CGEMV ", m, n, &c__0, &
			c__0);
		t1 = second_();
		cgemv_("C", m, n, &c_b2, &a[a_offset], lda, &b[b_offset], &
			c__1, &c_b1, &work[1], &c__1);
		t2 = second_();
		lstime_1.timng[gemv - 1] += t2 - t1;
		ccopy_(n, &work[1], &c__1, &b[b_offset], &c__1);
	    }
	}
    }

/*     Undo scaling */

    if (iascl == 1) {
	lstime_1.opcnt[gelss - 1] += (real) ((*n * *nrhs + minmn) * 6);
	clascl_("G", &c__0, &c__0, &anrm, &smlnum, n, nrhs, &b[b_offset], ldb,
		 info);
	slascl_("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], &
		minmn, info);
    } else if (iascl == 2) {
	lstime_1.opcnt[gelss - 1] += (real) ((*n * *nrhs + minmn) * 6);
	clascl_("G", &c__0, &c__0, &anrm, &bignum, n, nrhs, &b[b_offset], ldb,
		 info);
	slascl_("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], &
		minmn, info);
    }
    if (ibscl == 1) {
	lstime_1.opcnt[gelss - 1] += (real) (*n * 6 * *nrhs);
	clascl_("G", &c__0, &c__0, &smlnum, &bnrm, n, nrhs, &b[b_offset], ldb,
		 info);
    } else if (ibscl == 2) {
	lstime_1.opcnt[gelss - 1] += (real) (*n * 6 * *nrhs);
	clascl_("G", &c__0, &c__0, &bignum, &bnrm, n, nrhs, &b[b_offset], ldb,
		 info);
    }
L70:
    work[1].r = (real) maxwrk, work[1].i = 0.f;
    return 0;

/*     End of CGELSS */

} /* cgelss_ */

#undef b_ref
#undef b_subscr
#undef a_ref
#undef a_subscr