zgelsx.c

著名的LAPACK矩阵计算软件包, 是比较新的版本, 一般用到矩阵分解的朋友也许会用到

	lstime_1.opcnt[gelsx - 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[gelsx - 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);
	*rank = 0;
	goto L100;
    }

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

/*        Scale matrix norm up to SMLNUM */

	lstime_1.opcnt[gelsx - 1] += (doublereal) (*m * 6 * *nrhs);
	zlascl_("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[gelsx - 1] += (doublereal) (*m * 6 * *nrhs);
	zlascl_("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[geqpf - 1] += dopla_("ZGEQPF", m, n, &c__0, &c__0, &c__0);
    tim1 = dsecnd_();
    zgeqpf_(m, n, &a[a_offset], lda, &jpvt[1], &work[1], &work[mn + 1], &
	    rwork[1], info);
    tim2 = dsecnd_();
    lstime_1.timng[geqpf - 1] += tim2 - tim1;

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

       Determine RANK using incremental condition estimation */

    i__1 = ismin;
    work[i__1].r = 1., work[i__1].i = 0.;
    i__1 = ismax;
    work[i__1].r = 1., work[i__1].i = 0.;
    smax = z_abs(&a_ref(1, 1));
    smin = smax;
    if (z_abs(&a_ref(1, 1)) == 0.) {
	*rank = 0;
	i__1 = max(*m,*n);
	zlaset_("F", &i__1, nrhs, &c_b1, &c_b1, &b[b_offset], ldb);
	goto L100;
    } else {
	*rank = 1;
    }

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

	if (smaxpr * *rcond <= sminpr) {
	    lstime_1.opcnt[gelsx - 1] += (doublereal) (*rank * 6);
	    i__1 = *rank;
	    for (i__ = 1; i__ <= i__1; ++i__) {
		i__2 = ismin + i__ - 1;
		i__3 = ismin + i__ - 1;
		z__1.r = s1.r * work[i__3].r - s1.i * work[i__3].i, z__1.i = 
			s1.r * work[i__3].i + s1.i * work[i__3].r;
		work[i__2].r = z__1.r, work[i__2].i = z__1.i;
		i__2 = ismax + i__ - 1;
		i__3 = ismax + i__ - 1;
		z__1.r = s2.r * work[i__3].r - s2.i * work[i__3].i, z__1.i = 
			s2.r * work[i__3].i + s2.i * work[i__3].r;
		work[i__2].r = z__1.r, work[i__2].i = z__1.i;
/* L20: */
	    }
	    i__1 = ismin + *rank;
	    work[i__1].r = c1.r, work[i__1].i = c1.i;
	    i__1 = ismax + *rank;
	    work[i__1].r = c2.r, work[i__1].i = c2.i;
	    smin = sminpr;
	    smax = smaxpr;
	    ++(*rank);
	    goto L10;
	}
    }

/*     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[tzrqf - 1] += dopla_("ZTZRQF", rank, n, &c__0, &c__0, &
		c__0);
	tim1 = dsecnd_();
	ztzrqf_(rank, n, &a[a_offset], lda, &work[mn + 1], info);
	tim2 = dsecnd_();
	lstime_1.timng[tzrqf - 1] += tim2 - tim1;
    }

/*     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[unm2r - 1] += dopla_("ZUNMQR", m, nrhs, &mn, &c__0, &c__0);
    tim1 = dsecnd_();
    zunm2r_("Left", "Conjugate transpose", m, nrhs, &mn, &a[a_offset], lda, &
	    work[1], &b[b_offset], ldb, &work[(mn << 1) + 1], info);
    tim2 = dsecnd_();
    lstime_1.timng[unm2r - 1] += tim2 - tim1;

/*     workspace NRHS   

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

    lstime_1.opcnt[trsm - 1] += dopbl3_("ZTRSM ", rank, nrhs, &c__0);
    tim1 = dsecnd_();
    ztrsm_("Left", "Upper", "No transpose", "Non-unit", rank, nrhs, &c_b2, &a[
	    a_offset], lda, &b[b_offset], ldb);
    tim2 = dsecnd_();
    lstime_1.timng[trsm - 1] += tim2 - tim1;

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

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

    if (*rank < *n) {
	lstime_1.opcnt[latzm - 1] += (doublereal) (((*n - *rank) * *nrhs + *
		nrhs + (*n - *rank) * *nrhs << 3) * *rank);
	tim1 = dsecnd_();
	i__1 = *rank;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    i__2 = *n - *rank + 1;
	    d_cnjg(&z__1, &work[mn + i__]);
	    zlatzm_("Left", &i__2, nrhs, &a_ref(i__, *rank + 1), lda, &z__1, &
		    b_ref(i__, 1), &b_ref(*rank + 1, 1), ldb, &work[(mn << 1) 
		    + 1]);
/* L50: */
	}
	tim2 = dsecnd_();
	lstime_1.timng[latzm - 1] += tim2 - tim1;
    }

/*     workspace 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__) {
	    i__3 = (mn << 1) + i__;
	    work[i__3].r = 1., work[i__3].i = 0.;
/* L60: */
	}
	i__2 = *n;
	for (i__ = 1; i__ <= i__2; ++i__) {
	    i__3 = (mn << 1) + i__;
	    if (work[i__3].r == 1. && work[i__3].i == 0.) {
		if (jpvt[i__] != i__) {
		    k = i__;
		    i__3 = b_subscr(k, j);
		    t1.r = b[i__3].r, t1.i = b[i__3].i;
		    i__3 = b_subscr(jpvt[k], j);
		    t2.r = b[i__3].r, t2.i = b[i__3].i;
L70:
		    i__3 = b_subscr(jpvt[k], j);
		    b[i__3].r = t1.r, b[i__3].i = t1.i;
		    i__3 = (mn << 1) + k;
		    work[i__3].r = 0., work[i__3].i = 0.;
		    t1.r = t2.r, t1.i = t2.i;
		    k = jpvt[k];
		    i__3 = b_subscr(jpvt[k], j);
		    t2.r = b[i__3].r, t2.i = b[i__3].i;
		    if (jpvt[k] != i__) {
			goto L70;
		    }
		    i__3 = b_subscr(i__, j);
		    b[i__3].r = t1.r, b[i__3].i = t1.i;
		    i__3 = (mn << 1) + k;
		    work[i__3].r = 0., work[i__3].i = 0.;
		}
	    }
/* L80: */
	}
/* L90: */
    }

/*     Undo scaling */

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

L100:

    return 0;

/*     End of ZGELSX */

} /* zgelsx_ */

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