sgelsy.c

提供矩阵类的函数库

	    i__1 = 1, i__2 = mn + *n * 3 + 1, i__1 = max(i__1,i__2), i__2 = (
		    mn << 1) + *nrhs;
	    if (*lwork < max(i__1,i__2) && ! lquery) {
		*info = -12;
	    }
	}
    }

    if (*info != 0) {
	i__1 = -(*info);
	xerbla_("SGELSY", &i__1);
	return 0;
    } else if (lquery) {
	return 0;
    }

/*     Quick return if possible   

   Computing MIN */
    i__1 = min(*m,*n);
    if (min(i__1,*nrhs) == 0) {
	*rank = 0;
	return 0;
    }

/*     Get machine parameters */

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

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

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

/*        Scale matrix norm up to SMLNUM */

	lstime_1.opcnt[gelsy - 1] += (real) (*m * *n);
	slascl_("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[gelsy - 1] += (real) (*m * *n);
	slascl_("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);
	slaset_("F", &i__1, nrhs, &c_b31, &c_b31, &b[b_offset], ldb);
	*rank = 0;
	goto L70;
    }

    bnrm = slange_("M", m, nrhs, &b[b_offset], ldb, &work[1]);
    ibscl = 0;
    if (bnrm > 0.f && bnrm < smlnum) {

/*        Scale matrix norm up to SMLNUM */

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

/*        Scale matrix norm down to BIGNUM */

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

/*     Compute QR factorization with column pivoting of A:   
          A * P = Q * R */

    lstime_1.opcnt[geqp3 - 1] += sopla_("SGEQPF", m, n, &c__0, &c__0, &c__0);
    t1 = second_();
    i__1 = *lwork - mn;
    sgeqp3_(m, n, &a[a_offset], lda, &jpvt[1], &work[1], &work[mn + 1], &i__1,
	     info);
    t2 = second_();
    lstime_1.timng[geqp3 - 1] += t2 - t1;
    wsize = mn + work[mn + 1];

/*     workspace: MN+2*N+NB*(N+1).   
       Details of Householder rotations stored in WORK(1:MN).   

       Determine RANK using incremental condition estimation */

    work[ismin] = 1.f;
    work[ismax] = 1.f;
    smax = (r__1 = a_ref(1, 1), dabs(r__1));
    smin = smax;
    if ((r__1 = a_ref(1, 1), dabs(r__1)) == 0.f) {
	*rank = 0;
	i__1 = max(*m,*n);
	slaset_("F", &i__1, nrhs, &c_b31, &c_b31, &b[b_offset], ldb);
	goto L70;
    } else {
	*rank = 1;
    }

L10:
    if (*rank < mn) {
	i__ = *rank + 1;
	latime_1.ops = 0.f;
	slaic1_(&c__2, rank, &work[ismin], &smin, &a_ref(1, i__), &a_ref(i__, 
		i__), &sminpr, &s1, &c1);
	slaic1_(&c__1, rank, &work[ismax], &smax, &a_ref(1, i__), &a_ref(i__, 
		i__), &smaxpr, &s2, &c2);
	lstime_1.opcnt[gelsy - 1] = lstime_1.opcnt[gelsy - 1] + latime_1.ops 
		+ 1.f;

	if (smaxpr * *rcond <= sminpr) {
	    lstime_1.opcnt[gelsy - 1] += (real) (*rank << 1);
	    i__1 = *rank;
	    for (i__ = 1; i__ <= i__1; ++i__) {
		work[ismin + i__ - 1] = s1 * work[ismin + i__ - 1];
		work[ismax + i__ - 1] = s2 * work[ismax + i__ - 1];
/* L20: */
	    }
	    work[ismin + *rank] = c1;
	    work[ismax + *rank] = c2;
	    smin = sminpr;
	    smax = smaxpr;
	    ++(*rank);
	    goto L10;
	}
    }

/*     workspace: 3*MN.   

       Logically partition R = [ R11 R12 ]   
                               [  0  R22 ]   
       where R11 = R(1:RANK,1:RANK)   

       [R11,R12] = [ T11, 0 ] * Y */

    if (*rank < *n) {
	lstime_1.opcnt[tzrzf - 1] += sopla_("STZRQF", rank, n, &c__0, &c__0, &
		c__0);
	t1 = second_();
	i__1 = *lwork - (mn << 1);
	stzrzf_(rank, n, &a[a_offset], lda, &work[mn + 1], &work[(mn << 1) + 
		1], &i__1, info);
	t2 = second_();
	lstime_1.timng[tzrzf - 1] += t2 - t1;
    }

/*     workspace: 2*MN.   
       Details of Householder rotations stored in WORK(MN+1:2*MN)   

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

    lstime_1.opcnt[ormqr - 1] += sopla_("SORMQR", m, nrhs, &mn, &c__0, &c__0);
    t1 = second_();
    i__1 = *lwork - (mn << 1);
    sormqr_("Left", "Transpose", m, nrhs, &mn, &a[a_offset], lda, &work[1], &
	    b[b_offset], ldb, &work[(mn << 1) + 1], &i__1, info);
    t2 = second_();
    lstime_1.timng[ormqr - 1] += t2 - t1;
/* Computing MAX */
    r__1 = wsize, r__2 = (mn << 1) + work[(mn << 1) + 1];
    wsize = dmax(r__1,r__2);

/*     workspace: 2*MN+NB*NRHS.   

       B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS) */

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

    i__1 = *nrhs;
    for (j = 1; j <= i__1; ++j) {
	i__2 = *n;
	for (i__ = *rank + 1; i__ <= i__2; ++i__) {
	    b_ref(i__, j) = 0.f;
/* L30: */
	}
/* L40: */
    }

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

    if (*rank < *n) {
	nb = ilaenv_(&c__1, "SORMRQ", "LT", n, nrhs, rank, &c_n1, (ftnlen)6, (
		ftnlen)2);
	lstime_1.opcnt[ormrz - 1] += sopla_("SORMRQ", n, nrhs, rank, &c__0, &
		nb);
	t1 = second_();
	i__1 = *n - *rank;
	i__2 = *lwork - (mn << 1);
	sormrz_("Left", "Transpose", n, nrhs, rank, &i__1, &a[a_offset], lda, 
		&work[mn + 1], &b[b_offset], ldb, &work[(mn << 1) + 1], &i__2,
		 info);
	t2 = second_();
	lstime_1.timng[ormrz - 1] += t2 - t1;
    }

/*     workspace: 2*MN+NRHS.   

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

    i__1 = *nrhs;
    for (j = 1; j <= i__1; ++j) {
	i__2 = *n;
	for (i__ = 1; i__ <= i__2; ++i__) {
	    work[jpvt[i__]] = b_ref(i__, j);
/* L50: */
	}
	scopy_(n, &work[1], &c__1, &b_ref(1, j), &c__1);
/* L60: */
    }

