zgels.c

提供矩阵类的函数库

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

/*     Get machine parameters */

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

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

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

/*        Scale matrix norm up to SMLNUM */

	lstime_1.opcnt[gels - 1] += (doublereal) (*m * 6 * *n);
	zlascl_("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] += (doublereal) (*m * 6 * *n);
	zlascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, 
		info);
	iascl = 2;
    } else if (anrm == 0.) {

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

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

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

/*        Scale matrix norm up to SMLNUM */

	lstime_1.opcnt[gels - 1] += (doublereal) (brow * 6 * *nrhs);
	zlascl_("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] += (doublereal) (brow * 6 * *nrhs);
	zlascl_("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, "ZGEQRF", " ", m, n, &c_n1, &c_n1, (ftnlen)6, (
		ftnlen)1);
	lstime_1.opcnt[geqrf - 1] += dopla_("ZGEQRF", m, n, &c__0, &c__0, &nb);
	t1 = dsecnd_();
	i__1 = *lwork - mn;
	zgeqrf_(m, n, &a[a_offset], lda, &work[1], &work[mn + 1], &i__1, info)
		;
	t2 = dsecnd_();
	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, "ZUNMQR", "LC", m, nrhs, n, &c_n1, (ftnlen)6, 
		    (ftnlen)2);
	    lstime_1.opcnt[unmqr - 1] += dopla_("ZUNMQR", m, nrhs, n, &c__0, &
		    nb);
	    t1 = dsecnd_();
	    i__1 = *lwork - mn;
	    zunmqr_("Left", "Conjugate transpose", m, nrhs, n, &a[a_offset], 
		    lda, &work[1], &b[b_offset], ldb, &work[mn + 1], &i__1, 
		    info);
	    t2 = dsecnd_();
	    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] += dopbl3_("ZTRSM ", n, nrhs, &c__0);
	    t1 = dsecnd_();
	    ztrsm_("Left", "Upper", "No transpose", "Non-unit", n, nrhs, &
		    c_b2, &a[a_offset], lda, &b[b_offset], ldb);
	    t2 = dsecnd_();
	    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] += dopbl3_("ZTRSM ", n, nrhs, &c__0);
	    t1 = dsecnd_();
	    ztrsm_("Left", "Upper", "Conjugate transpose", "Non-unit", n, 
		    nrhs, &c_b2, &a[a_offset], lda, &b[b_offset], ldb);
	    t2 = dsecnd_();
	    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., b[i__3].i = 0.;
/* L10: */
		}
/* L20: */
	    }

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

	    nb = ilaenv_(&c__1, "ZUNMQR", "LN", m, nrhs, n, &c_n1, (ftnlen)6, 
		    (ftnlen)2);
	    lstime_1.opcnt[unmqr - 1] += dopla_("ZUNMQR", m, nrhs, n, &c__0, &
		    nb);
	    t1 = dsecnd_();
	    i__1 = *lwork - mn;
	    zunmqr_("Left", "No transpose", m, nrhs, n, &a[a_offset], lda, &
		    work[1], &b[b_offset], ldb, &work[mn + 1], &i__1, info);
	    t2 = dsecnd_();
	    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, "ZGELQF", " ", m, n, &c_n1, &c_n1, (ftnlen)6, (
		ftnlen)1);
	lstime_1.opcnt[gelqf - 1] += dopla_("ZGELQF", m, n, &c__0, &c__0, &nb);
	t1 = dsecnd_();
	i__1 = *lwork - mn;
	zgelqf_(m, n, &a[a_offset], lda, &work[1], &work[mn + 1], &i__1, info)
		;
	t2 = dsecnd_();
	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] += dopbl3_("ZTRSM ", m, nrhs, &c__0);
	    t1 = dsecnd_();
	    ztrsm_("Left", "Lower", "No transpose", "Non-unit", m, nrhs, &
		    c_b2, &a[a_offset], lda, &b[b_offset], ldb);
	    t2 = dsecnd_();
	    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., b[i__3].i = 0.;
/* L30: */
		}
/* L40: */
	    }

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

	    nb = ilaenv_(&c__1, "ZUNMLQ", "LC", n, nrhs, m, &c_n1, (ftnlen)6, 
		    (ftnlen)2);
	    lstime_1.opcnt[unmlq - 1] += dopla_("ZUNMLQ", n, nrhs, m, &c__0, &
		    nb);
	    t1 = dsecnd_();
	    i__1 = *lwork - mn;
	    zunmlq_("Left", "Conjugate transpose", n, nrhs, m, &a[a_offset], 
		    lda, &work[1], &b[b_offset], ldb, &work[mn + 1], &i__1, 
		    info);
	    t2 = dsecnd_();
	    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, "ZUNMLQ", "LN", n, nrhs, m, &c_n1, (ftnlen)6, 
		    (ftnlen)2);
	    lstime_1.opcnt[unmlq - 1] += dopla_("ZUNMLQ", n, nrhs, m, &c__0, &
		    nb);
	    t1 = dsecnd_();
	    i__1 = *lwork - mn;
	    zunmlq_("Left", "No transpose", n, nrhs, m, &a[a_offset], lda, &
		    work[1], &b[b_offset], ldb, &work[mn + 1], &i__1, info);
	    t2 = dsecnd_();
	    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] += dopbl3_("ZTRSM ", m, nrhs, &c__0);
	    t1 = dsecnd_();
	    ztrsm_("Left", "Lower", "Conjugate transpose", "Non-unit", m, 
		    nrhs, &c_b2, &a[a_offset], lda, &b[b_offset], ldb);
	    t2 = dsecnd_();
	    lstime_1.timng[trsm - 1] += t2 - t1;

	    scllen = *m;

	}

    }

/*     Undo scaling */

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

L50:
    d__1 = (doublereal) wsize;
    work[1].r = d__1, work[1].i = 0.;

    return 0;

/*     End of ZGELS */

} /* zgels_ */

#undef b_ref
#undef b_subscr