zgesdd.c

提供矩阵类的函数库

		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, &a[a_offset], lda, &work[
			itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
			ierr);

		if (*lwork >= *m * *n + *n * 3) {

/*              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)   
                Overwrite WORK(IU) by left singular vectors of A, copying   
                to A   
                (Cworkspace: need M*N+2*N, prefer M*N+N+N*NB)   
                (Rworkspace: need 0) */

		    zlaset_("F", m, n, &c_b1, &c_b1, &work[iu], &ldwrku);
		    zlacp2_("F", n, n, &rwork[iru], n, &work[iu], &ldwrku);
		    nb = ilaenv_(&c__1, "ZUNMBR", "QLN", m, n, n, &c_n1, (
			    ftnlen)6, (ftnlen)3);
		    latime_1.ops += dopla2_("ZUNMBR", "QLN", m, n, n, &c__0, &
			    nb);
		    i__1 = *lwork - nwork + 1;
		    zunmbr_("Q", "L", "N", m, n, n, &a[a_offset], lda, &work[
			    itauq], &work[iu], &ldwrku, &work[nwork], &i__1, &
			    ierr);
		    zlacpy_("F", m, n, &work[iu], &ldwrku, &a[a_offset], lda);
		} else {

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

		    nb = ilaenv_(&c__1, "ZUNGBR", "Q", m, n, n, &c_n1, (
			    ftnlen)6, (ftnlen)1);
		    latime_1.ops += dopla2_("ZUNGBR", "Q", m, n, n, &c__0, &
			    nb);
		    i__1 = *lwork - nwork + 1;
		    zungbr_("Q", m, n, n, &a[a_offset], lda, &work[itauq], &
			    work[nwork], &i__1, &ierr);

/*                 Multiply Q in A by real matrix RWORK(IRU), storing the   
                   result in WORK(IU), copying to A   
                   (CWorkspace: need N*N, prefer M*N)   
                   (Rworkspace: need 3*N*N, prefer N*N+2*M*N) */

		    nrwork = irvt;
		    i__1 = *m;
		    i__2 = ldwrku;
		    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,ldwrku);
			zlacrm_(&chunk, n, &a_ref(i__, 1), lda, &rwork[iru], 
				n, &work[iu], &ldwrku, &rwork[nrwork]);
			zlacpy_("F", &chunk, n, &work[iu], &ldwrku, &a_ref(
				i__, 1), lda);
/* L30: */
		    }
		}

	    } else if (wntqs) {

/*              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 = nrwork;
		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 A   
                (CWorkspace: need 3*N, prefer 2*N+N*NB)   
                (RWorkspace: 0) */

		zlaset_("F", m, n, &c_b1, &c_b1, &u[u_offset], ldu)
			;
		zlacp2_("F", n, n, &rwork[iru], n, &u[u_offset], ldu);
		nb = ilaenv_(&c__1, "ZUNMBR", "QLN", m, n, n, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += dopla2_("ZUNMBR", "QLN", m, n, n, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		zunmbr_("Q", "L", "N", m, n, n, &a[a_offset], lda, &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 A   
                (CWorkspace: need 3*N, prefer 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, &a[a_offset], lda, &work[
			itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
			ierr);
	    } else {

/*              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 = nrwork;
		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);

/*              Set the right corner of U to identity matrix */

		zlaset_("F", m, m, &c_b1, &c_b1, &u[u_offset], ldu)
			;
		i__2 = *m - *n;
		i__1 = *m - *n;
		zlaset_("F", &i__2, &i__1, &c_b1, &c_b2, &u_ref(*n + 1, *n + 
			1), ldu);

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

		zlacp2_("F", n, n, &rwork[iru], n, &u[u_offset], ldu);
		nb = ilaenv_(&c__1, "ZUNMBR", "QLN", m, m, n, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += dopla2_("ZUNMBR", "QLN", m, m, n, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		zunmbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &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 A   
                (CWorkspace: need 3*N, prefer 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, &a[a_offset], lda, &work[
			itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
			ierr);
	    }

	}

    } else {

/*        A has more columns than rows. If A has sufficiently more   
          columns than rows, first reduce using the LQ decomposition   
          (if sufficient workspace available) */

	if (*n >= mnthr1) {

	    if (wntqn) {

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

		itau = 1;
		nwork = itau + *m;

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

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

/*              Zero out above L */

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

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

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

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

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

	    } else if (wntqo) {

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

		ivt = 1;
		ldwkvt = *m;

/*              WORK(IVT) is M by M */

		il = ivt + ldwkvt * *m;
		if (*lwork >= *m * *n + *m * *m + *m * 3) {

/*                 WORK(IL) M by N */

		    ldwrkl = *m;
		    chunk = *n;
		} else {

/*                 WORK(IL) is M by CHUNK */

		    ldwrkl = *m;
		    chunk = (*lwork - *m * *m - *m * 3) / *m;
		}
		itau = il + ldwrkl * chunk;
		nwork = itau + *m;

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

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

/*              Copy L to WORK(IL), zeroing about above it */

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

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

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

/*              Bidiagonalize L in WORK(IL)   
                (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)   
                (RWorkspace: need M) */

		nb = ilaenv_(&c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (
			ftnlen)6, (ftnlen)1);
		latime_1.ops += dopla_("ZGEBRD", m, m, &c__0, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		zgebrd_(m, m, &work[il], &ldwrkl, &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 + *m;
		irvt = iru + *m * *m;
		nrwork = irvt + *m * *m;
		dbdsdc_("U", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
			rwork[irvt], m, 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 L   
                (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)   
                (RWorkspace: 0) */

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

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

		zlacp2_("F", m, m, &rwork[irvt], m, &work[ivt], &ldwkvt);
		nb = ilaenv_(&c__1, "ZUNMBR", "PRC", m, m, m, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += dopla2_("ZUNMBR", "PRC", m, m, m, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		zunmbr_("P", "R", "C", m, m, m, &work[il], &ldwrkl, &work[
			itaup], &work[ivt], &ldwkvt, &work[nwork], &i__2, &
			ierr);

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

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

	    } else if (wntqs) {

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

		il = 1