/*     workspace: N.   

       Undo scaling */

    if (iascl == 1) {
	lstime_1.opcnt[gelsy - 1] += (real) (*n * *nrhs + *rank * *rank);
	slascl_("G", &c__0, &c__0, &anrm, &smlnum, n, nrhs, &b[b_offset], ldb,
		 info);
	slascl_("U", &c__0, &c__0, &smlnum, &anrm, rank, rank, &a[a_offset], 
		lda, info);
    } else if (iascl == 2) {
	lstime_1.opcnt[gelsy - 1] += (real) (*n * *nrhs + *rank * *rank);
	slascl_("G", &c__0, &c__0, &anrm, &bignum, n, nrhs, &b[b_offset], ldb,
		 info);
	slascl_("U", &c__0, &c__0, &bignum, &anrm, rank, rank, &a[a_offset], 
		lda, info);
    }
    if (ibscl == 1) {
	lstime_1.opcnt[gelsy - 1] += (real) (*n * *nrhs);
	slascl_("G", &c__0, &c__0, &smlnum, &bnrm, n, nrhs, &b[b_offset], ldb,
		 info);
    } else if (ibscl == 2) {
	lstime_1.opcnt[gelsy - 1] += (real) (*n * *nrhs);
	slascl_("G", &c__0, &c__0, &bignum, &bnrm, n, nrhs, &b[b_offset], ldb,
		 info);
    }

L70:
    work[1] = (real) lwkopt;

    return 0;

/*     End of SGELSY */

} /* sgelsy_ */

#undef b_ref
#undef a_ref