cgels.c

提供矩阵类的函数库

   Computing MIN */
    i__1 = min(*m,*n);
    if (min(i__1,*nrhs) == 0) {
	i__1 = max(*m,*n);
	claset_("Full", &i__1, nrhs, &c_b1, &c_b1, &b[b_offset], ldb);
	return 0;
    }

/*     Get machine parameters */

    lstime_1.opcnt[gels - 1] += 2.f;
    smlnum = slamch_("S") / slamch_("P");
    bignum = 1.f / smlnum;
    slabad_(&smlnum, &bignum);

/*     Scale A, B if max element outside range [SMLNUM,BIGNUM] */

    anrm = clange_("M", m, n, &a[a_offset], lda, rwork);
    iascl = 0;
    if (anrm > 0.f && anrm < smlnum) {

/*        Scale matrix norm up to SMLNUM */

	lstime_1.opcnt[gels - 1] += (real) (*m * 6 * *n);
	clascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda, 
		info);
	iascl = 1;
    } else if (anrm > bignum) {

/*        Scale matrix norm down to BIGNUM */

	lstime_1.opcnt[gels - 1] += (real) (*m * 6 * *n);
	clascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, 
		info);
	iascl = 2;
    } else if (anrm == 0.f) {

/*        Matrix all zero. Return zero solution. */

	i__1 = max(*m,*n);
	claset_("F", &i__1, nrhs, &c_b1, &c_b1, &b[b_offset], ldb);
	goto L50;
    }

    brow = *m;
    if (tpsd) {
	brow = *n;
    }
    bnrm = clange_("M", &brow, nrhs, &b[b_offset], ldb, rwork);
    ibscl = 0;
    if (bnrm > 0.f && bnrm < smlnum) {

/*        Scale matrix norm up to SMLNUM */

	lstime_1.opcnt[gels - 1] += (real) (brow * 6 * *nrhs);
	clascl_("G", &c__0, &c__0, &bnrm, &smlnum, &brow, nrhs, &b[b_offset], 
		ldb, info);
	ibscl = 1;
    } else if (bnrm > bignum) {

/*        Scale matrix norm down to BIGNUM */

	lstime_1.opcnt[gels - 1] += (real) (brow * 6 * *nrhs);
	clascl_("G", &c__0, &c__0, &bnrm, &bignum, &brow, nrhs, &b[b_offset], 
		ldb, info);
	ibscl = 2;
    }

    if (*m >= *n) {

/*        compute QR factorization of A */

	nb = ilaenv_(&c__1, "CGEQRF", " ", m, n, &c_n1, &c_n1, (ftnlen)6, (
		ftnlen)1);
	lstime_1.opcnt[geqrf - 1] += sopla_("CGEQRF", m, n, &c__0, &c__0, &nb);
	t1 = second_();
	i__1 = *lwork - mn;
	cgeqrf_(m, n, &a[a_offset], lda, &work[1], &work[mn + 1], &i__1, info)
		;
	t2 = second_();
	lstime_1.timng[geqrf - 1] += t2 - t1;

/*        workspace at least N, optimally N*NB */

	if (! tpsd) {

/*           Least-Squares Problem min || A * X - B ||   

             B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS) */

	    nb = ilaenv_(&c__1, "CUNMQR", "LC", m, nrhs, n, &c_n1, (ftnlen)6, 
		    (ftnlen)2);
	    lstime_1.opcnt[unmqr - 1] += sopla_("CUNMQR", m, nrhs, n, &c__0, &
		    nb);
	    t1 = second_();
	    i__1 = *lwork - mn;
	    cunmqr_("Left", "Conjugate transpose", m, nrhs, n, &a[a_offset], 
		    lda, &work[1], &b[b_offset], ldb, &work[mn + 1], &i__1, 
		    info);
	    t2 = second_();
	    lstime_1.timng[unmqr - 1] += t2 - t1;

/*           workspace at least NRHS, optimally NRHS*NB   

             B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS) */

	    lstime_1.opcnt[trsm - 1] += sopbl3_("CTRSM ", n, nrhs, &c__0);
	    t1 = second_();
	    ctrsm_("Left", "Upper", "No transpose", "Non-unit", n, nrhs, &
		    c_b2, &a[a_offset], lda, &b[b_offset], ldb);
	    t2 = second_();
	    lstime_1.timng[trsm - 1] += t2 - t1;

	    scllen = *n;

	} else {

/*           Overdetermined system of equations A' * X = B   

             B(1:N,1:NRHS) := inv(R') * B(1:N,1:NRHS) */

	    lstime_1.opcnt[trsm - 1] += sopbl3_("CTRSM ", n, nrhs, &c__0);
	    t1 = second_();
	    ctrsm_("Left", "Upper", "Conjugate transpose", "Non-unit", n, 
		    nrhs, &c_b2, &a[a_offset], lda, &b[b_offset], ldb);
	    t2 = second_();
	    lstime_1.timng[trsm - 1] += t2 - t1;

/*           B(N+1:M,1:NRHS) = ZERO */

	    i__1 = *nrhs;
	    for (j = 1; j <= i__1; ++j) {
		i__2 = *m;
		for (i__ = *n + 1; i__ <= i__2; ++i__) {
		    i__3 = b_subscr(i__, j);
		    b[i__3].r = 0.f, b[i__3].i = 0.f;
/* L10: */
		}
/* L20: */
	    }

/*           B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS) */

	    nb = ilaenv_(&c__1, "CUNMQR", "LN", m, nrhs, n, &c_n1, (ftnlen)6, 
		    (ftnlen)2);
	    lstime_1.opcnt[unmqr - 1] += sopla_("CUNMQR", m, nrhs, n, &c__0, &
		    nb);
	    t1 = second_();
	    i__1 = *lwork - mn;
	    cunmqr_("Left", "No transpose", m, nrhs, n, &a[a_offset], lda, &
		    work[1], &b[b_offset], ldb, &work[mn + 1], &i__1, info);
	    t2 = second_();
	    lstime_1.timng[unmqr - 1] += t2 - t1;

/*           workspace at least NRHS, optimally NRHS*NB */

	    scllen = *m;

	}

    } else {

/*        Compute LQ factorization of A */

	nb = ilaenv_(&c__1, "CGELQF", " ", m, n, &c_n1, &c_n1, (ftnlen)6, (
		ftnlen)1);
	lstime_1.opcnt[gelqf - 1] += sopla_("CGELQF", m, n, &c__0, &c__0, &nb);
	t1 = second_();
	i__1 = *lwork - mn;
	cgelqf_(m, n, &a[a_offset], lda, &work[1], &work[mn + 1], &i__1, info)
		;
	t2 = second_();
	lstime_1.timng[gelqf - 1] += t2 - t1;

/*        workspace at least M, optimally M*NB. */

	if (! tpsd) {

/*           underdetermined system of equations A * X = B   

             B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS) */

	    lstime_1.opcnt[trsm - 1] += sopbl3_("CTRSM ", m, nrhs, &c__0);
	    t1 = second_();
	    ctrsm_("Left", "Lower", "No transpose", "Non-unit", m, nrhs, &
		    c_b2, &a[a_offset], lda, &b[b_offset], ldb);
	    t2 = second_();
	    lstime_1.timng[trsm - 1] += t2 - t1;

/*           B(M+1:N,1:NRHS) = 0 */

	    i__1 = *nrhs;
	    for (j = 1; j <= i__1; ++j) {
		i__2 = *n;
		for (i__ = *m + 1; i__ <= i__2; ++i__) {
		    i__3 = b_subscr(i__, j);
		    b[i__3].r = 0.f, b[i__3].i = 0.f;
/* L30: */
		}
/* L40: */
	    }

/*           B(1:N,1:NRHS) := Q(1:N,:)' * B(1:M,1:NRHS) */

	    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 - mn;
	    cunmlq_("Left", "Conjugate transpose", n, nrhs, m, &a[a_offset], 
		    lda, &work[1], &b[b_offset], ldb, &work[mn + 1], &i__1, 
		    info);
	    t2 = second_();
	    lstime_1.timng[unmlq - 1] += t2 - t1;

/*           workspace at least NRHS, optimally NRHS*NB */

	    scllen = *n;

	} else {

/*           overdetermined system min || A' * X - B ||   

             B(1:N,1:NRHS) := Q * B(1:N,1:NRHS) */

	    nb = ilaenv_(&c__1, "CUNMLQ", "LN", 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 - mn;
	    cunmlq_("Left", "No transpose", n, nrhs, m, &a[a_offset], lda, &
		    work[1], &b[b_offset], ldb, &work[mn + 1], &i__1, info);
	    t2 = second_();
	    lstime_1.timng[unmlq - 1] += t2 - t1;

/*           workspace at least NRHS, optimally NRHS*NB   

             B(1:M,1:NRHS) := inv(L') * B(1:M,1:NRHS) */

	    lstime_1.opcnt[trsm - 1] += sopbl3_("CTRSM ", m, nrhs, &c__0);
	    t1 = second_();
	    ctrsm_("Left", "Lower", "Conjugate transpose", "Non-unit", m, 
		    nrhs, &c_b2, &a[a_offset], lda, &b[b_offset], ldb);
	    t2 = second_();
	    lstime_1.timng[trsm - 1] += t2 - t1;

	    scllen = *m;

	}

    }

/*     Undo scaling */

    if (iascl == 1) {
	lstime_1.opcnt[gels - 1] += (real) (scllen * 6 * *nrhs);
	clascl_("G", &c__0, &c__0, &anrm, &smlnum, &scllen, nrhs, &b[b_offset]
		, ldb, info);
    } else if (iascl == 2) {
	lstime_1.opcnt[gels - 1] += (real) (scllen * 6 * *nrhs);
	clascl_("G", &c__0, &c__0, &anrm, &bignum, &scllen, nrhs, &b[b_offset]
		, ldb, info);
    }
    if (ibscl == 1) {
	lstime_1.opcnt[gels - 1] += (real) (scllen * 6 * *nrhs);
	clascl_("G", &c__0, &c__0, &smlnum, &bnrm, &scllen, nrhs, &b[b_offset]
		, ldb, info);
    } else if (ibscl == 2) {
	lstime_1.opcnt[gels - 1] += (real) (scllen * 6 * *nrhs);
	clascl_("G", &c__0, &c__0, &bignum, &bnrm, &scllen, nrhs, &b[b_offset]
		, ldb, info);
    }

L50:
    r__1 = (real) wsize;
    work[1].r = r__1, work[1].i = 0.f;

    return 0;

/*     End of CGELS */

} /* cgels_ */

#undef b_ref
#undef b_subscr