Better organizing SVD implementation and beginning to work on a REFER…

…ENCES file

AUTHORS

@@ -27,6 +27,3 @@ Jed Brown and Argonne National Laboratory:
 Jeff Hammond and Argonne National Laboratory:
   Several bug fixes and toolchain files for BlueGene/P and BlueGene/Q, and an 
   improvement for Elemental's threading support.
-
-Jack Poulson, Lexing Ying, and The University of Texas at Austin:
-  The parallel linear congruential generator.

REFERENCES

@@ -0,0 +1,28 @@
+NOTE: This list of references will hopefully grow very soon. I would like to
+      move towards making Elemental (and all of my open source efforts) closer
+      to a publication.
+
+[1] Tony F. Chan, "An improved algorithm for computing the Singular Value
+    Decomposition", ACM Transactions on Mathematical Software, Vol. 8, No. 1, 
+    pp. 72--83, 1982.
+    http://www.stat.uchicago.edu/~lekheng/courses/324/chan.pdf
+
+    NOTE: This paper introduced the idea of using a QR decomposition as a first
+          step in the SVD of a non-square matrix in order to accelerate the 
+          computation.
+
+[2] Gene H. Golub and Christian Reinsch, "Singular value decomposition and 
+    least squares solutions", Numerische Mathematik, Vol. 14, No. 5, 
+    pp. 403--420, 1970.
+
+    NOTE: This paper introduced the standard algorithm for computing the SVD.
+
+[3] Zhongxiao Jia, "Using cross-product matrices to compute the SVD",
+    Numerical Algorithms, 42:31--61, 2006.
+    http://faculty.math.tsinghua.edu.cn/~zjia/jianumalgo.pdf
+
+    NOTE: This paper could serve as a foundation for achieving high absolute
+          accuracy in a cross-product based algorithm for computing the SVD.
+          Such an approach should be more scalable than the current 
+          bidiagonalization-based approach.
+

doc/source/core/imports/lapack.rst

@@ -92,8 +92,8 @@ QR-based SVD
 Computes the singular value decomposition of a general matrix by running the 
 QR algorithm on the condensed bidiagonal matrix.
 
-.. cpp:function:: void lapack::SingularValues( int m, int n, R* A, int lda, R* s )
-.. cpp:function:: void lapack::SingularValues( int m, int n, Complex<R>* A, int lda, R* s )
+.. cpp:function:: void lapack::SVD( int m, int n, R* A, int lda, R* s )
+.. cpp:function:: void lapack::SVD( int m, int n, Complex<R>* A, int lda, R* s )
 
 Computes the singular values of a general matrix by running the QR algorithm
 on the condensed bidiagonal matrix.

doc/source/lapack-like/eigen_svd.rst

@@ -374,7 +374,7 @@ flips introduced by negative eigenvalues.
    singular vector matrices, :math:`U` and :math:`V`, such that 
    :math:`A=U \mathrm{diag}(s) V^H`.
 
-.. cpp:function:: void HermitianSingularValues( UpperOrLower uplo, DistMatrix<F>& A, DistMatrix<typename Base<F>::type,VR,STAR>& s )
+.. cpp:function:: void HermitianSVD( UpperOrLower uplo, DistMatrix<F>& A, DistMatrix<typename Base<F>::type,VR,STAR>& s )
 
    Return the singular values of :math:`A` in `s`. Note that the appropriate 
    triangle of `A` is overwritten during computation.
@@ -426,9 +426,9 @@ non-negative entries.
 
    Overwrites `A` with :math:`U`, `s` with the diagonal entries of :math:`\Sigma`, and `V` with :math:`V`. 
 
-.. cpp:function:: void SingularValues( Matrix<F>& A, Matrix<typename Base<F>::type>& s )
+.. cpp:function:: void SVD( Matrix<F>& A, Matrix<typename Base<F>::type>& s )
 
-.. cpp:function:: void SingularValues( DistMatrix<F>& A, DistMatrix<typename Base<F>::type,VR,STAR>& s )
+.. cpp:function:: void SVD( DistMatrix<F>& A, DistMatrix<typename Base<F>::type,VR,STAR>& s )
 
    Forms the singular values of :math:`A` in `s`. Note that `A` is overwritten in order to compute the singular values.
 

examples/lapack-like/SVD.cpp

@@ -54,7 +54,7 @@ main( int argc, char* argv[] )
         // Compute just the singular values 
         DistMatrix<R,VR,STAR> sOnly( g );
         DistMatrix<C> U( A );
