cgesdd.c

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

          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, "CGEQRF", " ", m, n, &c_n1, &c_n1, (
			ftnlen)6, (ftnlen)1);
		latime_1.ops += sopla_("CGEQRF", m, n, &c__0, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		cgeqrf_(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;
		claset_("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, "CGEBRD", " ", n, n, &c_n1, &c_n1, (
			ftnlen)6, (ftnlen)1);
		latime_1.ops += sopla_("CGEBRD", n, n, &c__0, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		cgebrd_(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) */

		sbdsdc_("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, "CGEQRF", " ", m, n, &c_n1, &c_n1, (
			ftnlen)6, (ftnlen)1);
		latime_1.ops += sopla_("CGEQRF", m, n, &c__0, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		cgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
			i__1, &ierr);

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

		clacpy_("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr);
		i__1 = *n - 1;
		i__2 = *n - 1;
		claset_("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, "CUNGQR", " ", m, n, n, &c_n1, (ftnlen)6, 
			(ftnlen)1);
		latime_1.ops += sopla_("CUNGQR", m, n, n, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		cungqr_(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, "CGEBRD", " ", n, n, &c_n1, &c_n1, (
			ftnlen)6, (ftnlen)1);
		latime_1.ops += sopla_("CGEBRD", n, n, &c__0, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		cgebrd_(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;
		sbdsdc_("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) */

		clacp2_("F", n, n, &rwork[iru], n, &work[iu], &ldwrku);
		nb = ilaenv_(&c__1, "CUNMBR", "QLN", n, n, n, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += sopla2_("CUNMBR", "QLN", n, n, n, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		cunmbr_("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) */

		clacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
		nb = ilaenv_(&c__1, "CUNMBR", "PRC", n, n, n, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += sopla2_("CUNMBR", "PRC", n, n, n, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		cunmbr_("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 += sopbl3_("CGEMM ", &chunk, n, n)
			    ;
		    cgemm_("N", "N", &chunk, n, n, &c_b2, &a_ref(i__, 1), lda,
			     &work[iu], &ldwrku, &c_b1, &work[ir], &ldwrkr);
		    clacpy_("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, "CGEQRF", " ", m, n, &c_n1, &c_n1, (
			ftnlen)6, (ftnlen)1);
		latime_1.ops += sopla_("CGEQRF", m, n, &c__0, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		cgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
			i__2, &ierr);

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

		clacpy_("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr);
		i__2 = *n - 1;
		i__1 = *n - 1;
		claset_("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, "CUNGQR", " ", m, n, n, &c_n1, (ftnlen)6, 
			(ftnlen)1);
		latime_1.ops += sopla_("CUNGQR", m, n, n, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		cungqr_(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, "CGEBRD", " ", n, n, &c_n1, &c_n1, (
			ftnlen)6, (ftnlen)1);
		latime_1.ops += sopla_("CGEBRD", n, n, &c__0, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		cgebrd_(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;
		sbdsdc_("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) */

		clacp2_("F", n, n, &rwork[iru], n, &u[u_offset], ldu);
		nb = ilaenv_(&c__1, "CUNMBR", "QLN", n, n, n, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += sopla2_("CUNMBR", "QLN", n, n, n, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		cunmbr_("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) */

		clacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
		nb = ilaenv_(&c__1, "CUNMBR", "PRC", n, n, n, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += sopla2_("CUNMBR", "PRC", n, n, n, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		cunmbr_("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) */

		clacpy_("F", n, n, &u[u_offset], ldu, &work[ir], &ldwrkr);
		latime_1.ops += sopbl3_("CGEMM ", m, n, n);
		cgemm_("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, "CGEQRF", " ", m, n, &c_n1, &c_n1, (
			ftnlen)6, (ftnlen)1);
		latime_1.ops += sopla_("CGEQRF", m, n, &c__0, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		cgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
			i__2, &ierr);
		clacpy_("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, "CUNGQR", " ", m, m, n, &c_n1, (ftnlen)6, 
			(ftnlen)1);
		latime_1.ops += sopla_("CUNGQR", m, m, n, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		cungqr_(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;
		claset_("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) */

		nb = ilaenv_(&c__1, "CGEBRD", " ", n, n, &c_n1, &c_n1, (
			ftnlen)6, (ftnlen)1);
		latime_1.ops += sopla_("CGEBRD", n, n, &c__0, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		cgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[
			itauq], &work[itaup], &work[nwork], &i__2, &ierr);
		iru = ie + *n;
		irvt = iru + *n * *n;
		nrwork = irvt + *n * *n;

/*              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) */

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