dgesdd.c

提供矩阵类的函数库

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

		ir = 1;

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

		if (*lwork >= *lda * *n + *n * *n + *n * 3 + bdspac) {
		    ldwrkr = *lda;
		} else {
		    ldwrkr = (*lwork - *n * *n - *n * 3 - bdspac) / *n;
		}
		itau = ir + ldwrkr * *n;
		nwork = itau + *n;

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

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

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

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

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

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

/*              Bidiagonalize R in VT, copying result to WORK(IR)   
                (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) */

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

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

		iu = nwork;
		nwork = iu + *n * *n;

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

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

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

		nb = ilaenv_(&c__1, "DORMBR", "QLN", n, n, n, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += dopla2_("DORMBR", "QLN", n, n, n, &c__0, &nb);
		i__1 = *lwork - nwork + 1;
		dormbr_("Q", "L", "N", n, n, n, &work[ir], &ldwrkr, &work[
			itauq], &work[iu], n, &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, &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   
                (Workspace: need 2*N*N, prefer N*N+M*N) */

		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_("DGEMM ", &chunk, n, n)
			    ;
		    dgemm_("N", "N", &chunk, n, n, &c_b301, &a_ref(i__, 1), 
			    lda, &work[iu], n, &c_b235, &work[ir], &ldwrkr);
		    dlacpy_("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   
                (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */

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

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

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

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

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

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

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

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

		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 R and VT   
                by right singular vectors of R   
                (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB) */

		nb = ilaenv_(&c__1, "DORMBR", "QLN", n, n, n, &c_n1, (ftnlen)
			6, (ftnlen)3);
		latime_1.ops += dopla2_("DORMBR", "QLN", n, n, n, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		dormbr_("Q", "L", "N", n, n, n, &work[ir], &ldwrkr, &work[
			itauq], &u[u_offset], ldu, &work[nwork], &i__2, &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__2 = *lwork - nwork + 1;
		dormbr_("P", "R", "T", 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   
                (Workspace: need N*N) */

		dlacpy_("F", n, n, &u[u_offset], ldu, &work[ir], &ldwrkr);
		latime_1.ops += dopbl3_("DGEMM ", m, n, n);
		dgemm_("N", "N", m, n, n, &c_b301, &a[a_offset], lda, &work[
			ir], &ldwrkr, &c_b235, &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   
                (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */

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

/*              Generate Q in U   
                (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
		nb = ilaenv_(&c__1, "DORGQR", " ", m, m, n, &c_n1, (ftnlen)6, 
			(ftnlen)1);
		latime_1.ops += dopla_("DORGQR", m, m, n, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		dorgqr_(m, m, n, &u[u_offset], ldu, &work[itau], &work[nwork],
			 &i__2, &ierr);

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

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

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

		nb = ilaenv_(&c__1, "DGEBRD", " ", n, n, &c_n1, &c_n1, (
			ftnlen)6, (ftnlen)1);
		latime_1.ops += dopla_("DGEBRD", n, n, &c__0, &c__0, &nb);
		i__2 = *lwork - nwork + 1;
		dgebrd_(n, n, &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 WORK(IU) and computing right   
                singular vectors of bidiagonal matrix in VT   
                (Workspace: need N+N*N+BDSPAC) */

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

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

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

/*              Multiply Q in U by left singular vectors of R in   
                WORK(IU), storing result in A   
                (Workspace: need N*N) */

		latime_1.ops += dopbl3_("DGEMM ", m, n, n);
		dgemm_("N", "N", m, n, n, &c_b301, &u[u_offset], ldu, &work[
			iu], &ldwrku, &c_b235, &a[a_offset], lda);

/*              Copy left singular vectors of A from A to U */

		dlacpy_("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);

	    }

	} else {

/*           M .LT. MNTHR   

             Path 5 (M at least N, but not much larger)   
             Reduce to bidiagonal form without QR decomposition */

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

/*           Bidiagonalize A   
             (Workspace: need 3*N+M, prefer 3*N+(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 N+BDSPAC) */

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

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

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

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

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

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

		    ir = nwork;
		    ldwrkr = (*lwork - *n * *n - *n * 3) / *n;
		}
		nwork = iu + ldwrku * *n;

/*              Perform bidiagonal SVD, computing left singular vectors   
                of bidiagonal matrix in WORK(IU) and computing right   
                singular vectors of bidiagonal matrix in VT