-        SingularValues( U, sOnly );
+        SVD( U, sOnly );
 
         // Compute the SVD of A 
         DistMatrix<C> V( g );

include/elemental/core/imports/lapack.hpp

@@ -127,13 +127,13 @@ void QRSVD
   double* s, dcomplex* U, int ldu, dcomplex* VAdj, int ldva );
 
 //
-// Compute the singular values of a general matrix using the QR algorithm
+// Compute the singular values of a general matrix (using the QR algorithm)
 //
 
-void SingularValues( int m, int n, float* A, int lda, float* s );
-void SingularValues( int m, int n, double* A, int lda, double* s );
-void SingularValues( int m, int n, scomplex* A, int lda, float* s );
-void SingularValues( int m, int n, dcomplex* A, int lda, double* s );
+void SVD( int m, int n, float* A, int lda, float* s );
+void SVD( int m, int n, double* A, int lda, double* s );
+void SVD( int m, int n, scomplex* A, int lda, float* s );
+void SVD( int m, int n, dcomplex* A, int lda, double* s );
 
 //
 // Compute the SVD of a bidiagonal matrix using the QR algorithm

include/elemental/lapack-like/ConditionNumber.hpp

@@ -25,7 +25,7 @@ ConditionNumber( const Matrix<F>& A )
 
     Matrix<F> B( A );
     Matrix<R> s;
-    SingularValues( B, s );
+    SVD( B, s );
 
     R cond = 1;
     const int numVals = s.Height();
@@ -48,7 +48,7 @@ ConditionNumber( const DistMatrix<F,U,V>& A )
 
     DistMatrix<F> B( A );
     DistMatrix<R,VR,STAR> s( A.Grid() );
-    SingularValues( B, s );
+    SVD( B, s );
 
     R cond = 1;
     const int numVals = s.Height();

include/elemental/lapack-like/HermitianNorm/KyFan.hpp

@@ -35,14 +35,14 @@ HermitianKyFanNorm( UpperOrLower uplo, const Matrix<F>& A, int k )
 // TODO: Enable sequential MRRR
 /*
 #ifdef HAVE_PMRRR
-    HermitianSingularValues( uplo, B, s );
+    HermitianSVD( uplo, B, s );
 #else
     MakeHermitian( uplo, B );
-    SingularValues( B, s );
+    SVD( B, s );
 #endif // ifdef HAVE_PMRRR
 */
     MakeHermitian( uplo, B );
-    SingularValues( B, s );
+    SVD( B, s );
 
     R norm = 0;
     for( int j=k-1; j>=0; --j )
@@ -67,10 +67,10 @@ HermitianKyFanNorm( UpperOrLower uplo, const DistMatrix<F,U,V>& A, int k )
     DistMatrix<F> B( A );
     DistMatrix<R,VR,STAR> s( A.Grid() );
 #ifdef HAVE_PMRRR
-    HermitianSingularValues( uplo, B, s );
+    HermitianSVD( uplo, B, s );
 #else
     MakeHermitian( uplo, B );
-    SingularValues( B, s );
+    SVD( B, s );
 #endif // ifdef HAVE_PMRRR
 
     R localNorm = 0;

include/elemental/lapack-like/HermitianNorm/Schatten.hpp

@@ -35,14 +35,14 @@ HermitianSchattenNorm
 // TODO: Enable sequential MRRR
 /*
 #ifdef HAVE_PMRRR
-    HermitianSingularValues( uplo, B, s );
+    HermitianSVD( uplo, B, s );
 #else
     MakeHermitian( uplo, B );
-    SingularValues( B, s );
+    SVD( B, s );
 #endif // ifdef HAVE_PMRRR
 */
     MakeHermitian( uplo, B );
-    SingularValues( B, s );
+    SVD( B, s );
 
     // TODO: Think of how to make this more stable
     const int k = s.Height();
@@ -69,10 +69,10 @@ HermitianSchattenNorm
     DistMatrix<F> B( A );
     DistMatrix<R,VR,STAR> s( A.Grid() );
 #ifdef HAVE_PMRRR
-    HermitianSingularValues( uplo, B, s );
+    HermitianSVD( uplo, B, s );
 #else
     MakeHermitian( uplo, B );
-    SingularValues( B, s );
+    SVD( B, s );
 #endif // ifdef HAVE_PMRRR
 
     // TODO: Think of how to make this more stable

