dgesdd.c

提供矩阵类的函数库

                (Workspace: need N+N*N+BDSPAC) */

		dbdsdc_("U", "I", n, &s[1], &work[ie], &work[iu], &ldwrku, &
			vt[vt_offset], ldvt, dum, idum, &work[nwork], &iwork[
			1], info);

/*              Overwrite VT by right singular vectors of A   
                (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */

		nb = ilaenv_(&c__1, "DORMBR", "PRT", n, n, n, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += dopla2_("DORMBR", "PRT", n, n, n, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		dormbr_("P", "R", "T", n, n, n, &a[a_offset], lda, &work[
			itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
			ierr);

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

/*                 Overwrite WORK(IU) by left singular vectors of A   
                   (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */

		    nb = ilaenv_(&c__1, "DORMBR", "QLN", m, n, n, &c_n1, (
			    ftnlen)6, (ftnlen)3);
		    latime_1.ops += dopla2_("DORMBR", "QLN", m, n, n, &c__0, &
			    nb);
		    i__2 = *lwork - nwork + 1;
		    dormbr_("Q", "L", "N", m, n, n, &a[a_offset], lda, &work[
			    itauq], &work[iu], &ldwrku, &work[nwork], &i__2, &
			    ierr);

/*                 Copy left singular vectors of A from WORK(IU) to A */

		    dlacpy_("F", m, n, &work[iu], &ldwrku, &a[a_offset], lda);
		} else {

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

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

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

		    i__2 = *m;
		    i__1 = ldwrkr;
		    for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ +=
			     i__1) {
/* Computing MIN */
			i__3 = *m - i__ + 1;
			chunk = min(i__3,ldwrkr);
			latime_1.ops += dopbl3_("DGEMM ", &chunk, n, n);
			dgemm_("N", "N", &chunk, n, n, &c_b301, &a_ref(i__, 1)
				, lda, &work[iu], &ldwrku, &c_b235, &work[ir],
				 &ldwrkr);
			dlacpy_("F", &chunk, n, &work[ir], &ldwrkr, &a_ref(
				i__, 1), lda);
/* L20: */
		    }
		}

	    } else if (wntqs) {

/*              Perform bidiagonal SVD, computing left singular vectors   
                of bidiagonal matrix in U and computing right singular   
                vectors of bidiagonal matrix in VT   
                (Workspace: need N+BDSPAC) */

		dlaset_("F", m, n, &c_b235, &c_b235, &u[u_offset], ldu);
		dbdsdc_("U", "I", n, &s[1], &work[ie], &u[u_offset], ldu, &vt[
			vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1], 
			info);

/*              Overwrite U by left singular vectors of A and VT   
                by right singular vectors of A   
                (Workspace: need 3*N, prefer 2*N+N*NB) */

		nb = ilaenv_(&c__1, "DORMBR", "QLN", m, n, n, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += dopla2_("DORMBR", "QLN", m, n, n, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		dormbr_("Q", "L", "N", m, n, n, &a[a_offset], lda, &work[
			itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
		nb = ilaenv_(&c__1, "DORMBR", "PRT", n, n, n, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += dopla2_("DORMBR", "PRT", n, n, n, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		dormbr_("P", "R", "T", n, n, n, &a[a_offset], lda, &work[
			itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
			ierr);
	    } else if (wntqa) {

/*              Perform bidiagonal SVD, computing left singular vectors   
                of bidiagonal matrix in U and computing right singular   
                vectors of bidiagonal matrix in VT   
                (Workspace: need N+BDSPAC) */

		dlaset_("F", m, m, &c_b235, &c_b235, &u[u_offset], ldu);
		dbdsdc_("U", "I", n, &s[1], &work[ie], &u[u_offset], ldu, &vt[
			vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1], 
			info);

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

		i__1 = *m - *n;
		i__2 = *m - *n;
		dlaset_("F", &i__1, &i__2, &c_b235, &c_b301, &u_ref(*n + 1, *
			n + 1), ldu);

/*              Overwrite U by left singular vectors of A and VT   
                by right singular vectors of A   
                (Workspace: need N*N+2*N+M, prefer N*N+2*N+M*NB) */

		nb = ilaenv_(&c__1, "DORMBR", "QLN", m, m, n, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += dopla2_("DORMBR", "QLN", m, m, n, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		dormbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
			itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
		nb = ilaenv_(&c__1, "DORMBR", "PRT", n, n, m, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += dopla2_("DORMBR", "PRT", n, n, m, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		dormbr_("P", "R", "T", n, n, m, &a[a_offset], lda, &work[
			itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
			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 >= mnthr) {

	    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   
                (Workspace: need 2*M, prefer M+M*NB) */

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

/*              Zero out above L */

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

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

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

/*              Perform bidiagonal SVD, computing singular values only   
                (Workspace: need M+BDSPAC) */

		dbdsdc_("U", "N", m, &s[1], &work[ie], dum, &c__1, dum, &c__1,
			 dum, idum, &work[nwork], &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;

/*              IVT is M by M */

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

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

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

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

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

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

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

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

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

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

		nb = ilaenv_(&c__1, "DGEBRD", " ", m, m, &c_n1, &c_n1, (
			ftnlen)6, (ftnlen)1);
		latime_1.ops += dopla_("DGEBRD", m, m, &c__0, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		dgebrd_(m, m, &work[il], &ldwrkl, &s[1], &work[ie], &work[
			itauq], &work[itaup], &work[nwork], &i__1, &ierr);

/*              Perform bidiagonal SVD, computing left singular vectors   
                of bidiagonal matrix in U, and computing right singular   
                vectors of bidiagonal matrix in WORK(IVT)   
                (Workspace: need M+M*M+BDSPAC) */

		dbdsdc_("U", "I", m, &s[1], &work[ie], &u[u_offset], ldu, &
			work[ivt], m, dum, idum, &work[nwork], &iwork[1], 
			info);

/*              Overwrite U by left singular vectors of L and WORK(IVT)   
                by right singular vectors of L   
                (Workspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB) */

		nb = ilaenv_(&c__1, "DORMBR", "QLN", m, m, m, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += dopla2_("DORMBR", "QLN", m, m, m, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		dormbr_("Q", "L", "N", m, m, m, &work[il], &ldwrkl, &work[
			itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
		nb = ilaenv_(&c__1, "DORMBR", "PRT", m, m, m, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += dopla2_("DORMBR", "PRT", m, m, m, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		dormbr_("P", "R", "T", m, m, m, &work[il], &ldwrkl, &work[
			itaup], &work[ivt], m, &work[nwork], &i__1, &ierr);

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

		i__1 = *n;
		i__2 = chunk;
		for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += 
			i__2) {
/* Computing MIN */
		    i__3 = *n - i__ + 1;
		    blk = min(i__3,chunk);
		    latime_1.ops += dopbl3_("DGEMM ", m, &blk, m);
		    dgemm_("N", "N", m, &blk, m, &c_b301, &work[ivt], m, &
			    a_ref(1, i__), lda, &c_b235, &work[il], &ldwrkl);
		    dlacpy_("F", m, &blk, &work[il], &ldwrkl, &a_ref(1, i__), 
			    lda);
/* L30: */
		}

	    } 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;

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

		ldwrkl = *m;
		itau = il + ldwrkl * *m;
		nwork = itau + *m;

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

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

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

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

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

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

/*              Bidiagonalize L in WORK(IU), copying result to U   
                (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)   

   Computing MAX */
		i__2 = 1, i__1 = ilaenv_(&c__1, "DGEBRD", " ", m, m, &c_n1, &
			c_n1, (ftnlen)6, (ftnlen)1);
		nb = max(i__2,i__1);
		latime_1.ops += dopla_("DGEBRD", m, m, &c__0, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		dgebrd_(m, m, &work[il], &ldwrkl, &s[1], &work[ie], &work[