support for compact SVD

nonclass.cc

@@ -801,11 +801,11 @@ void singular_decomposition(NRMat<double> &a, NRMat<double> *u, NRVec<double> &s
 	if(m<=0) m=(int)m0;
 	if(n<=0) n=(int)n0;
 	if(n>n0 || m>m0) laerror("bad dimension in singular_decomposition");
-	if (u) if (m > u->nrows() || m> u->ncols())
+	if (u) if (m > u->nrows() || s.size()> u->ncols())
 		laerror("inconsistent dimension of U Mat in singular_decomposition()");
 	if (s.size() < m && s.size() < n) 
 		laerror("inconsistent dimension of S Vec in singular_decomposition()");
-	if (v) if (n > v->nrows() || n > v->ncols())
+	if (v) if (s.size() > v->nrows() || n > v->ncols())
 		laerror("inconsistent dimension of V Mat in singular_decomposition()");
 
 	a.copyonwrite();
@@ -815,31 +815,31 @@ void singular_decomposition(NRMat<double> &a, NRMat<double> *u, NRVec<double> &s
 	
 	// C-order (transposed) input and swap u,v matrices,
 	// v should be transposed at the end
-	char jobu = u ? 'A' : 'N';
+        char jobu = u ? (u->nrows()==u->ncols() ? 'A' : 'S') : 'N';
-	char jobv = v ? 'A' : 'N';
+        char jobv = v ? (v->nrows()==v->ncols() ? 'A' : 'S') : 'N';
+
 	double work0;
 	FINT lwork = -1;
 	FINT r;
 
+	FINT lda=a.ncols();
+	FINT ldu= u ? u->ncols():0;
+	FINT ldv= v ? v->ncols():0;
 #ifdef FORINT
         FINT ntmp = n;
 	FINT mtmp = m;
-	FORNAME(dgesvd)(&jobv, &jobu, &ntmp, &mtmp, a, &n0, s, v?(*v)[0]:0, &n0,
+	FORNAME(dgesvd)(&jobv, &jobu, &ntmp, &mtmp, a, &lda, s, v?(*v)[0]:0, &ldv, u?(*u)[0]:0, &ldu, &work0, &lwork, &r);
-			u?(*u)[0]:0, &m0, &work0, &lwork, &r);
 #else
-	FORNAME(dgesvd)(&jobv, &jobu, &n, &m, a, &n0, s, v?(*v)[0]:0, &n0,
+	FORNAME(dgesvd)(&jobv, &jobu, &n, &m, a, &lda, s, v?(*v)[0]:0, &ldv, u?(*u)[0]:0, &ldu, &work0, &lwork, &r);
-			u?(*u)[0]:0, &m0, &work0, &lwork, &r);
 #endif
 
 	lwork = (FINT) work0;
 	double *work = new double[lwork];
 
 #ifdef FORINT
-	FORNAME(dgesvd)(&jobv, &jobu, &ntmp, &mtmp, a, &n0, s, v?(*v)[0]:0, &n0,
+	FORNAME(dgesvd)(&jobv, &jobu, &ntmp, &mtmp, a, &lda, s, v?(*v)[0]:0, &ldv, u?(*u)[0]:0, &ldu, work, &lwork, &r);
-			u?(*u)[0]:0, &m0, work, &lwork, &r);
 #else
-	FORNAME(dgesvd)(&jobv, &jobu, &n, &m, a, &n0, s, v?(*v)[0]:0, &n0,
+	FORNAME(dgesvd)(&jobv, &jobu, &n, &m, a, &lda, s, v?(*v)[0]:0, &ldv, u?(*u)[0]:0, &ldu, work, &lwork, &r);
-			u?(*u)[0]:0, &m0, work, &lwork, &r);
 #endif
 
 	delete[] work;
@@ -866,11 +866,11 @@ void singular_decomposition(NRMat<std::complex<double> > &a, NRMat<std::complex<
 	if(m<=0) m=(int)m0;
 	if(n<=0) n=(int)n0;
 	if(n>n0 || m>m0) laerror("bad dimension in singular_decomposition");
-	if (u) if (m > u->nrows() || m> u->ncols())
+	if (u) if (m > u->nrows() || s.size()> u->ncols())
 		laerror("inconsistent dimension of U Mat in singular_decomposition()");
 	if (s.size() < m && s.size() < n) 
 		laerror("inconsistent dimension of S Vec in singular_decomposition()");
-	if (v) if (n > v->nrows() || n > v->ncols())
+	if (v) if (s.size() > v->nrows() || n > v->ncols())
 		laerror("inconsistent dimension of V Mat in singular_decomposition()");
 
 	int nmin =  n<m?n:m;
@@ -881,32 +881,32 @@ void singular_decomposition(NRMat<std::complex<double> > &a, NRMat<std::complex<
 	
 	// C-order (transposed) input and swap u,v matrices,
 	// v should be transposed at the end
-	char jobu = u ? 'A' : 'N';
+	char jobu = u ? (u->nrows()==u->ncols() ? 'A' : 'S') : 'N';
-	char jobv = v ? 'A' : 'N';
+	char jobv = v ? (v->nrows()==v->ncols() ? 'A' : 'S') : 'N';
+
 	std::complex<double> work0;
 	FINT lwork = -1;
 	FINT r;
         double *rwork = new double[5*nmin];
 
+        FINT lda=a.ncols();
+        FINT ldu= u ? u->ncols():0;
+        FINT ldv= v ? v->ncols():0;
 #ifdef FORINT
         FINT ntmp = n;
 	FINT mtmp = m;
-	FORNAME(zgesvd)(&jobv, &jobu, &ntmp, &mtmp, a, &n0, s, v?(*v)[0]:0, &n0,
+	FORNAME(zgesvd)(&jobv, &jobu, &ntmp, &mtmp, a, &lda, s, v?(*v)[0]:0, &ldv, u?(*u)[0]:0, &ldu, &work0, &lwork, rwork, &r);
-			u?(*u)[0]:0, &m0, &work0, &lwork, rwork, &r);
 #else
-	FORNAME(zgesvd)(&jobv, &jobu, &n, &m, a, &n0, s, v?(*v)[0]:0, &n0,
+	FORNAME(zgesvd)(&jobv, &jobu, &n, &m, a, &lda, s, v?(*v)[0]:0, &ldv, u?(*u)[0]:0, &ldu, &work0, &lwork, rwork, &r);
-			u?(*u)[0]:0, &m0, &work0, &lwork, rwork, &r);
 #endif
 
 	lwork = (FINT) work0.real();
 	std::complex<double> *work = new std::complex<double>[lwork];
 
 #ifdef FORINT
-	FORNAME(zgesvd)(&jobv, &jobu, &ntmp, &mtmp, a, &n0, s, v?(*v)[0]:0, &n0,
+	FORNAME(zgesvd)(&jobv, &jobu, &ntmp, &mtmp, a, &lda, s, v?(*v)[0]:0, &ldv, u?(*u)[0]:0, &ldu, work, &lwork, rwork, &r);
-			u?(*u)[0]:0, &m0, work, &lwork, rwork, &r);
 #else
-	FORNAME(zgesvd)(&jobv, &jobu, &n, &m, a, &n0, s, v?(*v)[0]:0, &n0,
+	FORNAME(zgesvd)(&jobv, &jobu, &n, &m, a, &lda, s, v?(*v)[0]:0, &ldv, u?(*u)[0]:0, &ldu, work, &lwork, rwork, &r);
-			u?(*u)[0]:0, &m0, work, &lwork, rwork, &r);
 #endif
 
 	delete[] work;

t.cc

@@ -463,11 +463,11 @@ NRMat<double> u(a.nrows(),a.nrows()),v(a.ncols(),a.ncols());
 NRVec<double>s(a.ncols()<a.nrows()?a.ncols():a.nrows());
 singular_decomposition(a,&u,s,&v,0);
 cout <<u;
-NRMat<double> sdiag(0., u.ncols(),v.nrows());
+cout <<s;
-sdiag.diagonalset(s);
-cout <<sdiag;
 cout <<v;
+NRMat<double> sdiag(0., u.ncols(),v.nrows()); sdiag.diagonalset(s);
 cout << "Error "<<(u*sdiag*v-abak).norm()<<endl;
+
 NRMat<double> ai=calcinverse(abak2);
 cout <<"regular inverse "<<ai;
 NRVec<double>ss(s);ss.copyonwrite();
@@ -3481,7 +3481,7 @@ v.printsorted(cout,1,false);
 }
 
 
-if(1)
+if(0)
 {
 //grassmann product of n identical rank=2 tensors in m-dim space
 int n,m;
@@ -3498,6 +3498,29 @@ x.randomize(1);
 
 cout <<x;
 
+int sign;
+for(int j=0; j<m; ++j) 
+for(int i=0; i<m; ++i) 
+	{
+	FLATINDEX I(2);
+	I[0]=i; I[1]=j;
+	cout <<" test "<<i<<" "<<j<<" "<<x.index(&sign,I)<<endl;
+	}
+
+NRMat<double> xm=x.matrix();
+cout <<xm;
+
+Tensor<double> xu=x.unwind_index(0,0);
+NRMat<double> xum=xu.matrix();
+cout <<xum;
+
+Tensor<double> xut=x.unwind_index(0,1);
+NRMat<double> xutm=xut.matrix();
+cout <<xutm;
+
+if((xum-xutm.transpose()).norm()>1e-14) laerror("error in unwinding");
+
+
 //generate antisymmetrizer of even indices, with identity on odd indices
 NRVec<NRVec_from1<int> > indexclasses(1);
 indexclasses[0].resize(n);
@@ -3534,5 +3557,24 @@ y.apply_permutation_algebra(rhsvec,b,false,1.,0.);
 cout <<y;
 }
 
+if(1)
+{
+//compact SVD
+NRMat<double> a;
+cin >>a ;
+NRMat<double> abak=a;
+NRMat<double> abak2=a;
+int min = a.ncols()<a.nrows()?a.ncols():a.nrows();
+NRMat<double> u(a.nrows(),min),vt(min,a.ncols());
+NRVec<double>s(min);
+singular_decomposition(a,&u,s,&vt,0);
+cout <<u;
+cout <<s;
+cout <<vt;
+NRMat<double> sdiag(0., u.ncols(),vt.nrows()); sdiag.diagonalset(s);
+cout << "Error "<<(u*sdiag*vt-abak).norm()<<endl;
+
+}
+
 
 }//main