include/elemental/lapack-like/HermitianNorm/Two.hpp

@@ -33,14 +33,14 @@ HermitianTwoNorm( UpperOrLower uplo, const Matrix<F>& A )
 // TODO: Enable support for sequential MRRR
 /*
 #ifdef HAVE_PMRRR
-    HermitianSingularValues( uplo, B, s );
+    HermitianSVD( uplo, B, s );
 #else
     MakeHermitian( uplo, B );
-    SingularValues( B, s );
+    SVD( B, s );
 #endif
 */
     MakeHermitian( uplo, B );
-    SingularValues( B, s );
+    SVD( B, s );
     const R norm = InfinityNorm( s );
 #ifndef RELEASE
     PopCallStack();
@@ -59,10 +59,10 @@ HermitianTwoNorm( UpperOrLower uplo, const DistMatrix<F,U,V>& A )
     DistMatrix<F,U,V> B( A );
     DistMatrix<R,VR,STAR> s( A.Grid() );
 #ifdef HAVE_PMRRR
-    HermitianSingularValues( uplo, B, s );
+    HermitianSVD( uplo, B, s );
 #else
     MakeHermitian( uplo, B );
-    SingularValues( B, s );
+    SVD( B, s );
 #endif
     const R norm = InfinityNorm( s );
 #ifndef RELEASE

include/elemental/lapack-like/HermitianSVD.hpp

@@ -69,12 +69,12 @@ inline void HermitianSVD
 }
 
 template<typename F>
-inline void HermitianSingularValues
+inline void HermitianSVD
 ( UpperOrLower uplo, 
   DistMatrix<F>& A, DistMatrix<typename Base<F>::type,VR,STAR>& s )
 {
 #ifndef RELEASE
-    PushCallStack("HermitianSingularValues");
+    PushCallStack("HermitianSVD");
 #endif
     typedef typename Base<F>::type R;
 

include/elemental/lapack-like/Norm/KyFan.hpp

@@ -27,7 +27,7 @@ KyFanNorm( const Matrix<F>& A, int k )
     typedef typename Base<F>::type R;
     Matrix<F> B( A );
     Matrix<R> s;
-    SingularValues( B, s );
+    SVD( B, s );
 
     R norm = 0;
     for( int j=k-1; j>=0; --j )
@@ -51,7 +51,7 @@ KyFanNorm( const DistMatrix<F,U,V>& A, int k )
     typedef typename Base<F>::type R;
     DistMatrix<F> B( A );
     DistMatrix<R,VR,STAR> s( A.Grid() );
-    SingularValues( B, s );
+    SVD( B, s );
 
     R localNorm = 0;
     DistMatrix<R,VR,STAR> sTop( A.Grid() );

include/elemental/lapack-like/Norm/Schatten.hpp

@@ -24,7 +24,7 @@ SchattenNorm( const Matrix<F>& A, typename Base<F>::type p )
     typedef typename Base<F>::type R;
     Matrix<F> B( A );
     Matrix<R> s;
-    SingularValues( B, s );
+    SVD( B, s );
 
     // TODO: Think of how to make this more stable
     const int k = s.Height();
@@ -48,7 +48,7 @@ SchattenNorm( const DistMatrix<F,U,V>& A, typename Base<F>::type p )
     typedef typename Base<F>::type R;
     DistMatrix<F> B( A );
     DistMatrix<R,VR,STAR> s( A.Grid() );
-    SingularValues( B, s );
+    SVD( B, s );
 
     // TODO: Think of how to make this more stable
     const int kLocal = s.LocalHeight();

include/elemental/lapack-like/Norm/Two.hpp

@@ -25,7 +25,7 @@ TwoNorm( const Matrix<F>& A )
     typedef typename Base<F>::type R;
     Matrix<F> B( A );
     Matrix<R> s;
-    SingularValues( B, s );
+    SVD( B, s );
     const R norm = InfinityNorm( s );
 #ifndef RELEASE
     PopCallStack();
@@ -43,7 +43,7 @@ TwoNorm( const DistMatrix<F,U,V>& A )
     typedef typename Base<F>::type R;
     DistMatrix<F> B( A );
     DistMatrix<R,VR,STAR> s( A.Grid() );
-    SingularValues( B, s );
+    SVD( B, s );
     const R norm = InfinityNorm( s );
 #ifndef RELEASE
     PopCallStack();