zgesdd.c

提供矩阵类的函数库

    anrm = zlange_("M", m, n, &a[a_offset], lda, dum);
    iscl = 0;
    if (anrm > 0. && anrm < smlnum) {
	iscl = 1;
	latime_1.ops += (doublereal) (*m * 6 * *n);
	zlascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda, &
		ierr);
    } else if (anrm > bignum) {
	iscl = 1;
	latime_1.ops += (doublereal) (*m * 6 * *n);
	zlascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, &
		ierr);
    }

    if (*m >= *n) {

/*        A has at least as many rows as columns. If A has sufficiently   
          more rows than columns, first reduce using the QR   
          decomposition (if sufficient workspace available) */

	if (*m >= mnthr1) {

	    if (wntqn) {

/*              Path 1 (M much larger than N, JOBZ='N')   
                No singular vectors to be computed */

		itau = 1;
		nwork = itau + *n;

/*              Compute A=Q*R   
                (CWorkspace: need 2*N, prefer N+N*NB)   
                (RWorkspace: need 0) */

		nb = ilaenv_(&c__1, "ZGEQRF", " ", m, n, &c_n1, &c_n1, (
			ftnlen)6, (ftnlen)1);
		latime_1.ops += dopla_("ZGEQRF", m, n, &c__0, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
			i__1, &ierr);

/*              Zero out below R */

		i__1 = *n - 1;
		i__2 = *n - 1;
		zlaset_("L", &i__1, &i__2, &c_b1, &c_b1, &a_ref(2, 1), lda);
		ie = 1;
		itauq = 1;
		itaup = itauq + *n;
		nwork = itaup + *n;

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

		nb = ilaenv_(&c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (
			ftnlen)6, (ftnlen)1);
		latime_1.ops += dopla_("ZGEBRD", n, n, &c__0, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		zgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[
			itauq], &work[itaup], &work[nwork], &i__1, &ierr);
		nrwork = ie + *n;

/*              Perform bidiagonal SVD, compute singular values only   
                (CWorkspace: 0)   
                (RWorkspace: need BDSPAC) */

		dbdsdc_("U", "N", n, &s[1], &rwork[ie], dum, &c__1, dum, &
			c__1, dum, idum, &rwork[nrwork], &iwork[1], info);

	    } else if (wntqo) {

/*              Path 2 (M much larger than N, JOBZ='O')   
                N left singular vectors to be overwritten on A and   
                N right singular vectors to be computed in VT */

		iu = 1;

/*              WORK(IU) is N by N */

		ldwrku = *n;
		ir = iu + ldwrku * *n;
		if (*lwork >= *m * *n + *n * *n + *n * 3) {

/*                 WORK(IR) is M by N */

		    ldwrkr = *m;
		} else {
		    ldwrkr = (*lwork - *n * *n - *n * 3) / *n;
		}
		itau = ir + ldwrkr * *n;
		nwork = itau + *n;

/*              Compute A=Q*R   
                (CWorkspace: need N*N+2*N, prefer M*N+N+N*NB)   
                (RWorkspace: 0) */

		nb = ilaenv_(&c__1, "ZGEQRF", " ", m, n, &c_n1, &c_n1, (
			ftnlen)6, (ftnlen)1);
		latime_1.ops += dopla_("ZGEQRF", m, n, &c__0, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
			i__1, &ierr);

/*              Copy R to WORK( IR ), zeroing out below it */

		zlacpy_("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr);
		i__1 = *n - 1;
		i__2 = *n - 1;
		zlaset_("L", &i__1, &i__2, &c_b1, &c_b1, &work[ir + 1], &
			ldwrkr);

/*              Generate Q in A   
                (CWorkspace: need 2*N, prefer N+N*NB)   
                (RWorkspace: 0) */

		nb = ilaenv_(&c__1, "ZUNGQR", " ", m, n, n, &c_n1, (ftnlen)6, 
			(ftnlen)1);
		latime_1.ops += dopla_("ZUNGQR", m, n, n, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[nwork],
			 &i__1, &ierr);
		ie = 1;
		itauq = itau;
		itaup = itauq + *n;
		nwork = itaup + *n;

/*              Bidiagonalize R in WORK(IR)   
                (CWorkspace: need N*N+3*N, prefer M*N+2*N+2*N*NB)   
                (RWorkspace: need N) */

		nb = ilaenv_(&c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (
			ftnlen)6, (ftnlen)1);
		latime_1.ops += dopla_("ZGEBRD", n, n, &c__0, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		zgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &rwork[ie], &work[
			itauq], &work[itaup], &work[nwork], &i__1, &ierr);

/*              Perform bidiagonal SVD, computing left singular vectors   
                of R in WORK(IRU) and computing right singular vectors   
                of R in WORK(IRVT)   
                (CWorkspace: need 0)   
                (RWorkspace: need BDSPAC) */

		iru = ie + *n;
		irvt = iru + *n * *n;
		nrwork = irvt + *n * *n;
		dbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
			rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1], 
			info);

/*              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)   
                Overwrite WORK(IU) by the left singular vectors of R   
                (CWorkspace: need 2*N*N+3*N, prefer M*N+N*N+2*N+N*NB)   
                (RWorkspace: 0) */

		zlacp2_("F", n, n, &rwork[iru], n, &work[iu], &ldwrku);
		nb = ilaenv_(&c__1, "ZUNMBR", "QLN", n, n, n, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += dopla2_("ZUNMBR", "QLN", n, n, n, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		zunmbr_("Q", "L", "N", n, n, n, &work[ir], &ldwrkr, &work[
			itauq], &work[iu], &ldwrku, &work[nwork], &i__1, &
			ierr);

/*              Copy real matrix RWORK(IRVT) to complex matrix VT   
                Overwrite VT by the right singular vectors of R   
                (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)   
                (RWorkspace: 0) */

		zlacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
		nb = ilaenv_(&c__1, "ZUNMBR", "PRC", n, n, n, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += dopla2_("ZUNMBR", "PRC", n, n, n, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		zunmbr_("P", "R", "C", n, n, n, &work[ir], &ldwrkr, &work[
			itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
			ierr);

/*              Multiply Q in A by left singular vectors of R in   
                WORK(IU), storing result in WORK(IR) and copying to A   
                (CWorkspace: need 2*N*N, prefer N*N+M*N)   
                (RWorkspace: 0) */

		i__1 = *m;
		i__2 = ldwrkr;
		for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += 
			i__2) {
/* Computing MIN */
		    i__3 = *m - i__ + 1;
		    chunk = min(i__3,ldwrkr);
		    latime_1.ops += dopbl3_("ZGEMM ", &chunk, n, n)
			    ;
		    zgemm_("N", "N", &chunk, n, n, &c_b2, &a_ref(i__, 1), lda,
			     &work[iu], &ldwrku, &c_b1, &work[ir], &ldwrkr);
		    zlacpy_("F", &chunk, n, &work[ir], &ldwrkr, &a_ref(i__, 1)
			    , lda);
/* L10: */
		}

	    } else if (wntqs) {

/*              Path 3 (M much larger than N, JOBZ='S')   
                N left singular vectors to be computed in U and   
                N right singular vectors to be computed in VT */

		ir = 1;

/*              WORK(IR) is N by N */

		ldwrkr = *n;
		itau = ir + ldwrkr * *n;
		nwork = itau + *n;

/*              Compute A=Q*R   
                (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)   
                (RWorkspace: 0) */

		nb = ilaenv_(&c__1, "ZGEQRF", " ", m, n, &c_n1, &c_n1, (
			ftnlen)6, (ftnlen)1);
		latime_1.ops += dopla_("ZGEQRF", m, n, &c__0, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
			i__2, &ierr);

/*              Copy R to WORK(IR), zeroing out below it */

		zlacpy_("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr);
		i__2 = *n - 1;
		i__1 = *n - 1;
		zlaset_("L", &i__2, &i__1, &c_b1, &c_b1, &work[ir + 1], &
			ldwrkr);

/*              Generate Q in A   
                (CWorkspace: need 2*N, prefer N+N*NB)   
                (RWorkspace: 0) */

		nb = ilaenv_(&c__1, "ZUNGQR", " ", m, n, n, &c_n1, (ftnlen)6, 
			(ftnlen)1);
		latime_1.ops += dopla_("ZUNGQR", m, n, n, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[nwork],
			 &i__2, &ierr);
		ie = 1;
		itauq = itau;
		itaup = itauq + *n;
		nwork = itaup + *n;

/*              Bidiagonalize R in WORK(IR)   
                (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)   
                (RWorkspace: need N) */

		nb = ilaenv_(&c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (
			ftnlen)6, (ftnlen)1);
		latime_1.ops += dopla_("ZGEBRD", n, n, &c__0, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		zgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &rwork[ie], &work[
			itauq], &work[itaup], &work[nwork], &i__2, &ierr);

/*              Perform bidiagonal SVD, computing left singular vectors   
                of bidiagonal matrix in RWORK(IRU) and computing right   
                singular vectors of bidiagonal matrix in RWORK(IRVT)   
                (CWorkspace: need 0)   
                (RWorkspace: need BDSPAC) */

		iru = ie + *n;
		irvt = iru + *n * *n;
		nrwork = irvt + *n * *n;
		dbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
			rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1], 
			info);

/*              Copy real matrix RWORK(IRU) to complex matrix U   
                Overwrite U by left singular vectors of R   
                (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)   
                (RWorkspace: 0) */

		zlacp2_("F", n, n, &rwork[iru], n, &u[u_offset], ldu);
		nb = ilaenv_(&c__1, "ZUNMBR", "QLN", n, n, n, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += dopla2_("ZUNMBR", "QLN", n, n, n, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		zunmbr_("Q", "L", "N", n, n, n, &work[ir], &ldwrkr, &work[
			itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr);

/*              Copy real matrix RWORK(IRVT) to complex matrix VT   
                Overwrite VT by right singular vectors of R   
                (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)   
                (RWorkspace: 0) */

		zlacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
		nb = ilaenv_(&c__1, "ZUNMBR", "PRC", n, n, n, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += dopla2_("ZUNMBR", "PRC", n, n, n, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		zunmbr_("P", "R", "C", n, n, n, &work[ir], &ldwrkr, &work[
			itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
			ierr);

/*              Multiply Q in A by left singular vectors of R in   
                WORK(IR), storing result in U   
                (CWorkspace: need N*N)   
                (RWorkspace: 0) */

		zlacpy_("F", n, n, &u[u_offset], ldu, &work[ir], &ldwrkr);
		latime_1.ops += dopbl3_("ZGEMM ", m, n, n);
		zgemm_("N", "N", m, n, n, &c_b2, &a[a_offset], lda, &work[ir],
			 &ldwrkr, &c_b1, &u[u_offset], ldu);

	    } else if (wntqa) {

/*              Path 4 (M much larger than N, JOBZ='A')   
                M left singular vectors to be computed in U and   
                N right singular vectors to be computed in VT */

		iu = 1;

/*              WORK(IU) is N by N */

		ldwrku = *n;
		itau = iu + ldwrku * *n;
		nwork = itau + *n;

/*              Compute A=Q*R, copying result to U   
                (CWorkspace: need 2*N, prefer N+N*NB)   
                (RWorkspace: 0) */

		nb = ilaenv_(&c__1, "ZGEQRF", " ", m, n, &c_n1, &c_n1, (
			ftnlen)6, (ftnlen)1);
		latime_1.ops += dopla_("ZGEQRF", m, n, &c__0, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
			i__2, &ierr);
		zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);

/*              Generate Q in U   
                (CWorkspace: need N+M, prefer N+M*NB)   
                (RWorkspace: 0) */

		nb = ilaenv_(&c__1, "ZUNGQR", " ", m, m, n, &c_n1, (ftnlen)6, 
			(ftnlen)1);
		latime_1.ops += dopla_("ZUNGQR", m, m, n, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], &work[nwork],
			 &i__2, &ierr);

/*              Produce R in A, zeroing out below it */

		i__2 = *n - 1;
		i__1 = *n - 1;
		zlaset_("L", &i__2, &i__1, &c_b1, &c_b1, &a_ref(2, 1), lda);
		ie = 1;
		itauq = itau;
		itaup = itauq + *n;
		nwork = itaup + *n;

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