dgesdd.c

提供矩阵类的函数库

			itauq], &work[itaup], &work[nwork], &i__2, &ierr);

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

		dbdsdc_("U", "I", m, &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 L and VT   
                by right singular vectors of L   
                (Workspace: need M*M+3*M, prefer 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__2 = *lwork - nwork + 1;
		dormbr_("Q", "L", "N", m, m, m, &work[il], &ldwrkl, &work[
			itauq], &u[u_offset], ldu, &work[nwork], &i__2, &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__2 = *lwork - nwork + 1;
		dormbr_("P", "R", "T", m, m, m, &work[il], &ldwrkl, &work[
			itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
			ierr);

/*              Multiply right singular vectors of L in WORK(IL) by   
                Q in A, storing result in VT   
                (Workspace: need M*M) */

		dlacpy_("F", m, m, &vt[vt_offset], ldvt, &work[il], &ldwrkl);
		latime_1.ops += dopbl3_("DGEMM ", m, n, m);
		dgemm_("N", "N", m, n, m, &c_b301, &work[il], &ldwrkl, &a[
			a_offset], lda, &c_b235, &vt[vt_offset], ldvt);

	    } else if (wntqa) {

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

		ivt = 1;

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

		ldwkvt = *m;
		itau = ivt + ldwkvt * *m;
		nwork = itau + *m;

/*              Compute A=L*Q, copying result to VT   
                (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);
		dlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);

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

		nb = ilaenv_(&c__1, "DORGLQ", " ", n, n, m, &c_n1, (ftnlen)6, 
			(ftnlen)1);
		latime_1.ops += dopla_("DORGLQ", n, n, m, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		dorglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &work[
			nwork], &i__2, &ierr);

/*              Produce L in A, zeroing out other entries */

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

/*              Bidiagonalize L in A   
                (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__2 = *lwork - nwork + 1;
		dgebrd_(m, m, &a[a_offset], lda, &s[1], &work[ie], &work[
			itauq], &work[itaup], &work[nwork], &i__2, &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], &ldwkvt, 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 M*M+3*M, prefer 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__2 = *lwork - nwork + 1;
		dormbr_("Q", "L", "N", m, m, m, &a[a_offset], lda, &work[
			itauq], &u[u_offset], ldu, &work[nwork], &i__2, &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__2 = *lwork - nwork + 1;
		dormbr_("P", "R", "T", m, m, m, &a[a_offset], lda, &work[
			itaup], &work[ivt], &ldwkvt, &work[nwork], &i__2, &
			ierr);

/*              Multiply right singular vectors of L in WORK(IVT) by   
                Q in VT, storing result in A   
                (Workspace: need M*M) */

		latime_1.ops += dopbl3_("DGEMM ", m, n, m);
		dgemm_("N", "N", m, n, m, &c_b301, &work[ivt], &ldwkvt, &vt[
			vt_offset], ldvt, &c_b235, &a[a_offset], lda);

/*              Copy right singular vectors of A from A to VT */

		dlacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);

	    }

	} else {

/*           N .LT. MNTHR   

             Path 5t (N greater than M, but not much larger)   
             Reduce to bidiagonal form without LQ decomposition */

	    ie = 1;
	    itauq = ie + *m;
	    itaup = itauq + *m;
	    nwork = itaup + *m;

/*           Bidiagonalize A   
             (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)   

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

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

		dbdsdc_("L", "N", m, &s[1], &work[ie], dum, &c__1, dum, &c__1,
			 dum, idum, &work[nwork], &iwork[1], info);
	    } else if (wntqo) {
		ldwkvt = *m;
		ivt = nwork;
		if (*lwork >= *m * *n + *m * 3 + bdspac) {

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

		    dlaset_("F", m, n, &c_b235, &c_b235, &work[ivt], &ldwkvt);
		    nwork = ivt + ldwkvt * *n;
		} else {

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

		    nwork = ivt + ldwkvt * *m;
		    il = nwork;

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

		    chunk = (*lwork - *m * *m - *m * 3) / *m;
		}

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

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

/*              Overwrite U by left singular vectors of A   
                (Workspace: need M*M+2*M, prefer M*M+M+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__2 = *lwork - nwork + 1;
		dormbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
			itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr);

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

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

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

/*                 Copy right singular vectors of A from WORK(IVT) to A */

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

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

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

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

		    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_("DGEMM ", m, &blk, m);
			dgemm_("N", "N", m, &blk, m, &c_b301, &work[ivt], &
				ldwkvt, &a_ref(1, i__), lda, &c_b235, &work[
				il], m);
			dlacpy_("F", m, &blk, &work[il], m, &a_ref(1, i__), 
				lda);
/* L40: */
		    }
		}
	    } 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 M+BDSPAC) */

		dlaset_("F", m, n, &c_b235, &c_b235, &vt[vt_offset], ldvt);
		dbdsdc_("L", "I", m, &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*M, prefer 2*M+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", m, n, m, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += dopla2_("DORMBR", "PRT", m, n, m, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		dormbr_("P", "R", "T", m, n, m, &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 M+BDSPAC) */

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

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

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

/*              Overwrite U by left singular vectors of A and VT   
                by right singular vectors of A   
                (Workspace: need 2*M+N, prefer 2*M+N*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);
	    }

	}

    }

/*     Undo scaling if necessary */

    if (iscl == 1) {
	if (anrm > bignum) {
	    latime_1.ops += (doublereal) minmn;
	    dlascl_("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], &
		    minmn, &ierr);
	}
	if (anrm < smlnum) {
	    latime_1.ops += (doublereal) minmn;
	    dlascl_("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], &
		    minmn, &ierr);
	}
    }

/*     Return optimal workspace in WORK(1) */

    work[1] = (doublereal) maxwrk;

    return 0;

/*     End of DGESDD */

} /* dgesdd_ */

#undef vt_ref
#undef u_ref
#undef a_ref